boost/math/quaternion.hpp

0001 //  boost quaternion.hpp header file
0002
0003 //  (C) Copyright Hubert Holin 2001.
0004 //  Distributed under the Boost Software License, Version 1.0. (See
0005 //  accompanying file LICENSE_1_0.txt or copy at
0006 //  http://www.boost.org/LICENSE_1_0.txt)
0007
0008 // See http://www.boost.org for updates, documentation, and revision history.
0009
0010 #ifndef BOOST_QUATERNION_HPP
0011 #define BOOST_QUATERNION_HPP
0012
0013 #include <boost/math_fwd.hpp>
0014 #include <boost/math/tools/config.hpp>
0015 #include <locale>                                    // for the "<<" operator
0016
0017 #include <complex>
0018 #include <iosfwd>                                    // for the "<<" and ">>" operators
0019 #include <sstream>                                    // for the "<<" operator
0020
0021 #include <boost/math/special_functions/sinc.hpp>    // for the Sinus cardinal
0022 #include <boost/math/special_functions/sinhc.hpp>    // for the Hyperbolic Sinus cardinal
0023 #include <boost/math/tools/cxx03_warn.hpp>
0024
0025 #include <type_traits>
0026
0027 namespace boost
0028 {
0029    namespace math
0030    {
0031
0032       namespace detail {
0033
0034          template <class T>
0035          struct is_trivial_arithmetic_type_imp
0036          {
0037             typedef std::integral_constant<bool,
0038                noexcept(std::declval<T&>() += std::declval<T>())
0039                && noexcept(std::declval<T&>() -= std::declval<T>())
0040                && noexcept(std::declval<T&>() *= std::declval<T>())
0041                && noexcept(std::declval<T&>() /= std::declval<T>())
0042             > type;
0043          };
0044
0045          template <class T>
0046          struct is_trivial_arithmetic_type : public is_trivial_arithmetic_type_imp<T>::type {};
0047       }
0048
0049 #ifndef BOOST_NO_CXX14_CONSTEXPR
0050       namespace constexpr_detail
0051       {
0052          template <class T>
0053          constexpr void swap(T& a, T& b)
0054          {
0055             T t(a);
0056             a = b;
0057             b = t;
0058          }
0059        }
0060 #endif
0061
0062        template<typename T>
0063         class quaternion
0064         {
0065         public:
0066
0067             typedef T value_type;
0068
0069
0070             // constructor for H seen as R^4
0071             // (also default constructor)
0072
0073             constexpr explicit            quaternion( T const & requested_a = T(),
0074                                             T const & requested_b = T(),
0075                                             T const & requested_c = T(),
0076                                             T const & requested_d = T())
0077             :   a(requested_a),
0078                 b(requested_b),
0079                 c(requested_c),
0080                 d(requested_d)
0081             {
0082                 // nothing to do!
0083             }
0084
0085
0086             // constructor for H seen as C^2
0087
0088             constexpr explicit            quaternion( ::std::complex<T> const & z0,
0089                                             ::std::complex<T> const & z1 = ::std::complex<T>())
0090             :   a(z0.real()),
0091                 b(z0.imag()),
0092                 c(z1.real()),
0093                 d(z1.imag())
0094             {
0095                 // nothing to do!
0096             }
0097
0098
0099             // UNtemplated copy constructor
0100             constexpr quaternion(quaternion const & a_recopier)
0101                : a(a_recopier.R_component_1()),
0102                b(a_recopier.R_component_2()),
0103                c(a_recopier.R_component_3()),
0104                d(a_recopier.R_component_4()) {}
0105
0106             constexpr quaternion(quaternion && a_recopier)
0107                : a(std::move(a_recopier.R_component_1())),
0108                b(std::move(a_recopier.R_component_2())),
0109                c(std::move(a_recopier.R_component_3())),
0110                d(std::move(a_recopier.R_component_4())) {}
0111
0112             // templated copy constructor
0113
0114             template<typename X>
0115             constexpr explicit            quaternion(quaternion<X> const & a_recopier)
0116             :   a(static_cast<T>(a_recopier.R_component_1())),
0117                 b(static_cast<T>(a_recopier.R_component_2())),
0118                 c(static_cast<T>(a_recopier.R_component_3())),
0119                 d(static_cast<T>(a_recopier.R_component_4()))
0120             {
0121                 // nothing to do!
0122             }
0123
0124
0125             // destructor
0126             // (this is taken care of by the compiler itself)
0127
0128
0129             // accessors
0130             //
0131             // Note:    Like complex number, quaternions do have a meaningful notion of "real part",
0132             //            but unlike them there is no meaningful notion of "imaginary part".
0133             //            Instead there is an "unreal part" which itself is a quaternion, and usually
0134             //            nothing simpler (as opposed to the complex number case).
0135             //            However, for practicality, there are accessors for the other components
0136             //            (these are necessary for the templated copy constructor, for instance).
0137
0138             constexpr T real() const
0139             {
0140                return(a);
0141             }
0142
0143             constexpr quaternion<T> unreal() const
0144             {
0145                return(quaternion<T>(static_cast<T>(0), b, c, d));
0146             }
0147
0148             constexpr T R_component_1() const
0149             {
0150                return(a);
0151             }
0152
0153             constexpr T R_component_2() const
0154             {
0155                return(b);
0156             }
0157
0158             constexpr T R_component_3() const
0159             {
0160                return(c);
0161             }
0162
0163             constexpr T R_component_4() const
0164             {
0165                return(d);
0166             }
0167
0168             constexpr ::std::complex<T> C_component_1() const
0169             {
0170                return(::std::complex<T>(a, b));
0171             }
0172
0173             constexpr ::std::complex<T> C_component_2() const
0174             {
0175                return(::std::complex<T>(c, d));
0176             }
0177
0178             BOOST_CXX14_CONSTEXPR void swap(quaternion& o)
0179             {
0180 #ifndef BOOST_NO_CXX14_CONSTEXPR
0181                using constexpr_detail::swap;
0182 #else
0183                using std::swap;
0184 #endif
0185                swap(a, o.a);
0186                swap(b, o.b);
0187                swap(c, o.c);
0188                swap(d, o.d);
0189             }
0190
0191             // assignment operators
0192
0193             template<typename X>
0194             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<X> const  & a_affecter)
0195             {
0196                a = static_cast<T>(a_affecter.R_component_1());
0197                b = static_cast<T>(a_affecter.R_component_2());
0198                c = static_cast<T>(a_affecter.R_component_3());
0199                d = static_cast<T>(a_affecter.R_component_4());
0200
0201                return(*this);
0202             }
0203
0204             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> const & a_affecter)
0205             {
0206                a = a_affecter.a;
0207                b = a_affecter.b;
0208                c = a_affecter.c;
0209                d = a_affecter.d;
0210
0211                return(*this);
0212             }
0213
0214             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> && a_affecter)
0215             {
0216                a = std::move(a_affecter.a);
0217                b = std::move(a_affecter.b);
0218                c = std::move(a_affecter.c);
0219                d = std::move(a_affecter.d);
0220
0221                return(*this);
0222             }
0223
0224             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (T const & a_affecter)
0225             {
0226                a = a_affecter;
0227
0228                b = c = d = static_cast<T>(0);
0229
0230                return(*this);
0231             }
0232
0233             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (::std::complex<T> const & a_affecter)
0234             {
0235                a = a_affecter.real();
0236                b = a_affecter.imag();
0237
0238                c = d = static_cast<T>(0);
0239
0240                return(*this);
0241             }
0242
0243             // other assignment-related operators
0244             //
0245             // NOTE:    Quaternion multiplication is *NOT* commutative;
0246             //            symbolically, "q *= rhs;" means "q = q * rhs;"
0247             //            and "q /= rhs;" means "q = q * inverse_of(rhs);"
0248             //
0249             // Note2:   Each operator comes in 2 forms - one for the simple case where
0250             //          type T throws no exceptions, and one exception-safe version
0251             //          for the case where it might.
0252          private:
0253             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const std::true_type&)
0254             {
0255                a += rhs;
0256                return *this;
0257             }
0258             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const std::false_type&)
0259             {
0260                quaternion<T> result(a + rhs, b, c, d); // exception guard
0261                swap(result);
0262                return *this;
0263             }
0264             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const std::true_type&)
0265             {
0266                a += std::real(rhs);
0267                b += std::imag(rhs);
0268                return *this;
0269             }
0270             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const std::false_type&)
0271             {
0272                quaternion<T> result(a + std::real(rhs), b + std::imag(rhs), c, d); // exception guard
0273                swap(result);
0274                return *this;
0275             }
0276             template <class X>
0277             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const std::true_type&)
0278             {
0279                a += rhs.R_component_1();
0280                b += rhs.R_component_2();
0281                c += rhs.R_component_3();
0282                d += rhs.R_component_4();
0283                return *this;
0284             }
0285             template <class X>
0286             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const std::false_type&)
0287             {
0288                quaternion<T> result(a + rhs.R_component_1(), b + rhs.R_component_2(), c + rhs.R_component_3(), d + rhs.R_component_4()); // exception guard
0289                swap(result);
0290                return *this;
0291             }
0292
0293             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const std::true_type&)
0294             {
0295                a -= rhs;
0296                return *this;
0297             }
0298             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const std::false_type&)
0299             {
0300                quaternion<T> result(a - rhs, b, c, d); // exception guard
0301                swap(result);
0302                return *this;
0303             }
0304             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const std::true_type&)
0305             {
0306                a -= std::real(rhs);
0307                b -= std::imag(rhs);
0308                return *this;
0309             }
0310             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const std::false_type&)
0311             {
0312                quaternion<T> result(a - std::real(rhs), b - std::imag(rhs), c, d); // exception guard
0313                swap(result);
0314                return *this;
0315             }
0316             template <class X>
0317             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const std::true_type&)
0318             {
0319                a -= rhs.R_component_1();
0320                b -= rhs.R_component_2();
0321                c -= rhs.R_component_3();
0322                d -= rhs.R_component_4();
0323                return *this;
0324             }
0325             template <class X>
0326             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const std::false_type&)
0327             {
0328                quaternion<T> result(a - rhs.R_component_1(), b - rhs.R_component_2(), c - rhs.R_component_3(), d - rhs.R_component_4()); // exception guard
0329                swap(result);
0330                return *this;
0331             }
0332
0333             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const std::true_type&)
0334             {
0335                a *= rhs;
0336                b *= rhs;
0337                c *= rhs;
0338                d *= rhs;
0339                return *this;
0340             }
0341             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const std::false_type&)
0342             {
0343                quaternion<T> result(a * rhs, b * rhs, c * rhs, d * rhs); // exception guard
0344                swap(result);
0345                return *this;
0346             }
0347
0348             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const std::true_type&)
0349             {
0350                a /= rhs;
0351                b /= rhs;
0352                c /= rhs;
0353                d /= rhs;
0354                return *this;
0355             }
0356             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const std::false_type&)
0357             {
0358                quaternion<T> result(a / rhs, b / rhs, c / rhs, d / rhs); // exception guard
0359                swap(result);
0360                return *this;
0361             }
0362          public:
0363
0364             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (T const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
0365             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (::std::complex<T> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
0366             template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator += (quaternion<X> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
0367
0368             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (T const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
0369             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (::std::complex<T> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
0370             template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator -= (quaternion<X> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
0371
0372             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (T const & rhs) { return do_multiply(rhs, detail::is_trivial_arithmetic_type<T>()); }
0373
0374             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (::std::complex<T> const & rhs)
0375             {
0376                 T    ar = rhs.real();
0377                 T    br = rhs.imag();
0378                 quaternion<T> result(a*ar - b*br, a*br + b*ar, c*ar + d*br, -c*br+d*ar);
0379                 swap(result);
0380                 return(*this);
0381             }
0382
0383             template<typename X>
0384             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (quaternion<X> const & rhs)
0385             {
0386                 T    ar = static_cast<T>(rhs.R_component_1());
0387                 T    br = static_cast<T>(rhs.R_component_2());
0388                 T    cr = static_cast<T>(rhs.R_component_3());
0389                 T    dr = static_cast<T>(rhs.R_component_4());
0390
0391                 quaternion<T> result(a*ar - b*br - c*cr - d*dr, a*br + b*ar + c*dr - d*cr, a*cr - b*dr + c*ar + d*br, a*dr + b*cr - c*br + d*ar);
0392                 swap(result);
0393                 return(*this);
0394             }
0395
0396             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (T const & rhs) { return do_divide(rhs, detail::is_trivial_arithmetic_type<T>()); }
0397
0398             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (::std::complex<T> const & rhs)
0399             {
0400                 T    ar = rhs.real();
0401                 T    br = rhs.imag();
0402                 T    denominator = ar*ar+br*br;
0403                 quaternion<T> result((+a*ar + b*br) / denominator, (-a*br + b*ar) / denominator, (+c*ar - d*br) / denominator, (+c*br + d*ar) / denominator);
0404                 swap(result);
0405                 return(*this);
0406             }
0407
0408             template<typename X>
0409             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (quaternion<X> const & rhs)
0410             {
0411                 T    ar = static_cast<T>(rhs.R_component_1());
0412                 T    br = static_cast<T>(rhs.R_component_2());
0413                 T    cr = static_cast<T>(rhs.R_component_3());
0414                 T    dr = static_cast<T>(rhs.R_component_4());
0415
0416                 T    denominator = ar*ar+br*br+cr*cr+dr*dr;
0417                 quaternion<T> result((+a*ar+b*br+c*cr+d*dr)/denominator, (-a*br+b*ar-c*dr+d*cr)/denominator, (-a*cr+b*dr+c*ar-d*br)/denominator, (-a*dr-b*cr+c*br+d*ar)/denominator);
0418                 swap(result);
0419                 return(*this);
0420             }
0421         private:
0422            T a, b, c, d;
0423
0424         };
0425
0426 // swap:
0427 template <class T>
0428 BOOST_CXX14_CONSTEXPR void swap(quaternion<T>& a, quaternion<T>& b) { a.swap(b); }
0429
0430 // operator+
0431 template <class T1, class T2>
0432 inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0433 operator + (const quaternion<T1>& a, const T2& b)
0434 {
0435    return quaternion<T1>(static_cast<T1>(a.R_component_1() + b), a.R_component_2(), a.R_component_3(), a.R_component_4());
0436 }
0437 template <class T1, class T2>
0438 inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0439 operator + (const T1& a, const quaternion<T2>& b)
0440 {
0441    return quaternion<T2>(static_cast<T2>(b.R_component_1() + a), b.R_component_2(), b.R_component_3(), b.R_component_4());
0442 }
0443 template <class T1, class T2>
0444 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0445 operator + (const quaternion<T1>& a, const std::complex<T2>& b)
0446 {
0447    return quaternion<T1>(a.R_component_1() + std::real(b), a.R_component_2() + std::imag(b), a.R_component_3(), a.R_component_4());
0448 }
0449 template <class T1, class T2>
0450 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0451 operator + (const std::complex<T1>& a, const quaternion<T2>& b)
0452 {
0453    return quaternion<T1>(b.R_component_1() + std::real(a), b.R_component_2() + std::imag(a), b.R_component_3(), b.R_component_4());
0454 }
0455 template <class T>
0456 inline constexpr quaternion<T> operator + (const quaternion<T>& a, const quaternion<T>& b)
0457 {
0458    return quaternion<T>(a.R_component_1() + b.R_component_1(), a.R_component_2() + b.R_component_2(), a.R_component_3() + b.R_component_3(), a.R_component_4() + b.R_component_4());
0459 }
0460 // operator-
0461 template <class T1, class T2>
0462 inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0463 operator - (const quaternion<T1>& a, const T2& b)
0464 {
0465    return quaternion<T1>(static_cast<T1>(a.R_component_1() - b), a.R_component_2(), a.R_component_3(), a.R_component_4());
0466 }
0467 template <class T1, class T2>
0468 inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0469 operator - (const T1& a, const quaternion<T2>& b)
0470 {
0471    return quaternion<T2>(static_cast<T2>(a - b.R_component_1()), -b.R_component_2(), -b.R_component_3(), -b.R_component_4());
0472 }
0473 template <class T1, class T2>
0474 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0475 operator - (const quaternion<T1>& a, const std::complex<T2>& b)
0476 {
0477    return quaternion<T1>(a.R_component_1() - std::real(b), a.R_component_2() - std::imag(b), a.R_component_3(), a.R_component_4());
0478 }
0479 template <class T1, class T2>
0480 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0481 operator - (const std::complex<T1>& a, const quaternion<T2>& b)
0482 {
0483    return quaternion<T1>(std::real(a) - b.R_component_1(), std::imag(a) - b.R_component_2(), -b.R_component_3(), -b.R_component_4());
0484 }
0485 template <class T>
0486 inline constexpr quaternion<T> operator - (const quaternion<T>& a, const quaternion<T>& b)
0487 {
0488    return quaternion<T>(a.R_component_1() - b.R_component_1(), a.R_component_2() - b.R_component_2(), a.R_component_3() - b.R_component_3(), a.R_component_4() - b.R_component_4());
0489 }
0490
0491 // operator*
0492 template <class T1, class T2>
0493 inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0494 operator * (const quaternion<T1>& a, const T2& b)
0495 {
0496    return quaternion<T1>(static_cast<T1>(a.R_component_1() * b), a.R_component_2() * b, a.R_component_3() * b, a.R_component_4() * b);
0497 }
0498 template <class T1, class T2>
0499 inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0500 operator * (const T1& a, const quaternion<T2>& b)
0501 {
0502    return quaternion<T2>(static_cast<T2>(a * b.R_component_1()), a * b.R_component_2(), a * b.R_component_3(), a * b.R_component_4());
0503 }
0504 template <class T1, class T2>
0505 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0506 operator * (const quaternion<T1>& a, const std::complex<T2>& b)
0507 {
0508    quaternion<T1> result(a);
0509    result *= b;
0510    return result;
0511 }
0512 template <class T1, class T2>
0513 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0514 operator * (const std::complex<T1>& a, const quaternion<T2>& b)
0515 {
0516    quaternion<T1> result(a);
0517    result *= b;
0518    return result;
0519 }
0520 template <class T>
0521 inline BOOST_CXX14_CONSTEXPR quaternion<T> operator * (const quaternion<T>& a, const quaternion<T>& b)
0522 {
0523    quaternion<T> result(a);
0524    result *= b;
0525    return result;
0526 }
0527
0528 // operator/
0529 template <class T1, class T2>
0530 inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0531 operator / (const quaternion<T1>& a, const T2& b)
0532 {
0533    return quaternion<T1>(a.R_component_1() / b, a.R_component_2() / b, a.R_component_3() / b, a.R_component_4() / b);
0534 }
0535 template <class T1, class T2>
0536 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0537 operator / (const T1& a, const quaternion<T2>& b)
0538 {
0539    quaternion<T2> result(a);
0540    result /= b;
0541    return result;
0542 }
0543 template <class T1, class T2>
0544 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type
0545 operator / (const quaternion<T1>& a, const std::complex<T2>& b)
0546 {
0547    quaternion<T1> result(a);
0548    result /= b;
0549    return result;
0550 }
0551 template <class T1, class T2>
0552 inline BOOST_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type
0553 operator / (const std::complex<T1>& a, const quaternion<T2>& b)
0554 {
0555    quaternion<T2> result(a);
0556    result /= b;
0557    return result;
0558 }
0559 template <class T>
0560 inline BOOST_CXX14_CONSTEXPR quaternion<T> operator / (const quaternion<T>& a, const quaternion<T>& b)
0561 {
0562    quaternion<T> result(a);
0563    result /= b;
0564    return result;
0565 }
0566
0567
0568         template<typename T>
0569         inline constexpr const quaternion<T>&             operator + (quaternion<T> const & q)
0570         {
0571             return q;
0572         }
0573
0574
0575         template<typename T>
0576         inline constexpr quaternion<T>                    operator - (quaternion<T> const & q)
0577         {
0578             return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4()));
0579         }
0580
0581
0582         template<typename R, typename T>
0583         inline constexpr typename std::enable_if<std::is_convertible<R, T>::value, bool>::type operator == (R const & lhs, quaternion<T> const & rhs)
0584         {
0585             return    (
0586                         (rhs.R_component_1() == lhs)&&
0587                         (rhs.R_component_2() == static_cast<T>(0))&&
0588                         (rhs.R_component_3() == static_cast<T>(0))&&
0589                         (rhs.R_component_4() == static_cast<T>(0))
0590                     );
0591         }
0592
0593
0594         template<typename T, typename R>
0595         inline constexpr typename std::enable_if<std::is_convertible<R, T>::value, bool>::type operator == (quaternion<T> const & lhs, R const & rhs)
0596         {
0597            return rhs == lhs;
0598         }
0599
0600
0601         template<typename T>
0602         inline constexpr bool                                operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs)
0603         {
0604             return    (
0605                         (rhs.R_component_1() == lhs.real())&&
0606                         (rhs.R_component_2() == lhs.imag())&&
0607                         (rhs.R_component_3() == static_cast<T>(0))&&
0608                         (rhs.R_component_4() == static_cast<T>(0))
0609                     );
0610         }
0611
0612
0613         template<typename T>
0614         inline constexpr bool                                operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs)
0615         {
0616            return rhs == lhs;
0617         }
0618
0619
0620         template<typename T>
0621         inline constexpr bool                                operator == (quaternion<T> const & lhs, quaternion<T> const & rhs)
0622         {
0623             return    (
0624                         (rhs.R_component_1() == lhs.R_component_1())&&
0625                         (rhs.R_component_2() == lhs.R_component_2())&&
0626                         (rhs.R_component_3() == lhs.R_component_3())&&
0627                         (rhs.R_component_4() == lhs.R_component_4())
0628                     );
0629         }
0630
0631         template<typename R, typename T> inline constexpr bool operator != (R const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
0632         template<typename T, typename R> inline constexpr bool operator != (quaternion<T> const & lhs, R const & rhs) { return !(lhs == rhs); }
0633         template<typename T> inline constexpr bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
0634         template<typename T> inline constexpr bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs) { return !(lhs == rhs); }
0635         template<typename T> inline constexpr bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
0636
0637
0638         // Note:    we allow the following formats, with a, b, c, and d reals
0639         //            a
0640         //            (a), (a,b), (a,b,c), (a,b,c,d)
0641         //            (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))
0642         template<typename T, typename charT, class traits>
0643         ::std::basic_istream<charT,traits> &    operator >> (    ::std::basic_istream<charT,traits> & is,
0644                                                                 quaternion<T> & q)
0645         {
0646             const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc());
0647
0648             T    a = T();
0649             T    b = T();
0650             T    c = T();
0651             T    d = T();
0652
0653             ::std::complex<T>    u = ::std::complex<T>();
0654             ::std::complex<T>    v = ::std::complex<T>();
0655
0656             charT    ch = charT();
0657             char    cc;
0658
0659             is >> ch;                                        // get the first lexeme
0660
0661             if    (!is.good())    goto finish;
0662
0663             cc = ct.narrow(ch, char());
0664
0665             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,))
0666             {
0667                 is >> ch;                                    // get the second lexeme
0668
0669                 if    (!is.good())    goto finish;
0670
0671                 cc = ct.narrow(ch, char());
0672
0673                 if    (cc == '(')                        // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))
0674                 {
0675                     is.putback(ch);
0676
0677                     is >> u;                                // we extract the first and second components
0678                     a = u.real();
0679                     b = u.imag();
0680
0681                     if    (!is.good())    goto finish;
0682
0683                     is >> ch;                                // get the next lexeme
0684
0685                     if    (!is.good())    goto finish;
0686
0687                     cc = ct.narrow(ch, char());
0688
0689                     if        (cc == ')')                    // format: ((a)) or ((a,b))
0690                     {
0691                         q = quaternion<T>(a,b);
0692                     }
0693                     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,))
0694                     {
0695                         is >> v;                            // we extract the third and fourth components
0696                         c = v.real();
0697                         d = v.imag();
0698
0699                         if    (!is.good())    goto finish;
0700
0701                         is >> ch;                                // get the last lexeme
0702
0703                         if    (!is.good())    goto finish;
0704
0705                         cc = ct.narrow(ch, char());
0706
0707                         if    (cc == ')')                    // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,))
0708                         {
0709                             q = quaternion<T>(a,b,c,d);
0710                         }
0711                         else                            // error
0712                         {
0713                             is.setstate(::std::ios_base::failbit);
0714                         }
0715                     }
0716                     else                                // error
0717                     {
0718                         is.setstate(::std::ios_base::failbit);
0719                     }
0720                 }
0721                 else                                // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
0722                 {
0723                     is.putback(ch);
0724
0725                     is >> a;                                // we extract the first component
0726
0727                     if    (!is.good())    goto finish;
0728
0729                     is >> ch;                                // get the third lexeme
0730
0731                     if    (!is.good())    goto finish;
0732
0733                     cc = ct.narrow(ch, char());
0734
0735                     if        (cc == ')')                    // format: (a)
0736                     {
0737                         q = quaternion<T>(a);
0738                     }
0739                     else if    (cc == ',')                // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
0740                     {
0741                         is >> ch;                            // get the fourth lexeme
0742
0743                         if    (!is.good())    goto finish;
0744
0745                         cc = ct.narrow(ch, char());
0746
0747                         if    (cc == '(')                // read "(a,(", possible: (a,(c)), (a,(c,d))
0748                         {
0749                             is.putback(ch);
0750
0751                             is >> v;                        // we extract the third and fourth component
0752
0753                             c = v.real();
0754                             d = v.imag();
0755
0756                             if    (!is.good())    goto finish;
0757
0758                             is >> ch;                        // get the ninth lexeme
0759
0760                             if    (!is.good())    goto finish;
0761
0762                             cc = ct.narrow(ch, char());
0763
0764                             if    (cc == ')')                // format: (a,(c)) or (a,(c,d))
0765                             {
0766                                 q = quaternion<T>(a,b,c,d);
0767                             }
0768                             else                        // error
0769                             {
0770                                 is.setstate(::std::ios_base::failbit);
0771                             }
0772                         }
0773                         else                        // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d)
0774                         {
0775                             is.putback(ch);
0776
0777                             is >> b;                        // we extract the second component
0778
0779                             if    (!is.good())    goto finish;
0780
0781                             is >> ch;                        // get the fifth lexeme
0782
0783                             if    (!is.good())    goto finish;
0784
0785                             cc = ct.narrow(ch, char());
0786
0787                             if    (cc == ')')                // format: (a,b)
0788                             {
0789                                 q = quaternion<T>(a,b);
0790                             }
0791                             else if    (cc == ',')        // read "(a,b,", possible: (a,b,c), (a,b,c,d)
0792                             {
0793                                 is >> c;                    // we extract the third component
0794
0795                                 if    (!is.good())    goto finish;
0796
0797                                 is >> ch;                    // get the seventh lexeme
0798
0799                                 if    (!is.good())    goto finish;
0800
0801                                 cc = ct.narrow(ch, char());
0802
0803                                 if        (cc == ')')        // format: (a,b,c)
0804                                 {
0805                                     q = quaternion<T>(a,b,c);
0806                                 }
0807                                 else if    (cc == ',')    // read "(a,b,c,", possible: (a,b,c,d)
0808                                 {
0809                                     is >> d;                // we extract the fourth component
0810
0811                                     if    (!is.good())    goto finish;
0812
0813                                     is >> ch;                // get the ninth lexeme
0814
0815                                     if    (!is.good())    goto finish;
0816
0817                                     cc = ct.narrow(ch, char());
0818
0819                                     if    (cc == ')')        // format: (a,b,c,d)
0820                                     {
0821                                         q = quaternion<T>(a,b,c,d);
0822                                     }
0823                                     else                // error
0824                                     {
0825                                         is.setstate(::std::ios_base::failbit);
0826                                     }
0827                                 }
0828                                 else                    // error
0829                                 {
0830                                     is.setstate(::std::ios_base::failbit);
0831                                 }
0832                             }
0833                             else                        // error
0834                             {
0835                                 is.setstate(::std::ios_base::failbit);
0836                             }
0837                         }
0838                     }
0839                     else                                // error
0840                     {
0841                         is.setstate(::std::ios_base::failbit);
0842                     }
0843                 }
0844             }
0845             else                                        // format:    a
0846             {
0847                 is.putback(ch);
0848
0849                 is >> a;                                    // we extract the first component
0850
0851                 if    (!is.good())    goto finish;
0852
0853                 q = quaternion<T>(a);
0854             }
0855
0856             finish:
0857             return(is);
0858         }
0859
0860
0861         template<typename T, typename charT, class traits>
0862         ::std::basic_ostream<charT,traits> &    operator << (    ::std::basic_ostream<charT,traits> & os,
0863                                                                 quaternion<T> const & q)
0864         {
0865             ::std::basic_ostringstream<charT,traits>    s;
0866
0867             s.flags(os.flags());
0868             s.imbue(os.getloc());
0869             s.precision(os.precision());
0870
0871             s << '('    << q.R_component_1() << ','
0872                         << q.R_component_2() << ','
0873                         << q.R_component_3() << ','
0874                         << q.R_component_4() << ')';
0875
0876             return os << s.str();
0877         }
0878
0879
0880         // values
0881
0882         template<typename T>
0883         inline constexpr T real(quaternion<T> const & q)
0884         {
0885             return(q.real());
0886         }
0887
0888
0889         template<typename T>
0890         inline constexpr quaternion<T> unreal(quaternion<T> const & q)
0891         {
0892             return(q.unreal());
0893         }
0894
0895         template<typename T>
0896         inline T sup(quaternion<T> const & q)
0897         {
0898             using    ::std::abs;
0899             return (std::max)((std::max)(abs(q.R_component_1()), abs(q.R_component_2())), (std::max)(abs(q.R_component_3()), abs(q.R_component_4())));
0900         }
0901
0902
0903         template<typename T>
0904         inline T l1(quaternion<T> const & q)
0905         {
0906            using    ::std::abs;
0907            return abs(q.R_component_1()) + abs(q.R_component_2()) + abs(q.R_component_3()) + abs(q.R_component_4());
0908         }
0909
0910
0911         template<typename T>
0912         inline T abs(quaternion<T> const & q)
0913         {
0914             using    ::std::abs;
0915             using    ::std::sqrt;
0916
0917             T            maxim = sup(q);    // overflow protection
0918
0919             if    (maxim == static_cast<T>(0))
0920             {
0921                 return(maxim);
0922             }
0923             else
0924             {
0925                 T    mixam = static_cast<T>(1)/maxim;    // prefer multiplications over divisions
0926
0927                 T a = q.R_component_1() * mixam;
0928                 T b = q.R_component_2() * mixam;
0929                 T c = q.R_component_3() * mixam;
0930                 T d = q.R_component_4() * mixam;
0931
0932                 a *= a;
0933                 b *= b;
0934                 c *= c;
0935                 d *= d;
0936
0937                 return(maxim * sqrt(a + b + c + d));
0938             }
0939
0940             //return(sqrt(norm(q)));
0941         }
0942
0943
0944         // Note:    This is the Cayley norm, not the Euclidean norm...
0945
0946         template<typename T>
0947         inline BOOST_CXX14_CONSTEXPR T norm(quaternion<T>const  & q)
0948         {
0949             return(real(q*conj(q)));
0950         }
0951
0952
0953         template<typename T>
0954         inline constexpr quaternion<T> conj(quaternion<T> const & q)
0955         {
0956             return(quaternion<T>(   +q.R_component_1(),
0957                                     -q.R_component_2(),
0958                                     -q.R_component_3(),
0959                                     -q.R_component_4()));
0960         }
0961
0962
0963         template<typename T>
0964         inline quaternion<T>                    spherical(  T const & rho,
0965                                                             T const & theta,
0966                                                             T const & phi1,
0967                                                             T const & phi2)
0968         {
0969             using ::std::cos;
0970             using ::std::sin;
0971
0972             //T    a = cos(theta)*cos(phi1)*cos(phi2);
0973             //T    b = sin(theta)*cos(phi1)*cos(phi2);
0974             //T    c = sin(phi1)*cos(phi2);
0975             //T    d = sin(phi2);
0976
0977             T    courrant = static_cast<T>(1);
0978
0979             T    d = sin(phi2);
0980
0981             courrant *= cos(phi2);
0982
0983             T    c = sin(phi1)*courrant;
0984
0985             courrant *= cos(phi1);
0986
0987             T    b = sin(theta)*courrant;
0988             T    a = cos(theta)*courrant;
0989
0990             return(rho*quaternion<T>(a,b,c,d));
0991         }
0992
0993
0994         template<typename T>
0995         inline quaternion<T>                    semipolar(  T const & rho,
0996                                                             T const & alpha,
0997                                                             T const & theta1,
0998                                                             T const & theta2)
0999         {
1000             using ::std::cos;
1001             using ::std::sin;
1002
1003             T    a = cos(alpha)*cos(theta1);
1004             T    b = cos(alpha)*sin(theta1);
1005             T    c = sin(alpha)*cos(theta2);
1006             T    d = sin(alpha)*sin(theta2);
1007
1008             return(rho*quaternion<T>(a,b,c,d));
1009         }
1010
1011
1012         template<typename T>
1013         inline quaternion<T>                    multipolar( T const & rho1,
1014                                                             T const & theta1,
1015                                                             T const & rho2,
1016                                                             T const & theta2)
1017         {
1018             using ::std::cos;
1019             using ::std::sin;
1020
1021             T    a = rho1*cos(theta1);
1022             T    b = rho1*sin(theta1);
1023             T    c = rho2*cos(theta2);
1024             T    d = rho2*sin(theta2);
1025
1026             return(quaternion<T>(a,b,c,d));
1027         }
1028
1029
1030         template<typename T>
1031         inline quaternion<T>                    cylindrospherical(  T const & t,
1032                                                                     T const & radius,
1033                                                                     T const & longitude,
1034                                                                     T const & latitude)
1035         {
1036             using ::std::cos;
1037             using ::std::sin;
1038
1039
1040
1041             T    b = radius*cos(longitude)*cos(latitude);
1042             T    c = radius*sin(longitude)*cos(latitude);
1043             T    d = radius*sin(latitude);
1044
1045             return(quaternion<T>(t,b,c,d));
1046         }
1047
1048
1049         template<typename T>
1050         inline quaternion<T>                    cylindrical(T const & r,
1051                                                             T const & angle,
1052                                                             T const & h1,
1053                                                             T const & h2)
1054         {
1055             using ::std::cos;
1056             using ::std::sin;
1057
1058             T    a = r*cos(angle);
1059             T    b = r*sin(angle);
1060
1061             return(quaternion<T>(a,b,h1,h2));
1062         }
1063
1064
1065         // transcendentals
1066         // (please see the documentation)
1067
1068
1069         template<typename T>
1070         inline quaternion<T>                    exp(quaternion<T> const & q)
1071         {
1072             using    ::std::exp;
1073             using    ::std::cos;
1074
1075             using    ::boost::math::sinc_pi;
1076
1077             T    u = exp(real(q));
1078
1079             T    z = abs(unreal(q));
1080
1081             T    w = sinc_pi(z);
1082
1083             return(u*quaternion<T>(cos(z),
1084                 w*q.R_component_2(), w*q.R_component_3(),
1085                 w*q.R_component_4()));
1086         }
1087
1088
1089         template<typename T>
1090         inline quaternion<T>                    cos(quaternion<T> const & q)
1091         {
1092             using    ::std::sin;
1093             using    ::std::cos;
1094             using    ::std::cosh;
1095
1096             using    ::boost::math::sinhc_pi;
1097
1098             T    z = abs(unreal(q));
1099
1100             T    w = -sin(q.real())*sinhc_pi(z);
1101
1102             return(quaternion<T>(cos(q.real())*cosh(z),
1103                 w*q.R_component_2(), w*q.R_component_3(),
1104                 w*q.R_component_4()));
1105         }
1106
1107
1108         template<typename T>
1109         inline quaternion<T>                    sin(quaternion<T> const & q)
1110         {
1111             using    ::std::sin;
1112             using    ::std::cos;
1113             using    ::std::cosh;
1114
1115             using    ::boost::math::sinhc_pi;
1116
1117             T    z = abs(unreal(q));
1118
1119             T    w = +cos(q.real())*sinhc_pi(z);
1120
1121             return(quaternion<T>(sin(q.real())*cosh(z),
1122                 w*q.R_component_2(), w*q.R_component_3(),
1123                 w*q.R_component_4()));
1124         }
1125
1126
1127         template<typename T>
1128         inline quaternion<T>                    tan(quaternion<T> const & q)
1129         {
1130             return(sin(q)/cos(q));
1131         }
1132
1133
1134         template<typename T>
1135         inline quaternion<T>                    cosh(quaternion<T> const & q)
1136         {
1137             return((exp(+q)+exp(-q))/static_cast<T>(2));
1138         }
1139
1140
1141         template<typename T>
1142         inline quaternion<T>                    sinh(quaternion<T> const & q)
1143         {
1144             return((exp(+q)-exp(-q))/static_cast<T>(2));
1145         }
1146
1147
1148         template<typename T>
1149         inline quaternion<T>                    tanh(quaternion<T> const & q)
1150         {
1151             return(sinh(q)/cosh(q));
1152         }
1153
1154
1155         template<typename T>
1156         quaternion<T>                            pow(quaternion<T> const & q,
1157                                                     int n)
1158         {
1159             if        (n > 1)
1160             {
1161                 int    m = n>>1;
1162
1163                 quaternion<T>    result = pow(q, m);
1164
1165                 result *= result;
1166
1167                 if    (n != (m<<1))
1168                 {
1169                     result *= q; // n odd
1170                 }
1171
1172                 return(result);
1173             }
1174             else if    (n == 1)
1175             {
1176                 return(q);
1177             }
1178             else if    (n == 0)
1179             {
1180                 return(quaternion<T>(static_cast<T>(1)));
1181             }
1182             else    /* n < 0 */
1183             {
1184                 return(pow(quaternion<T>(static_cast<T>(1))/q,-n));
1185             }
1186         }
1187     }
1188 }
1189
1190 #endif /* BOOST_QUATERNION_HPP */