//  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 <boost/config.hpp> // for BOOST_NO_STD_LOCALE

#include <boost/detail/workaround.hpp>

#include <boost/type_traits/is_convertible.hpp>

#include <boost/utility/enable_if.hpp>

#ifndef    BOOST_NO_STD_LOCALE

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

#endif /* BOOST_NO_STD_LOCALE */

#include <complex>

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

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

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

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

#if defined(BOOST_NO_CXX11_NOEXCEPT) || defined(BOOST_NO_CXX11_RVALUE_REFERENCES) || defined(BOOST_NO_SFINAE_EXPR)

#include <boost/type_traits/is_pod.hpp>

#endif

namespace boost

{

   namespace math

   {

      namespace detail {

#if !defined(BOOST_NO_CXX11_NOEXCEPT) && !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && !defined(BOOST_NO_SFINAE_EXPR)

         template <class T>

         struct is_trivial_arithmetic_type_imp

         {

            typedef mpl::bool_<

               noexcept(std::declval<T&>() += std::declval<T>())

               && noexcept(std::declval<T&>() -= std::declval<T>())

               && noexcept(std::declval<T&>() *= std::declval<T>())

               && noexcept(std::declval<T&>() /= std::declval<T>())

            > type;

         };

         template <class T>

         struct is_trivial_arithmetic_type : public is_trivial_arithmetic_type_imp<T>::type {};

#else

         template <class T>

         struct is_trivial_arithmetic_type : public boost::is_pod<T> {};

#endif

      }

#ifndef BOOST_NO_CXX14_CONSTEXPR

      namespace constexpr_detail

      {

         template <class T>

         constexpr void swap(T& a, T& b)

         {

            T t(a);

            a = b;

            b = t;

         }

       }

#endif

       template<typename T>

        class quaternion

        {

        public:

            typedef T value_type;

            // constructor for H seen as R^4

            // (also default constructor)

            BOOST_CONSTEXPR 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

            BOOST_CONSTEXPR 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

            BOOST_CONSTEXPR quaternion(quaternion const & a_recopier)

               : a(a_recopier.R_component_1()),

               b(a_recopier.R_component_2()),

               c(a_recopier.R_component_3()),

               d(a_recopier.R_component_4()) {}

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES

            BOOST_CONSTEXPR quaternion(quaternion && a_recopier)

               : a(std::move(a_recopier.R_component_1())),

               b(std::move(a_recopier.R_component_2())),

               c(std::move(a_recopier.R_component_3())),

               d(std::move(a_recopier.R_component_4())) {}

#endif

            // templated copy constructor

            template<typename X>

            BOOST_CONSTEXPR 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_CONSTEXPR T real() const

            {

               return(a);

            }

            BOOST_CONSTEXPR quaternion<T> unreal() const

            {

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

            }

            BOOST_CONSTEXPR T R_component_1() const

            {

               return(a);

            }

            BOOST_CONSTEXPR T R_component_2() const

            {

               return(b);

            }

            BOOST_CONSTEXPR T R_component_3() const

            {

               return(c);

            }

            BOOST_CONSTEXPR T R_component_4() const

            {

               return(d);

            }

            BOOST_CONSTEXPR ::std::complex<T> C_component_1() const

            {

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

            }

            BOOST_CONSTEXPR ::std::complex<T> C_component_2() const

            {

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

            }

            BOOST_CXX14_CONSTEXPR void swap(quaternion& o)

            {

#ifndef BOOST_NO_CXX14_CONSTEXPR

               using constexpr_detail::swap;

#else

               using std::swap;

#endif

               swap(a, o.a);

               swap(b, o.b);

               swap(c, o.c);

               swap(d, o.d);

            }

            // assignment operators

            template<typename X>

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<X> const  & a_affecter)

            {

               a = static_cast<T>(a_affecter.R_component_1());

               b = static_cast<T>(a_affecter.R_component_2());

               c = static_cast<T>(a_affecter.R_component_3());

               d = static_cast<T>(a_affecter.R_component_4());

               return(*this);

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> const & a_affecter)

            {

               a = a_affecter.a;

               b = a_affecter.b;

               c = a_affecter.c;

               d = a_affecter.d;

               return(*this);

            }

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (quaternion<T> && a_affecter)

            {

               a = std::move(a_affecter.a);

               b = std::move(a_affecter.b);

               c = std::move(a_affecter.c);

               d = std::move(a_affecter.d);

               return(*this);

            }

#endif

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (T const & a_affecter)

