//  boost quaternion.hpp header file

//  (C) Copyright Hubert Holin 2001.

//  Distributed under the Boost Software License, Version 1.0. (See

//  accompanying file LICENSE_1_0.txt or copy at

//  http://www.boost.org/LICENSE_1_0.txt)

// See http://www.boost.org for updates, documentation, and revision history.

#ifndef BOOST_QUATERNION_HPP

#define BOOST_QUATERNION_HPP

#include <complex>

#include <iosfwd>                                    // for the "<<" and ">>" operators

#include <sstream>                                    // for the "<<" operator

#include <boost/config.hpp> // for BOOST_NO_STD_LOCALE

#include <boost/detail/workaround.hpp>

#ifndef    BOOST_NO_STD_LOCALE

    #include <locale>                                    // for the "<<" operator

#endif /* BOOST_NO_STD_LOCALE */

#include <valarray>

#include <boost/math/special_functions/sinc.hpp>    // for the Sinus cardinal

#include <boost/math/special_functions/sinhc.hpp>    // for the Hyperbolic Sinus cardinal

namespace boost

{

    namespace math

    {

#if     BOOST_WORKAROUND(__GNUC__, < 3)

        // gcc 2.95.x uses expression templates for valarray calculations, but

        // the result is not conforming. We need BOOST_GET_VALARRAY to get an

        // actual valarray result when we need to call a member function

    #define    BOOST_GET_VALARRAY(T,x)    ::std::valarray<T>(x)

        // gcc 2.95.x has an "std::ios" class that is similar to 

        // "std::ios_base", so we just use a #define

    #define    BOOST_IOS_BASE    ::std::ios

        // gcc 2.x ignores function scope using declarations,

        // put them in the scope of the enclosing namespace instead:

        using    ::std::valarray;

        using    ::std::sqrt;

        using    ::std::cos;

        using    ::std::sin;

        using    ::std::exp;

        using    ::std::cosh;

#endif /* BOOST_WORKAROUND(__GNUC__, < 3) */

#define    BOOST_QUATERNION_ACCESSOR_GENERATOR(type)                    \

            type                    real() const                        \

            {                                                           \

                return(a);                                              \

            }                                                           \

                                                                        \

            quaternion<type>        unreal() const                      \

            {                                                           \

                return(quaternion<type>(static_cast<type>(0),b,c,d));   \

            }                                                           \

                                                                        \

            type                    R_component_1() const               \

            {                                                           \

                return(a);                                              \

            }                                                           \

                                                                        \

            type                    R_component_2() const               \

            {                                                           \

                return(b);                                              \

            }                                                           \

                                                                        \

            type                    R_component_3() const               \

            {                                                           \

                return(c);                                              \

            }                                                           \

                                                                        \

            type                    R_component_4() const               \

            {                                                           \

                return(d);                                              \

            }                                                           \

                                                                        \

            ::std::complex<type>    C_component_1() const               \

            {                                                           \

                return(::std::complex<type>(a,b));                      \

            }                                                           \

                                                                        \

            ::std::complex<type>    C_component_2() const               \

            {                                                           \

                return(::std::complex<type>(c,d));                      \

            }

#define    BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(type)                               \

            template<typename X>                                                            \

            quaternion<type> &        operator = (quaternion<X> const  & a_affecter)        \

            {                                                                               \

                a = static_cast<type>(a_affecter.R_component_1());                          \

                b = static_cast<type>(a_affecter.R_component_2());                          \

                c = static_cast<type>(a_affecter.R_component_3());                          \

                d = static_cast<type>(a_affecter.R_component_4());                          \

                                                                                            \

                return(*this);                                                              \

            }                                                                               \

                                                                                            \

            quaternion<type> &        operator = (quaternion<type> const & a_affecter)      \

