--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ossrv_pub/boost_apis/boost/math/quaternion.hpp Tue Feb 02 02:01:42 2010 +0200
@@ -0,0 +1,1924 @@
+// 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 */