// 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 // // See for updates, documentation, and revision history. #ifndef BOOST_QUATERNION_HPP #define BOOST_QUATERNION_HPP #include #include // for the "<<" and ">>" operators #include // for the "<<" operator #include // for BOOST_NO_STD_LOCALE #include #ifndef BOOST_NO_STD_LOCALE #include // for the "<<" operator #endif /* BOOST_NO_STD_LOCALE */ #include // for the Sinus cardinal #include // for the Hyperbolic Sinus cardinal #include namespace boost { namespace math { namespace { template inline T arr_max (const T (&t)[N]) { #if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) using ::std::max; #endif //int assure_array_size_greater_than_zero[N>0?1:-1]; T _max = t[0]; for(size_t i=1;++i inline T arr_max_abs (const T (&t)[N]) { #if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) using ::std::max; #endif using ::std::abs; //int assure_array_size_greater_than_zero[N>0?1:-1]; T _max = abs(t[0]); for(size_t i=1;++i inline T arr_sum (const T (&t)[N]) { //int assure_array_size_greater_than_zero[N>0?1:-1]; T _sum = t[0]; for(size_t i=1;++i inline T arr_sum_abs (const T (&t)[N]) { #if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP) using ::std::abs; #endif //int assure_array_size_greater_than_zero[N>0?1:-1]; T _sum = abs(t[0]); for(size_t i=1;++i unreal() const \ { \ return(quaternion(static_cast(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 C_component_1() const \ { \ return(::std::complex(a,b)); \ } \ \ ::std::complex C_component_2() const \ { \ return(::std::complex(c,d)); \ } #define BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(type) \ template \ quaternion & operator = (quaternion const & a_affecter) \ { \ a = static_cast(a_affecter.R_component_1()); \ b = static_cast(a_affecter.R_component_2()); \ c = static_cast(a_affecter.R_component_3()); \ d = static_cast(a_affecter.R_component_4()); \ \ return(*this); \ } \ \ quaternion & operator = (quaternion const & a_affecter) \ { \ a = a_affecter.a; \ b = a_affecter.b; \ c = a_affecter.c; \ d = a_affecter.d; \ \ return(*this); \ } \ \ quaternion & operator = (type const & a_affecter) \ { \ a = a_affecter; \ \ b = c = d = static_cast(0); \ \ return(*this); \ } \ \ quaternion & operator = (::std::complex const & a_affecter) \ { \ a = a_affecter.real(); \ b = a_affecter.imag(); \ \ c = d = static_cast(0); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_DATA_GENERATOR(type) \ type a; \ type b; \ type c; \ type d; template 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 const & z0, ::std::complex const & z1 = ::std::complex()) : 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 explicit quaternion(quaternion const & a_recopier) : a(static_cast(a_recopier.R_component_1())), b(static_cast(a_recopier.R_component_2())), c(static_cast(a_recopier.R_component_3())), d(static_cast(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 & operator += (T const & rhs) { T at = a + rhs; // exception guard a = at; return(*this); } quaternion & operator += (::std::complex const & rhs) { T at = a + rhs.real(); // exception guard T bt = b + rhs.imag(); // exception guard a = at; b = bt; return(*this); } template quaternion & operator += (quaternion const & rhs) { T at = a + static_cast(rhs.R_component_1()); // exception guard T bt = b + static_cast(rhs.R_component_2()); // exception guard T ct = c + static_cast(rhs.R_component_3()); // exception guard T dt = d + static_cast(rhs.R_component_4()); // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion & operator -= (T const & rhs) { T at = a - rhs; // exception guard a = at; return(*this); } quaternion & operator -= (::std::complex const & rhs) { T at = a - rhs.real(); // exception guard T bt = b - rhs.imag(); // exception guard a = at; b = bt; return(*this); } template quaternion & operator -= (quaternion const & rhs) { T at = a - static_cast(rhs.R_component_1()); // exception guard T bt = b - static_cast(rhs.R_component_2()); // exception guard T ct = c - static_cast(rhs.R_component_3()); // exception guard T dt = d - static_cast(rhs.R_component_4()); // exception guard a = at; b = bt; c = ct; d = dt; return(*this); } quaternion & 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 & operator *= (::std::complex 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 quaternion & operator *= (quaternion const & rhs) { T ar = static_cast(rhs.R_component_1()); T br = static_cast(rhs.R_component_2()); T cr = static_cast(rhs.R_component_3()); T dr = static_cast(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 & 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 & operator /= (::std::complex 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 quaternion & operator /= (quaternion const & rhs) { T ar = static_cast(rhs.R_component_1()); T br = static_cast(rhs.R_component_2()); T cr = static_cast(rhs.R_component_3()); T dr = static_cast(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; template<> class quaternion; template<> class quaternion; // helper templates for converting copy constructors (declaration) namespace detail { template< typename T, typename U > quaternion quaternion_type_converter(quaternion const & rhs); } // implementation of quaternion specialization #define BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(type) \ explicit quaternion( type const & requested_a = static_cast(0), \ type const & requested_b = static_cast(0), \ type const & requested_c = static_cast(0), \ type const & requested_d = static_cast(0)) \ : a(requested_a), \ b(requested_b), \ c(requested_c), \ d(requested_d) \ { \ } \ \ explicit quaternion( ::std::complex const & z0, \ ::std::complex const & z1 = ::std::complex()) \ : a(z0.real()), \ b(z0.imag()), \ c(z1.real()), \ d(z1.imag()) \ { \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \ quaternion & operator += (type const & rhs) \ { \ a += rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \ quaternion & operator += (::std::complex const & rhs) \ { \ a += rhs.real(); \ b += rhs.imag(); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) \ template \ quaternion & operator += (quaternion const & rhs) \ { \ a += static_cast(rhs.R_component_1()); \ b += static_cast(rhs.R_component_2()); \ c += static_cast(rhs.R_component_3()); \ d += static_cast(rhs.R_component_4()); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \ quaternion & operator -= (type const & rhs) \ { \ a -= rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \ quaternion & operator -= (::std::complex const & rhs) \ { \ a -= rhs.real(); \ b -= rhs.imag(); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) \ template \ quaternion & operator -= (quaternion const & rhs) \ { \ a -= static_cast(rhs.R_component_1()); \ b -= static_cast(rhs.R_component_2()); \ c -= static_cast(rhs.R_component_3()); \ d -= static_cast(rhs.R_component_4()); \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \ quaternion & operator *= (type const & rhs) \ { \ a *= rhs; \ b *= rhs; \ c *= rhs; \ d *= rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \ quaternion & operator *= (::std::complex 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 \ quaternion & operator *= (quaternion const & rhs) \ { \ type ar = static_cast(rhs.R_component_1()); \ type br = static_cast(rhs.R_component_2()); \ type cr = static_cast(rhs.R_component_3()); \ type dr = static_cast(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 & operator /= (type const & rhs) \ { \ a /= rhs; \ b /= rhs; \ c /= rhs; \ d /= rhs; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \ quaternion & operator /= (::std::complex const & rhs) \ { \ \ type tr[2]; \ \ tr[0] = rhs.real(); \ tr[1] = rhs.imag(); \ \ type mixam = static_cast(1)/(arr_max_abs(tr)); \ \ tr[0] *= mixam; \ tr[1] *= mixam; \ \ 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[0] *= tr[0]; \ tr[1] *= tr[1]; \ \ type mixam_tr_sum = mixam/arr_sum(tr); \ tt[0] *= mixam_tr_sum; \ tt[1] *= mixam_tr_sum; \ tt[2] *= mixam_tr_sum; \ tt[3] *= mixam_tr_sum; \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \ template \ quaternion & operator /= (quaternion const & rhs) \ { \ \ type tr[4]; \ \ tr[0] = static_cast(rhs.R_component_1()); \ tr[1] = static_cast(rhs.R_component_2()); \ tr[2] = static_cast(rhs.R_component_3()); \ tr[3] = static_cast(rhs.R_component_4()); \ \ type mixam = static_cast(1)/(arr_max_abs(tr)); \ \ tr[0] *= mixam; \ tr[1] *= mixam; \ tr[2] *= mixam; \ tr[3] *= mixam; \ \ 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[0] *= tr[0]; \ tr[1] *= tr[1]; \ tr[2] *= tr[2]; \ tr[3] *= tr[3]; \ \ type mixam_tr_sum = mixam/arr_sum(tr); \ tt[0] *= mixam_tr_sum; \ tt[1] *= mixam_tr_sum; \ tt[2] *= mixam_tr_sum; \ tt[3] *= mixam_tr_sum; \ \ a = tt[0]; \ b = tt[1]; \ c = tt[2]; \ d = tt[3]; \ \ return(*this); \ } #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 { 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 const & a_recopier) { *this = detail::quaternion_type_converter(a_recopier); } explicit quaternion(quaternion const & a_recopier) { *this = detail::quaternion_type_converter(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 { 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 const & a_recopier) { *this = detail::quaternion_type_converter(a_recopier); } // explicit copy constructors (precision-loosing converters) explicit quaternion(quaternion const & a_recopier) { *this = detail::quaternion_type_converter(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 { 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 const & a_recopier) { *this = detail::quaternion_type_converter(a_recopier); } explicit quaternion(quaternion const & a_recopier) { *this = detail::quaternion_type_converter(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 res(lhs); \ res op##= rhs; \ return(res); \ } #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \ template \ inline quaternion operator op (T const & lhs, quaternion const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \ template \ inline quaternion operator op (quaternion const & lhs, T const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \ template \ inline quaternion operator op (::std::complex const & lhs, quaternion const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \ template \ inline quaternion operator op (quaternion const & lhs, ::std::complex const & rhs) \ BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) #define BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) \ template \ inline quaternion operator op (quaternion const & lhs, quaternion 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 inline quaternion operator + (quaternion const & q) { return(q); } template inline quaternion operator - (quaternion const & q) { return(quaternion(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4())); } template inline bool operator == (T const & lhs, quaternion const & rhs) { return ( (rhs.R_component_1() == lhs)&& (rhs.R_component_2() == static_cast(0))&& (rhs.R_component_3() == static_cast(0))&& (rhs.R_component_4() == static_cast(0)) ); } template inline bool operator == (quaternion const & lhs, T const & rhs) { return ( (lhs.R_component_1() == rhs)&& (lhs.R_component_2() == static_cast(0))&& (lhs.R_component_3() == static_cast(0))&& (lhs.R_component_4() == static_cast(0)) ); } template inline bool operator == (::std::complex const & lhs, quaternion const & rhs) { return ( (rhs.R_component_1() == lhs.real())&& (rhs.R_component_2() == lhs.imag())&& (rhs.R_component_3() == static_cast(0))&& (rhs.R_component_4() == static_cast(0)) ); } template inline bool operator == (quaternion const & lhs, ::std::complex const & rhs) { return ( (lhs.R_component_1() == rhs.real())&& (lhs.R_component_2() == rhs.imag())&& (lhs.R_component_3() == static_cast(0))&& (lhs.R_component_4() == static_cast(0)) ); } template inline bool operator == (quaternion const & lhs, quaternion 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 inline bool operator != (T const & lhs, quaternion const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template inline bool operator != (quaternion const & lhs, T const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template inline bool operator != (::std::complex const & lhs, quaternion const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template inline bool operator != (quaternion const & lhs, ::std::complex const & rhs) BOOST_QUATERNION_NOT_EQUAL_GENERATOR template inline bool operator != (quaternion const & lhs, quaternion 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)) template ::std::basic_istream & operator >> ( ::std::basic_istream & is, quaternion & q) { #ifdef BOOST_NO_STD_LOCALE #else const ::std::ctype & ct = ::std::use_facet< ::std::ctype >(is.getloc()); #endif /* BOOST_NO_STD_LOCALE */ T a = T(); T b = T(); T c = T(); T d = T(); ::std::complex u = ::std::complex(); ::std::complex v = ::std::complex(); 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(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(a,b,c,d); } else // error { is.setstate(::std::ios_base::failbit); } } else // error { is.setstate(::std::ios_base::failbit); } } 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(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(a,b,c,d); } else // error { is.setstate(::std::ios_base::failbit); } } 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(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(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(a,b,c,d); } else // error { is.setstate(::std::ios_base::failbit); } } else // error { is.setstate(::std::ios_base::failbit); } } else // error { is.setstate(::std::ios_base::failbit); } } } else // error { is.setstate(::std::ios_base::failbit); } } } else // format: a { is.putback(ch); is >> a; // we extract the first component if (!is.good()) goto finish; q = quaternion(a); } finish: return(is); } template ::std::basic_ostream & operator << ( ::std::basic_ostream & os, quaternion const & q) { ::std::basic_ostringstream s; 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 inline T real(quaternion const & q) { return(q.real()); } template inline quaternion unreal(quaternion const & q) { return(q.unreal()); } #define BOOST_QUATERNION_ARRAY_LOADER \ \ 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 inline T sup(quaternion const & q) { BOOST_QUATERNION_ARRAY_LOADER return(arr_max_abs(temp)); } template inline T l1(quaternion const & q) { BOOST_QUATERNION_ARRAY_LOADER return(sum_abs(temp)); } template inline T abs(quaternion const & q) { using ::std::sqrt; BOOST_QUATERNION_ARRAY_LOADER T maxim = (max_abs(temp)); // overflow protection if (maxim == static_cast(0)) { return(maxim); } else { T mixam = static_cast(1)/maxim; // prefer multiplications over divisions temp[0] *= mixam; temp[1] *= mixam; temp[2] *= mixam; temp[3] *= mixam; temp[0] *= temp[0]; temp[1] *= temp[1]; temp[2] *= temp[2]; temp[3] *= temp[3]; return(maxim*sqrt(sum(temp))); } } #undef BOOST_QUATERNION_ARRAY_LOADER // Note: This is the Cayley norm, not the Euclidian norm... template inline T norm(quaternionconst & q) { return(real(q*conj(q))); } template inline quaternion conj(quaternion const & q) { return(quaternion( +q.R_component_1(), -q.R_component_2(), -q.R_component_3(), -q.R_component_4())); } template inline quaternion 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(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(a,b,c,d)); } template inline quaternion 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(a,b,c,d)); } template inline quaternion 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(a,b,c,d)); } template inline quaternion 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,b,c,d)); } template inline quaternion 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(a,b,h1,h2)); } // transcendentals // (please see the documentation) template inline quaternion exp(quaternion 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(cos(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template inline quaternion cos(quaternion 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(cos(q.real())*cosh(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template inline quaternion sin(quaternion 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(sin(q.real())*cosh(z), w*q.R_component_2(), w*q.R_component_3(), w*q.R_component_4())); } template inline quaternion tan(quaternion const & q) { return(sin(q)/cos(q)); } template inline quaternion cosh(quaternion const & q) { return((exp(+q)+exp(-q))/static_cast(2)); } template inline quaternion sinh(quaternion const & q) { return((exp(+q)-exp(-q))/static_cast(2)); } template inline quaternion tanh(quaternion const & q) { return(sinh(q)/cosh(q)); } template quaternion pow(quaternion const & q, int n) { if (n > 1) { int m = n>>1; quaternion 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(static_cast(1))); } else /* n < 0 */ { return(pow(quaternion(static_cast(1))/q,-n)); } } // helper templates for converting copy constructors (definition) namespace detail { template< typename T, typename U > quaternion quaternion_type_converter(quaternion const & rhs) { return(quaternion( static_cast(rhs.R_component_1()), static_cast(rhs.R_component_2()), static_cast(rhs.R_component_3()), static_cast(rhs.R_component_4()))); } } } } #endif /* BOOST_QUATERNION_HPP */