            {                                                                               \

                a = a_affecter.a;                                                           \

                b = a_affecter.b;                                                           \

                c = a_affecter.c;                                                           \

                d = a_affecter.d;                                                           \

                                                                                            \

                return(*this);                                                              \

            }                                                                               \

                                                                                            \

            quaternion<type> &        operator = (type const & a_affecter)                  \

            {                                                                               \

                a = a_affecter;                                                             \

                                                                                            \

                b = c = d = static_cast<type>(0);                                           \

                                                                                            \

                return(*this);                                                              \

            }                                                                               \

                                                                                            \

            quaternion<type> &        operator = (::std::complex<type> const & a_affecter)  \

            {                                                                               \

                a = a_affecter.real();                                                      \

                b = a_affecter.imag();                                                      \

                                                                                            \

                c = d = static_cast<type>(0);                                               \

                                                                                            \

                return(*this);                                                              \

            }

#define    BOOST_QUATERNION_MEMBER_DATA_GENERATOR(type)       \

            type    a;                                        \

            type    b;                                        \

            type    c;                                        \

            type    d;

        template<typename T>

        class quaternion

        {

        public:

            typedef T value_type;

            // constructor for H seen as R^4

            // (also default constructor)

            explicit            quaternion( T const & requested_a = T(),

                                            T const & requested_b = T(),

                                            T const & requested_c = T(),

                                            T const & requested_d = T())

            :   a(requested_a),

                b(requested_b),

                c(requested_c),

                d(requested_d)

            {

                // nothing to do!

            }

            // constructor for H seen as C^2

            explicit            quaternion( ::std::complex<T> const & z0,

                                            ::std::complex<T> const & z1 = ::std::complex<T>())

            :   a(z0.real()),

                b(z0.imag()),

                c(z1.real()),

                d(z1.imag())

            {

                // nothing to do!

            }

            // UNtemplated copy constructor

            // (this is taken care of by the compiler itself)

            // templated copy constructor

            template<typename X>

            explicit            quaternion(quaternion<X> const & a_recopier)

            :   a(static_cast<T>(a_recopier.R_component_1())),

                b(static_cast<T>(a_recopier.R_component_2())),

                c(static_cast<T>(a_recopier.R_component_3())),

                d(static_cast<T>(a_recopier.R_component_4()))

            {

                // nothing to do!

            }

            // destructor

            // (this is taken care of by the compiler itself)

            // accessors

            //

            // Note:    Like complex number, quaternions do have a meaningful notion of "real part",

            //            but unlike them there is no meaningful notion of "imaginary part".

            //            Instead there is an "unreal part" which itself is a quaternion, and usually

            //            nothing simpler (as opposed to the complex number case).

            //            However, for practicallity, there are accessors for the other components

            //            (these are necessary for the templated copy constructor, for instance).

            BOOST_QUATERNION_ACCESSOR_GENERATOR(T)

            // assignment operators

            BOOST_QUATERNION_MEMBER_ASSIGNMENT_GENERATOR(T)

            // other assignment-related operators

            //

            // NOTE:    Quaternion multiplication is *NOT* commutative;

            //            symbolically, "q *= rhs;" means "q = q * rhs;"

            //            and "q /= rhs;" means "q = q * inverse_of(rhs);"

            quaternion<T> &        operator += (T const & rhs)

            {

                T    at = a + rhs;    // exception guard

                a = at;

                return(*this);

            }

            quaternion<T> &        operator += (::std::complex<T> const & rhs)

            {

                T    at = a + rhs.real();    // exception guard

                T    bt = b + rhs.imag();    // exception guard

                a = at; 

                b = bt;

                return(*this);

            }

            template<typename X>

            quaternion<T> &        operator += (quaternion<X> const & rhs)

