Ticket #11458: quaternion.hpp

File quaternion.hpp, 77.2 KB (added by Scott Macintire <scott.macintire@…>, 7 years ago)

editted header

Line 
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// <BOOST TRAC DOESNT LIKE LINKS>
7
8// See <BOOST TRAC DOESNT LIKE LINKS> for updates, documentation, and revision history.
9
10#ifndef BOOST_QUATERNION_HPP
11#define BOOST_QUATERNION_HPP
12
13
14#include <complex>
15#include <iosfwd> // for the "<<" and ">>" operators
16#include <sstream> // for the "<<" operator
17
18#include <boost/config.hpp> // for BOOST_NO_STD_LOCALE
19#include <boost/detail/workaround.hpp>
20#ifndef BOOST_NO_STD_LOCALE
21 #include <locale> // for the "<<" operator
22#endif /* BOOST_NO_STD_LOCALE */
23
24
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 <cmath>
29
30namespace boost
31{
32 namespace math
33 {
34 namespace {
35 template<typename T, size_t N> inline
36 T arr_max (const T (&t)[N])
37 {
38#if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP)
39 using ::std::max;
40#endif
41 //int assure_array_size_greater_than_zero[N>0?1:-1];
42 T _max = t[0];
43 for(size_t i=1;++i<N; )
44 {
45 _max = max(_max,t[i]);
46 }
47 return _max;
48 }
49
50 template<typename T, size_t N> inline
51 T arr_max_abs (const T (&t)[N])
52 {
53#if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP)
54 using ::std::max;
55#endif
56 using ::std::abs;
57 //int assure_array_size_greater_than_zero[N>0?1:-1];
58 T _max = abs(t[0]);
59 for(size_t i=1;++i<N; )
60 {
61 _max = max(_max,abs(t[i]));
62 }
63 return _max;
64 }
65
66 template<typename T, size_t N> inline
67 T arr_sum (const T (&t)[N])
68 {
69 //int assure_array_size_greater_than_zero[N>0?1:-1];
70 T _sum = t[0];
71 for(size_t i=1;++i<N; )
72 {
73 _sum+=t[i];
74 }
75 return _sum;
76 }
77
78 template<typename T, size_t N> inline
79 T arr_sum_abs (const T (&t)[N])
80 {
81#if defined(BOOST_NO_ARGUMENT_DEPENDENT_LOOKUP)
82 using ::std::abs;
83#endif
84 //int assure_array_size_greater_than_zero[N>0?1:-1];
85 T _sum = abs(t[0]);
86 for(size_t i=1;++i<N; )
87 {
88 _sum+=abs(t[i]);
89 }
90 return _sum;
91 }
92 }//namespace detail
93
94#define BOOST_QUATERNION_ACCESSOR_GENERATOR(type) \
95 type real() const \
96 { \
97 return(a); \
98 } \
99 \
100 quaternion<type> unreal() const \
101 { \
102 return(quaternion<type>(static_cast<type>(0),b,c,d)); \
103 } \
104 \
105 type R_component_1() const \
106 { \
107 return(a); \
108 } \
109 \
110 type R_component_2() const \
111 { \
112 return(b); \
113 } \
114 \
115 type R_component_3() const \
116 { \
117 return(c); \
118 } \
119 \
120 type R_component_4() const \
121 { \
122 return(d); \
123 } \
124 \
125 ::std::complex<type> C_component_1() const \
126 { \
127 return(::std::complex<type>(a,b)); \
128 } \
129 \
130 ::std::complex<type> C_component_2() const \
131 { \
132 return(::std::complex<type>(c,d)); \
133 }
134
135
136#define BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(type) \
137 template<typename X> \
138 quaternion<type> & operator = (quaternion<X> const & a_affecter) \
139 { \
140 a = static_cast<type>(a_affecter.R_component_1()); \
141 b = static_cast<type>(a_affecter.R_component_2()); \
142 c = static_cast<type>(a_affecter.R_component_3()); \
143 d = static_cast<type>(a_affecter.R_component_4()); \
144 \
145 return(*this); \
146 } \
147 \
148 quaternion<type> & operator = (quaternion<type> const & a_affecter) \
149 { \
150 a = a_affecter.a; \
151 b = a_affecter.b; \
152 c = a_affecter.c; \
153 d = a_affecter.d; \
154 \
155 return(*this); \
156 } \
157 \
158 quaternion<type> & operator = (type const & a_affecter) \
159 { \
160 a = a_affecter; \
161 \
162 b = c = d = static_cast<type>(0); \
163 \
164 return(*this); \
165 } \
166 \
167 quaternion<type> & operator = (::std::complex<type> const & a_affecter) \
168 { \
169 a = a_affecter.real(); \
170 b = a_affecter.imag(); \
171 \
172 c = d = static_cast<type>(0); \
173 \
174 return(*this); \
175 }
176
177
178#define BOOST_QUATERNION_MEMBER_DATA_GENERATOR(type) \
179 type a; \
180 type b; \
181 type c; \
182 type d;
183
184
185 template<typename T>
186 class quaternion
187 {
188 public:
189
190 typedef T value_type;
191
192
193 // constructor for H seen as R^4
194 // (also default constructor)
195
196 explicit quaternion( T const & requested_a = T(),
197 T const & requested_b = T(),
198 T const & requested_c = T(),
199 T const & requested_d = T())
200 : a(requested_a),
201 b(requested_b),
202 c(requested_c),
203 d(requested_d)
204 {
205 // nothing to do!
206 }
207
208
209 // constructor for H seen as C^2
210
211 explicit quaternion( ::std::complex<T> const & z0,
212 ::std::complex<T> const & z1 = ::std::complex<T>())
213 : a(z0.real()),
214 b(z0.imag()),
215 c(z1.real()),
216 d(z1.imag())
217 {
218 // nothing to do!
219 }
220
221
222 // UNtemplated copy constructor
223 // (this is taken care of by the compiler itself)
224
225
226 // templated copy constructor
227
228 template<typename X>
229 explicit quaternion(quaternion<X> const & a_recopier)
230 : a(static_cast<T>(a_recopier.R_component_1())),
231 b(static_cast<T>(a_recopier.R_component_2())),
232 c(static_cast<T>(a_recopier.R_component_3())),
233 d(static_cast<T>(a_recopier.R_component_4()))
234 {
235 // nothing to do!
236 }
237
238
239 // destructor
240 // (this is taken care of by the compiler itself)
241
242
243 // accessors
244 //
245 // Note: Like complex number, quaternions do have a meaningful notion of "real part",
246 // but unlike them there is no meaningful notion of "imaginary part".
247 // Instead there is an "unreal part" which itself is a quaternion, and usually
248 // nothing simpler (as opposed to the complex number case).
249 // However, for practicallity, there are accessors for the other components
250 // (these are necessary for the templated copy constructor, for instance).
251
252 BOOST_QUATERNION_ACCESSOR_GENERATOR(T)
253
254 // assignment operators
255
256 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(T)
257
258 // other assignment-related operators
259 //
260 // NOTE: Quaternion multiplication is *NOT* commutative;
261 // symbolically, "q *= rhs;" means "q = q * rhs;"
262 // and "q /= rhs;" means "q = q * inverse_of(rhs);"
263
264 quaternion<T> & operator += (T const & rhs)
265 {
266 T at = a + rhs; // exception guard
267
268 a = at;
269
270 return(*this);
271 }
272
273
274 quaternion<T> & operator += (::std::complex<T> const & rhs)
275 {
276 T at = a + rhs.real(); // exception guard
277 T bt = b + rhs.imag(); // exception guard
278
279 a = at;
280 b = bt;
281
282 return(*this);
283 }
284
285
286 template<typename X>
287 quaternion<T> & operator += (quaternion<X> const & rhs)
288 {
289 T at = a + static_cast<T>(rhs.R_component_1()); // exception guard
290 T bt = b + static_cast<T>(rhs.R_component_2()); // exception guard
291 T ct = c + static_cast<T>(rhs.R_component_3()); // exception guard
292 T dt = d + static_cast<T>(rhs.R_component_4()); // exception guard
293
294 a = at;
295 b = bt;
296 c = ct;
297 d = dt;
298
299 return(*this);
300 }
301
302
303
304 quaternion<T> & operator -= (T const & rhs)
305 {
306 T at = a - rhs; // exception guard
307
308 a = at;
309
310 return(*this);
311 }
312
313
314 quaternion<T> & operator -= (::std::complex<T> const & rhs)
315 {
316 T at = a - rhs.real(); // exception guard
317 T bt = b - rhs.imag(); // exception guard
318
319 a = at;
320 b = bt;
321
322 return(*this);
323 }
324
325
326 template<typename X>
327 quaternion<T> & operator -= (quaternion<X> const & rhs)
328 {
329 T at = a - static_cast<T>(rhs.R_component_1()); // exception guard
330 T bt = b - static_cast<T>(rhs.R_component_2()); // exception guard
331 T ct = c - static_cast<T>(rhs.R_component_3()); // exception guard
332 T dt = d - static_cast<T>(rhs.R_component_4()); // exception guard
333
334 a = at;
335 b = bt;
336 c = ct;
337 d = dt;
338
339 return(*this);
340 }
341
342
343 quaternion<T> & operator *= (T const & rhs)
344 {
345 T at = a * rhs; // exception guard
346 T bt = b * rhs; // exception guard
347 T ct = c * rhs; // exception guard
348 T dt = d * rhs; // exception guard
349
350 a = at;
351 b = bt;
352 c = ct;
353 d = dt;
354
355 return(*this);
356 }
357
358
359 quaternion<T> & operator *= (::std::complex<T> const & rhs)
360 {
361 T ar = rhs.real();
362 T br = rhs.imag();
363
364 T at = +a*ar-b*br;
365 T bt = +a*br+b*ar;
366 T ct = +c*ar+d*br;
367 T dt = -c*br+d*ar;
368
369 a = at;
370 b = bt;
371 c = ct;
372 d = dt;
373
374 return(*this);
375 }
376
377
378 template<typename X>
379 quaternion<T> & operator *= (quaternion<X> const & rhs)
380 {
381 T ar = static_cast<T>(rhs.R_component_1());
382 T br = static_cast<T>(rhs.R_component_2());
383 T cr = static_cast<T>(rhs.R_component_3());
384 T dr = static_cast<T>(rhs.R_component_4());
385
386 T at = +a*ar-b*br-c*cr-d*dr;
387 T bt = +a*br+b*ar+c*dr-d*cr; //(a*br+ar*b)+(c*dr-cr*d);
388 T ct = +a*cr-b*dr+c*ar+d*br; //(a*cr+ar*c)+(d*br-dr*b);
389 T dt = +a*dr+b*cr-c*br+d*ar; //(a*dr+ar*d)+(b*cr-br*c);
390
391 a = at;
392 b = bt;
393 c = ct;
394 d = dt;
395
396 return(*this);
397 }
398
399
400
401 quaternion<T> & operator /= (T const & rhs)
402 {
403 T at = a / rhs; // exception guard
404 T bt = b / rhs; // exception guard
405 T ct = c / rhs; // exception guard
406 T dt = d / rhs; // exception guard
407
408 a = at;
409 b = bt;
410 c = ct;
411 d = dt;
412
413 return(*this);
414 }
415
416
417 quaternion<T> & operator /= (::std::complex<T> const & rhs)
418 {
419 T ar = rhs.real();
420 T br = rhs.imag();
421
422 T denominator = ar*ar+br*br;
423
424 T at = (+a*ar+b*br)/denominator; //(a*ar+b*br)/denominator;
425 T bt = (-a*br+b*ar)/denominator; //(ar*b-a*br)/denominator;
426 T ct = (+c*ar-d*br)/denominator; //(ar*c-d*br)/denominator;
427 T dt = (+c*br+d*ar)/denominator; //(ar*d+br*c)/denominator;
428
429 a = at;
430 b = bt;
431 c = ct;
432 d = dt;
433
434 return(*this);
435 }
436
437
438 template<typename X>
439 quaternion<T> & operator /= (quaternion<X> const & rhs)
440 {
441 T ar = static_cast<T>(rhs.R_component_1());
442 T br = static_cast<T>(rhs.R_component_2());
443 T cr = static_cast<T>(rhs.R_component_3());
444 T dr = static_cast<T>(rhs.R_component_4());
445
446 T denominator = ar*ar+br*br+cr*cr+dr*dr;
447
448 T at = (+a*ar+b*br+c*cr+d*dr)/denominator; //(a*ar+b*br+c*cr+d*dr)/denominator;
449 T bt = (-a*br+b*ar-c*dr+d*cr)/denominator; //((ar*b-a*br)+(cr*d-c*dr))/denominator;
450 T ct = (-a*cr+b*dr+c*ar-d*br)/denominator; //((ar*c-a*cr)+(dr*b-d*br))/denominator;
451 T dt = (-a*dr-b*cr+c*br+d*ar)/denominator; //((ar*d-a*dr)+(br*c-b*cr))/denominator;
452
453 a = at;
454 b = bt;
455 c = ct;
456 d = dt;
457
458 return(*this);
459 }
460
461
462 protected:
463
464 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(T)
465
466
467 private:
468
469 };
470
471
472 // declaration of quaternion specialization
473
474 template<> class quaternion<float>;
475 template<> class quaternion<double>;
476 template<> class quaternion<long double>;
477
478
479 // helper templates for converting copy constructors (declaration)
480
481 namespace detail
482 {
483
484 template< typename T,
485 typename U
486 >
487 quaternion<T> quaternion_type_converter(quaternion<U> const & rhs);
488 }
489
490
491 // implementation of quaternion specialization
492
493
494#define BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(type) \
495 explicit quaternion( type const & requested_a = static_cast<type>(0), \
496 type const & requested_b = static_cast<type>(0), \
497 type const & requested_c = static_cast<type>(0), \
498 type const & requested_d = static_cast<type>(0)) \
499 : a(requested_a), \
500 b(requested_b), \
501 c(requested_c), \
502 d(requested_d) \
503 { \
504 } \
505 \
506 explicit quaternion( ::std::complex<type> const & z0, \
507 ::std::complex<type> const & z1 = ::std::complex<type>()) \
508 : a(z0.real()), \
509 b(z0.imag()), \
510 c(z1.real()), \
511 d(z1.imag()) \
512 { \
513 }
514
515
516#define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \
517 quaternion<type> & operator += (type const & rhs) \
518 { \
519 a += rhs; \
520 \
521 return(*this); \
522 }
523
524#define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \
525 quaternion<type> & operator += (::std::complex<type> const & rhs) \
526 { \
527 a += rhs.real(); \
528 b += rhs.imag(); \
529 \
530 return(*this); \
531 }
532
533#define BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type) \
534 template<typename X> \
535 quaternion<type> & operator += (quaternion<X> const & rhs) \
536 { \
537 a += static_cast<type>(rhs.R_component_1()); \
538 b += static_cast<type>(rhs.R_component_2()); \
539 c += static_cast<type>(rhs.R_component_3()); \
540 d += static_cast<type>(rhs.R_component_4()); \
541 \
542 return(*this); \
543 }
544
545#define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \
546 quaternion<type> & operator -= (type const & rhs) \
547 { \
548 a -= rhs; \
549 \
550 return(*this); \
551 }
552
553#define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \
554 quaternion<type> & operator -= (::std::complex<type> const & rhs) \
555 { \
556 a -= rhs.real(); \
557 b -= rhs.imag(); \
558 \
559 return(*this); \
560 }
561
562#define BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type) \
563 template<typename X> \
564 quaternion<type> & operator -= (quaternion<X> const & rhs) \
565 { \
566 a -= static_cast<type>(rhs.R_component_1()); \
567 b -= static_cast<type>(rhs.R_component_2()); \
568 c -= static_cast<type>(rhs.R_component_3()); \
569 d -= static_cast<type>(rhs.R_component_4()); \
570 \
571 return(*this); \
572 }
573
574#define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \
575 quaternion<type> & operator *= (type const & rhs) \
576 { \
577 a *= rhs; \
578 b *= rhs; \
579 c *= rhs; \
580 d *= rhs; \
581 \
582 return(*this); \
583 }
584
585#define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \
586 quaternion<type> & operator *= (::std::complex<type> const & rhs) \
587 { \
588 type ar = rhs.real(); \
589 type br = rhs.imag(); \
590 \
591 type at = +a*ar-b*br; \
592 type bt = +a*br+b*ar; \
593 type ct = +c*ar+d*br; \
594 type dt = -c*br+d*ar; \
595 \
596 a = at; \
597 b = bt; \
598 c = ct; \
599 d = dt; \
600 \
601 return(*this); \
602 }
603
604#define BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type) \
605 template<typename X> \
606 quaternion<type> & operator *= (quaternion<X> const & rhs) \
607 { \
608 type ar = static_cast<type>(rhs.R_component_1()); \
609 type br = static_cast<type>(rhs.R_component_2()); \
610 type cr = static_cast<type>(rhs.R_component_3()); \
611 type dr = static_cast<type>(rhs.R_component_4()); \
612 \
613 type at = +a*ar-b*br-c*cr-d*dr; \
614 type bt = +a*br+b*ar+c*dr-d*cr; \
615 type ct = +a*cr-b*dr+c*ar+d*br; \
616 type dt = +a*dr+b*cr-c*br+d*ar; \
617 \
618 a = at; \
619 b = bt; \
620 c = ct; \
621 d = dt; \
622 \
623 return(*this); \
624 }
625
626// There is quite a lot of repetition in the code below. This is intentional.
627// The last conditional block is the normal form, and the others merely
628// consist of workarounds for various compiler deficiencies. Hopefuly, when
629// more compilers are conformant and we can retire support for those that are
630// not, we will be able to remove the clutter. This is makes the situation
631// (painfully) explicit.
632
633#define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \
634 quaternion<type> & operator /= (type const & rhs) \
635 { \
636 a /= rhs; \
637 b /= rhs; \
638 c /= rhs; \
639 d /= rhs; \
640 \
641 return(*this); \
642 }
643
644#define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \
645 quaternion<type> & operator /= (::std::complex<type> const & rhs) \
646 { \
647 \
648 type tr[2]; \
649 \
650 tr[0] = rhs.real(); \
651 tr[1] = rhs.imag(); \
652 \
653 type mixam = static_cast<type>(1)/(arr_max_abs(tr)); \
654 \
655 tr[0] *= mixam; \
656 tr[1] *= mixam; \
657 \
658 type tt[4]; \
659 \
660 tt[0] = +a*tr[0]+b*tr[1]; \
661 tt[1] = -a*tr[1]+b*tr[0]; \
662 tt[2] = +c*tr[0]-d*tr[1]; \
663 tt[3] = +c*tr[1]+d*tr[0]; \
664 \
665 tr[0] *= tr[0]; \
666 tr[1] *= tr[1]; \
667 \
668 type mixam_tr_sum = mixam/arr_sum(tr); \
669 tt[0] *= mixam_tr_sum; \
670 tt[1] *= mixam_tr_sum; \
671 tt[2] *= mixam_tr_sum; \
672 tt[3] *= mixam_tr_sum; \
673 \
674 a = tt[0]; \
675 b = tt[1]; \
676 c = tt[2]; \
677 d = tt[3]; \
678 \
679 return(*this); \
680 }
681
682
683#define BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type) \
684 template<typename X> \
685 quaternion<type> & operator /= (quaternion<X> const & rhs) \
686 { \
687 \
688 type tr[4]; \
689 \
690 tr[0] = static_cast<type>(rhs.R_component_1()); \
691 tr[1] = static_cast<type>(rhs.R_component_2()); \
692 tr[2] = static_cast<type>(rhs.R_component_3()); \
693 tr[3] = static_cast<type>(rhs.R_component_4()); \
694 \
695 type mixam = static_cast<type>(1)/(arr_max_abs(tr)); \
696 \
697 tr[0] *= mixam; \
698 tr[1] *= mixam; \
699 tr[2] *= mixam; \
700 tr[3] *= mixam; \
701 \
702 type tt[4]; \
703 \
704 tt[0] = +a*tr[0]+b*tr[1]+c*tr[2]+d*tr[3]; \
705 tt[1] = -a*tr[1]+b*tr[0]-c*tr[3]+d*tr[2]; \
706 tt[2] = -a*tr[2]+b*tr[3]+c*tr[0]-d*tr[1]; \
707 tt[3] = -a*tr[3]-b*tr[2]+c*tr[1]+d*tr[0]; \
708 \
709 tr[0] *= tr[0]; \
710 tr[1] *= tr[1]; \
711 tr[2] *= tr[2]; \
712 tr[3] *= tr[3]; \
713 \
714 type mixam_tr_sum = mixam/arr_sum(tr); \
715 tt[0] *= mixam_tr_sum; \
716 tt[1] *= mixam_tr_sum; \
717 tt[2] *= mixam_tr_sum; \
718 tt[3] *= mixam_tr_sum; \
719 \
720 a = tt[0]; \
721 b = tt[1]; \
722 c = tt[2]; \
723 d = tt[3]; \
724 \
725 return(*this); \
726 }
727
728#define BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \
729 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type) \
730 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type) \
731 BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type)
732
733#define BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \
734 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type) \
735 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type) \
736 BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type)
737
738#define BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \
739 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1(type) \
740 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2(type) \
741 BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3(type)
742
743#define BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type) \
744 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1(type) \
745 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2(type) \
746 BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3(type)
747
748#define BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(type) \
749 BOOST_QUATERNION_MEMBER_ADD_GENERATOR(type) \
750 BOOST_QUATERNION_MEMBER_SUB_GENERATOR(type) \
751 BOOST_QUATERNION_MEMBER_MUL_GENERATOR(type) \
752 BOOST_QUATERNION_MEMBER_DIV_GENERATOR(type)
753
754
755 template<>
756 class quaternion<float>
757 {
758 public:
759
760 typedef float value_type;
761
762 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(float)
763
764 // UNtemplated copy constructor
765 // (this is taken care of by the compiler itself)
766
767 // explicit copy constructors (precision-loosing converters)
768
769 explicit quaternion(quaternion<double> const & a_recopier)
770 {
771 *this = detail::quaternion_type_converter<float, double>(a_recopier);
772 }
773
774 explicit quaternion(quaternion<long double> const & a_recopier)
775 {
776 *this = detail::quaternion_type_converter<float, long double>(a_recopier);
777 }
778
779 // destructor
780 // (this is taken care of by the compiler itself)
781
782 // accessors
783 //
784 // Note: Like complex number, quaternions do have a meaningful notion of "real part",
785 // but unlike them there is no meaningful notion of "imaginary part".
786 // Instead there is an "unreal part" which itself is a quaternion, and usually
787 // nothing simpler (as opposed to the complex number case).
788 // However, for practicallity, there are accessors for the other components
789 // (these are necessary for the templated copy constructor, for instance).
790
791 BOOST_QUATERNION_ACCESSOR_GENERATOR(float)
792
793 // assignment operators
794
795 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(float)
796
797 // other assignment-related operators
798 //
799 // NOTE: Quaternion multiplication is *NOT* commutative;
800 // symbolically, "q *= rhs;" means "q = q * rhs;"
801 // and "q /= rhs;" means "q = q * inverse_of(rhs);"
802
803 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(float)
804
805
806 protected:
807
808 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(float)
809
810
811 private:
812
813 };
814
815
816 template<>
817 class quaternion<double>
818 {
819 public:
820
821 typedef double value_type;
822
823 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(double)
824
825 // UNtemplated copy constructor
826 // (this is taken care of by the compiler itself)
827
828 // converting copy constructor
829
830 explicit quaternion(quaternion<float> const & a_recopier)
831 {
832 *this = detail::quaternion_type_converter<double, float>(a_recopier);
833 }
834
835 // explicit copy constructors (precision-loosing converters)
836
837 explicit quaternion(quaternion<long double> const & a_recopier)
838 {
839 *this = detail::quaternion_type_converter<double, long double>(a_recopier);
840 }
841
842 // destructor
843 // (this is taken care of by the compiler itself)
844
845 // accessors
846 //
847 // Note: Like complex number, quaternions do have a meaningful notion of "real part",
848 // but unlike them there is no meaningful notion of "imaginary part".
849 // Instead there is an "unreal part" which itself is a quaternion, and usually
850 // nothing simpler (as opposed to the complex number case).
851 // However, for practicallity, there are accessors for the other components
852 // (these are necessary for the templated copy constructor, for instance).
853
854 BOOST_QUATERNION_ACCESSOR_GENERATOR(double)
855
856 // assignment operators
857
858 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(double)
859
860 // other assignment-related operators
861 //
862 // NOTE: Quaternion multiplication is *NOT* commutative;
863 // symbolically, "q *= rhs;" means "q = q * rhs;"
864 // and "q /= rhs;" means "q = q * inverse_of(rhs);"
865
866 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(double)
867
868
869 protected:
870
871 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(double)
872
873
874 private:
875
876 };
877
878
879 template<>
880 class quaternion<long double>
881 {
882 public:
883
884 typedef long double value_type;
885
886 BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(long double)
887
888 // UNtemplated copy constructor
889 // (this is taken care of by the compiler itself)
890
891 // converting copy constructors
892
893 explicit quaternion(quaternion<float> const & a_recopier)
894 {
895 *this = detail::quaternion_type_converter<long double, float>(a_recopier);
896 }
897
898 explicit quaternion(quaternion<double> const & a_recopier)
899 {
900 *this = detail::quaternion_type_converter<long double, double>(a_recopier);
901 }
902
903 // destructor
904 // (this is taken care of by the compiler itself)
905
906 // accessors
907 //
908 // Note: Like complex number, quaternions do have a meaningful notion of "real part",
909 // but unlike them there is no meaningful notion of "imaginary part".
910 // Instead there is an "unreal part" which itself is a quaternion, and usually
911 // nothing simpler (as opposed to the complex number case).
912 // However, for practicallity, there are accessors for the other components
913 // (these are necessary for the templated copy constructor, for instance).
914
915 BOOST_QUATERNION_ACCESSOR_GENERATOR(long double)
916
917 // assignment operators
918
919 BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(long double)
920
921 // other assignment-related operators
922 //
923 // NOTE: Quaternion multiplication is *NOT* commutative;
924 // symbolically, "q *= rhs;" means "q = q * rhs;"
925 // and "q /= rhs;" means "q = q * inverse_of(rhs);"
926
927 BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR(long double)
928
929
930 protected:
931
932 BOOST_QUATERNION_MEMBER_DATA_GENERATOR(long double)
933
934
935 private:
936
937 };
938
939
940#undef BOOST_QUATERNION_MEMBER_ALGEBRAIC_GENERATOR
941#undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR
942#undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR
943#undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR
944#undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR
945#undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1
946#undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2
947#undef BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3
948#undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1
949#undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2
950#undef BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3
951#undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_1
952#undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_2
953#undef BOOST_QUATERNION_MEMBER_MUL_GENERATOR_3
954#undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_1
955#undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_2
956#undef BOOST_QUATERNION_MEMBER_DIV_GENERATOR_3
957
958#undef BOOST_QUATERNION_CONSTRUCTOR_GENERATOR
959
960
961#undef BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR
962
963#undef BOOST_QUATERNION_MEMBER_DATA_GENERATOR
964
965#undef BOOST_QUATERNION_ACCESSOR_GENERATOR
966
967
968 // operators
969
970#define BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op) \
971 { \
972 quaternion<T> res(lhs); \
973 res op##= rhs; \
974 return(res); \
975 }
976
977#define BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \
978 template<typename T> \
979 inline quaternion<T> operator op (T const & lhs, quaternion<T> const & rhs) \
980 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op)
981
982#define BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \
983 template<typename T> \
984 inline quaternion<T> operator op (quaternion<T> const & lhs, T const & rhs) \
985 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op)
986
987#define BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \
988 template<typename T> \
989 inline quaternion<T> operator op (::std::complex<T> const & lhs, quaternion<T> const & rhs) \
990 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op)
991
992#define BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \
993 template<typename T> \
994 inline quaternion<T> operator op (quaternion<T> const & lhs, ::std::complex<T> const & rhs) \
995 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op)
996
997#define BOOST_QUATERNION_OPERATOR_GENERATOR_3(op) \
998 template<typename T> \
999 inline quaternion<T> operator op (quaternion<T> const & lhs, quaternion<T> const & rhs) \
1000 BOOST_QUATERNION_OPERATOR_GENERATOR_BODY(op)
1001
1002#define BOOST_QUATERNION_OPERATOR_GENERATOR(op) \
1003 BOOST_QUATERNION_OPERATOR_GENERATOR_1_L(op) \
1004 BOOST_QUATERNION_OPERATOR_GENERATOR_1_R(op) \
1005 BOOST_QUATERNION_OPERATOR_GENERATOR_2_L(op) \
1006 BOOST_QUATERNION_OPERATOR_GENERATOR_2_R(op) \
1007 BOOST_QUATERNION_OPERATOR_GENERATOR_3(op)
1008
1009
1010 BOOST_QUATERNION_OPERATOR_GENERATOR(+)
1011 BOOST_QUATERNION_OPERATOR_GENERATOR(-)
1012 BOOST_QUATERNION_OPERATOR_GENERATOR(*)
1013 BOOST_QUATERNION_OPERATOR_GENERATOR(/)
1014
1015
1016#undef BOOST_QUATERNION_OPERATOR_GENERATOR
1017
1018#undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_L
1019#undef BOOST_QUATERNION_OPERATOR_GENERATOR_1_R
1020#undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_L
1021#undef BOOST_QUATERNION_OPERATOR_GENERATOR_2_R
1022#undef BOOST_QUATERNION_OPERATOR_GENERATOR_3
1023
1024#undef BOOST_QUATERNION_OPERATOR_GENERATOR_BODY
1025
1026
1027 template<typename T>
1028 inline quaternion<T> operator + (quaternion<T> const & q)
1029 {
1030 return(q);
1031 }
1032
1033
1034 template<typename T>
1035 inline quaternion<T> operator - (quaternion<T> const & q)
1036 {
1037 return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4()));
1038 }
1039
1040
1041 template<typename T>
1042 inline bool operator == (T const & lhs, quaternion<T> const & rhs)
1043 {
1044 return (
1045 (rhs.R_component_1() == lhs)&&
1046 (rhs.R_component_2() == static_cast<T>(0))&&
1047 (rhs.R_component_3() == static_cast<T>(0))&&
1048 (rhs.R_component_4() == static_cast<T>(0))
1049 );
1050 }
1051
1052
1053 template<typename T>
1054 inline bool operator == (quaternion<T> const & lhs, T const & rhs)
1055 {
1056 return (
1057 (lhs.R_component_1() == rhs)&&
1058 (lhs.R_component_2() == static_cast<T>(0))&&
1059 (lhs.R_component_3() == static_cast<T>(0))&&
1060 (lhs.R_component_4() == static_cast<T>(0))
1061 );
1062 }
1063
1064
1065 template<typename T>
1066 inline bool operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs)
1067 {
1068 return (
1069 (rhs.R_component_1() == lhs.real())&&
1070 (rhs.R_component_2() == lhs.imag())&&
1071 (rhs.R_component_3() == static_cast<T>(0))&&
1072 (rhs.R_component_4() == static_cast<T>(0))
1073 );
1074 }
1075
1076
1077 template<typename T>
1078 inline bool operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs)
1079 {
1080 return (
1081 (lhs.R_component_1() == rhs.real())&&
1082 (lhs.R_component_2() == rhs.imag())&&
1083 (lhs.R_component_3() == static_cast<T>(0))&&
1084 (lhs.R_component_4() == static_cast<T>(0))
1085 );
1086 }
1087
1088
1089 template<typename T>
1090 inline bool operator == (quaternion<T> const & lhs, quaternion<T> const & rhs)
1091 {
1092 return (
1093 (rhs.R_component_1() == lhs.R_component_1())&&
1094 (rhs.R_component_2() == lhs.R_component_2())&&
1095 (rhs.R_component_3() == lhs.R_component_3())&&
1096 (rhs.R_component_4() == lhs.R_component_4())
1097 );
1098 }
1099
1100
1101#define BOOST_QUATERNION_NOT_EQUAL_GENERATOR \
1102 { \
1103 return(!(lhs == rhs)); \
1104 }
1105
1106 template<typename T>
1107 inline bool operator != (T const & lhs, quaternion<T> const & rhs)
1108 BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1109
1110 template<typename T>
1111 inline bool operator != (quaternion<T> const & lhs, T const & rhs)
1112 BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1113
1114 template<typename T>
1115 inline bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs)
1116 BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1117
1118 template<typename T>
1119 inline bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs)
1120 BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1121
1122 template<typename T>
1123 inline bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs)
1124 BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1125
1126#undef BOOST_QUATERNION_NOT_EQUAL_GENERATOR
1127
1128
1129 // Note: we allow the following formats, whith a, b, c, and d reals
1130 // a
1131 // (a), (a,b), (a,b,c), (a,b,c,d)
1132 // (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))
1133 template<typename T, typename charT, class traits>
1134 ::std::basic_istream<charT,traits> & operator >> ( ::std::basic_istream<charT,traits> & is,
1135 quaternion<T> & q)
1136 {
1137
1138#ifdef BOOST_NO_STD_LOCALE
1139#else
1140 const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc());
1141#endif /* BOOST_NO_STD_LOCALE */
1142
1143 T a = T();
1144 T b = T();
1145 T c = T();
1146 T d = T();
1147
1148 ::std::complex<T> u = ::std::complex<T>();
1149 ::std::complex<T> v = ::std::complex<T>();
1150
1151 charT ch = charT();
1152 char cc;
1153
1154 is >> ch; // get the first lexeme
1155
1156 if (!is.good()) goto finish;
1157
1158#ifdef BOOST_NO_STD_LOCALE
1159 cc = ch;
1160#else
1161 cc = ct.narrow(ch, char());
1162#endif /* BOOST_NO_STD_LOCALE */
1163
1164 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,))
1165 {
1166 is >> ch; // get the second lexeme
1167
1168 if (!is.good()) goto finish;
1169
1170#ifdef BOOST_NO_STD_LOCALE
1171 cc = ch;
1172#else
1173 cc = ct.narrow(ch, char());
1174#endif /* BOOST_NO_STD_LOCALE */
1175
1176 if (cc == '(') // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))
1177 {
1178 is.putback(ch);
1179
1180 is >> u; // we extract the first and second components
1181 a = u.real();
1182 b = u.imag();
1183
1184 if (!is.good()) goto finish;
1185
1186 is >> ch; // get the next lexeme
1187
1188 if (!is.good()) goto finish;
1189
1190#ifdef BOOST_NO_STD_LOCALE
1191 cc = ch;
1192#else
1193 cc = ct.narrow(ch, char());
1194#endif /* BOOST_NO_STD_LOCALE */
1195
1196 if (cc == ')') // format: ((a)) or ((a,b))
1197 {
1198 q = quaternion<T>(a,b);
1199 }
1200 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,))
1201 {
1202 is >> v; // we extract the third and fourth components
1203 c = v.real();
1204 d = v.imag();
1205
1206 if (!is.good()) goto finish;
1207
1208 is >> ch; // get the last lexeme
1209
1210 if (!is.good()) goto finish;
1211
1212#ifdef BOOST_NO_STD_LOCALE
1213 cc = ch;
1214#else
1215 cc = ct.narrow(ch, char());
1216#endif /* BOOST_NO_STD_LOCALE */
1217
1218 if (cc == ')') // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,))
1219 {
1220 q = quaternion<T>(a,b,c,d);
1221 }
1222 else // error
1223 {
1224 is.setstate(::std::ios_base::failbit);
1225 }
1226 }
1227 else // error
1228 {
1229 is.setstate(::std::ios_base::failbit);
1230 }
1231 }
1232 else // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
1233 {
1234 is.putback(ch);
1235
1236 is >> a; // we extract the first component
1237
1238 if (!is.good()) goto finish;
1239
1240 is >> ch; // get the third lexeme
1241
1242 if (!is.good()) goto finish;
1243
1244#ifdef BOOST_NO_STD_LOCALE
1245 cc = ch;
1246#else
1247 cc = ct.narrow(ch, char());
1248#endif /* BOOST_NO_STD_LOCALE */
1249
1250 if (cc == ')') // format: (a)
1251 {
1252 q = quaternion<T>(a);
1253 }
1254 else if (cc == ',') // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))
1255 {
1256 is >> ch; // get the fourth lexeme
1257
1258 if (!is.good()) goto finish;
1259
1260#ifdef BOOST_NO_STD_LOCALE
1261 cc = ch;
1262#else
1263 cc = ct.narrow(ch, char());
1264#endif /* BOOST_NO_STD_LOCALE */
1265
1266 if (cc == '(') // read "(a,(", possible: (a,(c)), (a,(c,d))
1267 {
1268 is.putback(ch);
1269
1270 is >> v; // we extract the third and fourth component
1271
1272 c = v.real();
1273 d = v.imag();
1274
1275 if (!is.good()) goto finish;
1276
1277 is >> ch; // get the ninth lexeme
1278
1279 if (!is.good()) goto finish;
1280
1281#ifdef BOOST_NO_STD_LOCALE
1282 cc = ch;
1283#else
1284 cc = ct.narrow(ch, char());
1285#endif /* BOOST_NO_STD_LOCALE */
1286
1287 if (cc == ')') // format: (a,(c)) or (a,(c,d))
1288 {
1289 q = quaternion<T>(a,b,c,d);
1290 }
1291 else // error
1292 {
1293 is.setstate(::std::ios_base::failbit);
1294 }
1295 }
1296 else // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d)
1297 {
1298 is.putback(ch);
1299
1300 is >> b; // we extract the second component
1301
1302 if (!is.good()) goto finish;
1303
1304 is >> ch; // get the fifth lexeme
1305
1306 if (!is.good()) goto finish;
1307
1308#ifdef BOOST_NO_STD_LOCALE
1309 cc = ch;
1310#else
1311 cc = ct.narrow(ch, char());
1312#endif /* BOOST_NO_STD_LOCALE */
1313
1314 if (cc == ')') // format: (a,b)
1315 {
1316 q = quaternion<T>(a,b);
1317 }
1318 else if (cc == ',') // read "(a,b,", possible: (a,b,c), (a,b,c,d)
1319 {
1320 is >> c; // we extract the third component
1321
1322 if (!is.good()) goto finish;
1323
1324 is >> ch; // get the seventh lexeme
1325
1326 if (!is.good()) goto finish;
1327
1328#ifdef BOOST_NO_STD_LOCALE
1329 cc = ch;
1330#else
1331 cc = ct.narrow(ch, char());
1332#endif /* BOOST_NO_STD_LOCALE */
1333
1334 if (cc == ')') // format: (a,b,c)
1335 {
1336 q = quaternion<T>(a,b,c);
1337 }
1338 else if (cc == ',') // read "(a,b,c,", possible: (a,b,c,d)
1339 {
1340 is >> d; // we extract the fourth component
1341
1342 if (!is.good()) goto finish;
1343
1344 is >> ch; // get the ninth lexeme
1345
1346 if (!is.good()) goto finish;
1347
1348#ifdef BOOST_NO_STD_LOCALE
1349 cc = ch;
1350#else
1351 cc = ct.narrow(ch, char());
1352#endif /* BOOST_NO_STD_LOCALE */
1353
1354 if (cc == ')') // format: (a,b,c,d)
1355 {
1356 q = quaternion<T>(a,b,c,d);
1357 }
1358 else // error
1359 {
1360 is.setstate(::std::ios_base::failbit);
1361 }
1362 }
1363 else // error
1364 {
1365 is.setstate(::std::ios_base::failbit);
1366 }
1367 }
1368 else // error
1369 {
1370 is.setstate(::std::ios_base::failbit);
1371 }
1372 }
1373 }
1374 else // error
1375 {
1376 is.setstate(::std::ios_base::failbit);
1377 }
1378 }
1379 }
1380 else // format: a
1381 {
1382 is.putback(ch);
1383
1384 is >> a; // we extract the first component
1385
1386 if (!is.good()) goto finish;
1387
1388 q = quaternion<T>(a);
1389 }
1390
1391 finish:
1392 return(is);
1393 }
1394
1395
1396 template<typename T, typename charT, class traits>
1397 ::std::basic_ostream<charT,traits> & operator << ( ::std::basic_ostream<charT,traits> & os,
1398 quaternion<T> const & q)
1399 {
1400 ::std::basic_ostringstream<charT,traits> s;
1401
1402 s.flags(os.flags());
1403#ifdef BOOST_NO_STD_LOCALE
1404#else
1405 s.imbue(os.getloc());
1406#endif /* BOOST_NO_STD_LOCALE */
1407 s.precision(os.precision());
1408
1409 s << '(' << q.R_component_1() << ','
1410 << q.R_component_2() << ','
1411 << q.R_component_3() << ','
1412 << q.R_component_4() << ')';
1413
1414 return os << s.str();
1415 }
1416
1417
1418 // values
1419
1420 template<typename T>
1421 inline T real(quaternion<T> const & q)
1422 {
1423 return(q.real());
1424 }
1425
1426
1427 template<typename T>
1428 inline quaternion<T> unreal(quaternion<T> const & q)
1429 {
1430 return(q.unreal());
1431 }
1432
1433
1434#define BOOST_QUATERNION_ARRAY_LOADER \
1435 \
1436 T temp[4]; \
1437 \
1438 temp[0] = q.R_component_1(); \
1439 temp[1] = q.R_component_2(); \
1440 temp[2] = q.R_component_3(); \
1441 temp[3] = q.R_component_4();
1442
1443
1444 template<typename T>
1445 inline T sup(quaternion<T> const & q)
1446 {
1447 BOOST_QUATERNION_ARRAY_LOADER
1448
1449 return(arr_max_abs(temp));
1450 }
1451
1452
1453 template<typename T>
1454 inline T l1(quaternion<T> const & q)
1455 {
1456 BOOST_QUATERNION_ARRAY_LOADER
1457
1458 return(sum_abs(temp));
1459 }
1460
1461
1462 template<typename T>
1463 inline T abs(quaternion<T> const & q)
1464 {
1465 using ::std::sqrt;
1466
1467 BOOST_QUATERNION_ARRAY_LOADER
1468
1469 T maxim = (max_abs(temp)); // overflow protection
1470
1471 if (maxim == static_cast<T>(0))
1472 {
1473 return(maxim);
1474 }
1475 else
1476 {
1477 T mixam = static_cast<T>(1)/maxim; // prefer multiplications over divisions
1478
1479 temp[0] *= mixam;
1480 temp[1] *= mixam;
1481 temp[2] *= mixam;
1482 temp[3] *= mixam;
1483
1484 temp[0] *= temp[0];
1485 temp[1] *= temp[1];
1486 temp[2] *= temp[2];
1487 temp[3] *= temp[3];
1488
1489 return(maxim*sqrt(sum(temp)));
1490 }
1491 }
1492
1493
1494#undef BOOST_QUATERNION_ARRAY_LOADER
1495
1496
1497 // Note: This is the Cayley norm, not the Euclidian norm...
1498
1499 template<typename T>
1500 inline T norm(quaternion<T>const & q)
1501 {
1502 return(real(q*conj(q)));
1503 }
1504
1505
1506 template<typename T>
1507 inline quaternion<T> conj(quaternion<T> const & q)
1508 {
1509 return(quaternion<T>( +q.R_component_1(),
1510 -q.R_component_2(),
1511 -q.R_component_3(),
1512 -q.R_component_4()));
1513 }
1514
1515
1516 template<typename T>
1517 inline quaternion<T> spherical( T const & rho,
1518 T const & theta,
1519 T const & phi1,
1520 T const & phi2)
1521 {
1522 using ::std::cos;
1523 using ::std::sin;
1524
1525 //T a = cos(theta)*cos(phi1)*cos(phi2);
1526 //T b = sin(theta)*cos(phi1)*cos(phi2);
1527 //T c = sin(phi1)*cos(phi2);
1528 //T d = sin(phi2);
1529
1530 T courrant = static_cast<T>(1);
1531
1532 T d = sin(phi2);
1533
1534 courrant *= cos(phi2);
1535
1536 T c = sin(phi1)*courrant;
1537
1538 courrant *= cos(phi1);
1539
1540 T b = sin(theta)*courrant;
1541 T a = cos(theta)*courrant;
1542
1543 return(rho*quaternion<T>(a,b,c,d));
1544 }
1545
1546
1547 template<typename T>
1548 inline quaternion<T> semipolar( T const & rho,
1549 T const & alpha,
1550 T const & theta1,
1551 T const & theta2)
1552 {
1553 using ::std::cos;
1554 using ::std::sin;
1555
1556 T a = cos(alpha)*cos(theta1);
1557 T b = cos(alpha)*sin(theta1);
1558 T c = sin(alpha)*cos(theta2);
1559 T d = sin(alpha)*sin(theta2);
1560
1561 return(rho*quaternion<T>(a,b,c,d));
1562 }
1563
1564
1565 template<typename T>
1566 inline quaternion<T> multipolar( T const & rho1,
1567 T const & theta1,
1568 T const & rho2,
1569 T const & theta2)
1570 {
1571 using ::std::cos;
1572 using ::std::sin;
1573
1574 T a = rho1*cos(theta1);
1575 T b = rho1*sin(theta1);
1576 T c = rho2*cos(theta2);
1577 T d = rho2*sin(theta2);
1578
1579 return(quaternion<T>(a,b,c,d));
1580 }
1581
1582
1583 template<typename T>
1584 inline quaternion<T> cylindrospherical( T const & t,
1585 T const & radius,
1586 T const & longitude,
1587 T const & latitude)
1588 {
1589 using ::std::cos;
1590 using ::std::sin;
1591
1592
1593
1594 T b = radius*cos(longitude)*cos(latitude);
1595 T c = radius*sin(longitude)*cos(latitude);
1596 T d = radius*sin(latitude);
1597
1598 return(quaternion<T>(t,b,c,d));
1599 }
1600
1601
1602 template<typename T>
1603 inline quaternion<T> cylindrical(T const & r,
1604 T const & angle,
1605 T const & h1,
1606 T const & h2)
1607 {
1608 using ::std::cos;
1609 using ::std::sin;
1610
1611 T a = r*cos(angle);
1612 T b = r*sin(angle);
1613
1614 return(quaternion<T>(a,b,h1,h2));
1615 }
1616
1617
1618 // transcendentals
1619 // (please see the documentation)
1620
1621
1622 template<typename T>
1623 inline quaternion<T> exp(quaternion<T> const & q)
1624 {
1625 using ::std::exp;
1626 using ::std::cos;
1627
1628 using ::boost::math::sinc_pi;
1629
1630 T u = exp(real(q));
1631
1632 T z = abs(unreal(q));
1633
1634 T w = sinc_pi(z);
1635
1636 return(u*quaternion<T>(cos(z),
1637 w*q.R_component_2(), w*q.R_component_3(),
1638 w*q.R_component_4()));
1639 }
1640
1641
1642 template<typename T>
1643 inline quaternion<T> cos(quaternion<T> const & q)
1644 {
1645 using ::std::sin;
1646 using ::std::cos;
1647 using ::std::cosh;
1648
1649 using ::boost::math::sinhc_pi;
1650
1651 T z = abs(unreal(q));
1652
1653 T w = -sin(q.real())*sinhc_pi(z);
1654
1655 return(quaternion<T>(cos(q.real())*cosh(z),
1656 w*q.R_component_2(), w*q.R_component_3(),
1657 w*q.R_component_4()));
1658 }
1659
1660
1661 template<typename T>
1662 inline quaternion<T> sin(quaternion<T> const & q)
1663 {
1664 using ::std::sin;
1665 using ::std::cos;
1666 using ::std::cosh;
1667
1668 using ::boost::math::sinhc_pi;
1669
1670 T z = abs(unreal(q));
1671
1672 T w = +cos(q.real())*sinhc_pi(z);
1673
1674 return(quaternion<T>(sin(q.real())*cosh(z),
1675 w*q.R_component_2(), w*q.R_component_3(),
1676 w*q.R_component_4()));
1677 }
1678
1679
1680 template<typename T>
1681 inline quaternion<T> tan(quaternion<T> const & q)
1682 {
1683 return(sin(q)/cos(q));
1684 }
1685
1686
1687 template<typename T>
1688 inline quaternion<T> cosh(quaternion<T> const & q)
1689 {
1690 return((exp(+q)+exp(-q))/static_cast<T>(2));
1691 }
1692
1693
1694 template<typename T>
1695 inline quaternion<T> sinh(quaternion<T> const & q)
1696 {
1697 return((exp(+q)-exp(-q))/static_cast<T>(2));
1698 }
1699
1700
1701 template<typename T>
1702 inline quaternion<T> tanh(quaternion<T> const & q)
1703 {
1704 return(sinh(q)/cosh(q));
1705 }
1706
1707
1708 template<typename T>
1709 quaternion<T> pow(quaternion<T> const & q,
1710 int n)
1711 {
1712 if (n > 1)
1713 {
1714 int m = n>>1;
1715
1716 quaternion<T> result = pow(q, m);
1717
1718 result *= result;
1719
1720 if (n != (m<<1))
1721 {
1722 result *= q; // n odd
1723 }
1724
1725 return(result);
1726 }
1727 else if (n == 1)
1728 {
1729 return(q);
1730 }
1731 else if (n == 0)
1732 {
1733 return(quaternion<T>(static_cast<T>(1)));
1734 }
1735 else /* n < 0 */
1736 {
1737 return(pow(quaternion<T>(static_cast<T>(1))/q,-n));
1738 }
1739 }
1740
1741
1742 // helper templates for converting copy constructors (definition)
1743
1744 namespace detail
1745 {
1746
1747 template< typename T,
1748 typename U
1749 >
1750 quaternion<T> quaternion_type_converter(quaternion<U> const & rhs)
1751 {
1752 return(quaternion<T>( static_cast<T>(rhs.R_component_1()),
1753 static_cast<T>(rhs.R_component_2()),
1754 static_cast<T>(rhs.R_component_3()),
1755 static_cast<T>(rhs.R_component_4())));
1756 }
1757 }
1758 }
1759}
1760
1761#endif /* BOOST_QUATERNION_HPP */