1191 lines · plain
1// boost quaternion.hpp header file2 3// (C) Copyright Hubert Holin 2001.4// Distributed under the Boost Software License, Version 1.0. (See5// accompanying file LICENSE_1_0.txt or copy at6// 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_HPP11#define BOOST_QUATERNION_HPP12 13#include <boost/math_fwd.hpp>14#include <boost/math/tools/config.hpp>15#include <locale> // for the "<<" operator16 17#include <complex>18#include <iosfwd> // for the "<<" and ">>" operators19#include <sstream> // for the "<<" operator20 21#include <boost/math/special_functions/sinc.hpp> // for the Sinus cardinal22#include <boost/math/special_functions/sinhc.hpp> // for the Hyperbolic Sinus cardinal23#include <boost/math/tools/cxx03_warn.hpp>24 25#include <type_traits>26 27namespace boost28{29 namespace math30 {31 32 namespace detail {33 34 template <class T>35 struct is_trivial_arithmetic_type_imp36 {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_MATH_NO_CXX14_CONSTEXPR50 namespace constexpr_detail51 {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#endif61 62 template<typename T>63 class quaternion64 {65 public:66 67 typedef T value_type;68 69 70 // constructor for H seen as R^471 // (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^287 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 constructor100 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 constructor113 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 // destructor126 // (this is taken care of by the compiler itself)127 128 129 // accessors130 //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 usually134 // nothing simpler (as opposed to the complex number case).135 // However, for practicality, there are accessors for the other components136 // (these are necessary for the templated copy constructor, for instance).137 138 constexpr T real() const139 {140 return(a);141 }142 143 constexpr quaternion<T> unreal() const144 {145 return(quaternion<T>(static_cast<T>(0), b, c, d));146 }147 148 constexpr T R_component_1() const149 {150 return(a);151 }152 153 constexpr T R_component_2() const154 {155 return(b);156 }157 158 constexpr T R_component_3() const159 {160 return(c);161 }162 163 constexpr T R_component_4() const164 {165 return(d);166 }167 168 constexpr ::std::complex<T> C_component_1() const169 {170 return(::std::complex<T>(a, b));171 }172 173 constexpr ::std::complex<T> C_component_2() const174 {175 return(::std::complex<T>(c, d));176 }177 178 BOOST_MATH_CXX14_CONSTEXPR void swap(quaternion& o)179 {180#ifndef BOOST_MATH_NO_CXX14_CONSTEXPR181 using constexpr_detail::swap;182#else183 using std::swap;184#endif185 swap(a, o.a);186 swap(b, o.b);187 swap(c, o.c);188 swap(d, o.d);189 }190 191 // assignment operators192 193 template<typename X>194 BOOST_MATH_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_MATH_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_MATH_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_MATH_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_MATH_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 operators244 //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 where250 // type T throws no exceptions, and one exception-safe version251 // for the case where it might.252 private:253 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & do_add(T const & rhs, const std::true_type&)254 {255 a += rhs;256 return *this;257 }258 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & do_add(T const & rhs, const std::false_type&)259 {260 quaternion<T> result(a + rhs, b, c, d); // exception guard261 swap(result);262 return *this;263 }264 BOOST_MATH_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_MATH_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 guard273 swap(result);274 return *this;275 }276 template <class X>277 BOOST_MATH_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_MATH_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 guard289 swap(result);290 return *this;291 }292 293 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & do_subtract(T const & rhs, const std::true_type&)294 {295 a -= rhs;296 return *this;297 }298 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & do_subtract(T const & rhs, const std::false_type&)299 {300 quaternion<T> result(a - rhs, b, c, d); // exception guard301 swap(result);302 return *this;303 }304 BOOST_MATH_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_MATH_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 guard313 swap(result);314 return *this;315 }316 template <class X>317 BOOST_MATH_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_MATH_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 guard329 swap(result);330 return *this;331 }332 333 BOOST_MATH_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_MATH_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 guard344 swap(result);345 return *this;346 }347 348 BOOST_MATH_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_MATH_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 guard359 swap(result);360 return *this;361 }362 public:363 364 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & operator += (T const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }365 BOOST_MATH_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_MATH_CXX14_CONSTEXPR quaternion<T> & operator += (quaternion<X> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }367 368 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & operator -= (T const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }369 BOOST_MATH_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_MATH_CXX14_CONSTEXPR quaternion<T> & operator -= (quaternion<X> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }371 372 BOOST_MATH_CXX14_CONSTEXPR quaternion<T> & operator *= (T const & rhs) { return do_multiply(rhs, detail::is_trivial_arithmetic_type<T>()); }373 374 BOOST_MATH_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_MATH_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_MATH_CXX14_CONSTEXPR quaternion<T> & operator /= (T const & rhs) { return do_divide(rhs, detail::is_trivial_arithmetic_type<T>()); }397 398 BOOST_MATH_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_MATH_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:427template <class T>428BOOST_MATH_CXX14_CONSTEXPR void swap(quaternion<T>& a, quaternion<T>& b) { a.swap(b); }429 430// operator+431template <class T1, class T2>432inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type433operator + (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}437template <class T1, class T2>438inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type439operator + (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}443template <class T1, class T2>444inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type445operator + (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}449template <class T1, class T2>450inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type451operator + (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}455template <class T>456inline 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-461template <class T1, class T2>462inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type463operator - (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}467template <class T1, class T2>468inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type469operator - (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}473template <class T1, class T2>474inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type475operator - (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}479template <class T1, class T2>480inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type481operator - (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}485template <class T>486inline 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*492template <class T1, class T2>493inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type494operator * (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}498template <class T1, class T2>499inline constexpr typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type500operator * (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}504template <class T1, class T2>505inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type506operator * (const quaternion<T1>& a, const std::complex<T2>& b)507{508 quaternion<T1> result(a);509 result *= b;510 return result;511}512template <class T1, class T2>513inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type514operator * (const std::complex<T1>& a, const quaternion<T2>& b)515{516 quaternion<T1> result(a);517 result *= b;518 return result;519}520template <class T>521inline BOOST_MATH_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/529template <class T1, class T2>530inline constexpr typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type531operator / (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}535template <class T1, class T2>536inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type537operator / (const T1& a, const quaternion<T2>& b)538{539 quaternion<T2> result(a);540 result /= b;541 return result;542}543template <class T1, class T2>544inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T2, T1>::value, quaternion<T1> >::type545operator / (const quaternion<T1>& a, const std::complex<T2>& b)546{547 quaternion<T1> result(a);548 result /= b;549 return result;550}551template <class T1, class T2>552inline BOOST_MATH_CXX14_CONSTEXPR typename std::enable_if<std::is_convertible<T1, T2>::value, quaternion<T2> >::type553operator / (const std::complex<T1>& a, const quaternion<T2>& b)554{555 quaternion<T2> result(a);556 result /= b;557 return result;558}559template <class T>560inline BOOST_MATH_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 reals639 // a640 // (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 lexeme660 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 lexeme668 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 components678 a = u.real();679 b = u.imag();680 681 if (!is.good()) goto finish;682 683 is >> ch; // get the next lexeme684 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 components696 c = v.real();697 d = v.imag();698 699 if (!is.good()) goto finish;700 701 is >> ch; // get the last lexeme702 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 // error712 {713 is.setstate(::std::ios_base::failbit);714 }715 }716 else // error717 {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 component726 727 if (!is.good()) goto finish;728 729 is >> ch; // get the third lexeme730 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 lexeme742 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 component752 753 c = v.real();754 d = v.imag();755 756 if (!is.good()) goto finish;757 758 is >> ch; // get the ninth lexeme759 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 // error769 {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 component778 779 if (!is.good()) goto finish;780 781 is >> ch; // get the fifth lexeme782 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 component794 795 if (!is.good()) goto finish;796 797 is >> ch; // get the seventh lexeme798 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 component810 811 if (!is.good()) goto finish;812 813 is >> ch; // get the ninth lexeme814 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 // error824 {825 is.setstate(::std::ios_base::failbit);826 }827 }828 else // error829 {830 is.setstate(::std::ios_base::failbit);831 }832 }833 else // error834 {835 is.setstate(::std::ios_base::failbit);836 }837 }838 }839 else // error840 {841 is.setstate(::std::ios_base::failbit);842 }843 }844 }845 else // format: a846 {847 is.putback(ch);848 849 is >> a; // we extract the first component850 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 // values881 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 protection918 919 if (maxim == static_cast<T>(0))920 {921 return(maxim);922 }923 else924 {925 T mixam = static_cast<T>(1)/maxim; // prefer multiplications over divisions926 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_MATH_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 // transcendentals1066 // (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 odd1170 }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 */1191