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