//    boost octonion.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_OCTONION_HPP

#define BOOST_OCTONION_HPP

#include <boost/math/quaternion.hpp>

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_OCTONION_ACCESSOR_GENERATOR(type)                      \

            type                        real() const                    \

            {                                                           \

                return(a);                                              \

            }                                                           \

                                                                        \

            octonion<type>                unreal() const                \

            {                                                           \

                return( octonion<type>(static_cast<type>(0),b,c,d,e,f,g,h));   \

            }                                                           \

                                                                        \

            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);                                              \

            }                                                           \

                                                                        \

            type                            R_component_5() const       \

            {                                                           \

                return(e);                                              \

            }                                                           \

                                                                        \

            type                            R_component_6() const       \

            {                                                           \

                return(f);                                              \

            }                                                           \

                                                                        \

            type                            R_component_7() const       \

            {                                                           \

                return(g);                                              \

            }                                                           \

                                                                        \

            type                            R_component_8() const       \

            {                                                           \

                return(h);                                              \

            }                                                           \

                                                                        \

            ::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));                      \

            }                                                           \

                                                                        \

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

            {                                                           \

                return(::std::complex<type>(e,f));                      \

            }                                                           \

                                                                        \

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

            {                                                           \

                return(::std::complex<type>(g,h));                      \

            }                                                           \

                                                                        \

            ::boost::math::quaternion<type>    H_component_1() const    \

            {                                                           \

                return(::boost::math::quaternion<type>(a,b,c,d));       \

            }                                                           \

                                                                        \

            ::boost::math::quaternion<type>    H_component_2() const    \

            {                                                           \

                return(::boost::math::quaternion<type>(e,f,g,h));       \

            }

#define    BOOST_OCTONION_MEMBER_ASSIGNMENT_GENERATOR(type)                                         \

            template<typename X>                                                                    \

            octonion<type> &        operator = (octonion<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());                                  \

                e = static_cast<type>(a_affecter.R_component_5());                                  \

                f = static_cast<type>(a_affecter.R_component_6());                                  \

                g = static_cast<type>(a_affecter.R_component_7());                                  \

                h = static_cast<type>(a_affecter.R_component_8());                                  \

                                                                                                    \

                return(*this);                                                                      \

            }                                                                                       \

                                                                                                    \

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

            {                                                                                       \

                a = a_affecter.a;                                                                   \

                b = a_affecter.b;                                                                   \

                c = a_affecter.c;                                                                   \

                d = a_affecter.d;                                                                   \

                e = a_affecter.e;                                                                   \

                f = a_affecter.f;                                                                   \

                g = a_affecter.g;                                                                   \

                h = a_affecter.h;                                                                   \

                                                                                                    \

                return(*this);                                                                      \

            }                                                                                       \

                                                                                                    \

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

            {                                                                                       \

                a = a_affecter;                                                                     \

                                                                                                    \

                b = c = d = e = f= g = h = static_cast<type>(0);                                    \

                                                                                                    \

                return(*this);                                                                      \

            }                                                                                       \

                                                                                                    \

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

            {                                                                                       \

                a = a_affecter.real();                                                              \

                b = a_affecter.imag();                                                              \

                                                                                                    \

                c = d = e = f = g = h = static_cast<type>(0);                                       \

                                                                                                    \

                return(*this);                                                                      \

            }                                                                                       \

                                                                                                    \

            octonion<type> &        operator = (::boost::math::quaternion<type> const & a_affecter) \

            {                                                                                       \

                a = a_affecter.R_component_1();                                                     \

                b = a_affecter.R_component_2();                                                     \

                c = a_affecter.R_component_3();                                                     \

                d = a_affecter.R_component_4();                                                     \

                                                                                                    \

                e = f = g = h = static_cast<type>(0);                                               \

                                                                                                    \

                return(*this);                                                                      \

            }

