boost/math/quaternion.hpp
// boost quaternion.hpp header file // (C) Copyright Hubert Holin 2001. // Distributed under the Boost Software License, Version 1.0. (See // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // See http://www.boost.org for updates, documentation, and revision history. #ifndef BOOST_QUATERNION_HPP #define BOOST_QUATERNION_HPP #include <complex> #include <iosfwd> // for the "<<" and ">>" operators #include <sstream> // for the "<<" operator #include <boost/config.hpp> // for BOOST_NO_STD_LOCALE #include <boost/detail/workaround.hpp> #ifndef BOOST_NO_STD_LOCALE #include <locale> // for the "<<" operator #endif /* BOOST_NO_STD_LOCALE */ #include <valarray> #include <boost/math/special_functions/sinc.hpp> // for the Sinus cardinal #include <boost/math/special_functions/sinhc.hpp> // for the Hyperbolic Sinus cardinal namespace boost { namespace math { #if BOOST_WORKAROUND(__GNUC__, < 3) // gcc 2.95.x uses expression templates for valarray calculations, but // the result is not conforming. We need BOOST_GET_VALARRAY to get an // actual valarray result when we need to call a member function #define BOOST_GET_VALARRAY(T,x) ::std::valarray<T>(x) // gcc 2.95.x has an "std::ios" class that is similar to // "std::ios_base", so we just use a #define #define BOOST_IOS_BASE ::std::ios // gcc 2.x ignores function scope using declarations, // put them in the scope of the enclosing namespace instead: using ::std::valarray; using ::std::sqrt; using ::std::cos; using ::std::sin; using ::std::exp; using ::std::cosh; #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ #define BOOST_QUATERNION_ACCESSOR_GENERATOR(type) \ type real() const \ { \ return(a); \ } \ \ quaternion<type> unreal() const \ { \ return(quaternion<type>(static_cast<type>(0),b,c,d)); \ } \ \ type R_component_1() const \ { \ return(a); \ } \ \ type R_component_2() const \ { \ return(b); \ } \ \ type R_component_3() const \ { \ return(c); \ } \ \ type R_component_4() const \ { \ return(d); \ } \ \ ::std::complex<type> C_component_1() const \ { \ return(::std::complex<type>(a,b)); \ } \ \ ::std::complex<type> C_component_2() const \ { \ return(::std::complex<type>(c,d)); \ } #define BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(type) \ template<typename X> \ quaternion<type> & operator = (quaternion<X> const & a_affecter) \ { \ a = static_cast<type>(a_affecter.R_component_1()); \ b = static_cast<type>(a_affecter.R_component_2()); \ c = static_cast<type>(a_affecter.R_component_3()); \ d = static_cast<type>(a_affecter.R_component_4()); \ \ return(*this); \ } \ \ quaternion<type> & operator = (quaternion<type> const & a_affecter) \ { \ a = a_affecter.a; \ b = a_affecter.b; \ c = a_affecter.c; \ d = a_affecter.d; \ \ return(*this); \ } \ \ quaternion<type> & operator = (type const & a_affecter) \ { \ a = a_affecter; \ \ b = c = d = static_cast<type>(0); \ \ return(*this); \ } \ \ quaternion<type> & operator = (::std::complex<type> const & a_affecter) \ { \ a = a_affecter.real(); \ b = a_affecter.imag(); \ \ c = d = static_cast<type>(0); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_DATA_GENERATOR(type) \ type a; \ type b; \ type c; \ type d; template<typename T> class quaternion { public: typedef T value_type; // constructor for H seen as R^4 // (also default constructor) explicit quaternion( T const & requested_a = T(), T const & requested_b = T(), T const & requested_c = T(), T const & requested_d = T()) : a(requested_a), b(requested_b), c(requested_c), d(requested_d) { // nothing to do! } // constructor for H seen as C^2 explicit quaternion( ::std::complex<T> const & z0, ::std::complex<T> const & z1 = ::std::complex<T>()) : a(z0.real()), b(z0.imag()), c(z1.real()), d(z1.imag()) { // nothing to do! } // UNtemplated copy constructor // (this is taken care of by the compiler itself) // templated copy constructor template<typename X> explicit quaternion(quaternion<X> const & a_recopier) : a(static_cast<T>(a_recopier.R_component_1())), b(static_cast<T>(a_recopier.R_component_2())), c(static_cast<T>(a_recopier.R_component_3())), d(static_cast<T>(a_recopier.R_component_4())) { // nothing to do! } // destructor // (this is taken care of by the compiler itself) // accessors // // Note: Like complex number, quaternions do have a meaningful notion of "real part", // but unlike them there is no meaningful notion of "imaginary part". // Instead there is an "unreal part" which itself is a quaternion, and usually // nothing simpler (as opposed to the complex number case). // However, for practicallity, there are accessors for the other components // (these are necessary for the templated copy constructor, for instance). BOOST_QUATERNION_ACCESSOR_GENERATOR(T) // assignment operators BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(T) // other assignment-related operators // // NOTE: Quaternion multiplication is *NOT* commutative; // symbolically, "q *= rhs;" means "q = q * rhs;" // and "q /= rhs;" means "q = q * inverse_of(rhs);" quaternion<T> & operator += (T const & rhs) { T at = a + rhs; // exception guard a = at; return(*this); } quaternion<T> & operator += (::std::complex<T> const & rhs) { T at = a + rhs.real(); // exception guard T bt = b + rhs.imag(); // exception guard a = at; b = bt; return(*this); } template<typename X> quaternion<T> & operator += (quaternion<X> const & rhs) { T at = a + static_cast<T>(rhs.R_component_1()); // exception guard T bt = b + static_cast<T>(rhs.R_component_2()); // exception guard T ct = c + static_cast<T>(rhs.R_component_3()); // exception guard T dt = d + static_cast<T>(rhs.R_component_4()); // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion<T> & operator -= (T const & rhs) { T at = a - rhs; // exception guard a = at; return(*this); } quaternion<T> & operator -= (::std::complex<T> const & rhs) { T at = a - rhs.real(); // exception guard T bt = b - rhs.imag(); // exception guard a = at; b = bt; return(*this); } template<typename X> quaternion<T> & operator -= (quaternion<X> const & rhs) { T at = a - static_cast<T>(rhs.R_component_1()); // exception guard T bt = b - static_cast<T>(rhs.R_component_2()); // exception guard T ct = c - static_cast<T>(rhs.R_component_3()); // exception guard T dt = d - static_cast<T>(rhs.R_component_4()); // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion<T> & operator *= (T const & rhs) { T at = a * rhs; // exception guard T bt = b * rhs; // exception guard T ct = c * rhs; // exception guard T dt = d * rhs; // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion<T> & operator *= (::std::complex<T> const & rhs) { T ar = rhs.real(); T br = rhs.imag(); T at = +a*ar-b*br; T bt = +a*br+b*ar; T ct = +c*ar+d*br; T dt = -c*br+d*ar; a = at; b = bt; c = ct; d = dt; return(*this); } template<typename X> quaternion<T> & operator *= (quaternion<X> const & rhs) { T ar = static_cast<T>(rhs.R_component_1()); T br = static_cast<T>(rhs.R_component_2()); T cr = static_cast<T>(rhs.R_component_3()); T dr = static_cast<T>(rhs.R_component_4()); T at = +a*ar-b*br-c*cr-d*dr; T bt = +a*br+b*ar+c*dr-d*cr; //(a*br+ar*b)+(c*dr-cr*d); T ct = +a*cr-b*dr+c*ar+d*br; //(a*cr+ar*c)+(d*br-dr*b); T dt = +a*dr+b*cr-c*br+d*ar; //(a*dr+ar*d)+(b*cr-br*c); a = at; b = bt; c = ct; d = dt; return(*this); } quaternion<T> & operator /= (T const & rhs) { T at = a / rhs; // exception guard T bt = b / rhs; // exception guard T ct = c / rhs; // exception guard T dt = d / rhs; // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion<T> & operator /= (::std::complex<T> const & rhs) { T ar = rhs.real(); T br = rhs.imag(); T denominator = ar*ar+br*br; T at = (+a*ar+b*br)/denominator; //(a*ar+b*br)/denominator; T bt = (-a*br+b*ar)/denominator; //(ar*b-a*br)/denominator; T ct = (+c*ar-d*br)/denominator; //(ar*c-d*br)/denominator; T dt = (+c*br+d*ar)/denominator; //(ar*d+br*c)/denominator; a = at; b = bt; c = ct; d = dt; return(*this); } template<typename X> quaternion<T> & operator /= (quaternion<X> const & rhs) { T ar = static_cast<T>(rhs.R_component_1()); T br = static_cast<T>(rhs.R_component_2()); T cr = static_cast<T>(rhs.R_component_3()); T dr = static_cast<T>(rhs.R_component_4()); T denominator = ar*ar+br*br+cr*cr+dr*dr; T at = (+a*ar+b*br+c*cr+d*dr)/denominator; //(a*ar+b*br+c*cr+d*dr)/denominator; T bt = (-a*br+b*ar-c*dr+d*cr)/denominator; //((ar*b-a*br)+(cr*d-c*dr))/denominator; T ct = (-a*cr+b*dr+c*ar-d*br)/denominator; //((ar*c-a*cr)+(dr*b-d*br))/denominator; T dt = (-a*dr-b*cr+c*br+d*ar)/denominator; //((ar*d-a*dr)+(br*c-b*cr))/denominator; a = at; b = bt; c = ct; d = dt; return(*this); } protected: BOOST_QUATERNION_MEMBER_DATA_GENERATOR(T) private: }; // declaration of quaternion specialization template<> class quaternion<float>; template<> class quaternion<double>; template<> class quaternion<long double>; // helper templates for converting copy constructors (declaration) namespace detail { template< typename T, typename U > quaternion<T> quaternion_type_converter(quaternion<U> const & rhs); } // implementation of quaternion specialization #define BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(type) \ explicit quaternion( type const & requested_a = static_cast<type>(0), \ type const & requested_b = static_cast<type>(0), \ type const & requested_c = static_cast<type>(0), \ type const & requested_d = static_cast<type>(0)) \ : a(requested_a), \ b(requested_b), \ c(requested_c), \ d(requested_d) \ { \ } \ \ explicit quaternion( ::std::complex<type> const & z0, \ ::std::complex<type> const & z1 = ::std::complex<type>()) \ : a(z0.real()), \ b(z0.imag()), \ c(z1.real()), \ d(z1.imag()) \ { \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \ quaternion<type> & operator += (type const & rhs) \ { \ a += rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \ quaternion<type> & operator += (::std::complex<type> const & rhs) \ { \ a += rhs.real(); \ b += rhs.imag(); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator += (quaternion<X> const & rhs) \ { \ a += static_cast<type>(rhs.R_component_1()); \ b += static_cast<type>(rhs.R_component_2()); \ c += static_cast<type>(rhs.R_component_3()); \ d += static_cast<type>(rhs.R_component_4()); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \ quaternion<type> & operator -= (type const & rhs) \ { \ a -= rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \ quaternion<type> & operator -= (::std::complex<type> const & rhs) \ { \ a -= rhs.real(); \ b -= rhs.imag(); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator -= (quaternion<X> const & rhs) \ { \ a -= static_cast<type>(rhs.R_component_1()); \ b -= static_cast<type>(rhs.R_component_2()); \ c -= static_cast<type>(rhs.R_component_3()); \ d -= static_cast<type>(rhs.R_component_4()); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \ quaternion<type> & operator *= (type const & rhs) \ { \ a *= rhs; \ b *= rhs; \ c *= rhs; \ d *= rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \ quaternion<type> & operator *= (::std::complex<type> const & rhs) \ { \ type ar = rhs.real(); \ type br = rhs.imag(); \ \ type at = +a*ar-b*br; \ type bt = +a*br+b*ar; \ type ct = +c*ar+d*br; \ type dt = -c*br+d*ar; \ \ a = at; \ b = bt; \ c = ct; \ d = dt; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator *= (quaternion<X> const & rhs) \ { \ type ar = static_cast<type>(rhs.R_component_1()); \ type br = static_cast<type>(rhs.R_component_2()); \ type cr = static_cast<type>(rhs.R_component_3()); \ type dr = static_cast<type>(rhs.R_component_4()); \ \ type at = +a*ar-b*br-c*cr-d*dr; \ type bt = +a*br+b*ar+c*dr-d*cr; \ type ct = +a*cr-b*dr+c*ar+d*br; \ type dt = +a*dr+b*cr-c*br+d*ar; \ \ a = at; \ b = bt; \ c = ct; \ d = dt; \ \ return(*this); \ } // There is quite a lot of repetition in the code below. This is intentional. // The last conditional block is the normal form, and the others merely // consist of workarounds for various compiler deficiencies. Hopefuly, when // more compilers are conformant and we can retire support for those that are // not, we will be able to remove the clutter. This is makes the situation // (painfully) explicit. #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \ quaternion<type> & operator /= (type const & rhs) \ { \ a /= rhs; \ b /= rhs; \ c /= rhs; \ d /= rhs; \ \ return(*this); \ } #if defined(__GNUC__) && (__GNUC__ < 3) #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ quaternion<type> & operator /= (::std::complex<type> const & rhs) \ { \ using ::std::valarray; \ \ valarray<type> tr(2); \ \ tr[0] = rhs.real(); \ tr[1] = rhs.imag(); \ \ type mixam = (BOOST_GET_VALARRAY(type,static_cast<type>(1)/abs(tr)).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]; \ tt[1] = -a*tr[1]+b*tr[0]; \ tt[2] = +c*tr[0]-d*tr[1]; \ tt[3] = +c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #elif defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ quaternion<type> & operator /= (::std::complex<type> const & rhs) \ { \ using ::std::valarray; \ using ::std::abs; \ \ valarray<type> tr(2); \ \ tr[0] = rhs.real(); \ tr[1] = rhs.imag(); \ \ type mixam = static_cast<type>(1)/(abs(tr).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]; \ tt[1] = -a*tr[1]+b*tr[0]; \ tt[2] = +c*tr[0]-d*tr[1]; \ tt[3] = +c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #else #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ quaternion<type> & operator /= (::std::complex<type> const & rhs) \ { \ using ::std::valarray; \ \ valarray<type> tr(2); \ \ tr[0] = rhs.real(); \ tr[1] = rhs.imag(); \ \ type mixam = static_cast<type>(1)/(abs(tr).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]; \ tt[1] = -a*tr[1]+b*tr[0]; \ tt[2] = +c*tr[0]-d*tr[1]; \ tt[3] = +c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #endif /* defined(__GNUC__) && (__GNUC__ < 3) */ /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ #if defined(__GNUC__) && (__GNUC__ < 3) #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator /= (quaternion<X> const & rhs) \ { \ using ::std::valarray; \ \ valarray<type> tr(4); \ \ tr[0] = static_cast<type>(rhs.R_component_1()); \ tr[1] = static_cast<type>(rhs.R_component_2()); \ tr[2] = static_cast<type>(rhs.R_component_3()); \ tr[3] = static_cast<type>(rhs.R_component_4()); \ \ type mixam = (BOOST_GET_VALARRAY(type,static_cast<type>(1)/abs(tr)).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \ tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \ tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \ tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #elif defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator /= (quaternion<X> const & rhs) \ { \ using ::std::valarray; \ using ::std::abs; \ \ valarray<type> tr(4); \ \ tr[0] = static_cast<type>(rhs.R_component_1()); \ tr[1] = static_cast<type>(rhs.R_component_2()); \ tr[2] = static_cast<type>(rhs.R_component_3()); \ tr[3] = static_cast<type>(rhs.R_component_4()); \ \ type mixam = static_cast<type>(1)/(abs(tr).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \ tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \ tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \ tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #else #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ template<typename X> \ quaternion<type> & operator /= (quaternion<X> const & rhs) \ { \ using ::std::valarray; \ \ valarray<type> tr(4); \ \ tr[0] = static_cast<type>(rhs.R_component_1()); \ tr[1] = static_cast<type>(rhs.R_component_2()); \ tr[2] = static_cast<type>(rhs.R_component_3()); \ tr[3] = static_cast<type>(rhs.R_component_4()); \ \ type mixam = static_cast<type>(1)/(abs(tr).max)(); \ \ tr *= mixam; \ \ valarray<type> tt(4); \ \ tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \ tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \ tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \ tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \ \ tr *= tr; \ \ tt *= (mixam/tr.sum()); \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #endif /* defined(__GNUC__) && (__GNUC__ < 3) */ /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \ BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \ BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \ BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \ BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \ BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \ BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type) #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \ BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) #define BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \ BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type) template<> class quaternion<float> { public: typedef float value_type; BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(float) // UNtemplated copy constructor // (this is taken care of by the compiler itself) // explicit copy constructors (precision-loosing converters) explicit quaternion(quaternion<double> const & a_recopier) { *this = detail::quaternion_type_converter<float, double>(a_recopier); } explicit quaternion(quaternion<long double> const & a_recopier) { *this = detail::quaternion_type_converter<float, long double>(a_recopier); } // destructor // (this is taken care of by the compiler itself) // accessors // // Note: Like complex number, quaternions do have a meaningful notion of "real part", // but unlike them there is no meaningful notion of "imaginary part". // Instead there is an "unreal part" which itself is a quaternion, and usually // nothing simpler (as opposed to the complex number case). // However, for practicallity, there are accessors for the other components // (these are necessary for the templated copy constructor, for instance). BOOST_QUATERNION_ACCESSOR_GENERATOR(float) // assignment operators BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(float) // other assignment-related operators // // NOTE: Quaternion multiplication is *NOT* commutative; // symbolically, "q *= rhs;" means "q = q * rhs;" // and "q /= rhs;" means "q = q * inverse_of(rhs);" BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(float) protected: BOOST_QUATERNION_MEMBER_DATA_GENERATOR(float) private: }; template<> class quaternion<double> { public: typedef double value_type; BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(double) // UNtemplated copy constructor // (this is taken care of by the compiler itself) // converting copy constructor explicit quaternion(quaternion<float> const & a_recopier) { *this = detail::quaternion_type_converter<double, float>(a_recopier); } // explicit copy constructors (precision-loosing converters) explicit quaternion(quaternion<long double> const & a_recopier) { *this = detail::quaternion_type_converter<double, long double>(a_recopier); } // destructor // (this is taken care of by the compiler itself) // accessors // // Note: Like complex number, quaternions do have a meaningful notion of "real part", // but unlike them there is no meaningful notion of "imaginary part". // Instead there is an "unreal part" which itself is a quaternion, and usually // nothing simpler (as opposed to the complex number case). // However, for practicallity, there are accessors for the other components // (these are necessary for the templated copy constructor, for instance). BOOST_QUATERNION_ACCESSOR_GENERATOR(double) // assignment operators BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(double) // other assignment-related operators // // NOTE: Quaternion multiplication is *NOT* commutative; // symbolically, "q *= rhs;" means "q = q * rhs;" // and "q /= rhs;" means "q = q * inverse_of(rhs);" BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(double) protected: BOOST_QUATERNION_MEMBER_DATA_GENERATOR(double) private: }; template<> class quaternion<long double> { public: typedef long double value_type; BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(long double) // UNtemplated copy constructor // (this is taken care of by the compiler itself) // converting copy constructors explicit quaternion(quaternion<float> const & a_recopier) { *this = detail::quaternion_type_converter<long double, float>(a_recopier); } explicit quaternion(quaternion<double> const & a_recopier) { *this = detail::quaternion_type_converter<long double, double>(a_recopier); } // destructor // (this is taken care of by the compiler itself) // accessors // // Note: Like complex number, quaternions do have a meaningful notion of "real part", // but unlike them there is no meaningful notion of "imaginary part". // Instead there is an "unreal part" which itself is a quaternion, and usually // nothing simpler (as opposed to the complex number case). // However, for practicallity, there are accessors for the other components // (these are necessary for the templated copy constructor, for instance). BOOST_QUATERNION_ACCESSOR_GENERATOR(long double) // assignment operators BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(long double) // other assignment-related operators // // NOTE: Quaternion multiplication is *NOT* commutative; // symbolically, "q *= rhs;" means "q = q * rhs;" // and "q /= rhs;" means "q = q * inverse_of(rhs);" BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(long double) protected: BOOST_QUATERNION_MEMBER_DATA_GENERATOR(long double) private: }; #undef BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2 #undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2 #undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2 #undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2 #undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3 #undef BOOST_QUATERNION_CONSTRUCTOR_GENERATOR #undef BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR #undef BOOST_QUATERNION_MEMBER_DATA_GENERATOR #undef BOOST_QUATERNION_ACCESSOR_GENERATOR // operators #define BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) \ { \ quaternion<T> res(lhs); \ res op##= rhs; \ return(res); \ } #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \ template<typename T> \ inline quaternion<T> operator op (T const & lhs, quaternion<T> const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \ template<typename T> \ inline quaternion<T> operator op (quaternion<T> const & lhs, T const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \ template<typename T> \ inline quaternion<T> operator op (::std::complex<T> const & lhs, quaternion<T> const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \ template<typename T> \ inline quaternion<T> operator op (quaternion<T> const & lhs, ::std::complex<T> const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) \ template<typename T> \ inline quaternion<T> operator op (quaternion<T> const & lhs, quaternion<T> const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR(op) \ BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \ BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \ BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \ BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \ BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) BOOST_QUATERNION_OPERATOR_GENERATOR(+) BOOST_QUATERNION_OPERATOR_GENERATOR(-) BOOST_QUATERNION_OPERATOR_GENERATOR(*) BOOST_QUATERNION_OPERATOR_GENERATOR(/) #undef BOOST_QUATERNION_OPERATOR_GENERATOR #undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_L #undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_R #undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_L #undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_R #undef BOOST_QUATERNION_OPERATOR_GENERATOR_3 #undef BOOST_QUATERNION_OPERATOR_GENERATOR_BODY template<typename T> inline quaternion<T> operator + (quaternion<T> const & q) { return(q); } template<typename T> inline quaternion<T> operator - (quaternion<T> const & q) { return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4())); } template<typename T> inline bool operator == (T const & lhs, quaternion<T> const & rhs) { return ( (rhs.R_component_1() == lhs)&& (rhs.R_component_2() == static_cast<T>(0))&& (rhs.R_component_3() == static_cast<T>(0))&& (rhs.R_component_4() == static_cast<T>(0)) ); } template<typename T> inline bool operator == (quaternion<T> const & lhs, T const & rhs) { return ( (lhs.R_component_1() == rhs)&& (lhs.R_component_2() == static_cast<T>(0))&& (lhs.R_component_3() == static_cast<T>(0))&& (lhs.R_component_4() == static_cast<T>(0)) ); } template<typename T> inline bool operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs) { return ( (rhs.R_component_1() == lhs.real())&& (rhs.R_component_2() == lhs.imag())&& (rhs.R_component_3() == static_cast<T>(0))&& (rhs.R_component_4() == static_cast<T>(0)) ); } template<typename T> inline bool operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs) { return ( (lhs.R_component_1() == rhs.real())&& (lhs.R_component_2() == rhs.imag())&& (lhs.R_component_3() == static_cast<T>(0))&& (lhs.R_component_4() == static_cast<T>(0)) ); } template<typename T> inline bool operator == (quaternion<T> const & lhs, quaternion<T> const & rhs) { return ( (rhs.R_component_1() == lhs.R_component_1())&& (rhs.R_component_2() == lhs.R_component_2())&& (rhs.R_component_3() == lhs.R_component_3())&& (rhs.R_component_4() == lhs.R_component_4()) ); } #define BOOST_QUATERNION_NOT_EQUAL_GENERATOR \ { \ return(!(lhs == rhs)); \ } template<typename T> inline bool operator != (T const & lhs, quaternion<T> const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template<typename T> inline bool operator != (quaternion<T> const & lhs, T const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template<typename T> inline bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template<typename T> inline bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template<typename T> inline bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR #undef BOOST_QUATERNION_NOT_EQUAL_GENERATOR // Note: we allow the following formats, whith a, b, c, and d reals // a // (a), (a,b), (a,b,c), (a,b,c,d) // (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b),(c,d)) #if BOOST_WORKAROUND(__GNUC__, < 3) template<typename T> std::istream & operator >> ( ::std::istream & is, quaternion<T> & q) #else template<typename T, typename charT, class traits> ::std::basic_istream<charT,traits> & operator >> ( ::std::basic_istream<charT,traits> & is, quaternion<T> & q) #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ { #if BOOST_WORKAROUND(__GNUC__, < 3) typedef char charT; #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ #ifdef BOOST_NO_STD_LOCALE #else const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc()); #endif /* BOOST_NO_STD_LOCALE */ T a = T(); T b = T(); T c = T(); T d = T(); ::std::complex<T> u = ::std::complex<T>(); ::std::complex<T> v = ::std::complex<T>(); charT ch = charT(); char cc; is >> ch; // get the first lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == '(') // read "(", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) { is >> ch; // get the second lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == '(') // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) { is.putback(ch); is >> u; // we extract the first and second components a = u.real(); b = u.imag(); if (!is.good()) goto finish; is >> ch; // get the next lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: ((a)) or ((a,b)) { q = quaternion<T>(a,b); } else if (cc == ',') // read "((a)," or "((a,b),", possible: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,)) { is >> v; // we extract the third and fourth components c = v.real(); d = v.imag(); if (!is.good()) goto finish; is >> ch; // get the last lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,)) { q = quaternion<T>(a,b,c,d); } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } else // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)) { is.putback(ch); is >> a; // we extract the first component if (!is.good()) goto finish; is >> ch; // get the third lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: (a) { q = quaternion<T>(a); } else if (cc == ',') // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)) { is >> ch; // get the fourth lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == '(') // read "(a,(", possible: (a,(c)), (a,(c,d)) { is.putback(ch); is >> v; // we extract the third and fourth component c = v.real(); d = v.imag(); if (!is.good()) goto finish; is >> ch; // get the ninth lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: (a,(c)) or (a,(c,d)) { q = quaternion<T>(a,b,c,d); } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } else // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d) { is.putback(ch); is >> b; // we extract the second component if (!is.good()) goto finish; is >> ch; // get the fifth lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: (a,b) { q = quaternion<T>(a,b); } else if (cc == ',') // read "(a,b,", possible: (a,b,c), (a,b,c,d) { is >> c; // we extract the third component if (!is.good()) goto finish; is >> ch; // get the seventh lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: (a,b,c) { q = quaternion<T>(a,b,c); } else if (cc == ',') // read "(a,b,c,", possible: (a,b,c,d) { is >> d; // we extract the fourth component if (!is.good()) goto finish; is >> ch; // get the ninth lexeme if (!is.good()) goto finish; #ifdef BOOST_NO_STD_LOCALE cc = ch; #else cc = ct.narrow(ch, char()); #endif /* BOOST_NO_STD_LOCALE */ if (cc == ')') // format: (a,b,c,d) { q = quaternion<T>(a,b,c,d); } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } } else // error { #if BOOST_WORKAROUND(__GNUC__, < 3) is.setstate(::std::ios::failbit); #else is.setstate(::std::ios_base::failbit); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } } } else // format: a { is.putback(ch); is >> a; // we extract the first component if (!is.good()) goto finish; q = quaternion<T>(a); } finish: return(is); } #if BOOST_WORKAROUND(__GNUC__, < 3) template<typename T> ::std::ostream & operator << ( ::std::ostream & os, quaternion<T> const & q) #else template<typename T, typename charT, class traits> ::std::basic_ostream<charT,traits> & operator << ( ::std::basic_ostream<charT,traits> & os, quaternion<T> const & q) #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ { #if BOOST_WORKAROUND(__GNUC__, < 3) ::std::ostringstream s; #else ::std::basic_ostringstream<charT,traits> s; #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ s.flags(os.flags()); #ifdef BOOST_NO_STD_LOCALE #else s.imbue(os.getloc()); #endif /* BOOST_NO_STD_LOCALE */ s.precision(os.precision()); s << '(' << q.R_component_1() << ',' << q.R_component_2() << ',' << q.R_component_3() << ',' << q.R_component_4() << ')'; return os << s.str(); } // values template<typename T> inline T real(quaternion<T> const & q) { return(q.real()); } template<typename T> inline quaternion<T> unreal(quaternion<T> const & q) { return(q.unreal()); } #define BOOST_QUATERNION_VALARRAY_LOADER \ using ::std::valarray; \ \ valarray<T> temp(4); \ \ temp[0] = q.R_component_1(); \ temp[1] = q.R_component_2(); \ temp[2] = q.R_component_3(); \ temp[3] = q.R_component_4(); template<typename T> inline T sup(quaternion<T> const & q) { #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP using ::std::abs; #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ BOOST_QUATERNION_VALARRAY_LOADER #if BOOST_WORKAROUND(__GNUC__, < 3) return((BOOST_GET_VALARRAY(T, abs(temp)).max)()); #else return((abs(temp).max)()); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } template<typename T> inline T l1(quaternion<T> const & q) { #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP using ::std::abs; #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ BOOST_QUATERNION_VALARRAY_LOADER #if BOOST_WORKAROUND(__GNUC__, < 3) return(BOOST_GET_VALARRAY(T, abs(temp)).sum()); #else return(abs(temp).sum()); #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ } template<typename T> inline T abs(quaternion<T> const & q) { #ifdef BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP using ::std::abs; #endif /* BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP */ using ::std::sqrt; BOOST_QUATERNION_VALARRAY_LOADER #if BOOST_WORKAROUND(__GNUC__, < 3) T maxim = (BOOST_GET_VALARRAY(T, abs(temp)).max)(); // overflow protection #else T maxim = (abs(temp).max)(); // overflow protection #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ if (maxim == static_cast<T>(0)) { return(maxim); } else { T mixam = static_cast<T>(1)/maxim; // prefer multiplications over divisions temp *= mixam; temp *= temp; return(maxim*sqrt(temp.sum())); } //return(sqrt(norm(q))); } #undef BOOST_QUATERNION_VALARRAY_LOADER // Note: This is the Cayley norm, not the Euclidian norm... template<typename T> inline T norm(quaternion<T>const & q) { return(real(q*conj(q))); } template<typename T> inline quaternion<T> conj(quaternion<T> const & q) { return(quaternion<T>( +q.R_component_1(), -q.R_component_2(), -q.R_component_3(), -q.R_component_4())); } template<typename T> inline quaternion<T> spherical( T const & rho, T const & theta, T const & phi1, T const & phi2) { using ::std::cos; using ::std::sin; //T a = cos(theta)*cos(phi1)*cos(phi2); //T b = sin(theta)*cos(phi1)*cos(phi2); //T c = sin(phi1)*cos(phi2); //T d = sin(phi2); T courrant = static_cast<T>(1); T d = sin(phi2); courrant *= cos(phi2); T c = sin(phi1)*courrant; courrant *= cos(phi1); T b = sin(theta)*courrant; T a = cos(theta)*courrant; return(rho*quaternion<T>(a,b,c,d)); } template<typename T> inline quaternion<T> semipolar( T const & rho, T const & alpha, T const & theta1, T const & theta2) { using ::std::cos; using ::std::sin; T a = cos(alpha)*cos(theta1); T b = cos(alpha)*sin(theta1); T c = sin(alpha)*cos(theta2); T d = sin(alpha)*sin(theta2); return(rho*quaternion<T>(a,b,c,d)); } template<typename T> inline quaternion<T> multipolar( T const & rho1, T const & theta1, T const & rho2, T const & theta2) { using ::std::cos; using ::std::sin; T a = rho1*cos(theta1); T b = rho1*sin(theta1); T c = rho2*cos(theta2); T d = rho2*sin(theta2); return(quaternion<T>(a,b,c,d)); } template<typename T> inline quaternion<T> cylindrospherical( T const & t, T const & radius, T const & longitude, T const & latitude) { using ::std::cos; using ::std::sin; T b = radius*cos(longitude)*cos(latitude); T c = radius*sin(longitude)*cos(latitude); T d = radius*sin(latitude); return(quaternion<T>(t,b,c,d)); } template<typename T> inline quaternion<T> cylindrical(T const & r, T const & angle, T const & h1, T const & h2) { using ::std::cos; using ::std::sin; T a = r*cos(angle); T b = r*sin(angle); return(quaternion<T>(a,b,h1,h2)); } // transcendentals // (please see the documentation) template<typename T> inline quaternion<T> exp(quaternion<T> const & q) { using ::std::exp; using ::std::cos; using ::boost::math::sinc_pi; T u = exp(real(q)); T z = abs(unreal(q)); T w = sinc_pi(z); return(u*quaternion<T>(cos(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template<typename T> inline quaternion<T> cos(quaternion<T> const & q) { using ::std::sin; using ::std::cos; using ::std::cosh; using ::boost::math::sinhc_pi; T z = abs(unreal(q)); T w = -sin(q.real())*sinhc_pi(z); return(quaternion<T>(cos(q.real())*cosh(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template<typename T> inline quaternion<T> sin(quaternion<T> const & q) { using ::std::sin; using ::std::cos; using ::std::cosh; using ::boost::math::sinhc_pi; T z = abs(unreal(q)); T w = +cos(q.real())*sinhc_pi(z); return(quaternion<T>(sin(q.real())*cosh(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template<typename T> inline quaternion<T> tan(quaternion<T> const & q) { return(sin(q)/cos(q)); } template<typename T> inline quaternion<T> cosh(quaternion<T> const & q) { return((exp(+q)+exp(-q))/static_cast<T>(2)); } template<typename T> inline quaternion<T> sinh(quaternion<T> const & q) { return((exp(+q)-exp(-q))/static_cast<T>(2)); } template<typename T> inline quaternion<T> tanh(quaternion<T> const & q) { return(sinh(q)/cosh(q)); } template<typename T> quaternion<T> pow(quaternion<T> const & q, int n) { if (n > 1) { int m = n>>1; quaternion<T> result = pow(q, m); result *= result; if (n != (m<<1)) { result *= q; // n odd } return(result); } else if (n == 1) { return(q); } else if (n == 0) { return(quaternion<T>(1)); } else /* n < 0 */ { return(pow(quaternion<T>(1)/q,-n)); } } // helper templates for converting copy constructors (definition) namespace detail { template< typename T, typename U > quaternion<T> quaternion_type_converter(quaternion<U> const & rhs) { return(quaternion<T>( static_cast<T>(rhs.R_component_1()), static_cast<T>(rhs.R_component_2()), static_cast<T>(rhs.R_component_3()), static_cast<T>(rhs.R_component_4()))); } } } } #if BOOST_WORKAROUND(__GNUC__, < 3) #undef BOOST_GET_VALARRAY #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ #endif /* BOOST_QUATERNION_HPP */