ceph/src/boost/boost/math/quaternion.hpp

   1 //  boost quaternion.hpp header file
   2
   3 //  (C) Copyright Hubert Holin 2001.
   4 //  Distributed under the Boost Software License, Version 1.0. (See
   5 //  accompanying file LICENSE_1_0.txt or copy at
   6 //  http://www.boost.org/LICENSE_1_0.txt)
   7
   8 // See http://www.boost.org for updates, documentation, and revision history.
   9
  10 #ifndef BOOST_QUATERNION_HPP
  11 #define BOOST_QUATERNION_HPP
  12
  13 #include <boost/config.hpp> // for BOOST_NO_STD_LOCALE
  14 #include <boost/math_fwd.hpp>
  15 #include <boost/detail/workaround.hpp>
  16 #include <boost/type_traits/is_convertible.hpp>
  17 #include <boost/utility/enable_if.hpp>
  18 #ifndef    BOOST_NO_STD_LOCALE
  19 #  include <locale>                                    // for the "<<" operator
  20 #endif /* BOOST_NO_STD_LOCALE */
  21
  22 #include <complex>
  23 #include <iosfwd>                                    // for the "<<" and ">>" operators
  24 #include <sstream>                                    // for the "<<" operator
  25
  26 #include <boost/math/special_functions/sinc.hpp>    // for the Sinus cardinal
  27 #include <boost/math/special_functions/sinhc.hpp>    // for the Hyperbolic Sinus cardinal
  28
  29 #if defined(BOOST_NO_CXX11_NOEXCEPT) || defined(BOOST_NO_CXX11_RVALUE_REFERENCES) || defined(BOOST_NO_SFINAE_EXPR)
  30 #include <boost/type_traits/is_pod.hpp>
  31 #endif
  32
  33 namespace boost
  34 {
  35    namespace math
  36    {
  37
  38       namespace detail {
  39
  40 #if !defined(BOOST_NO_CXX11_NOEXCEPT) && !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && !defined(BOOST_NO_SFINAE_EXPR)
  41
  42          template <class T>
  43          struct is_trivial_arithmetic_type_imp
  44          {
  45             typedef mpl::bool_<
  46                noexcept(std::declval<T&>() += std::declval<T>())
  47                && noexcept(std::declval<T&>() -= std::declval<T>())
  48                && noexcept(std::declval<T&>() *= std::declval<T>())
  49                && noexcept(std::declval<T&>() /= std::declval<T>())
  50             > type;
  51          };
  52
  53          template <class T>
  54          struct is_trivial_arithmetic_type : public is_trivial_arithmetic_type_imp<T>::type {};
  55 #else
  56
  57          template <class T>
  58          struct is_trivial_arithmetic_type : public boost::is_pod<T> {};
  59
  60 #endif
  61
  62       }
  63
  64 #ifndef BOOST_NO_CXX14_CONSTEXPR
  65       namespace constexpr_detail
  66       {
  67          template <class T>
  68          constexpr void swap(T& a, T& b)
  69          {
  70             T t(a);
  71             a = b;
  72             b = t;
  73          }
  74        }
  75 #endif
  76
  77        template<typename T>
  78         class quaternion
  79         {
  80         public:
  81
  82             typedef T value_type;
  83
  84
  85             // constructor for H seen as R^4
  86             // (also default constructor)
  87
  88             BOOST_CONSTEXPR explicit            quaternion( T const & requested_a = T(),
  89                                             T const & requested_b = T(),
  90                                             T const & requested_c = T(),
  91                                             T const & requested_d = T())
  92             :   a(requested_a),
  93                 b(requested_b),
  94                 c(requested_c),
  95                 d(requested_d)
  96             {
  97                 // nothing to do!
  98             }
  99
 100
 101             // constructor for H seen as C^2
 102
 103             BOOST_CONSTEXPR explicit            quaternion( ::std::complex<T> const & z0,
 104                                             ::std::complex<T> const & z1 = ::std::complex<T>())
 105             :   a(z0.real()),
 106                 b(z0.imag()),
 107                 c(z1.real()),
 108                 d(z1.imag())
 109             {
 110                 // nothing to do!
 111             }
 112
 113
 114             // UNtemplated copy constructor
 115             BOOST_CONSTEXPR quaternion(quaternion const & a_recopier)
 116                : a(a_recopier.R_component_1()),
 117                b(a_recopier.R_component_2()),
 118                c(a_recopier.R_component_3()),
 119                d(a_recopier.R_component_4()) {}
 120 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 121             BOOST_CONSTEXPR quaternion(quaternion && a_recopier)
 122                : a(std::move(a_recopier.R_component_1())),
 123                b(std::move(a_recopier.R_component_2())),
 124                c(std::move(a_recopier.R_component_3())),
 125                d(std::move(a_recopier.R_component_4())) {}
 126 #endif
 127
 128             // templated copy constructor
 129
 130             template<typename X>
 131             BOOST_CONSTEXPR explicit            quaternion(quaternion<X> const & a_recopier)
 132             :   a(static_cast<T>(a_recopier.R_component_1())),
 133                 b(static_cast<T>(a_recopier.R_component_2())),
 134                 c(static_cast<T>(a_recopier.R_component_3())),
 135                 d(static_cast<T>(a_recopier.R_component_4()))
 136             {
 137                 // nothing to do!
 138             }
 139
 140
 141             // destructor
 142             // (this is taken care of by the compiler itself)
 143
 144
 145             // accessors
 146             //
 147             // Note:    Like complex number, quaternions do have a meaningful notion of "real part",
 148             //            but unlike them there is no meaningful notion of "imaginary part".
 149             //            Instead there is an "unreal part" which itself is a quaternion, and usually
 150             //            nothing simpler (as opposed to the complex number case).
 151             //            However, for practicallity, there are accessors for the other components
 152             //            (these are necessary for the templated copy constructor, for instance).
 153
 154             BOOST_CONSTEXPR T real() const
 155             {
 156                return(a);
 157             }
 158
 159             BOOST_CONSTEXPR quaternion<T> unreal() const
 160             {
 161                return(quaternion<T>(static_cast<T>(0), b, c, d));
 162             }
 163
 164             BOOST_CONSTEXPR T R_component_1() const
 165             {
 166                return(a);
 167             }
 168
 169             BOOST_CONSTEXPR T R_component_2() const
 170             {
 171                return(b);
 172             }
 173
 174             BOOST_CONSTEXPR T R_component_3() const
 175             {
 176                return(c);
 177             }
 178
 179             BOOST_CONSTEXPR T R_component_4() const
 180             {
 181                return(d);
 182             }
 183
 184             BOOST_CONSTEXPR ::std::complex<T> C_component_1() const
 185             {
 186                return(::std::complex<T>(a, b));
 187             }
 188
 189             BOOST_CONSTEXPR ::std::complex<T> C_component_2() const
 190             {
 191                return(::std::complex<T>(c, d));
 192             }
 193
 194             BOOST_CXX14_CONSTEXPR void swap(quaternion& o)
 195             {
 196 #ifndef BOOST_NO_CXX14_CONSTEXPR
 197                using constexpr_detail::swap;
 198 #else
 199                using std::swap;
 200 #endif
 201                swap(a, o.a);
 202                swap(b, o.b);
 203                swap(c, o.c);
 204                swap(d, o.d);
 205             }
 206
 207             // assignment operators
 208
 209             template<typename X>
 210             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<X> const  & a_affecter)
 211             {
 212                a = static_cast<T>(a_affecter.R_component_1());
 213                b = static_cast<T>(a_affecter.R_component_2());
 214                c = static_cast<T>(a_affecter.R_component_3());
 215                d = static_cast<T>(a_affecter.R_component_4());
 216
 217                return(*this);
 218             }
 219
 220             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> const & a_affecter)
 221             {
 222                a = a_affecter.a;
 223                b = a_affecter.b;
 224                c = a_affecter.c;
 225                d = a_affecter.d;
 226
 227                return(*this);
 228             }
 229 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 230             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> && a_affecter)
 231             {
 232                a = std::move(a_affecter.a);
 233                b = std::move(a_affecter.b);
 234                c = std::move(a_affecter.c);
 235                d = std::move(a_affecter.d);
 236
 237                return(*this);
 238             }
 239 #endif
 240             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (T const & a_affecter)
 241             {
 242                a = a_affecter;
 243
 244                b = c = d = static_cast<T>(0);
 245
 246                return(*this);
 247             }
 248
 249             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (::std::complex<T> const & a_affecter)
 250             {
 251                a = a_affecter.real();
 252                b = a_affecter.imag();
 253
 254                c = d = static_cast<T>(0);
 255
 256                return(*this);
 257             }
 258
 259             // other assignment-related operators
 260             //
 261             // NOTE:    Quaternion multiplication is *NOT* commutative;
 262             //            symbolically, "q *= rhs;" means "q = q * rhs;"
 263             //            and "q /= rhs;" means "q = q * inverse_of(rhs);"
 264             //
 265             // Note2:   Each operator comes in 2 forms - one for the simple case where
 266             //          type T throws no exceptions, and one exception-safe version
 267             //          for the case where it might.
 268          private:
 269             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const mpl::true_&)
 270             {
 271                a += rhs;
 272                return *this;
 273             }
 274             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const mpl::false_&)
 275             {
 276                quaternion<T> result(a + rhs, b, c, d); // exception guard
 277                swap(result);
 278                return *this;
 279             }
 280             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const mpl::true_&)
 281             {
 282                a += std::real(rhs);
 283                b += std::imag(rhs);
 284                return *this;
 285             }
 286             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const mpl::false_&)
 287             {
 288                quaternion<T> result(a + std::real(rhs), b + std::imag(rhs), c, d); // exception guard
 289                swap(result);
 290                return *this;
 291             }
 292             template <class X>
 293             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const mpl::true_&)
 294             {
 295                a += rhs.R_component_1();
 296                b += rhs.R_component_2();
 297                c += rhs.R_component_3();
 298                d += rhs.R_component_4();
 299                return *this;
 300             }
 301             template <class X>
 302             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const mpl::false_&)
 303             {
 304                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
 305                swap(result);
 306                return *this;
 307             }
 308
 309             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const mpl::true_&)
 310             {
 311                a -= rhs;
 312                return *this;
 313             }
 314             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const mpl::false_&)
 315             {
 316                quaternion<T> result(a - rhs, b, c, d); // exception guard
 317                swap(result);
 318                return *this;
 319             }
 320             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const mpl::true_&)
 321             {
 322                a -= std::real(rhs);
 323                b -= std::imag(rhs);
 324                return *this;
 325             }
 326             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const mpl::false_&)
 327             {
 328                quaternion<T> result(a - std::real(rhs), b - std::imag(rhs), c, d); // exception guard
 329                swap(result);
 330                return *this;
 331             }
 332             template <class X>
 333             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const mpl::true_&)
 334             {
 335                a -= rhs.R_component_1();
 336                b -= rhs.R_component_2();
 337                c -= rhs.R_component_3();
 338                d -= rhs.R_component_4();
 339                return *this;
 340             }
 341             template <class X>
 342             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const mpl::false_&)
 343             {
 344                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
 345                swap(result);
 346                return *this;
 347             }
 348
 349             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const mpl::true_&)
 350             {
 351                a *= rhs;
 352                b *= rhs;
 353                c *= rhs;
 354                d *= rhs;
 355                return *this;
 356             }
 357             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const mpl::false_&)
 358             {
 359                quaternion<T> result(a * rhs, b * rhs, c * rhs, d * rhs); // exception guard
 360                swap(result);
 361                return *this;
 362             }
 363
 364             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const mpl::true_&)
 365             {
 366                a /= rhs;
 367                b /= rhs;
 368                c /= rhs;
 369                d /= rhs;
 370                return *this;
 371             }
 372             BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const mpl::false_&)
 373             {
 374                quaternion<T> result(a / rhs, b / rhs, c / rhs, d / rhs); // exception guard
 375                swap(result);
 376                return *this;
 377             }
 378          public:
 379
 380             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (T const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
 381             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (::std::complex<T> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
 382             template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator += (quaternion<X> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }
 383
 384             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (T const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
 385             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (::std::complex<T> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
 386             template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator -= (quaternion<X> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }
 387
 388             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (T const & rhs) { return do_multiply(rhs, detail::is_trivial_arithmetic_type<T>()); }
 389
 390             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (::std::complex<T> const & rhs)
 391             {
 392                 T    ar = rhs.real();
 393                 T    br = rhs.imag();
 394                 quaternion<T> result(a*ar - b*br, a*br + b*ar, c*ar + d*br, -c*br+d*ar);
 395                 swap(result);
 396                 return(*this);
 397             }
 398
 399             template<typename X>
 400             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (quaternion<X> const & rhs)
 401             {
 402                 T    ar = static_cast<T>(rhs.R_component_1());
 403                 T    br = static_cast<T>(rhs.R_component_2());
 404                 T    cr = static_cast<T>(rhs.R_component_3());
 405                 T    dr = static_cast<T>(rhs.R_component_4());
 406
 407                 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);
 408                 swap(result);
 409                 return(*this);
 410             }
 411
 412             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (T const & rhs) { return do_divide(rhs, detail::is_trivial_arithmetic_type<T>()); }
 413
 414             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (::std::complex<T> const & rhs)
 415             {
 416                 T    ar = rhs.real();
 417                 T    br = rhs.imag();
 418                 T    denominator = ar*ar+br*br;
 419                 quaternion<T> result((+a*ar + b*br) / denominator, (-a*br + b*ar) / denominator, (+c*ar - d*br) / denominator, (+c*br + d*ar) / denominator);
 420                 swap(result);
 421                 return(*this);
 422             }
 423
 424             template<typename X>
 425             BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (quaternion<X> const & rhs)
 426             {
 427                 T    ar = static_cast<T>(rhs.R_component_1());
 428                 T    br = static_cast<T>(rhs.R_component_2());
 429                 T    cr = static_cast<T>(rhs.R_component_3());
 430                 T    dr = static_cast<T>(rhs.R_component_4());
 431
 432                 T    denominator = ar*ar+br*br+cr*cr+dr*dr;
 433                 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);
 434                 swap(result);
 435                 return(*this);
 436             }
 437         private:
 438            T a, b, c, d;
 439
 440         };
 441
 442 // swap:
 443 template <class T>
 444 BOOST_CXX14_CONSTEXPR void swap(quaternion<T>& a, quaternion<T>& b) { a.swap(b); }
 445
 446 // operator+
 447 template <class T1, class T2>
 448 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 449 operator + (const quaternion<T1>& a, const T2& b)
 450 {
 451    return quaternion<T1>(static_cast<T1>(a.R_component_1() + b), a.R_component_2(), a.R_component_3(), a.R_component_4());
 452 }
 453 template <class T1, class T2>
 454 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 455 operator + (const T1& a, const quaternion<T2>& b)
 456 {
 457    return quaternion<T2>(static_cast<T2>(b.R_component_1() + a), b.R_component_2(), b.R_component_3(), b.R_component_4());
 458 }
 459 template <class T1, class T2>
 460 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 461 operator + (const quaternion<T1>& a, const std::complex<T2>& b)
 462 {
 463    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());
 464 }
 465 template <class T1, class T2>
 466 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 467 operator + (const std::complex<T1>& a, const quaternion<T2>& b)
 468 {
 469    return quaternion<T1>(b.R_component_1() + real(a), b.R_component_2() + imag(a), b.R_component_3(), b.R_component_4());
 470 }
 471 template <class T>
 472 inline BOOST_CONSTEXPR quaternion<T> operator + (const quaternion<T>& a, const quaternion<T>& b)
 473 {
 474    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());
 475 }
 476 // operator-
 477 template <class T1, class T2>
 478 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 479 operator - (const quaternion<T1>& a, const T2& b)
 480 {
 481    return quaternion<T1>(static_cast<T1>(a.R_component_1() - b), a.R_component_2(), a.R_component_3(), a.R_component_4());
 482 }
 483 template <class T1, class T2>
 484 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 485 operator - (const T1& a, const quaternion<T2>& b)
 486 {
 487    return quaternion<T2>(static_cast<T2>(a - b.R_component_1()), -b.R_component_2(), -b.R_component_3(), -b.R_component_4());
 488 }
 489 template <class T1, class T2>
 490 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 491 operator - (const quaternion<T1>& a, const std::complex<T2>& b)
 492 {
 493    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());
 494 }
 495 template <class T1, class T2>
 496 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 497 operator - (const std::complex<T1>& a, const quaternion<T2>& b)
 498 {
 499    return quaternion<T1>(real(a) - b.R_component_1(), imag(a) - b.R_component_2(), -b.R_component_3(), -b.R_component_4());
 500 }
 501 template <class T>
 502 inline BOOST_CONSTEXPR quaternion<T> operator - (const quaternion<T>& a, const quaternion<T>& b)
 503 {
 504    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());
 505 }
 506
 507 // operator*
 508 template <class T1, class T2>
 509 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 510 operator * (const quaternion<T1>& a, const T2& b)
 511 {
 512    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);
 513 }
 514 template <class T1, class T2>
 515 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 516 operator * (const T1& a, const quaternion<T2>& b)
 517 {
 518    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());
 519 }
 520 template <class T1, class T2>
 521 inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 522 operator * (const quaternion<T1>& a, const std::complex<T2>& b)
 523 {
 524    quaternion<T1> result(a);
 525    result *= b;
 526    return result;
 527 }
 528 template <class T1, class T2>
 529 inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 530 operator * (const std::complex<T1>& a, const quaternion<T2>& b)
 531 {
 532    quaternion<T1> result(a);
 533    result *= b;
 534    return result;
 535 }
 536 template <class T>
 537 inline BOOST_CXX14_CONSTEXPR quaternion<T> operator * (const quaternion<T>& a, const quaternion<T>& b)
 538 {
 539    quaternion<T> result(a);
 540    result *= b;
 541    return result;
 542 }
 543
 544 // operator/
 545 template <class T1, class T2>
 546 inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 547 operator / (const quaternion<T1>& a, const T2& b)
 548 {
 549    return quaternion<T1>(a.R_component_1() / b, a.R_component_2() / b, a.R_component_3() / b, a.R_component_4() / b);
 550 }
 551 template <class T1, class T2>
 552 inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 553 operator / (const T1& a, const quaternion<T2>& b)
 554 {
 555    quaternion<T2> result(a);
 556    result /= b;
 557    return result;
 558 }
 559 template <class T1, class T2>
 560 inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type
 561 operator / (const quaternion<T1>& a, const std::complex<T2>& b)
 562 {
 563    quaternion<T1> result(a);
 564    result /= b;
 565    return result;
 566 }
 567 template <class T1, class T2>
 568 inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type
 569 operator / (const std::complex<T1>& a, const quaternion<T2>& b)
 570 {
 571    quaternion<T2> result(a);
 572    result /= b;
 573    return result;
 574 }
 575 template <class T>
 576 inline BOOST_CXX14_CONSTEXPR quaternion<T> operator / (const quaternion<T>& a, const quaternion<T>& b)
 577 {
 578    quaternion<T> result(a);
 579    result /= b;
 580    return result;
 581 }
 582
 583
 584         template<typename T>
 585         inline BOOST_CONSTEXPR const quaternion<T>&             operator + (quaternion<T> const & q)
 586         {
 587             return q;
 588         }
 589
 590
 591         template<typename T>
 592         inline BOOST_CONSTEXPR quaternion<T>                    operator - (quaternion<T> const & q)
 593         {
 594             return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4()));
 595         }
 596
 597
 598         template<typename R, typename T>
 599         inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<R, T>::value, bool>::type operator == (R const & lhs, quaternion<T> const & rhs)
 600         {
 601             return    (
 602                         (rhs.R_component_1() == lhs)&&
 603                         (rhs.R_component_2() == static_cast<T>(0))&&
 604                         (rhs.R_component_3() == static_cast<T>(0))&&
 605                         (rhs.R_component_4() == static_cast<T>(0))
 606                     );
 607         }
 608
 609
 610         template<typename T, typename R>
 611         inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<R, T>::value, bool>::type operator == (quaternion<T> const & lhs, R const & rhs)
 612         {
 613            return rhs == lhs;
 614         }
 615
 616
 617         template<typename T>
 618         inline BOOST_CONSTEXPR bool                                operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs)
 619         {
 620             return    (
 621                         (rhs.R_component_1() == lhs.real())&&
 622                         (rhs.R_component_2() == lhs.imag())&&
 623                         (rhs.R_component_3() == static_cast<T>(0))&&
 624                         (rhs.R_component_4() == static_cast<T>(0))
 625                     );
 626         }
 627
 628
 629         template<typename T>
 630         inline BOOST_CONSTEXPR bool                                operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs)
 631         {
 632            return rhs == lhs;
 633         }
 634
 635
 636         template<typename T>
 637         inline BOOST_CONSTEXPR bool                                operator == (quaternion<T> const & lhs, quaternion<T> const & rhs)
 638         {
 639             return    (
 640                         (rhs.R_component_1() == lhs.R_component_1())&&
 641                         (rhs.R_component_2() == lhs.R_component_2())&&
 642                         (rhs.R_component_3() == lhs.R_component_3())&&
 643                         (rhs.R_component_4() == lhs.R_component_4())
 644                     );
 645         }
 646
 647         template<typename R, typename T> inline BOOST_CONSTEXPR bool operator != (R const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
 648         template<typename T, typename R> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, R const & rhs) { return !(lhs == rhs); }
 649         template<typename T> inline BOOST_CONSTEXPR bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
 650         template<typename T> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs) { return !(lhs == rhs); }
 651         template<typename T> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }
 652
 653
 654         // Note:    we allow the following formats, whith a, b, c, and d reals
 655         //            a
 656         //            (a), (a,b), (a,b,c), (a,b,c,d)
 657         //            (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))
 658         template<typename T, typename charT, class traits>
 659         ::std::basic_istream<charT,traits> &    operator >> (    ::std::basic_istream<charT,traits> & is,
 660                                                                 quaternion<T> & q)
 661         {
 662
 663 #ifdef    BOOST_NO_STD_LOCALE
 664 #else
 665             const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc());
 666 #endif /* BOOST_NO_STD_LOCALE */
 667
 668             T    a = T();
 669             T    b = T();
 670             T    c = T();
 671             T    d = T();
 672
 673             ::std::complex<T>    u = ::std::complex<T>();
 674             ::std::complex<T>    v = ::std::complex<T>();
 675
 676             charT    ch = charT();
 677             char    cc;
 678
 679             is >> ch;                                        // get the first lexeme
 680
 681             if    (!is.good())    goto finish;
 682
 683 #ifdef    BOOST_NO_STD_LOCALE
 684             cc = ch;
 685 #else
 686             cc = ct.narrow(ch, char());
 687 #endif /* BOOST_NO_STD_LOCALE */
 688
 689             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,))
 690             {
 691                 is >> ch;                                    // get the second lexeme
 692
 693                 if    (!is.good())    goto finish;
 694
 695 #ifdef    BOOST_NO_STD_LOCALE
 696                 cc = ch;
 697 #else
 698                 cc = ct.narrow(ch, char());
 699 #endif /* BOOST_NO_STD_LOCALE */
 700
 701                 if    (cc == '(')                        // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))
 702                 {
 703                     is.putback(ch);
 704
 705                     is >> u;                                // we extract the first and second components
 706                     a = u.real();
 707                     b = u.imag();
 708
 709                     if    (!is.good())    goto finish;
 710
 711                     is >> ch;                                // get the next lexeme
 712
 713                     if    (!is.good())    goto finish;
 714
 715 #ifdef    BOOST_NO_STD_LOCALE
 716                     cc = ch;
 717 #else
 718                     cc = ct.narrow(ch, char());
 719 #endif /* BOOST_NO_STD_LOCALE */
 720
 721                     if        (cc == ')')                    // format: ((a)) or ((a,b))
 722                     {
 723                         q = quaternion<T>(a,b);
 724                     }
 725                     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,))
 726                     {
 727                         is >> v;                            // we extract the third and fourth components
 728                         c = v.real();
 729                         d = v.imag();
 730
 731                         if    (!is.good())    goto finish;
 732
 733                         is >> ch;                                // get the last lexeme
 734
 735                         if    (!is.good())    goto finish;
 736
 737 #ifdef    BOOST_NO_STD_LOCALE
 738                         cc = ch;
 739 #else
 740                         cc = ct.narrow(ch, char());
 741 #endif /* BOOST_NO_STD_LOCALE */
 742
 743                         if    (cc == ')')                    // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,))
 744                         {
 745                             q = quaternion<T>(a,b,c,d);
 746                         }
 747                         else                            // error
 748                         {
 749                             is.setstate(::std::ios_base::failbit);
 750                         }
 751                     }
 752                     else                                // error
 753                     {
 754                         is.setstate(::std::ios_base::failbit);
 755                     }
 756                 }
 757                 else                                // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
 758                 {
 759                     is.putback(ch);
 760
 761                     is >> a;                                // we extract the first component
 762
 763                     if    (!is.good())    goto finish;
 764
 765                     is >> ch;                                // get the third lexeme
 766
 767                     if    (!is.good())    goto finish;
 768
 769 #ifdef    BOOST_NO_STD_LOCALE
 770                     cc = ch;
 771 #else
 772                     cc = ct.narrow(ch, char());
 773 #endif /* BOOST_NO_STD_LOCALE */
 774
 775                     if        (cc == ')')                    // format: (a)
 776                     {
 777                         q = quaternion<T>(a);
 778                     }
 779                     else if    (cc == ',')                // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
 780                     {
 781                         is >> ch;                            // get the fourth lexeme
 782
 783                         if    (!is.good())    goto finish;
 784
 785 #ifdef    BOOST_NO_STD_LOCALE
 786                         cc = ch;
 787 #else
 788                         cc = ct.narrow(ch, char());
 789 #endif /* BOOST_NO_STD_LOCALE */
 790
 791                         if    (cc == '(')                // read "(a,(", possible: (a,(c)), (a,(c,d))
 792                         {
 793                             is.putback(ch);
 794
 795                             is >> v;                        // we extract the third and fourth component
 796
 797                             c = v.real();
 798                             d = v.imag();
 799
 800                             if    (!is.good())    goto finish;
 801
 802                             is >> ch;                        // get the ninth lexeme
 803
 804                             if    (!is.good())    goto finish;
 805
 806 #ifdef    BOOST_NO_STD_LOCALE
 807                             cc = ch;
 808 #else
 809                             cc = ct.narrow(ch, char());
 810 #endif /* BOOST_NO_STD_LOCALE */
 811
 812                             if    (cc == ')')                // format: (a,(c)) or (a,(c,d))
 813                             {
 814                                 q = quaternion<T>(a,b,c,d);
 815                             }
 816                             else                        // error
 817                             {
 818                                 is.setstate(::std::ios_base::failbit);
 819                             }
 820                         }
 821                         else                        // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d)
 822                         {
 823                             is.putback(ch);
 824
 825                             is >> b;                        // we extract the second component
 826
 827                             if    (!is.good())    goto finish;
 828
 829                             is >> ch;                        // get the fifth lexeme
 830
 831                             if    (!is.good())    goto finish;
 832
 833 #ifdef    BOOST_NO_STD_LOCALE
 834                             cc = ch;
 835 #else
 836                             cc = ct.narrow(ch, char());
 837 #endif /* BOOST_NO_STD_LOCALE */
 838
 839                             if    (cc == ')')                // format: (a,b)
 840                             {
 841                                 q = quaternion<T>(a,b);
 842                             }
 843                             else if    (cc == ',')        // read "(a,b,", possible: (a,b,c), (a,b,c,d)
 844                             {
 845                                 is >> c;                    // we extract the third component
 846
 847                                 if    (!is.good())    goto finish;
 848
 849                                 is >> ch;                    // get the seventh lexeme
 850
 851                                 if    (!is.good())    goto finish;
 852
 853 #ifdef    BOOST_NO_STD_LOCALE
 854                                 cc = ch;
 855 #else
 856                                 cc = ct.narrow(ch, char());
 857 #endif /* BOOST_NO_STD_LOCALE */
 858
 859                                 if        (cc == ')')        // format: (a,b,c)
 860                                 {
 861                                     q = quaternion<T>(a,b,c);
 862                                 }
 863                                 else if    (cc == ',')    // read "(a,b,c,", possible: (a,b,c,d)
 864                                 {
 865                                     is >> d;                // we extract the fourth component
 866
 867                                     if    (!is.good())    goto finish;
 868
 869                                     is >> ch;                // get the ninth lexeme
 870
 871                                     if    (!is.good())    goto finish;
 872
 873 #ifdef    BOOST_NO_STD_LOCALE
 874                                     cc = ch;
 875 #else
 876                                     cc = ct.narrow(ch, char());
 877 #endif /* BOOST_NO_STD_LOCALE */
 878
 879                                     if    (cc == ')')        // format: (a,b,c,d)
 880                                     {
 881                                         q = quaternion<T>(a,b,c,d);
 882                                     }
 883                                     else                // error
 884                                     {
 885                                         is.setstate(::std::ios_base::failbit);
 886                                     }
 887                                 }
 888                                 else                    // error
 889                                 {
 890                                     is.setstate(::std::ios_base::failbit);
 891                                 }
 892                             }
 893                             else                        // error
 894                             {
 895                                 is.setstate(::std::ios_base::failbit);
 896                             }
 897                         }
 898                     }
 899                     else                                // error
 900                     {
 901                         is.setstate(::std::ios_base::failbit);
 902                     }
 903                 }
 904             }
 905             else                                        // format:    a
 906             {
 907                 is.putback(ch);
 908
 909                 is >> a;                                    // we extract the first component
 910
 911                 if    (!is.good())    goto finish;
 912
 913                 q = quaternion<T>(a);
 914             }
 915
 916             finish:
 917             return(is);
 918         }
 919
 920
 921         template<typename T, typename charT, class traits>
 922         ::std::basic_ostream<charT,traits> &    operator << (    ::std::basic_ostream<charT,traits> & os,
 923                                                                 quaternion<T> const & q)
 924         {
 925             ::std::basic_ostringstream<charT,traits>    s;
 926
 927             s.flags(os.flags());
 928 #ifdef    BOOST_NO_STD_LOCALE
 929 #else
 930             s.imbue(os.getloc());
 931 #endif /* BOOST_NO_STD_LOCALE */
 932             s.precision(os.precision());
 933
 934             s << '('    << q.R_component_1() << ','
 935                         << q.R_component_2() << ','
 936                         << q.R_component_3() << ','
 937                         << q.R_component_4() << ')';
 938
 939             return os << s.str();
 940         }
 941
 942
 943         // values
 944
 945         template<typename T>
 946         inline BOOST_CONSTEXPR T real(quaternion<T> const & q)
 947         {
 948             return(q.real());
 949         }
 950
 951
 952         template<typename T>
 953         inline BOOST_CONSTEXPR quaternion<T> unreal(quaternion<T> const & q)
 954         {
 955             return(q.unreal());
 956         }
 957
 958         template<typename T>
 959         inline T sup(quaternion<T> const & q)
 960         {
 961             using    ::std::abs;
 962             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())));
 963         }
 964
 965
 966         template<typename T>
 967         inline T l1(quaternion<T> const & q)
 968         {
 969            using    ::std::abs;
 970            return abs(q.R_component_1()) + abs(q.R_component_2()) + abs(q.R_component_3()) + abs(q.R_component_4());
 971         }
 972
 973
 974         template<typename T>
 975         inline T abs(quaternion<T> const & q)
 976         {
 977             using    ::std::abs;
 978             using    ::std::sqrt;
 979
 980             T            maxim = sup(q);    // overflow protection
 981
 982             if    (maxim == static_cast<T>(0))
 983             {
 984                 return(maxim);
 985             }
 986             else
 987             {
 988                 T    mixam = static_cast<T>(1)/maxim;    // prefer multiplications over divisions
 989
 990                 T a = q.R_component_1() * mixam;
 991                 T b = q.R_component_2() * mixam;
 992                 T c = q.R_component_3() * mixam;
 993                 T d = q.R_component_4() * mixam;
 994
 995                 a *= a;
 996                 b *= b;
 997                 c *= c;
 998                 d *= d;
 999
1000                 return(maxim * sqrt(a + b + c + d));
1001             }
1002
1003             //return(sqrt(norm(q)));
1004         }
1005
1006
1007         // Note:    This is the Cayley norm, not the Euclidian norm...
1008
1009         template<typename T>
1010         inline BOOST_CXX14_CONSTEXPR T norm(quaternion<T>const  & q)
1011         {
1012             return(real(q*conj(q)));
1013         }
1014
1015
1016         template<typename T>
1017         inline BOOST_CONSTEXPR quaternion<T> conj(quaternion<T> const & q)
1018         {
1019             return(quaternion<T>(   +q.R_component_1(),
1020                                     -q.R_component_2(),
1021                                     -q.R_component_3(),
1022                                     -q.R_component_4()));
1023         }
1024
1025
1026         template<typename T>
1027         inline quaternion<T>                    spherical(  T const & rho,
1028                                                             T const & theta,
1029                                                             T const & phi1,
1030                                                             T const & phi2)
1031         {
1032             using ::std::cos;
1033             using ::std::sin;
1034
1035             //T    a = cos(theta)*cos(phi1)*cos(phi2);
1036             //T    b = sin(theta)*cos(phi1)*cos(phi2);
1037             //T    c = sin(phi1)*cos(phi2);
1038             //T    d = sin(phi2);
1039
1040             T    courrant = static_cast<T>(1);
1041
1042             T    d = sin(phi2);
1043
1044             courrant *= cos(phi2);
1045
1046             T    c = sin(phi1)*courrant;
1047
1048             courrant *= cos(phi1);
1049
1050             T    b = sin(theta)*courrant;
1051             T    a = cos(theta)*courrant;
1052
1053             return(rho*quaternion<T>(a,b,c,d));
1054         }
1055
1056
1057         template<typename T>
1058         inline quaternion<T>                    semipolar(  T const & rho,
1059                                                             T const & alpha,
1060                                                             T const & theta1,
1061                                                             T const & theta2)
1062         {
1063             using ::std::cos;
1064             using ::std::sin;
1065
1066             T    a = cos(alpha)*cos(theta1);
1067             T    b = cos(alpha)*sin(theta1);
1068             T    c = sin(alpha)*cos(theta2);
1069             T    d = sin(alpha)*sin(theta2);
1070
1071             return(rho*quaternion<T>(a,b,c,d));
1072         }
1073
1074
1075         template<typename T>
1076         inline quaternion<T>                    multipolar( T const & rho1,
1077                                                             T const & theta1,
1078                                                             T const & rho2,
1079                                                             T const & theta2)
1080         {
1081             using ::std::cos;
1082             using ::std::sin;
1083
1084             T    a = rho1*cos(theta1);
1085             T    b = rho1*sin(theta1);
1086             T    c = rho2*cos(theta2);
1087             T    d = rho2*sin(theta2);
1088
1089             return(quaternion<T>(a,b,c,d));
1090         }
1091
1092
1093         template<typename T>
1094         inline quaternion<T>                    cylindrospherical(  T const & t,
1095                                                                     T const & radius,
1096                                                                     T const & longitude,
1097                                                                     T const & latitude)
1098         {
1099             using ::std::cos;
1100             using ::std::sin;
1101
1102
1103
1104             T    b = radius*cos(longitude)*cos(latitude);
1105             T    c = radius*sin(longitude)*cos(latitude);
1106             T    d = radius*sin(latitude);
1107
1108             return(quaternion<T>(t,b,c,d));
1109         }
1110
1111
1112         template<typename T>
1113         inline quaternion<T>                    cylindrical(T const & r,
1114                                                             T const & angle,
1115                                                             T const & h1,
1116                                                             T const & h2)
1117         {
1118             using ::std::cos;
1119             using ::std::sin;
1120
1121             T    a = r*cos(angle);
1122             T    b = r*sin(angle);
1123
1124             return(quaternion<T>(a,b,h1,h2));
1125         }
1126
1127
1128         // transcendentals
1129         // (please see the documentation)
1130
1131
1132         template<typename T>
1133         inline quaternion<T>                    exp(quaternion<T> const & q)
1134         {
1135             using    ::std::exp;
1136             using    ::std::cos;
1137
1138             using    ::boost::math::sinc_pi;
1139
1140             T    u = exp(real(q));
1141
1142             T    z = abs(unreal(q));
1143
1144             T    w = sinc_pi(z);
1145
1146             return(u*quaternion<T>(cos(z),
1147                 w*q.R_component_2(), w*q.R_component_3(),
1148                 w*q.R_component_4()));
1149         }
1150
1151
1152         template<typename T>
1153         inline quaternion<T>                    cos(quaternion<T> const & q)
1154         {
1155             using    ::std::sin;
1156             using    ::std::cos;
1157             using    ::std::cosh;
1158
1159             using    ::boost::math::sinhc_pi;
1160
1161             T    z = abs(unreal(q));
1162
1163             T    w = -sin(q.real())*sinhc_pi(z);
1164
1165             return(quaternion<T>(cos(q.real())*cosh(z),
1166                 w*q.R_component_2(), w*q.R_component_3(),
1167                 w*q.R_component_4()));
1168         }
1169
1170
1171         template<typename T>
1172         inline quaternion<T>                    sin(quaternion<T> const & q)
1173         {
1174             using    ::std::sin;
1175             using    ::std::cos;
1176             using    ::std::cosh;
1177
1178             using    ::boost::math::sinhc_pi;
1179
1180             T    z = abs(unreal(q));
1181
1182             T    w = +cos(q.real())*sinhc_pi(z);
1183
1184             return(quaternion<T>(sin(q.real())*cosh(z),
1185                 w*q.R_component_2(), w*q.R_component_3(),
1186                 w*q.R_component_4()));
1187         }
1188
1189
1190         template<typename T>
1191         inline quaternion<T>                    tan(quaternion<T> const & q)
1192         {
1193             return(sin(q)/cos(q));
1194         }
1195
1196
1197         template<typename T>
1198         inline quaternion<T>                    cosh(quaternion<T> const & q)
1199         {
1200             return((exp(+q)+exp(-q))/static_cast<T>(2));
1201         }
1202
1203
1204         template<typename T>
1205         inline quaternion<T>                    sinh(quaternion<T> const & q)
1206         {
1207             return((exp(+q)-exp(-q))/static_cast<T>(2));
1208         }
1209
1210
1211         template<typename T>
1212         inline quaternion<T>                    tanh(quaternion<T> const & q)
1213         {
1214             return(sinh(q)/cosh(q));
1215         }
1216
1217
1218         template<typename T>
1219         quaternion<T>                            pow(quaternion<T> const & q,
1220                                                     int n)
1221         {
1222             if        (n > 1)
1223             {
1224                 int    m = n>>1;
1225
1226                 quaternion<T>    result = pow(q, m);
1227
1228                 result *= result;
1229
1230                 if    (n != (m<<1))
1231                 {
1232                     result *= q; // n odd
1233                 }
1234
1235                 return(result);
1236             }
1237             else if    (n == 1)
1238             {
1239                 return(q);
1240             }
1241             else if    (n == 0)
1242             {
1243                 return(quaternion<T>(static_cast<T>(1)));
1244             }
1245             else    /* n < 0 */
1246             {
1247                 return(pow(quaternion<T>(static_cast<T>(1))/q,-n));
1248             }
1249         }
1250     }
1251 }
1252
1253 #endif /* BOOST_QUATERNION_HPP */