Mercurial > hg > vamp-build-and-test

comparison DEPENDENCIES/generic/include/boost/math/quaternion.hpp @ 16:2665513ce2d3

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