quaternion.cpp

/******************************************
 * C++ Quaternions
 * Version: 1.0.9
 * Author:  Douglas Wilhelm Harder
 * Date:    2008/03/03
 *
 * Copyright (c) 2006-2008 by Douglas Wilhelm Harder.
 * All rights reserved.
 ******************************************/

#include "Complex.h"
#include "Quaternion.h"
#include "Support.h"
#include <iostream>
#include <cmath>
#include <string>

/******************************************
 * Constructors
 ******************************************/

template <typename T> Quaternion<T>::Quaternion( T re, T im, T jm, T km ):r(re), i(im), j(jm), k(km) {
}

template <typename T> Quaternion<T>::Quaternion( T real ):r(real), i(0.0), j(0.0), k(0.0) {
}


/******************************************
 * Assignment Operator
 ******************************************/

template <typename T> const Quaternion<T> & Quaternion<T>::operator = ( const T & x ) {
    r = x;
    i = 0.0;
    j = 0.0;
    k = 0.0;

    return *this;
}

/******************************************
 * Mutating Arithmetic Operators
 ******************************************/

template <typename T> Quaternion<T> & Quaternion<T>::operator += ( const Quaternion<T> & q ) {
    r += q.r;
    i += q.i;
    j += q.j;
    k += q.k;

    return *this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator -= ( const Quaternion<T> & q ) {
    r -= q.r;
    i -= q.i;
    j -= q.j;
    k -= q.k;

    return *this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator *= ( const Quaternion<T> & q ) {
    T RE = r, I = i, J = j;

    r = RE*q.r - I*q.i - J*q.j - k*q.k;
    i = RE*q.i + I*q.r + J*q.k - k*q.j;
    j = RE*q.j - I*q.k + J*q.r + k*q.i;
    k = RE*q.k + I*q.j - J*q.i + k*q.r;

    return * this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator /= ( const Quaternion<T> & q ) {
    T denom = q.norm();
    T RE = r, I = i, J = j;

    r = ( RE*q.r + I*q.i + J*q.j + k*q.k)/denom;
    i = (-RE*q.i + I*q.r - J*q.k + k*q.j)/denom;
    j = (-RE*q.j + I*q.k + J*q.r - k*q.i)/denom;
    k = (-RE*q.k - I*q.j + J*q.i + k*q.r)/denom;

    return * this;
}


template <typename T> Quaternion<T> & Quaternion<T>::operator += ( T x ) {
    r += x;

    return *this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator -= ( T x ) {
    r -= x;

    return *this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator *= ( T x ) {
    if ( Support<T>::is_inf( x ) && norm() > 0 ) {
        r *= ( r == 0 )?Support<T>::sign( x ):x;
        i *= ( i == 0 )?Support<T>::sign( x ):x;
        j *= ( j == 0 )?Support<T>::sign( x ):x;
        k *= ( k == 0 )?Support<T>::sign( x ):x;
    } else {
        r *= x;
        i *= x;
        j *= x;
        k *= x;
    }

    return * this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator /= ( T x ) {
    if ( x == 0.0 && norm() > 0 ) {
        r /= ( r == 0 ) ? Support<T>::sign( x ):x;
        i /= ( i == 0 ) ? Support<T>::sign( x ):x;
        j /= ( j == 0 ) ? Support<T>::sign( x ):x;
        k /= ( k == 0 ) ? Support<T>::sign( x ):x;
    } else {
        r /= x;
        i /= x;
        j /= x;
        k /= x;
    }

    return * this;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator ++() {
    ++r;

    return *this;
}

template <typename T> Quaternion<T> Quaternion<T>::operator ++( int ) {
    Quaternion copy = *this;
    ++r;

    return copy;
}

template <typename T> Quaternion<T> & Quaternion<T>::operator --() {
    --r;

    return *this;
}

template <typename T> Quaternion<T> Quaternion<T>::operator --( int ) {
    Quaternion copy = *this;
    --r;

    return copy;
}

/******************************************
 * Real-valued Functions
 ******************************************/

template <typename T> T Quaternion<T>::real() const {
    return r;
}

template <typename T> T Quaternion<T>::operator []( int n ) const {
    return reinterpret_cast<const T *>( this )[n];
}

template <typename T> T &Quaternion<T>::operator []( int n ) {
    return reinterpret_cast<T *>( this )[n];
}

template <typename T> T Quaternion<T>::imag_i() const {
    return i;
}

template <typename T> T Quaternion<T>::imag_j() const {
    return j;
}

template <typename T> T Quaternion<T>::imag_k() const {
    return k;
}

template <typename T> T Quaternion<T>::csgn() const {
    return is_zero() ? 0.0 : Support<T>::sign( r );
}

template <typename T> T Quaternion<T>::abs() const {
    return is_inf()? Support<T>::POS_INF : std::sqrt( r*r + i*i + j*j + k*k);
}

template <typename T> T Quaternion<T>::norm() const {
    return is_inf() ? Support<T>::POS_INF : r*r + i*i + j*j + k*k;
}

template <typename T> T Quaternion<T>::abs_imag() const {
    return is_inf() ? Support<T>::POS_INF : std::sqrt( i*i + j*j + k*k );
}

template <typename T> T Quaternion<T>::norm_imag() const {
    return is_inf() ? Support<T>::POS_INF : i*i + j*j + k*k;
}

template <typename T> T Quaternion<T>::arg() const {
    return std::atan2( abs_imag(), r );
}

/******************************************
 * Quaternion<T>-valued Functions
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::imag() const {
    return Quaternion<T>( 0.0, i, j, k );
}

template <typename T> Quaternion<T> Quaternion<T>::conj() const {
    return Quaternion<T>( r, -i, -j, -k );
}

template <typename T> Quaternion<T> Quaternion<T>::operator * () const {
    return Quaternion<T>( r, -i, -j, -k );
}

template <typename T> Quaternion<T> Quaternion<T>::signum() const {
    T absq = abs();

    if ( absq == 0.0 || Support<T>::is_nan( absq ) || Support<T>::is_inf( absq ) ) {
        return *this;
    } else {
        return Quaternion<T>( r/absq, i/absq, j/absq, k/absq );
    }
}

template <typename T> Quaternion<T> Quaternion<T>::sqr() const {
    return Quaternion<T>( r*r - i*i - j*j - k*k, 2*r*i, 2*r*j, 2*r*k );
}

template <typename T> Quaternion<T> Quaternion<T>::sqrt() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).sqrt();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::rotate( const Quaternion<T> & q ) const {
    // the assumption is that |q| = 1
    // in this case, q.inverse() == q.conj()

    T rr = q.r*q.r;
    T ii = q.i*q.i;
    T jj = q.j*q.j;
    T kk = q.k*q.k;

    T ri = q.r*q.i;
    T ij = q.i*q.j;
    T ik = q.i*q.k;

    T rj = q.r*q.j;
    T jk = q.j*q.k;

    T rk = q.r*q.k;

    return Quaternion<T>(
        ( ( r == 0 ) ? r : r*(rr + ii + jj + kk) ),
        i*(rr + ii - jj - kk) + 2*(               j*(-rk + ij) + k*( rj + ik)),
        j*(rr - ii + jj - kk) + 2*(i*( rk + ij)                + k*(-ri + jk)),
        k*(rr - ii - jj + kk) + 2*(i*(-rj + ik) + j*( ri + jk)               )
    );
}

/******************************************
 * Boolean-valued Functions
 ******************************************/


template <typename T> bool Quaternion<T>::is_imaginary() const {
    return ( r == 0.0 );
}

template <typename T> bool Quaternion<T>::is_inf() const {
    return
        ( r == Support<T>::POS_INF ) ||
        ( r == Support<T>::NEG_INF ) ||
        ( i == Support<T>::POS_INF ) ||
        ( i == Support<T>::NEG_INF ) ||
        ( j == Support<T>::POS_INF ) ||
        ( j == Support<T>::NEG_INF ) ||
        ( k == Support<T>::POS_INF ) ||
        ( k == Support<T>::NEG_INF );
}

template <typename T> bool Quaternion<T>::is_nan() const {
    return ( r != r ) || ( i != i ) || ( j != j ) || ( k != k );
}

template <typename T> bool Quaternion<T>::is_neg_inf() const {
    return ( r == Support<T>::NEG_INF ) && ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

template <typename T> bool Quaternion<T>::is_pos_inf() const {
    return ( r == Support<T>::POS_INF ) && ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

template <typename T> bool Quaternion<T>::is_real() const {
    return ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

template <typename T> bool Quaternion<T>::is_real_inf() const {
    return ( r == Support<T>::POS_INF || r == Support<T>::NEG_INF ) && ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

template <typename T> bool Quaternion<T>::is_zero() const {
    return ( r == 0.0 ) && ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

/******************************************
 * Multiplier Function
 ******************************************/

template <typename T> inline Quaternion<T> Quaternion<T>::multiplier( T r, T mltplr, const Quaternion<T> & q ) {
    if ( Support<T>::is_nan( mltplr ) || Support<T>::is_inf( mltplr ) ) {
        if ( q.i == 0 && q.j == 0 && q.k == 0 ) {
            return Quaternion<T>( r, mltplr*q.i, mltplr*q.j, mltplr*q.k );
        } else {
            return Quaternion<T>(
                r,
                (q.i == 0) ? Support<T>::sign( mltplr )*q.i : mltplr*q.i,
                (q.j == 0) ? Support<T>::sign( mltplr )*q.j : mltplr*q.j,
                (q.k == 0) ? Support<T>::sign( mltplr )*q.k : mltplr*q.k
            );
        }
    } else {
        return Quaternion<T>( r, mltplr*q.i, mltplr*q.j, mltplr*q.k );
    }
}

template <typename T> inline Quaternion<T> Quaternion<T>::make_inf( T r, T i ) {
    return Quaternion<T>( r, i, i, i );
}

template <typename T> inline Quaternion<T> Quaternion<T>::make_i( T r, T i ) {
    return Quaternion<T>( r, i, 0.0, 0.0 );
}

/******************************************
 * Exponential and Logarithmic Functions
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::exp() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).exp();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::log() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).log();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::log10() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).log10();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::pow( const Quaternion<T> & q ) const {
    return ( log() * q ).exp();
}

template <typename T> Quaternion<T> Quaternion<T>::pow(T x) const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).pow( x );

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::inverse() const {
    if ( is_zero() ) {
        return Quaternion(
             Support<T>::sign(  r )*Support<T>::POS_INF,
            -Support<T>::sign(  i )*Support<T>::POS_INF,
            -Support<T>::sign(  j )*Support<T>::POS_INF,
            -Support<T>::sign(  k )*Support<T>::POS_INF
        );
    } else if ( is_inf() ) {
        return Quaternion(
             Support<T>::sign(  r )*0.0,
            -Support<T>::sign(  i )*0.0,
            -Support<T>::sign(  j )*0.0,
            -Support<T>::sign(  k )*0.0
        );
    } else if ( is_nan() ) {
        return Quaternion( Support<T>::NaN, Support<T>::NaN, Support<T>::NaN, Support<T>::NaN );
    } else {
        T denom = norm();

        return Quaternion( r/denom, -i/denom, -j/denom, -k/denom );
    }
}

template <typename T> Quaternion<T> Quaternion<T>::cross(const Quaternion<T>& other) const {
    T q1, q2, q3;
    q1 = j * other[3] - k * other[2];
    q2 = k * other[1] - i * other[3];
    q3 = i * other[2] - j * other[1];
    return Quaternion(0, q1, q2, q3);
}

/********************************************************************
 * Trigonometric and Hyperbolic Functions
 *
 * For each function f:H -> H, we define:
 *
 *            ~                       Im(q)     ~
 *  f(q) = Re(f(Re(q) + i|Im(q)|)) + ------- Im(f(Re(q) + i|Im(q)|))
 *                                   |Im(q)|
 *       ~
 * where f:C -> C is the complex equivalent of the
 * function f.
 ********************************************************************/

/**********************************************************
 * Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::sin() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).sin();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Complementary Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::cos() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).cos();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Tangent Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::tan() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).tan();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::sec() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).sec();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Complementary Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::csc() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).csc();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Complementary Tangent Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::cot() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).cot();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::sinh() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).sinh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Complementary Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::cosh() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).cosh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Tangent Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::tanh() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).tanh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::sech() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).sech();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Complementary Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::csch() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).csch();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Hyperbolic Complementary Tangent Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::coth() const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).coth();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

// Real Branch Cut:    (-oo, -1) U (1, oo)

/**********************************************************
 * Inverse Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::asin( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    // Branch Cuts:   (-oo, -1) U (1, oo)

    if ( absIm == 0 ) {
        if ( r > 1 ) {
            // Branch cut (1, oo)

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( Support<T>::PI2, std::log( r + std::sqrt( r*r - 1 ) ) );
            } else {
                return multiplier( Support<T>::PI2, std::log( r + std::sqrt( r*r - 1 ) )/absq, q );
            }
        } else if ( r < -1 ) {
            // Branch cut (-oo, -1)

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( -Support<T>::PI2, std::log( -r + std::sqrt( r*r - 1 ) ) );
            } else {
                return multiplier( -Support<T>::PI2, std::log( -r + std::sqrt( r*r - 1 ) )/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).asin();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );

}

// Real Branch Cut:    (-oo, -1) U (1, oo)

/**********************************************************
 * Inverse Complementary Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acos( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    // Branch Cuts:   (-oo, -1) U (1, oo)

    if ( absIm == 0 ) {
        if ( r > 1 ) {
            // Branch cut (1, oo)

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, -std::log( r + std::sqrt( r*r - 1 ) ) );
            } else {
                return multiplier( 0.0, -std::log( r + std::sqrt( r*r - 1 ) )/absq, q );
            }
        } else if ( r < -1 ) {
            // Branch cut (-oo, -1)

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( Support<T>::PI, -std::log( -r + std::sqrt( r*r - 1 ) ) );
            } else {
                return multiplier( Support<T>::PI, -std::log( -r + std::sqrt( r*r - 1 ) )/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).acos();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

// Complex Branch Cut:    (-ooi, -i] U [i, ooi)

/**********************************************************
 * Inverse Tangent Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::atan() const {
    T absIm = abs_imag();

    if ( r == 0 ) {
        if ( absIm == 1 ) {
                return multiplier( Support<T>::NaN, Support<T>::POS_INF, *this );
        } else if ( absIm > 1 ) {
            // Branch cut [ui, Inf)
            //  - ui is a unit purely-imaginary quaternion

            T p = absIm + 1;
            T m = absIm - 1;

            T mltplr = 0.25*std::log( (p*p)/(m*m) )/absIm;

            return multiplier( Support<T>::sign( r )*Support<T>::PI2, mltplr, *this );
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).atan();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Inverse Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::asec( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    // Branch Cut: (-1, 1)

    if ( absIm == 0 ) {
        if ( r == 0.0 ) {
            return multiplier( Support<T>::POS_INF, Support<T>::POS_INF, q );
        } else if ( r > 0.0 && r < 1.0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, std::log( 1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) ) );
            } else {
                return multiplier( 0.0, std::log( 1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) )/absq, q );
            }
        } else if ( r > -1.0 && r < 0.0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( Support<T>::PI, std::log( -1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) ) );
            } else {
                return multiplier( Support<T>::PI, std::log( -1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) )/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).asec();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

// Real Branch Cut:    (-1, 1)
/**********************************************************
 * Inverse Complementary Secant Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acsc( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    if ( absIm == 0.0 ) {
        if ( r == 0.0 ) {
            return multiplier( Support<T>::POS_INF, Support<T>::POS_INF, q );
        } else if ( r > 0.0 && r <  1.0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( Support<T>::PI2, -std::log( 1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) ) );
            } else {
                return multiplier( Support<T>::PI2, -std::log( 1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) )/absq, q );
            }
        } else if ( r > -1.0 && r < 0.0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( -Support<T>::PI2, -std::log( -1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) ) );
            } else {
                return multiplier( -Support<T>::PI2, -std::log( -1.0/r + std::sqrt( 1.0/(r*r) - 1.0 ) )/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).acsc();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Complementary Tangent Function
 * Complex Branch Cut:    (-i, i)
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acot() const {
    T absIm = abs_imag();

    if ( r == 0 ) {
        if ( absIm == 0 ) {
            return Quaternion<T>( Support<T>::PI2, -i, -j, -k );
        } else if ( absIm < 1 ) {
            // Branch cut [ui, Inf)
            //  - ui is a unit purely-imaginary quaternion

            T p = absIm + 1;
            T m = absIm - 1;

            T mltplr = -0.25*std::log( (p*p)/(m*m) )/absIm;

            return multiplier( Support<T>::PI2, mltplr, *this );
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).acot();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

// Complex Branch Cut:    (-ooi, -i) U (i, ooi)

/**********************************************************
 * Inverse Hyperbolic Sine Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::asinh() const {
    T absIm = abs_imag();

    if ( r == 0 ) {
        if ( absIm > 1 ) {
            return multiplier( Support<T>::sign( r )*std::log( absIm + std::sqrt(absIm*absIm - 1) ), Support<T>::PI2/absIm, *this );
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).asinh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Hyperbolic Complementary Sine Function
 * Real Branch Cut:    (-oo, 1)
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acosh( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    if ( absIm == 0 ) {
        if ( r < -1 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( std::log(-r + std::sqrt(r*r - 1)), Support<T>::PI*Support<T>::sign(i) );
            } else {
                return multiplier( std::log(-r + std::sqrt(r*r - 1)), Support<T>::PI/absq, q );
            }
        } else if ( r == -1 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, Support<T>::PI*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, Support<T>::PI/absq, q );
            }
        } else if ( r < 0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, std::acos(r)*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, std::acos(r)/absq, q );
            }
        } else if ( r == 0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, Support<T>::PI2*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, Support<T>::PI2/absq, q );
            }
        } else if ( r < 1 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, std::acos( r )*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, std::acos(r)/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).acosh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Hyperbolic Tangent Function
 * Real Branch Cut:    (-oo, -1] U [1, oo)
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::atanh( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    if ( absIm == 0 ) {
        if ( r == -1 ) {
            return make_inf( Support<T>::NEG_INF, Support<T>::NaN );
        } else if ( r == 1 ) {
            return make_inf( Support<T>::POS_INF, Support<T>::NaN );
        } else if ( r < -1 || r > 1 ) {
            T p = r + 1;
            T m = r - 1;

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.25*std::log( (p*p)/(m*m) ), -Support<T>::sign(r)*Support<T>::PI2 );
            } else {
                return multiplier( 0.25*std::log( (p*p)/(m*m) ), -Support<T>::sign(r)*Support<T>::PI2/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).atanh();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Hyperbolic Secant Function
 * Real Branch Cut:    (-oo, 0] U (1, oo)
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::asech( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    if ( absIm == 0 ) {
        if ( r < -1 || r > 1 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, -std::acos( 1/r )*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, -std::acos( 1/r )/absq, q );
            }
        } else if ( r == -1 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, Support<T>::PI*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, Support<T>::PI/absq, q );
            }
        } else if ( r < 0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( std::log( -1/r + std::sqrt( 1/(r*r) - 1 ) ), -Support<T>::PI );
            } else {
                return multiplier( std::log( -1/r + std::sqrt( 1/(r*r) - 1 ) ), -Support<T>::PI/absq, q );
            }
        } else if ( r == 0 ) {
            return make_inf( Support<T>::POS_INF, Support<T>::NaN );
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).asech();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Hyperbolic Complementary Secant Function
 * Complex Branch Cut:    (-i, i)
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acsch() const {
    T absIm = abs_imag();

    if ( r == 0 ) {
            if ( absIm == 0 ) {
                return make_inf( Support<T>::NEG_INF, Support<T>::NaN );
            } else if ( absIm < 1 ) {
                return multiplier( Support<T>::sign( r )*std::log( 1/absIm + std::sqrt(1/(absIm*absIm) - 1) ), -Support<T>::PI2/absIm, *this );
            }
    }

    Complex<T> z = Complex<T>( r, absIm ).acsch();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}


/**********************************************************
 * Inverse Hyperbolic Complementary Tangent Function
 * Real Branch Cut:    [-1, 1]
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::acoth( const Quaternion<T> & q ) const {
    T absIm = abs_imag();

    if ( absIm == 0 ) {
        if ( r == -1 ) {
            return make_inf( Support<T>::NEG_INF, Support<T>::NaN );
        } else if ( r == 1 ) {
            return make_inf( Support<T>::POS_INF, Support<T>::NaN );
        } else if ( r == 0 ) {
            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.0, -Support<T>::PI2*Support<T>::sign(i) );
            } else {
                return multiplier( 0.0, -Support<T>::PI2/absq, q );
            }
        } else if ( r > -1 && r < 0 ) {
            T p = r + 1;
            T m = r - 1;

            T absq = q.abs_imag();

            if ( q == I || absq == 0 ) {
                return make_i( 0.25*std::log( (p*p)/(m*m) ), -Support<T>::sign(r)*Support<T>::PI2 );
            } else {
                return multiplier( 0.25*std::log( (p*p)/(m*m) ), -Support<T>::sign(r)*Support<T>::PI2/absq, q );
            }
        }
    }

    Complex<T> z = Complex<T>( r, absIm ).acoth();

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/**********************************************************
 * Bessel J Function
 **********************************************************/

template <typename T> Quaternion<T> Quaternion<T>::bessel_J( int n ) const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).bessel_J( n );

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/******************************************
 * Integer Functions
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::ceil() const {
    return Quaternion<T>(
        std::ceil(r), std::ceil(i), std::ceil(j), std::ceil(k)
    );
}

template <typename T> Quaternion<T> Quaternion<T>::floor() const {
    return Quaternion<T>(
        std::floor(r), std::floor(i), std::floor(j), std::floor(k)
    );
}


/******************************************
 * Horner's Rule
 *
 *   The polynomial is defined by giving the highest
 *   coefficient first:
 *
 *            n - 1         n - 2
 *      v[0]*q      + v[1]*q      + ... + v[n-2]*q + v[n-1]
 *
 *   This is the same as with Matlab.  Because quaternions are
 *   not commutative, this only makes sense if the coefficients
 *   and offsets are real.
 *
 *        Re(q) + i |Imag(q)|
 *
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::horner( T * v, unsigned int n ) const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).horner( v, n );

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

template <typename T> Quaternion<T> Quaternion<T>::horner( T * v, T * c, unsigned int n ) const {
    T absIm = abs_imag();
    Complex<T> z = Complex<T>( r, absIm ).horner( v, c, n );

    T mltplr;

    if ( absIm == 0.0 ) {
        mltplr = z.imag_i();
    } else {
        mltplr = z.imag_i() / absIm;
    }

    return multiplier( z.real(), mltplr, *this );
}

/******************************************
 * Random Factories
 ******************************************/


template <typename T> Quaternion<T> Quaternion<T>::random() {
    return Quaternion<T>(
        (static_cast<T>( rand() ))/RAND_MAX,
        (static_cast<T>( rand() ))/RAND_MAX,
        (static_cast<T>( rand() ))/RAND_MAX,
        (static_cast<T>( rand() ))/RAND_MAX
    );
}

template <typename T> Quaternion<T> Quaternion<T>::random_imag() {
    return Quaternion<T>( 0.0, (static_cast<T>( rand() ))/RAND_MAX, (static_cast<T>( rand() ))/RAND_MAX, (static_cast<T>( rand() ))/RAND_MAX );
}

template <typename T> Quaternion<T> Quaternion<T>::random_real() {
    return Quaternion<T>( (static_cast<T>( rand() ))/RAND_MAX, 0.0, 0.0, 0.0 );
}

/******************************************
 * Binary Arithmetic Operators
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::operator + ( const Quaternion<T> & z ) const {
    return Quaternion<T>( r + z.r, i + z.i, j + z.j, k + z.k );
}

template <typename T> Quaternion<T> Quaternion<T>::operator + ( T x ) const {
    return Quaternion<T>( r + x, i, j, k );
}

template <typename T> Quaternion<T> operator + ( T x, const Quaternion<T> & z ) {
    return Quaternion<T>( x + z.real(), z.imag_i(), z.imag_j(), z.imag_k() );
}

template <typename T> Quaternion<T> operator + ( long x, const Quaternion<T> & z ) {
    return Quaternion<T>( static_cast<T>( x ) + z.real(), z.imag_i(), z.imag_j(), z.imag_k() );
}

template <typename T> Quaternion<T> Quaternion<T>::operator - ( const Quaternion<T> & z ) const {
    return Quaternion<T>( r - z.r, i - z.i, j - z.j, k - z.k );
}

template <typename T> Quaternion<T> Quaternion<T>::operator - ( T x ) const {
    return Quaternion<T>( r - x, i, j, k );
}

template <typename T> Quaternion<T> operator - ( T x, const Quaternion<T> & z ) {
    return Quaternion<T>( x - z.real(), -z.imag_i(), -z.imag_j(), -z.imag_k() );
}

template <typename T> Quaternion<T> operator - ( long x, const Quaternion<T> & z ) {
    return Quaternion<T>( static_cast<T>( x ) - z.real(), -z.imag_i(), -z.imag_j(), -z.imag_k() );
}

template <typename T> Quaternion<T> Quaternion<T>::operator * ( const Quaternion<T> & q ) const {
    return Quaternion<T>(
        r*q.r - i*q.i - j*q.j - k*q.k,
        r*q.i + i*q.r + j*q.k - k*q.j,
        r*q.j - i*q.k + j*q.r + k*q.i,
        r*q.k + i*q.j - j*q.i + k*q.r
    );
}

template <typename T> Quaternion<T> Quaternion<T>::operator * ( T x ) const {
    if ( Support<T>::is_inf( x ) && norm() > 0 ) {
        return Quaternion<T>(
            ( (  r == 0 )?Support<T>::sign(x):x )*r,
            ( (  i == 0 )?Support<T>::sign(x):x )*i,
            ( (  j == 0 )?Support<T>::sign(x):x )*j,
            ( (  k == 0 )?Support<T>::sign(x):x )*k
        );
    } else {
        return Quaternion<T>( x*r, x*i, x*j, x*k );
    }
}

template <typename T> Quaternion<T> operator * ( T x, const Quaternion<T> & q ) {
    return q.operator * ( x );
}

template <typename T> Quaternion<T> operator * ( long x, const Quaternion<T> & q ) {
    return q.operator * ( static_cast<T>( x ) );
}

template <typename T> Quaternion<T> Quaternion<T>::operator / ( const Quaternion<T> & q ) const {
    T denom = q.norm();

    return Quaternion<T>(
        ( r*q.r + i*q.i + j*q.j + k*q.k)/denom,
        (-r*q.i + i*q.r - j*q.k + k*q.j)/denom,
        (-r*q.j + i*q.k + j*q.r - k*q.i)/denom,
        (-r*q.k - i*q.j + j*q.i + k*q.r)/denom
    );
}

template <typename T> Quaternion<T> Quaternion<T>::operator / ( T x ) const {
    if ( x == 0.0 && norm() > 0 ) {
        return Quaternion<T>(
            r / ( ( r == 0 )?Support<T>::sign( x ):x ),
            i / ( ( i == 0 )?Support<T>::sign( x ):x ),
            j / ( ( j == 0 )?Support<T>::sign( x ):x ),
            k / ( ( k == 0 )?Support<T>::sign( x ):x )
        );
    } else {
        return Quaternion<T>( r/x, i/x, j/x, k/x );
    }
}

template <typename T> Quaternion<T> operator / ( T x, const Quaternion<T> & q ) {
    T mltplr = x/q.norm();

    return Quaternion<T>( mltplr*q.real(), -mltplr*q.imag_i(), -mltplr*q.imag_j(), -mltplr*q.imag_k() );
}

template <typename T> Quaternion<T> operator / ( long x, const Quaternion<T> & q ) {
    T mltplr = static_cast<T>( x )/q.norm();

    return Quaternion<T>( mltplr*q.real(), -mltplr*q.imag_i(), -mltplr*q.imag_j(), -mltplr*q.imag_k() );
}

/******************************************
 * Unary Arithmetic Operators
 ******************************************/

template <typename T> Quaternion<T> Quaternion<T>::operator - () const {
    return Quaternion<T>( -r, -i, -j, -k );
}

/******************************************
 * Binary Boolean Operators
 ******************************************/

template <typename T> bool Quaternion<T>::operator == ( const Quaternion<T> & q ) const {
    return ( r == q.r ) && ( i == q.i ) && ( j == q.j ) && ( k == q.k );
}

template <typename T> bool Quaternion<T>::operator == ( T x ) const {
    return ( r == x ) && ( i == 0.0 ) && ( j == 0.0 ) && ( k == 0.0 );
}

template <typename T> bool operator == ( T x, const Quaternion<T> & q ) {
    return q.operator == ( x );
}

template <typename T> bool operator == ( long x, const Quaternion<T> & q ) {
    return q.operator == ( static_cast<T>( x ) );
}

template <typename T> bool Quaternion<T>::operator != ( const Quaternion<T> & q ) const {
    return ( r != q.r ) || ( i != q.i ) || ( j != q.j ) || ( k != q.k );
}

template <typename T> bool Quaternion<T>::operator != ( T x ) const {
    return ( r != x ) || ( i != 0.0 ) || ( j != 0.0 ) || ( k != 0.0 );
}

template <typename T> bool operator != ( T x, const Quaternion<T> & q ) {
    return q.operator != ( x );
}

template <typename T> bool operator != ( long x, const Quaternion<T> & q ) {
    return q.operator != ( static_cast<T>( x ) );
}

/******************************************
 * IO Stream Operators
 ******************************************/

template <typename T> std::ostream & operator << ( std::ostream & out, const Quaternion<T> & z ) {
    Support<T>::print_real( z.real(), out );
    Support<T>::print_imaginary( z.imag_i(), 'i', out );
    Support<T>::print_imaginary( z.imag_j(), 'j', out );
    Support<T>::print_imaginary( z.imag_k(), 'k', out );

    return out;
}

/******************************************
 * ************************************** *
 * *                                    * *
 * *      Procedural Functions          * *
 * *                                    * *
 * ************************************** *
 ******************************************/

/******************************************
 * Real-valued Functions
 ******************************************/

template <typename T> T real( const Quaternion<T> & q ) {
    return q.real();
}

template <typename T> T imag_i( const Quaternion<T> & q ) {
    return q.imag_i();
}

template <typename T> T imag_j( const Quaternion<T> & q ) {
    return q.imag_j();
}

template <typename T> T imag_k( const Quaternion<T> & q ) {
    return q.imag_k();
}

template <typename T> T csgn( const Quaternion<T> & q ) {
    return q.csgn();
}

template <typename T> T abs( const Quaternion<T> & q ) {
    return q.abs();
}

template <typename T> T norm( const Quaternion<T> & q ) {
    return q.norm();
}

template <typename T> T abs_imag( const Quaternion<T> & z ) {
    return z.abs_imag();
}

template <typename T> T norm_imag( const Quaternion<T> & z ) {
    return z.norm_imag();
}

/******************************************
 * Quaternion-valued Functions
 ******************************************/

template <typename T> Quaternion<T> imag( const Quaternion<T> & q ) {
    return q.imag();
}

template <typename T> Quaternion<T> conj( const Quaternion<T> & q ) {
    return q.conj();
}

template <typename T> Quaternion<T> signum( const Quaternion<T> & q ) {
    return q.signum();
}

template <typename T> Quaternion<T> sqr( const Quaternion<T> & q ) {
    return q.sqr();
}

template <typename T> Quaternion<T> sqrt( const Quaternion<T> & q ) {
    return q.sqrt();
}

template <typename T> Quaternion<T> rotate( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.rotate( p );
}

/******************************************
 * Exponential and Logarithmic Functions
 ******************************************/

template <typename T> Quaternion<T> exp( const Quaternion<T> & q ) {
    return q.exp();
}

template <typename T> Quaternion<T> log( const Quaternion<T> & q ) {
    return q.log();
}

template <typename T> Quaternion<T> log10( const Quaternion<T> & q ) {
    return q.log10();
}

template <typename T> Quaternion<T> pow( const Quaternion<T> & q, const Quaternion<T> & w ) {
    return q.pow( w );
}

template <typename T> Quaternion<T> pow( const Quaternion<T> & q, T x ) {
    return q.pow( x );
}

template <typename T> Quaternion<T> inverse( const Quaternion<T> & q ) {
    return q.inverse();
}

/******************************************
 * Trigonometric and Hyperbolic Functions
 ******************************************/

template <typename T> Quaternion<T> sin( const Quaternion<T> & q ) {
    return q.sin();
}

template <typename T> Quaternion<T> cos( const Quaternion<T> & q ) {
    return q.cos();
}

template <typename T> Quaternion<T> tan( const Quaternion<T> & q ) {
    return q.tan();
}

template <typename T> Quaternion<T> sec( const Quaternion<T> & q ) {
    return q.sec();
}

template <typename T> Quaternion<T> csc( const Quaternion<T> & q ) {
    return q.csc();
}

template <typename T> Quaternion<T> cot( const Quaternion<T> & q ) {
    return q.cot();
}

template <typename T> Quaternion<T> sinh( const Quaternion<T> & q ) {
    return q.sinh();
}

template <typename T> Quaternion<T> cosh( const Quaternion<T> & q ) {
    return q.cosh();
}

template <typename T> Quaternion<T> tanh( const Quaternion<T> & q ) {
    return q.tanh();
}

template <typename T> Quaternion<T> sech( const Quaternion<T> & q ) {
    return q.sech();
}

template <typename T> Quaternion<T> csch( const Quaternion<T> & q ) {
    return q.csch();
}

template <typename T> Quaternion<T> coth( const Quaternion<T> & q ) {
    return q.coth();
}

template <typename T> Quaternion<T> asin( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.asin( p );
}

template <typename T> Quaternion<T> acos( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.acos( p );
}

template <typename T> Quaternion<T> atan( const Quaternion<T> & q ) {
    return q.atan();
}

template <typename T> Quaternion<T> asec( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.asec( p );
}

template <typename T> Quaternion<T> acsc( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return p.acsc( p );
}

template <typename T> Quaternion<T> acot( const Quaternion<T> & q ) {
    return q.acot();
}

template <typename T> Quaternion<T> asinh( const Quaternion<T> & q ) {
    return q.asinh();
}

template <typename T> Quaternion<T> acosh( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.acosh( p );
}

template <typename T> Quaternion<T> atanh( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.atanh( p );
}

template <typename T> Quaternion<T> asech( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.asech( p );
}

template <typename T> Quaternion<T> acsch( const Quaternion<T> & q ) {
    return q.acsch();
}

template <typename T> Quaternion<T> acoth( const Quaternion<T> & q, const Quaternion<T> & p ) {
    return q.acoth( p );
}

template <typename T> Quaternion<T> bessel_J( int n, const Quaternion<T> & z ) {
    return z.bessel_J( n );
}

/******************************************
 * Integer Functions
 ******************************************/

template <typename T> Quaternion<T> floor( const Quaternion<T> & q ) {
    return q.floor();
}

template <typename T> Quaternion<T> ceil( const Quaternion<T> & q ) {
    return q.ceil();
}

/******************************************
 * Horner's Rule
 ******************************************/

template <typename T> Quaternion<T> horner( const Quaternion<T> & q, T * v, unsigned int n ) {
    return q.horner( v, n );
}

template <typename T> Quaternion<T> horner( const Quaternion<T> & q, T * v, T * c, unsigned int n ) {
    return q.horner( v, c, n );
}

/**************************************************
 * ********************************************** *
 * *                                            * *
 * *    Double-precision Floating-point         * *
 * *    Instance of Template                    * *
 * *                                            * *
 * ********************************************** *
 **************************************************/

template class Quaternion<double>;
template std::ostream & operator << ( std::ostream &, const Quaternion<double> & );

template Quaternion<double> operator + ( double, const Quaternion<double> & );
template Quaternion<double> operator + ( long, const Quaternion<double> & );
template Quaternion<double> operator - ( double, const Quaternion<double> & );
template Quaternion<double> operator - ( long, const Quaternion<double> & );
template Quaternion<double> operator * ( double, const Quaternion<double> & );
template Quaternion<double> operator * ( long, const Quaternion<double> & );
template Quaternion<double> operator / ( double, const Quaternion<double> & );
template Quaternion<double> operator / ( long, const Quaternion<double> & );

template bool operator == ( double, const Quaternion<double> & );
template bool operator == ( long, const Quaternion<double> & );
template bool operator != ( double, const Quaternion<double> & );
template bool operator != ( long, const Quaternion<double> & );

template <> const Quaternion<double> Quaternion<double>::ZERO = Quaternion<double>( 0, 0, 0, 0 );
template <> const Quaternion<double>  Quaternion<double>::ONE = Quaternion<double>( 1, 0, 0, 0 );
template <> const Quaternion<double>    Quaternion<double>::I = Quaternion<double>( 0, 1, 0, 0 );
template <> const Quaternion<double>    Quaternion<double>::J = Quaternion<double>( 0, 0, 1, 0 );
template <> const Quaternion<double>    Quaternion<double>::K = Quaternion<double>( 0, 0, 0, 1 );

template <> const Quaternion<double> Quaternion<double>::UNITS[4] = {
    Quaternion<double>::ONE,
    Quaternion<double>::I,
    Quaternion<double>::J,
    Quaternion<double>::K
};

template double real( const Quaternion<double> & );
template double imag_i( const Quaternion<double> & );
template double imag_j( const Quaternion<double> & );
template double imag_k( const Quaternion<double> & );
template double csgn( const Quaternion<double> & );
template double abs( const Quaternion<double> & );
template double norm( const Quaternion<double> & );
template double abs_imag( const Quaternion<double> & );
template double norm_imag( const Quaternion<double> & );
template Quaternion<double> imag( const Quaternion<double> & );
template Quaternion<double> conj( const Quaternion<double> & );
template Quaternion<double> signum( const Quaternion<double> & );
template Quaternion<double> sqr( const Quaternion<double> & );
template Quaternion<double> sqrt( const Quaternion<double> & );
template Quaternion<double> rotate( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> exp( const Quaternion<double> & );
template Quaternion<double> log( const Quaternion<double> & );
template Quaternion<double> log10( const Quaternion<double> & );
template Quaternion<double> pow( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> pow( const Quaternion<double> &, double );
template Quaternion<double> inverse( const Quaternion<double> & );
template Quaternion<double> sin( const Quaternion<double> & );
template Quaternion<double> cos( const Quaternion<double> & );
template Quaternion<double> tan( const Quaternion<double> & );
template Quaternion<double> sec( const Quaternion<double> & );
template Quaternion<double> csc( const Quaternion<double> & );
template Quaternion<double> cot( const Quaternion<double> & );
template Quaternion<double> sinh( const Quaternion<double> & );
template Quaternion<double> cosh( const Quaternion<double> & );
template Quaternion<double> tanh( const Quaternion<double> & );
template Quaternion<double> sech( const Quaternion<double> & );
template Quaternion<double> csch( const Quaternion<double> & );
template Quaternion<double> coth( const Quaternion<double> & );
template Quaternion<double> asin( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> acos( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> atan( const Quaternion<double> & );
template Quaternion<double> asec( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> acsc( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> acot( const Quaternion<double> & );
template Quaternion<double> asinh( const Quaternion<double> & );
template Quaternion<double> acosh( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> atanh( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> asech( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> acsch( const Quaternion<double> & );
template Quaternion<double> acoth( const Quaternion<double> &, const Quaternion<double> & );
template Quaternion<double> bessel_J( int, const Quaternion<double> & );
template Quaternion<double> floor( const Quaternion<double> & );
template Quaternion<double> ceil( const Quaternion<double> & );
template Quaternion<double> horner( const Quaternion<double> &, double *, unsigned int );
template Quaternion<double> horner( const Quaternion<double> &, double *, double *, unsigned int );

/**************************************************
 * ********************************************** *
 * *                                            * *
 * *    Floating-point Instance of Template     * *
 * *                                            * *
 * ********************************************** *
 **************************************************/

template class Quaternion<float>;
template std::ostream & operator << ( std::ostream &, const Quaternion<float> & );

template Quaternion<float> operator + ( float, const Quaternion<float> & );
template Quaternion<float> operator + ( long, const Quaternion<float> & );
template Quaternion<float> operator - ( float, const Quaternion<float> & );
template Quaternion<float> operator - ( long, const Quaternion<float> & );
template Quaternion<float> operator * ( float, const Quaternion<float> & );
template Quaternion<float> operator * ( long, const Quaternion<float> & );
template Quaternion<float> operator / ( float, const Quaternion<float> & );
template Quaternion<float> operator / ( long, const Quaternion<float> & );

template bool operator == ( float, const Quaternion<float> & );
template bool operator == ( long, const Quaternion<float> & );
template bool operator != ( float, const Quaternion<float> & );
template bool operator != ( long, const Quaternion<float> & );

template <> const Quaternion<float> Quaternion<float>::ZERO = Quaternion<float>( 0, 0, 0, 0 );
template <> const Quaternion<float>  Quaternion<float>::ONE = Quaternion<float>( 1, 0, 0, 0 );
template <> const Quaternion<float>    Quaternion<float>::I = Quaternion<float>( 0, 1, 0, 0 );
template <> const Quaternion<float>    Quaternion<float>::J = Quaternion<float>( 0, 0, 1, 0 );
template <> const Quaternion<float>    Quaternion<float>::K = Quaternion<float>( 0, 0, 0, 1 );

template <> const Quaternion<float> Quaternion<float>::UNITS[4] = {
    Quaternion<float>::ONE,
    Quaternion<float>::I,
    Quaternion<float>::J,
    Quaternion<float>::K
};

template float real( const Quaternion<float> & );
template float imag_i( const Quaternion<float> & );
template float imag_j( const Quaternion<float> & );
template float imag_k( const Quaternion<float> & );
template float csgn( const Quaternion<float> & );
template float abs( const Quaternion<float> & );
template float norm( const Quaternion<float> & );
template float abs_imag( const Quaternion<float> & );
template float norm_imag( const Quaternion<float> & );
template Quaternion<float> imag( const Quaternion<float> & );
template Quaternion<float> conj( const Quaternion<float> & );
template Quaternion<float> signum( const Quaternion<float> & );
template Quaternion<float> sqr( const Quaternion<float> & );
template Quaternion<float> sqrt( const Quaternion<float> & );
template Quaternion<float> rotate( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> exp( const Quaternion<float> & );
template Quaternion<float> log( const Quaternion<float> & );
template Quaternion<float> log10( const Quaternion<float> & );
template Quaternion<float> pow( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> pow( const Quaternion<float> &, float );
template Quaternion<float> inverse( const Quaternion<float> & );
template Quaternion<float> sin( const Quaternion<float> & );
template Quaternion<float> cos( const Quaternion<float> & );
template Quaternion<float> tan( const Quaternion<float> & );
template Quaternion<float> sec( const Quaternion<float> & );
template Quaternion<float> csc( const Quaternion<float> & );
template Quaternion<float> cot( const Quaternion<float> & );
template Quaternion<float> sinh( const Quaternion<float> & );
template Quaternion<float> cosh( const Quaternion<float> & );
template Quaternion<float> tanh( const Quaternion<float> & );
template Quaternion<float> sech( const Quaternion<float> & );
template Quaternion<float> csch( const Quaternion<float> & );
template Quaternion<float> coth( const Quaternion<float> & );
template Quaternion<float> asin( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> acos( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> atan( const Quaternion<float> & );
template Quaternion<float> asec( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> acsc( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> acot( const Quaternion<float> & );
template Quaternion<float> asinh( const Quaternion<float> & );
template Quaternion<float> acosh( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> atanh( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> asech( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> acsch( const Quaternion<float> & );
template Quaternion<float> acoth( const Quaternion<float> &, const Quaternion<float> & );
template Quaternion<float> bessel_J( int, const Quaternion<float> & );
template Quaternion<float> floor( const Quaternion<float> & );
template Quaternion<float> ceil( const Quaternion<float> & );
template Quaternion<float> horner( const Quaternion<float> &, float *, unsigned int );
template Quaternion<float> horner( const Quaternion<float> &, float *, float *, unsigned int );

/************************************************************************
 * ******************************************************************** *
 * *                                                                  * *
 * *    Long Double-precision Floating-point Instance of Template     * *
 * *                                                                  * *
 * ******************************************************************** *
 ************************************************************************/

template class Quaternion<long double>;
template std::ostream & operator << ( std::ostream &, const Quaternion<long double> & );

template Quaternion<long double> operator + ( long double, const Quaternion<long double> & );
template Quaternion<long double> operator + ( long, const Quaternion<long double> & );
template Quaternion<long double> operator - ( long double, const Quaternion<long double> & );
template Quaternion<long double> operator - ( long, const Quaternion<long double> & );
template Quaternion<long double> operator * ( long double, const Quaternion<long double> & );
template Quaternion<long double> operator * ( long, const Quaternion<long double> & );
template Quaternion<long double> operator / ( long double, const Quaternion<long double> & );
template Quaternion<long double> operator / ( long, const Quaternion<long double> & );

template bool operator == ( long double, const Quaternion<long double> & );
template bool operator == ( long, const Quaternion<long double> & );
template bool operator != ( long double, const Quaternion<long double> & );
template bool operator != ( long, const Quaternion<long double> & );

template <> const Quaternion<long double> Quaternion<long double>::ZERO = Quaternion<long double>( 0, 0, 0, 0 );
template <> const Quaternion<long double>  Quaternion<long double>::ONE = Quaternion<long double>( 1, 0, 0, 0 );
template <> const Quaternion<long double>    Quaternion<long double>::I = Quaternion<long double>( 0, 1, 0, 0 );
template <> const Quaternion<long double>    Quaternion<long double>::J = Quaternion<long double>( 0, 0, 1, 0 );
template <> const Quaternion<long double>    Quaternion<long double>::K = Quaternion<long double>( 0, 0, 0, 1 );

template <> const Quaternion<long double> Quaternion<long double>::UNITS[4] = {
    Quaternion<long double>::ONE,
    Quaternion<long double>::I,
    Quaternion<long double>::J,
    Quaternion<long double>::K
};

template long double real( const Quaternion<long double> & );
template long double imag_i( const Quaternion<long double> & );
template long double imag_j( const Quaternion<long double> & );
template long double imag_k( const Quaternion<long double> & );
template long double csgn( const Quaternion<long double> & );
template long double abs( const Quaternion<long double> & );
template long double norm( const Quaternion<long double> & );
template long double abs_imag( const Quaternion<long double> & );
template long double norm_imag( const Quaternion<long double> & );
template Quaternion<long double> imag( const Quaternion<long double> & );
template Quaternion<long double> conj( const Quaternion<long double> & );
template Quaternion<long double> signum( const Quaternion<long double> & );
template Quaternion<long double> sqr( const Quaternion<long double> & );
template Quaternion<long double> sqrt( const Quaternion<long double> & );
template Quaternion<long double> rotate( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> exp( const Quaternion<long double> & );
template Quaternion<long double> log( const Quaternion<long double> & );
template Quaternion<long double> log10( const Quaternion<long double> & );
template Quaternion<long double> pow( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> pow( const Quaternion<long double> &, long double );
template Quaternion<long double> inverse( const Quaternion<long double> & );
template Quaternion<long double> sin( const Quaternion<long double> & );
template Quaternion<long double> cos( const Quaternion<long double> & );
template Quaternion<long double> tan( const Quaternion<long double> & );
template Quaternion<long double> sec( const Quaternion<long double> & );
template Quaternion<long double> csc( const Quaternion<long double> & );
template Quaternion<long double> cot( const Quaternion<long double> & );
template Quaternion<long double> sinh( const Quaternion<long double> & );
template Quaternion<long double> cosh( const Quaternion<long double> & );
template Quaternion<long double> tanh( const Quaternion<long double> & );
template Quaternion<long double> sech( const Quaternion<long double> & );
template Quaternion<long double> csch( const Quaternion<long double> & );
template Quaternion<long double> coth( const Quaternion<long double> & );
template Quaternion<long double> asin( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> acos( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> atan( const Quaternion<long double> & );
template Quaternion<long double> asec( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> acsc( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> acot( const Quaternion<long double> & );
template Quaternion<long double> asinh( const Quaternion<long double> & );
template Quaternion<long double> acosh( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> atanh( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> asech( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> acsch( const Quaternion<long double> & );
template Quaternion<long double> acoth( const Quaternion<long double> &, const Quaternion<long double> & );
template Quaternion<long double> bessel_J( int, const Quaternion<long double> & );
template Quaternion<long double> floor( const Quaternion<long double> & );
template Quaternion<long double> ceil( const Quaternion<long double> & );
template Quaternion<long double> horner( const Quaternion<long double> &, long double *, unsigned int );
template Quaternion<long double> horner( const Quaternion<long double> &, long double *, long double *, unsigned int );