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