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