            {

               a = a_affecter;

               b = c = d = static_cast<T>(0);

               return(*this);

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator = (::std::complex<T> const & a_affecter)

            {

               a = a_affecter.real();

               b = a_affecter.imag();

               c = d = static_cast<T>(0);

               return(*this);

            }

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

            //

            // Note2:   Each operator comes in 2 forms - one for the simple case where

            //          type T throws no exceptions, and one exception-safe version

            //          for the case where it might.

         private:

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const mpl::true_&)

            {

               a += rhs;

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(T const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a + rhs, b, c, d); // exception guard

               swap(result);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const mpl::true_&)

            {

               a += std::real(rhs);

               b += std::imag(rhs);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(std::complex<T> const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a + std::real(rhs), b + std::imag(rhs), c, d); // exception guard

               swap(result);

               return *this;

            }

            template <class X>

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const mpl::true_&)

            {

               a += rhs.R_component_1();

               b += rhs.R_component_2();

               c += rhs.R_component_3();

               d += rhs.R_component_4();

               return *this;

            }

            template <class X>

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_add(quaternion<X> const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a + rhs.R_component_1(), b + rhs.R_component_2(), c + rhs.R_component_3(), d + rhs.R_component_4()); // exception guard

               swap(result);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const mpl::true_&)

            {

               a -= rhs;

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(T const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a - rhs, b, c, d); // exception guard

               swap(result);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const mpl::true_&)

            {

               a -= std::real(rhs);

               b -= std::imag(rhs);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(std::complex<T> const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a - std::real(rhs), b - std::imag(rhs), c, d); // exception guard

               swap(result);

               return *this;

            }

            template <class X>

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const mpl::true_&)

            {

               a -= rhs.R_component_1();

               b -= rhs.R_component_2();

               c -= rhs.R_component_3();

               d -= rhs.R_component_4();

               return *this;

            }

            template <class X>

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_subtract(quaternion<X> const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a - rhs.R_component_1(), b - rhs.R_component_2(), c - rhs.R_component_3(), d - rhs.R_component_4()); // exception guard

               swap(result);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const mpl::true_&)

            {

               a *= rhs;

               b *= rhs;

               c *= rhs;

               d *= rhs;

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_multiply(T const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a * rhs, b * rhs, c * rhs, d * rhs); // exception guard

               swap(result);

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const mpl::true_&)

            {

               a /= rhs;

               b /= rhs;

               c /= rhs;

               d /= rhs;

               return *this;

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        do_divide(T const & rhs, const mpl::false_&)

            {

               quaternion<T> result(a / rhs, b / rhs, c / rhs, d / rhs); // exception guard

               swap(result);

               return *this;

            }

         public:

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (T const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator += (::std::complex<T> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }

            template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator += (quaternion<X> const & rhs) { return do_add(rhs, detail::is_trivial_arithmetic_type<T>()); }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (T const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator -= (::std::complex<T> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }

            template<typename X> BOOST_CXX14_CONSTEXPR quaternion<T> & operator -= (quaternion<X> const & rhs) { return do_subtract(rhs, detail::is_trivial_arithmetic_type<T>()); }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator *= (T const & rhs) { return do_multiply(rhs, detail::is_trivial_arithmetic_type<T>()); }

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

            {

                T    ar = rhs.real();

                T    br = rhs.imag();

                quaternion<T> result(a*ar - b*br, a*br + b*ar, c*ar + d*br, -c*br+d*ar);

                swap(result);

                return(*this);

            }

            template<typename X>

            BOOST_CXX14_CONSTEXPR 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());

                quaternion<T> result(a*ar - b*br - c*cr - d*dr, a*br + b*ar + c*dr - d*cr, a*cr - b*dr + c*ar + d*br, a*dr + b*cr - c*br + d*ar);

                swap(result);

                return(*this);

            }

            BOOST_CXX14_CONSTEXPR quaternion<T> &        operator /= (T const & rhs) { return do_divide(rhs, detail::is_trivial_arithmetic_type<T>()); }

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

            {

                T    ar = rhs.real();

                T    br = rhs.imag();

                T    denominator = ar*ar+br*br;

                quaternion<T> result((+a*ar + b*br) / denominator, (-a*br + b*ar) / denominator, (+c*ar - d*br) / denominator, (+c*br + d*ar) / denominator);

                swap(result);

                return(*this);

            }

            template<typename X>

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

                quaternion<T> result((+a*ar+b*br+c*cr+d*dr)/denominator, (-a*br+b*ar-c*dr+d*cr)/denominator, (-a*cr+b*dr+c*ar-d*br)/denominator, (-a*dr-b*cr+c*br+d*ar)/denominator);

                swap(result);

                return(*this);

            }

        private:

           T a, b, c, d;

        };

