Ugly, not-to-be-pushed sucking in of all of Boost to get windows to work.
[dyninst.git] / external / boost / math / quaternion.hpp
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>(1));
1891             }
1892             else    /* n < 0 */
1893             {
1894                 return(pow(quaternion<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 */