#define    BOOST_OCTONION_MEMBER_DATA_GENERATOR(type) \

            type    a;                                \

            type    b;                                \

            type    c;                                \

            type    d;                                \

            type    e;                                \

            type    f;                                \

            type    g;                                \

            type    h;                                \

        template<typename T>

        class octonion

        {

        public:

            typedef T value_type;

            // constructor for O seen as R^8

            // (also default constructor)

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

                                                T const & requested_b = T(),

                                                T const & requested_c = T(),

                                                T const & requested_d = T(),

                                                T const & requested_e = T(),

                                                T const & requested_f = T(),

                                                T const & requested_g = T(),

                                                T const & requested_h = T())

            :   a(requested_a),

                b(requested_b),

                c(requested_c),

                d(requested_d),

                e(requested_e),

                f(requested_f),

                g(requested_g),

                h(requested_h)

            {

                // nothing to do!

            }

            // constructor for H seen as C^4

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

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

                                                ::std::complex<T> const & z2 = ::std::complex<T>(),

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

            :   a(z0.real()),

                b(z0.imag()),

                c(z1.real()),

                d(z1.imag()),

                e(z2.real()),

                f(z2.imag()),

                g(z3.real()),

                h(z3.imag())

            {

                // nothing to do!

            }

            // constructor for O seen as H^2

            explicit                octonion(   ::boost::math::quaternion<T> const & q0,

                                                ::boost::math::quaternion<T> const & q1 = ::boost::math::quaternion<T>())

            :   a(q0.R_component_1()),

                b(q0.R_component_2()),

                c(q0.R_component_3()),

                d(q0.R_component_4()),

                e(q1.R_component_1()),

                f(q1.R_component_2()),

                g(q1.R_component_3()),

                h(q1.R_component_4())

            {

                // nothing to do!

            }

            // UNtemplated copy constructor

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

            // templated copy constructor

            template<typename X>

            explicit                octonion(octonion<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())),

                e(static_cast<T>(a_recopier.R_component_5())),

                f(static_cast<T>(a_recopier.R_component_6())),

                g(static_cast<T>(a_recopier.R_component_7())),

                h(static_cast<T>(a_recopier.R_component_8()))

            {

                // nothing to do!

            }

            // destructor

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

            // accessors

            //

            // Note:    Like complex number, octonions 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 an octonion, 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_OCTONION_ACCESSOR_GENERATOR(T)

            // assignment operators

            BOOST_OCTONION_MEMBER_ASSIGNMENT_GENERATOR(T)

            // other assignment-related operators

            //

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

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

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

            //            octonion multiplication is also *NOT* associative

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

            {

                T    at = a + rhs;    // exception guard

                a = at;

                return(*this);

            }

            octonion<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);

            }

            octonion<T> &            operator += (::boost::math::quaternion<T> const & rhs)

            {

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

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

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

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

                a = at; 

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            template<typename X>

            octonion<T> &            operator += (octonion<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

                T    et = e + static_cast<T>(rhs.R_component_5());    // exception guard

                T    ft = f + static_cast<T>(rhs.R_component_6());    // exception guard

                T    gt = g + static_cast<T>(rhs.R_component_7());    // exception guard

                T    ht = h + static_cast<T>(rhs.R_component_8());    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

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

            {

                T    at = a - rhs;    // exception guard

                a = at;

                return(*this);

            }

            octonion<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);

            }

            octonion<T> &            operator -= (::boost::math::quaternion<T> const & rhs)

            {

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

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

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

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

                a = at; 

                b = bt;

                c = ct;

                d = dt;

                return(*this);

            }

            template<typename X>

            octonion<T> &            operator -= (octonion<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

                T    et = e - static_cast<T>(rhs.R_component_5());    // exception guard

                T    ft = f - static_cast<T>(rhs.R_component_6());    // exception guard

                T    gt = g - static_cast<T>(rhs.R_component_7());    // exception guard

                T    ht = h - static_cast<T>(rhs.R_component_8());    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

            octonion<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

                T    et = e * rhs;    // exception guard

                T    ft = f * rhs;    // exception guard

                T    gt = g * rhs;    // exception guard

                T    ht = h * rhs;    // exception guard

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

            octonion<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;

                T    et = +e*ar+f*br;

                T    ft = -e*br+f*ar;

                T    gt = +g*ar-h*br;

                T    ht = +g*br+h*ar;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

            octonion<T> &            operator *= (::boost::math::quaternion<T> const & rhs)

            {

                T    ar = rhs.R_component_1();

                T    br = rhs.R_component_2();

                T    cr = rhs.R_component_2();

                T    dr = rhs.R_component_2();

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

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

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

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

                T    et = +e*ar+f*br+g*cr+h*dr;

                T    ft = -e*br+f*ar-g*dr+h*cr;

                T    gt = -e*cr+f*dr+g*ar-h*br;

                T    ht = -e*dr-f*cr+g*br+h*ar;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

            template<typename X>

            octonion<T> &            operator *= (octonion<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    er = static_cast<T>(rhs.R_component_5());

                T    fr = static_cast<T>(rhs.R_component_6());

                T    gr = static_cast<T>(rhs.R_component_7());

                T    hr = static_cast<T>(rhs.R_component_8());

                T    at = +a*ar-b*br-c*cr-d*dr-e*er-f*fr-g*gr-h*hr;

                T    bt = +a*br+b*ar+c*dr-d*cr+e*fr-f*er-g*hr+h*gr;