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