Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

quaternion.hpp @ 857

Revision 857, 91.3 KB checked in by igarcia, 19 years ago (diff)

Line
1	// boost quaternion.hpp header file
2
3	// (C) Copyright Hubert Holin 2001.
4	// Distributed under the Boost Software License, Version 1.0. (See
5	// accompanying file LICENSE_1_0.txt or copy at
6	// 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 = +aar-bbr;
325	T bt = +abr+bar;
326	T ct = +car+dbr;
327	T dt = -cbr+dar;
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 = +aar-bbr-ccr-ddr;
347	T bt = +abr+bar+cdr-dcr; //(abr+arb)+(cdr-crd);
348	T ct = +acr-bdr+car+dbr; //(acr+arc)+(dbr-drb);
349	T dt = +adr+bcr-cbr+dar; //(adr+ard)+(bcr-brc);
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 = arar+brbr;
383
384	T at = (+aar+bbr)/denominator; //(aar+bbr)/denominator;
385	T bt = (-abr+bar)/denominator; //(arb-abr)/denominator;
386	T ct = (+car-dbr)/denominator; //(arc-dbr)/denominator;
387	T dt = (+cbr+dar)/denominator; //(ard+brc)/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 = arar+brbr+crcr+drdr;
407
408	T at = (+aar+bbr+ccr+ddr)/denominator; //(aar+bbr+ccr+ddr)/denominator;
409	T bt = (-abr+bar-cdr+dcr)/denominator; //((arb-abr)+(crd-cdr))/denominator;
410	T ct = (-acr+bdr+car-dbr)/denominator; //((arc-acr)+(drb-dbr))/denominator;
411	T dt = (-adr-bcr+cbr+dar)/denominator; //((ard-adr)+(brc-bcr))/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 = +aar-bbr; \
552	type bt = +abr+bar; \
553	type ct = +car+dbr; \
554	type dt = -cbr+dar; \
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 = +aar-bbr-ccr-ddr; \
574	type bt = +abr+bar+cdr-dcr; \
575	type ct = +acr-bdr+car+dbr; \
576	type dt = +adr+bcr-cbr+dar; \
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] = +atr[0]+btr[1]; \
622	tt[1] = -atr[1]+btr[0]; \
623	tt[2] = +ctr[0]-dtr[1]; \
624	tt[3] = +ctr[1]+dtr[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] = +atr[0]+btr[1]; \
656	tt[1] = -atr[1]+btr[0]; \
657	tt[2] = +ctr[0]-dtr[1]; \
658	tt[3] = +ctr[1]+dtr[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] = +atr[0]+btr[1]; \
689	tt[1] = -atr[1]+btr[0]; \
690	tt[2] = +ctr[0]-dtr[1]; \
691	tt[3] = +ctr[1]+dtr[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] = +atr[0]+btr[1]+ctr[2]+dtr[3]; \
727	tt[1] = -atr[1]+btr[0]-ctr[3]+dtr[2]; \
728	tt[2] = -atr[2]+btr[3]+ctr[0]-dtr[1]; \
729	tt[3] = -atr[3]-btr[2]+ctr[1]+dtr[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] = +atr[0]+btr[1]+ctr[2]+dtr[3]; \
764	tt[1] = -atr[1]+btr[0]-ctr[3]+dtr[2]; \
765	tt[2] = -atr[2]+btr[3]+ctr[0]-dtr[1]; \
766	tt[3] = -atr[3]-btr[2]+ctr[1]+dtr[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] = +atr[0]+btr[1]+ctr[2]+dtr[3]; \
800	tt[1] = -atr[1]+btr[0]-ctr[3]+dtr[2]; \
801	tt[2] = -atr[2]+btr[3]+ctr[0]-dtr[1]; \
802	tt[3] = -atr[3]-btr[2]+ctr[1]+dtr[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 = radiuscos(longitude)cos(latitude);
1752	T c = radiussin(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	wq.R_component_2(), wq.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	wq.R_component_2(), wq.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	wq.R_component_2(), wq.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 */

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: NonGTP/Boost/boost/math/quaternion.hpp @ 857

Download in other formats: