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

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

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

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

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

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

                T    ht = +a*hr-b*gr+c*fr+d*er-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);

            }

            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    denominator = ar*ar+br*br;

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

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

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

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

                T    et = (+e*ar-f*br)/denominator;

                T    ft = (+e*br+f*ar)/denominator;

                T    gt = (+g*ar+h*br)/denominator;

                T    ht = (+g*br+h*ar)/denominator;

                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    denominator = ar*ar+br*br+cr*cr+dr*dr;

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

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

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

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

                T    et = (+e*ar-f*br-g*cr-h*dr)/denominator;

                T    ft = (+e*br+f*ar+g*dr-h*cr)/denominator;

                T    gt = (+e*cr-f*dr+g*ar+h*br)/denominator;

                T    ht = (+e*dr+f*cr-g*br+h*ar)/denominator;

                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    denominator = ar*ar+br*br+cr*cr+dr*dr+er*er+fr*fr+gr*gr+hr*hr;

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

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

                T    ct = (-a*cr+b*dr+c*ar-d*br-e*gr-f*hr+g*er+h*fr)/denominator;

                T    dt = (-a*dr-b*cr+c*br+d*ar-e*hr+f*gr-g*fr+h*er)/denominator;

                T    et = (-a*er+b*fr+c*gr+d*hr+e*ar-f*br-g*cr-h*dr)/denominator;

                T    ft = (-a*fr-b*er+c*hr-d*gr+e*br+f*ar+g*dr-h*cr)/denominator;

                T    gt = (-a*gr-b*hr-c*er+d*fr+e*cr-f*dr+g*ar+h*br)/denominator;

                T    ht = (-a*hr+b*gr-c*fr-d*er+e*dr+f*cr-g*br+h*ar)/denominator;

                a = at;

                b = bt;

                c = ct;

                d = dt;

                e = et;

                f = ft;

                g = gt;

                h = ht;

                return(*this);

            }

        protected:

            BOOST_OCTONION_MEMBER_DATA_GENERATOR(T)

        private:

        };

        // declaration of octonion specialization

        template<>    class octonion<float>;

        template<>    class octonion<double>;

        template<>    class octonion<long double>;

        // helper templates for converting copy constructors (declaration)

        namespace detail

        {

            template<   typename T,

                        typename U

                    >

            octonion<T>    octonion_type_converter(octonion<U> const & rhs);

        }

        // implementation of octonion specialization

#define    BOOST_OCTONION_CONSTRUCTOR_GENERATOR(type)                                                                               \

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

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

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

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

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

            :   a(requested_a),                                                                                                     \

                b(requested_b),                                                                                                     \

                c(requested_c),                                                                                                     \

                d(requested_d),                                                                                                     \

                e(requested_e),                                                                                                     \

                f(requested_f),                                                                                                     \

                g(requested_g),                                                                                                     \

                h(requested_h)                                                                                                      \

            {                                                                                                                       \

            }                                                                                                                       \

                                                                                                                                    \

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

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

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

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

            :   a(z0.real()),                                                                                                       \

                b(z0.imag()),                                                                                                       \

                c(z1.real()),                                                                                                       \

                d(z1.imag()),                                                                                                       \

                e(z2.real()),                                                                                                       \

                f(z2.imag()),                                                                                                       \

                g(z3.real()),                                                                                                       \

                h(z3.imag())                                                                                                        \

            {                                                                                                                       \

            }                                                                                                                       \

                                                                                                                                    \

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

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

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

            {                                                                                                                       \

            }

#define    BOOST_OCTONION_MEMBER_ADD_GENERATOR_1(type)                  \

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

            {                                                           \

                a += rhs;                                               \

                                                                        \

                return(*this);                                          \

            }

#define    BOOST_OCTONION_MEMBER_ADD_GENERATOR_2(type)                                  \

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

            {                                                                           \

                a += rhs.real();                                                        \

                b += rhs.imag();                                                        \

                                                                                        \

                return(*this);                                                          \

            }

#define    BOOST_OCTONION_MEMBER_ADD_GENERATOR_3(type)                                              \

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

            {                                                                                       \

                a += rhs.R_component_1();                                                           \

                b += rhs.R_component_2();                                                           \

                c += rhs.R_component_3();                                                           \

                d += rhs.R_component_4();                                                           \

                                                                                                    \

                return(*this);                                                                      \

            }

#define    BOOST_OCTONION_MEMBER_ADD_GENERATOR_4(type)                          \

            template<typename X>                                                \

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

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

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

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

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

                                                                                \

                return(*this);                                                  \

            }

#define    BOOST_OCTONION_MEMBER_SUB_GENERATOR_1(type)                  \

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

            {                                                           \

                a -= rhs;                                               \

                                                                        \

                return(*this);                                          \

            }

#define    BOOST_OCTONION_MEMBER_SUB_GENERATOR_2(type)                                  \

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

            {                                                                           \

                a -= rhs.real();                                                        \

                b -= rhs.imag();                                                        \

                                                                                        \

                return(*this);                                                          \

            }

#define    BOOST_OCTONION_MEMBER_SUB_GENERATOR_3(type)                                              \

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

            {                                                                                       \

                a -= rhs.R_component_1();                                                           \

                b -= rhs.R_component_2();                                                           \

                c -= rhs.R_component_3();                                                           \

                d -= rhs.R_component_4();                                                           \

                                                                                                    \

                return(*this);                                                                      \

            }

#define    BOOST_OCTONION_MEMBER_SUB_GENERATOR_4(type)                        \

            template<typename X>                                              \

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

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

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

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

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

                                                                              \

                return(*this);                                                \

            }

#define    BOOST_OCTONION_MEMBER_MUL_GENERATOR_1(type)                   \

            octonion<type> &            operator *= (type const & rhs)   \

            {                                                            \

                a *= rhs;                                                \

                b *= rhs;                                                \

                c *= rhs;                                                \

                d *= rhs;                                                \

                e *= rhs;                                                \

                f *= rhs;                                                \

                g *= rhs;                                                \

                h *= rhs;                                                \

                                                                         \

                return(*this);                                           \

            }

#define    BOOST_OCTONION_MEMBER_MUL_GENERATOR_2(type)                                  \

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

            {                                                                           \

                type    ar = rhs.real();                                                \

                type    br = rhs.imag();                                                \

                                                                                        \

                type    at = +a*ar-b*br;                                                \

                type    bt = +a*br+b*ar;                                                \

                type    ct = +c*ar+d*br;                                                \

                type    dt = -c*br+d*ar;                                                \

                type    et = +e*ar+f*br;                                                \

                type    ft = -e*br+f*ar;                                                \

                type    gt = +g*ar-h*br;                                                \

                type    ht = +g*br+h*ar;                                                \

                                                                                        \

                a = at;                                                                 \

                b = bt;                                                                 \

                c = ct;                                                                 \

                d = dt;                                                                 \

                e = et;                                                                 \

                f = ft;                                                                 \

                g = gt;                                                                 \

                h = ht;                                                                 \

                                                                                        \

                return(*this);                                                          \

            }

#define    BOOST_OCTONION_MEMBER_MUL_GENERATOR_3(type)                                                    \

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

            {                                                                                       \

                type    ar = rhs.R_component_1();                                                   \

                type    br = rhs.R_component_2();                                                   \

                type    cr = rhs.R_component_2();                                                   \

                type    dr = rhs.R_component_2();                                                   \

                                                                                                    \

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

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

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

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

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

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

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

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

            }

#define    BOOST_OCTONION_MEMBER_MUL_GENERATOR_4(type)                          \

            template<typename X>                                                \

            octonion<type> &            operator *= (octonion<X> const & rhs)   \

            {                                                                   \

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

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

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

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

                type    er = static_cast<type>(rhs.R_component_5());            \

                type    fr = static_cast<type>(rhs.R_component_6());            \

                type    gr = static_cast<type>(rhs.R_component_7());            \

                type    hr = static_cast<type>(rhs.R_component_8());            \

                                                                                \

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

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

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

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

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

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

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

                type    ht = +a*hr-b*gr+c*fr+d*er-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);                                                  \

            }

// There is quite a lot of repetition in the code below. This is intentional.

// The last conditional block is the normal form, and the others merely

// consist of workarounds for various compiler deficiencies. Hopefuly, when

// more compilers are conformant and we can retire support for those that are

// not, we will be able to remove the clutter. This is makes the situation

// (painfully) explicit.

#define    BOOST_OCTONION_MEMBER_DIV_GENERATOR_1(type)                  \

            octonion<type> &            operator /= (type const & rhs)  \

            {                                                           \

                a /= rhs;                                               \

                b /= rhs;                                               \

                c /= rhs;                                               \

                d /= rhs;                                               \

                                                                        \

                return(*this);                                          \

            }

#if defined(__GNUC__) && (__GNUC__ < 3)

    #define    BOOST_OCTONION_MEMBER_DIV_GENERATOR_2(type)                                              \

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

            {                                                                                           \

                using    ::std::valarray;                                                               \

                                                                                                        \

                valarray<type>    tr(2);                                                                \

                                                                                                        \

                tr[0] = rhs.real();                                                                     \

                tr[1] = rhs.imag();                                                                     \

                                                                                                        \

                type            mixam = (BOOST_GET_VALARRAY(type,static_cast<type>(1)/abs(tr)).max)();  \

                                                                                                        \

                tr *= mixam;                                                                            \

                                                                                                        \

                valarray<type>    tt(8);                                                                \

                                                                                                        \

                tt[0] = +a*tr[0]-b*tr[1];                                                               \

                tt[1] = -a*tr[1]+b*tr[0];                                                               \

                tt[2] = +c*tr[0]-d*tr[1];                                                               \

                tt[3] = +c*tr[1]+d*tr[0];                                                               \

                tt[4] = +e*tr[0]-f*tr[1];                                                               \

                tt[5] = +e*tr[1]+f*tr[0];                                                               \

                tt[6] = +g*tr[0]+h*tr[1];                                                               \

                tt[7] = +g*tr[1]+h*tr[0];                                                               \

                                                                                                        \

                tr *= tr;                                                                               \

                                                                                                        \

                tt *= (mixam/tr.sum());                                                                 \

                                                                                                        \

                a = tt[0];                                                                              \

                b = tt[1];                                                                              \

                c = tt[2];                                                                              \

                d = tt[3];                                                                              \