// swap:

template <class T>

BOOST_CXX14_CONSTEXPR void swap(quaternion<T>& a, quaternion<T>& b) { a.swap(b); }

// operator+

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator + (const quaternion<T1>& a, const T2& b)

{

   return quaternion<T1>(static_cast<T1>(a.R_component_1() + b), a.R_component_2(), a.R_component_3(), a.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator + (const T1& a, const quaternion<T2>& b)

{

   return quaternion<T2>(static_cast<T2>(b.R_component_1() + a), b.R_component_2(), b.R_component_3(), b.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator + (const quaternion<T1>& a, const std::complex<T2>& b)

{

   return quaternion<T1>(a.R_component_1() + std::real(b), a.R_component_2() + std::imag(b), a.R_component_3(), a.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator + (const std::complex<T1>& a, const quaternion<T2>& b)

{

   return quaternion<T1>(b.R_component_1() + real(a), b.R_component_2() + imag(a), b.R_component_3(), b.R_component_4());

}

template <class T>

inline BOOST_CONSTEXPR quaternion<T> operator + (const quaternion<T>& a, const quaternion<T>& b)

{

   return quaternion<T>(a.R_component_1() + b.R_component_1(), a.R_component_2() + b.R_component_2(), a.R_component_3() + b.R_component_3(), a.R_component_4() + b.R_component_4());

}

// operator-

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator - (const quaternion<T1>& a, const T2& b)

{

   return quaternion<T1>(static_cast<T1>(a.R_component_1() - b), a.R_component_2(), a.R_component_3(), a.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator - (const T1& a, const quaternion<T2>& b)

{

   return quaternion<T2>(static_cast<T2>(a - b.R_component_1()), -b.R_component_2(), -b.R_component_3(), -b.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator - (const quaternion<T1>& a, const std::complex<T2>& b)

{

   return quaternion<T1>(a.R_component_1() - std::real(b), a.R_component_2() - std::imag(b), a.R_component_3(), a.R_component_4());

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator - (const std::complex<T1>& a, const quaternion<T2>& b)

{

   return quaternion<T1>(real(a) - b.R_component_1(), imag(a) - b.R_component_2(), -b.R_component_3(), -b.R_component_4());

}

template <class T>

inline BOOST_CONSTEXPR quaternion<T> operator - (const quaternion<T>& a, const quaternion<T>& b)

{

   return quaternion<T>(a.R_component_1() - b.R_component_1(), a.R_component_2() - b.R_component_2(), a.R_component_3() - b.R_component_3(), a.R_component_4() - b.R_component_4());

}

// operator*

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator * (const quaternion<T1>& a, const T2& b)

{

   return quaternion<T1>(static_cast<T1>(a.R_component_1() * b), a.R_component_2() * b, a.R_component_3() * b, a.R_component_4() * b);

}

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator * (const T1& a, const quaternion<T2>& b)

{

   return quaternion<T2>(static_cast<T2>(a * b.R_component_1()), a * b.R_component_2(), a * b.R_component_3(), a * b.R_component_4());

}

template <class T1, class T2>

inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator * (const quaternion<T1>& a, const std::complex<T2>& b)

{

   quaternion<T1> result(a);

   result *= b;

   return result;

}

template <class T1, class T2>

inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator * (const std::complex<T1>& a, const quaternion<T2>& b)

{

   quaternion<T1> result(a);

   result *= b;

   return result;

}

template <class T>

inline BOOST_CXX14_CONSTEXPR quaternion<T> operator * (const quaternion<T>& a, const quaternion<T>& b)

{

   quaternion<T> result(a);

   result *= b;

   return result;

}

// operator/

template <class T1, class T2>

inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator / (const quaternion<T1>& a, const T2& b)

{

   return quaternion<T1>(a.R_component_1() / b, a.R_component_2() / b, a.R_component_3() / b, a.R_component_4() / b);

}

template <class T1, class T2>

inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator / (const T1& a, const quaternion<T2>& b)

{

   quaternion<T2> result(a);

   result /= b;

   return result;

}

template <class T1, class T2>

inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T2, T1>::value, quaternion<T1> >::type

operator / (const quaternion<T1>& a, const std::complex<T2>& b)

{

   quaternion<T1> result(a);

   result /= b;

   return result;

}

template <class T1, class T2>

inline BOOST_CXX14_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<T1, T2>::value, quaternion<T2> >::type

operator / (const std::complex<T1>& a, const quaternion<T2>& b)

{

   quaternion<T2> result(a);

   result /= b;

   return result;

}

template <class T>

inline BOOST_CXX14_CONSTEXPR quaternion<T> operator / (const quaternion<T>& a, const quaternion<T>& b)

{

   quaternion<T> result(a);

   result /= b;

   return result;

}

        template<typename T>

        inline BOOST_CONSTEXPR const quaternion<T>&             operator + (quaternion<T> const & q)

        {

            return q;

        }

        template<typename T>

        inline BOOST_CONSTEXPR quaternion<T>                    operator - (quaternion<T> const & q)

        {

            return(quaternion<T>(-q.R_component_1(),-q.R_component_2(),-q.R_component_3(),-q.R_component_4()));

        }

        template<typename R, typename T>

        inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<R, T>::value, bool>::type operator == (R const & lhs, quaternion<T> const & rhs)

        {

            return    (

                        (rhs.R_component_1() == lhs)&&

                        (rhs.R_component_2() == static_cast<T>(0))&&

                        (rhs.R_component_3() == static_cast<T>(0))&&

                        (rhs.R_component_4() == static_cast<T>(0))

                    );

        }

        template<typename T, typename R>

        inline BOOST_CONSTEXPR typename boost::enable_if_c<boost::is_convertible<R, T>::value, bool>::type operator == (quaternion<T> const & lhs, R const & rhs)

        {

           return rhs == lhs;

        }

        template<typename T>

        inline BOOST_CONSTEXPR bool                                operator == (::std::complex<T> const & lhs, quaternion<T> const & rhs)

        {

            return    (

                        (rhs.R_component_1() == lhs.real())&&

                        (rhs.R_component_2() == lhs.imag())&&

                        (rhs.R_component_3() == static_cast<T>(0))&&

                        (rhs.R_component_4() == static_cast<T>(0))

                    );

        }

        template<typename T>

        inline BOOST_CONSTEXPR bool                                operator == (quaternion<T> const & lhs, ::std::complex<T> const & rhs)

        {

           return rhs == lhs;

        }

        template<typename T>

        inline BOOST_CONSTEXPR bool                                operator == (quaternion<T> const & lhs, quaternion<T> const & rhs)

        {

            return    (

                        (rhs.R_component_1() == lhs.R_component_1())&&

                        (rhs.R_component_2() == lhs.R_component_2())&&

                        (rhs.R_component_3() == lhs.R_component_3())&&

                        (rhs.R_component_4() == lhs.R_component_4())

                    );

        }

        template<typename R, typename T> inline BOOST_CONSTEXPR bool operator != (R const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }

        template<typename T, typename R> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, R const & rhs) { return !(lhs == rhs); }

        template<typename T> inline BOOST_CONSTEXPR bool operator != (::std::complex<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }

        template<typename T> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, ::std::complex<T> const & rhs) { return !(lhs == rhs); }

        template<typename T> inline BOOST_CONSTEXPR bool operator != (quaternion<T> const & lhs, quaternion<T> const & rhs) { return !(lhs == rhs); }

        // Note:    we allow the following formats, whith a, b, c, and d reals

        //            a

        //            (a), (a,b), (a,b,c), (a,b,c,d)

        //            (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b),(c,d))

        template<typename T, typename charT, class traits>

        ::std::basic_istream<charT,traits> &    operator >> (    ::std::basic_istream<charT,traits> & is,

                                                                quaternion<T> & q)

        {

#ifdef    BOOST_NO_STD_LOCALE

#else

            const ::std::ctype<charT> & ct = ::std::use_facet< ::std::ctype<charT> >(is.getloc());

#endif /* BOOST_NO_STD_LOCALE */

            T    a = T();

            T    b = T();

            T    c = T();

            T    d = T();

            ::std::complex<T>    u = ::std::complex<T>();

            ::std::complex<T>    v = ::std::complex<T>();

            charT    ch = charT();

            char    cc;

            is >> ch;                                        // get the first lexeme

            if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

            cc = ch;

#else

            cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

            if    (cc == '(')                            // read "(", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d)), ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))

            {

                is >> ch;                                    // get the second lexeme

                if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                cc = ch;

#else

                cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                if    (cc == '(')                        // read "((", possible: ((a)), ((a),c), ((a),(c)), ((a),(c,d)), ((a,b)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))

                {

                    is.putback(ch);

                    is >> u;                                // we extract the first and second components

                    a = u.real();

                    b = u.imag();

                    if    (!is.good())    goto finish;

                    is >> ch;                                // get the next lexeme

                    if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                    cc = ch;

#else

                    cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                    if        (cc == ')')                    // format: ((a)) or ((a,b))

                    {

                        q = quaternion<T>(a,b);

                    }

                    else if    (cc == ',')                // read "((a)," or "((a,b),", possible: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)), ((a,b,),(c,d,))

                    {

                        is >> v;                            // we extract the third and fourth components

                        c = v.real();

                        d = v.imag();

                        if    (!is.good())    goto finish;

                        is >> ch;                                // get the last lexeme

                        if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                        cc = ch;

#else

                        cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                        if    (cc == ')')                    // format: ((a),c), ((a),(c)), ((a),(c,d)), ((a,b),c), ((a,b),(c)) or ((a,b,),(c,d,))

                        {

                            q = quaternion<T>(a,b,c,d);

                        }

                        else                            // error

                        {

                            is.setstate(::std::ios_base::failbit);

                        }

                    }

                    else                                // error

                    {

                        is.setstate(::std::ios_base::failbit);

                    }

                }

                else                                // read "(a", possible: (a), (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))

                {

                    is.putback(ch);

                    is >> a;                                // we extract the first component

                    if    (!is.good())    goto finish;

                    is >> ch;                                // get the third lexeme

                    if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                    cc = ch;

#else

                    cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                    if        (cc == ')')                    // format: (a)

                    {

                        q = quaternion<T>(a);

                    }

                    else if    (cc == ',')                // read "(a,", possible: (a,b), (a,b,c), (a,b,c,d), (a,(c)), (a,(c,d))

                    {

                        is >> ch;                            // get the fourth lexeme

                        if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                        cc = ch;

#else

                        cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                        if    (cc == '(')                // read "(a,(", possible: (a,(c)), (a,(c,d))

                        {

                            is.putback(ch);

                            is >> v;                        // we extract the third and fourth component

                            c = v.real();

                            d = v.imag();

                            if    (!is.good())    goto finish;

                            is >> ch;                        // get the ninth lexeme

                            if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                            cc = ch;

#else

                            cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                            if    (cc == ')')                // format: (a,(c)) or (a,(c,d))

                            {

                                q = quaternion<T>(a,b,c,d);

                            }

                            else                        // error

                            {

                                is.setstate(::std::ios_base::failbit);

                            }

                        }

                        else                        // read "(a,b", possible: (a,b), (a,b,c), (a,b,c,d)

                        {

                            is.putback(ch);

                            is >> b;                        // we extract the second component

                            if    (!is.good())    goto finish;

                            is >> ch;                        // get the fifth lexeme

                            if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                            cc = ch;

#else

                            cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                            if    (cc == ')')                // format: (a,b)

                            {

                                q = quaternion<T>(a,b);

                            }

                            else if    (cc == ',')        // read "(a,b,", possible: (a,b,c), (a,b,c,d)

                            {

                                is >> c;                    // we extract the third component

                                if    (!is.good())    goto finish;

                                is >> ch;                    // get the seventh lexeme

                                if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                                cc = ch;

#else

                                cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                                if        (cc == ')')        // format: (a,b,c)

                                {

                                    q = quaternion<T>(a,b,c);

                                }

                                else if    (cc == ',')    // read "(a,b,c,", possible: (a,b,c,d)

                                {

                                    is >> d;                // we extract the fourth component

                                    if    (!is.good())    goto finish;

                                    is >> ch;                // get the ninth lexeme

                                    if    (!is.good())    goto finish;

#ifdef    BOOST_NO_STD_LOCALE

                                    cc = ch;

#else

                                    cc = ct.narrow(ch, char());

#endif /* BOOST_NO_STD_LOCALE */

                                    if    (cc == ')')        // format: (a,b,c,d)

                                    {

                                        q = quaternion<T>(a,b,c,d);