            {

                T    at = a + static_cast<T>(rhs.R_component_1());    // exception guard

                T    bt = b + static_cast<T>(rhs.R_component_2());    // exception guard

                T    ct = c + static_cast<T>(rhs.R_component_3());    // exception guard

                T    dt = d + static_cast<T>(rhs.R_component_4());    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            quaternion<T> &        operator -= (T const & rhs)

            {

                T    at = a - rhs;    // exception guard

                a = at;

                return(*this);

            }

            quaternion<T> &        operator -= (::std::complex<T> const & rhs)

            {

                T    at = a - rhs.real();    // exception guard

                T    bt = b - rhs.imag();    // exception guard

                a = at;

                b = bt;

                return(*this);

            }

            template<typename X>

            quaternion<T> &        operator -= (quaternion<X> const & rhs)

            {

                T    at = a - static_cast<T>(rhs.R_component_1());    // exception guard

                T    bt = b - static_cast<T>(rhs.R_component_2());    // exception guard

                T    ct = c - static_cast<T>(rhs.R_component_3());    // exception guard

                T    dt = d - static_cast<T>(rhs.R_component_4());    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            quaternion<T> &        operator *= (T const & rhs)

            {

                T    at = a * rhs;    // exception guard

                T    bt = b * rhs;    // exception guard

                T    ct = c * rhs;    // exception guard

                T    dt = d * rhs;    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            quaternion<T> &        operator *= (::std::complex<T> const & rhs)

            {

                T    ar = rhs.real();

                T    br = rhs.imag();

                T    at = +a*ar-b*br;

                T    bt = +a*br+b*ar;

                T    ct = +c*ar+d*br;

                T    dt = -c*br+d*ar;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            template<typename X>

            quaternion<T> &        operator *= (quaternion<X> const & rhs)

            {

                T    ar = static_cast<T>(rhs.R_component_1());

                T    br = static_cast<T>(rhs.R_component_2());

                T    cr = static_cast<T>(rhs.R_component_3());

                T    dr = static_cast<T>(rhs.R_component_4());

                T    at = +a*ar-b*br-c*cr-d*dr;

                T    bt = +a*br+b*ar+c*dr-d*cr;    //(a*br+ar*b)+(c*dr-cr*d);

                T    ct = +a*cr-b*dr+c*ar+d*br;    //(a*cr+ar*c)+(d*br-dr*b);

                T    dt = +a*dr+b*cr-c*br+d*ar;    //(a*dr+ar*d)+(b*cr-br*c);

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            quaternion<T> &        operator /= (T const & rhs)

            {

                T    at = a / rhs;    // exception guard

                T    bt = b / rhs;    // exception guard

                T    ct = c / rhs;    // exception guard

                T    dt = d / rhs;    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            quaternion<T> &        operator /= (::std::complex<T> const & rhs)

            {

                T    ar = rhs.real();

                T    br = rhs.imag();

                T    denominator = ar*ar+br*br;

                T    at = (+a*ar+b*br)/denominator;    //(a*ar+b*br)/denominator;

                T    bt = (-a*br+b*ar)/denominator;    //(ar*b-a*br)/denominator;

                T    ct = (+c*ar-d*br)/denominator;    //(ar*c-d*br)/denominator;

                T    dt = (+c*br+d*ar)/denominator;    //(ar*d+br*c)/denominator;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            template<typename X>

            quaternion<T> &        operator /= (quaternion<X> const & rhs)

            {

                T    ar = static_cast<T>(rhs.R_component_1());

                T    br = static_cast<T>(rhs.R_component_2());

                T    cr = static_cast<T>(rhs.R_component_3());

                T    dr = static_cast<T>(rhs.R_component_4());

                T    denominator = ar*ar+br*br+cr*cr+dr*dr;

                T    at = (+a*ar+b*br+c*cr+d*dr)/denominator;    //(a*ar+b*br+c*cr+d*dr)/denominator;

                T    bt = (-a*br+b*ar-c*dr+d*cr)/denominator;    //((ar*b-a*br)+(cr*d-c*dr))/denominator;

                T    ct = (-a*cr+b*dr+c*ar-d*br)/denominator;    //((ar*c-a*cr)+(dr*b-d*br))/denominator;

                T    dt = (-a*dr-b*cr+c*br+d*ar)/denominator;    //((ar*d-a*dr)+(br*c-b*cr))/denominator;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

        protected:

            BOOST_QUATERNION_MEMBER_DATA_GENERATOR(T)

        private:

        };

        // declaration of quaternion specialization

        template<>    class quaternion<float>;

        template<>    class quaternion<double>;

        template<>    class quaternion<long double>;

        // helper templates for converting copy constructors (declaration)

        namespace detail

        {

            template<   typename T,

                        typename U

                    >

            quaternion<T>    quaternion_type_converter(quaternion<U> const & rhs);

        }

        // implementation of quaternion specialization

#define    BOOST_QUATERNION_CONSTRUCTOR_GENERATOR(type)                                                 \

            explicit            quaternion( type const & requested_a = static_cast<type>(0),            \

                                            type const & requested_b = static_cast<type>(0),            \

                                            type const & requested_c = static_cast<type>(0),            \

                                            type const & requested_d = static_cast<type>(0))            \

            :   a(requested_a),                                                                         \

                b(requested_b),                                                                         \

                c(requested_c),                                                                         \

                d(requested_d)                                                                          \

            {                                                                                           \

            }                                                                                           \

                                                                                                        \

            explicit            quaternion( ::std::complex<type> const & z0,                            \

                                            ::std::complex<type> const & z1 = ::std::complex<type>())   \

            :   a(z0.real()),                                                                           \

                b(z0.imag()),                                                                           \

                c(z1.real()),                                                                           \

                d(z1.imag())                                                                            \

            {                                                                                           \

            }

#define    BOOST_QUATERNION_MEMBER_ADD_GENERATOR_1(type)             \

            quaternion<type> &        operator += (type const & rhs) \

            {                                                        \

                a += rhs;                                            \

                                                                     \

                return(*this);                                       \

            }

#define    BOOST_QUATERNION_MEMBER_ADD_GENERATOR_2(type)                             \

            quaternion<type> &        operator += (::std::complex<type> const & rhs) \

            {                                                                        \

                a += rhs.real();                                                     \

                b += rhs.imag();                                                     \

                                                                                     \

                return(*this);                                                       \

            }

#define    BOOST_QUATERNION_MEMBER_ADD_GENERATOR_3(type)                      \

            template<typename X>                                              \

            quaternion<type> &        operator += (quaternion<X> const & rhs) \

            {                                                                 \

                a += static_cast<type>(rhs.R_component_1());                  \

                b += static_cast<type>(rhs.R_component_2());                  \

                c += static_cast<type>(rhs.R_component_3());                  \

                d += static_cast<type>(rhs.R_component_4());                  \

                                                                              \

                return(*this);                                                \

            }

#define    BOOST_QUATERNION_MEMBER_SUB_GENERATOR_1(type)             \

            quaternion<type> &        operator -= (type const & rhs) \

            {                                                        \

                a -= rhs;                                            \

                                                                     \

                return(*this);                                       \

            }

#define    BOOST_QUATERNION_MEMBER_SUB_GENERATOR_2(type)                             \

            quaternion<type> &        operator -= (::std::complex<type> const & rhs) \

            {                                                                        \

                a -= rhs.real();                                                     \

                b -= rhs.imag();                                                     \

                                                                                     \

                return(*this);                                                       \

            }

#define    BOOST_QUATERNION_MEMBER_SUB_GENERATOR_3(type)                      \

            template<typename X>                                              \

            quaternion<type> &        operator -= (quaternion<X> const & rhs) \

            {                                                                 \

                a -= static_cast<type>(rhs.R_component_1());                  \

                b -= static_cast<type>(rhs.R_component_2());                  \

                c -= static_cast<type>(rhs.R_component_3());                  \