Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

quaternion.hpp @ 857

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

Rev	Line
[857]	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: