/*
   vecmat3.h
   version 2

   Ramses van Zon
   April/May 2009

   Fast three dimensional "Vector" and "Matrix" classes using
   expression templates, with elements of any type.  

   Read 'vecmat3v2.pdf' for information on how to use these
   classes.

   For the curious: At the end of this header file, this file's 
   main principles are explained.
*/

#ifndef _VECMAT3_
#define _VECMAT3_

#include <cmath>

/* 
 *  Handle difference between cxx (alpha) and g++ (e.g. linux)
 *  regarding ostream.
 */
#ifdef __alpha
#include <iostream.h>
#define OSTREAM ostream
#else
#include <iostream>
#define OSTREAM std::ostream
#endif

/*
 *  Some useful short-hand notational macros (while macros are evil,
 *  so is c++'s template notation, and a few abbreviations will
 *  alleviate our notational pain)
 */
#define ENODE(A,B,C)                    class A,int B,class C
#define EXPRESSION_TEMPLATE             template <typename T,ENODE(X,Y,Z)>
#define EXPRESSION_TEMPLATE_PAIR        template <typename T,ENODE(A,B,C),ENODE(D,E,F)> 
#define EXPRESSION_TEMPLATE_TRIPLET     template <typename T,ENODE(A,B,C),ENODE(D,E,F),ENODE(G,H,I)>
#define EXPRESSION_TEMPLATE_MEMBER      template <ENODE(X,Y,Z)>
#define ANOTHER_EXPRESSION_TEMPLATE     template <ENODE(U,V,W)>
#define CONVERTIBLE_TEMPLATE            template <typename CONVERT>
#define CONVERTIBLE_EXPRESSION_TEMPLATE template <typename T,ENODE(X,Y,Z),typename CONVERT>
#define CONVERTIBLE_ANOTHER_EXPRESSION_TEMPLATE template <ENODE(U,V,W),typename CONVERT>
#define VECTOR          Vector<T,X,Y,Z>
#define ANOTHER_VECTOR  Vector<T,U,V,W>
#define VECTOR1         Vector<T,A,B,C>
#define VECTOR2         Vector<T,D,E,F>
#define VECTOR3         Vector<T,G,H,I>
#define MATRIX          Matrix<T,X,Y,Z>
#define ANOTHER_MATRIX  Matrix<T,U,V,W>
#define MATRIX1         Matrix<T,A,B,C>
#define MATRIX2         Matrix<T,D,E,F>

/*
 *  A dedicated namespace for vector and matrix classes and
 *  operations.
 */
namespace vecmat3 {

    /*
     *  Enumeration type for possible templated operations for Vectors
     *  and Matrices.
     */
    enum {
        NoOp,        ///< tag to indicate no operation
        RowOp,       ///< tag to indicate row operation
        ColOp,       ///< tag to indicate columns operation
        PlusOp,      ///< tag for addition operator    
        MinusOp,     ///< tag for subtraction operator
        TimesOp,     ///< tag for multiplication operator
        NegativeOp,  ///< tag for negation operator
        TransposeOp, ///< tag to indicate transpose
        DyadicOp,    ///< tag for dyadic operator
        USER         ///< in case the user defines templated operations
    };
   
    /*
     *  Forward declaration for definitions Vector<T> & Matrix<T> and
     *  CommaOp.
     */
    class Base;  
    template <typename T,class A=Base,int B=NoOp,class C=Base> class Vector;
    template <typename T,class A=Base,int B=NoOp,class C=Base> class Matrix;
    template <typename T> class CommaOp;
    
    /*
     *  The Vector class
     */
    template <typename T>
    class Vector<T> 
    {
      public:
        T x;   ///< x element of the vector
        T y;   ///< y element of the vector
        T z;   ///< z element of the vector

        /*
         *  Constructors
         */
        EXPRESSION_TEMPLATE_MEMBER
        inline            Vector<T> ( const VECTOR& v );
        inline            Vector<T> () {}  
        inline explicit   Vector<T> ( T x_, T y_ = 0, T z_ = 0 ); 
        
        /*
         *  Assignment operators
         */
        EXPRESSION_TEMPLATE_MEMBER 
        inline const Vector<T>& operator= ( const VECTOR & v );
        inline const Vector<T>& operator= ( const Vector<T> & v );
        inline CommaOp<T>       operator= ( T e );
        
        /*
         *  Methods
         */
        inline T          nrm()  const;
        inline T          nrm2() const;
        inline void       zero();  
        
        /*
         *  Template evaluation
         */
        template <int I> 
        inline T          eval() const;

        /*
         *  Operators
         */
        inline const T&   operator() ( const int i ) const;
        inline T&         operator() ( const int i );
        
        /*
         *  Templated operators
         */
        CONVERTIBLE_TEMPLATE
        inline Vector<T>& operator*= ( const CONVERT s ); 

        CONVERTIBLE_TEMPLATE
        inline Vector<T>& operator/= ( const CONVERT s ); 
        
        EXPRESSION_TEMPLATE_MEMBER 
        inline Vector<T>& operator+= ( const VECTOR& v );

        EXPRESSION_TEMPLATE_MEMBER
        inline Vector<T>& operator-= ( const VECTOR& v );
    };

    /*
     *  The Matrix class
     */
    template <typename T>
    class Matrix<T> 
    {
      public:
        T xx;   ///< xx components of the matrix  
        T xy;   ///< xy components of the matrix  
        T xz;   ///< xz components of the matrix
        T yx;   ///< yx components of the matrix  
        T yy;   ///< yy components of the matrix  
        T yz;   ///< yz components of the matrix
        T zx;   ///< zx components of the matrix  
        T zy;   ///< zy components of the matrix  
        T zz;   ///< zz components of the matrix

        /*
         *  Constructors
         */
        EXPRESSION_TEMPLATE_MEMBER
        inline            Matrix<T>( const MATRIX& m );
        inline            Matrix<T>() {}
        inline explicit   Matrix<T>( T a,     T b = 0, T c = 0,
                                     T d = 0, T e = 0, T f = 0,
                                     T g = 0, T h = 0, T i = 0 );
   
        /*
         *  Assignment operators
         */
        EXPRESSION_TEMPLATE_MEMBER
        inline const Matrix<T>& operator= ( const MATRIX& m );
        inline const Matrix<T>& operator= ( const Matrix<T>& m );
        inline CommaOp<T> operator= ( T e );
   
        /*
         *  Evaluate components
         */
        template <int I, int J>
        inline T          eval() const; 
   
        /*
         *  Methods
         */
        inline Vector<T,Matrix<T>,RowOp> row   ( const int i ) const;
        inline Vector<T,Matrix<T>,ColOp> column( const int j ) const; 
        inline T          nrm2() const;                
        inline T          nrm()  const;                 
        inline T          tr()   const;
        inline T          det()  const;
        inline void       zero();
        inline void       one();
        inline void       reorthogonalize();

        /*
         *  Templated methods
         */
        EXPRESSION_TEMPLATE_MEMBER
        inline void       setRow    ( const int i, const VECTOR & v );
        EXPRESSION_TEMPLATE_MEMBER
        inline void       setColumn ( const int j, const VECTOR & v ); 
   
        /*
         *  Operators
         */
        inline const T&   operator() ( const int i, const int j ) const;
        inline T&         operator() ( const int i, const int j );

        /*
         *  Templated operators
         */
        CONVERTIBLE_TEMPLATE
        inline Matrix<T>& operator*= ( const CONVERT s );

        CONVERTIBLE_TEMPLATE
        inline Matrix<T>& operator/= ( const CONVERT s );
   
        EXPRESSION_TEMPLATE_MEMBER
        inline Matrix<T>& 
        operator+= ( const MATRIX & m );
   
        EXPRESSION_TEMPLATE_MEMBER
        inline Matrix<T>& 
        operator-= ( const MATRIX & m );
   
        EXPRESSION_TEMPLATE_MEMBER
        inline Matrix<T>& 
        operator*= ( const MATRIX & m );
    };

    /*
     * Definitions of vector-matrix operations
     */

    /*
     *  The operator+ takes two vector expressions and returns the
     *  appropriate vector expression.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,PlusOp,VECTOR2>  
    operator+ ( const VECTOR1 & v1, 
                const VECTOR2 & v2 );

    /* 
     *  The operator- takes two vector expressions and returning an
     *  appropriate vector expression.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,MinusOp,VECTOR2> 
    operator- ( const VECTOR1 & v1, 
                const VECTOR2 & v2 );

    /* 
     *  The operator^ Takes two vector expressions and returns the
     *  appropriate vector expression.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,TimesOp,VECTOR2> 
    operator^ ( const VECTOR1 & v1, 
                const VECTOR2 & v2 );

    /*
     *  The operator/ takes a vector expression and a T, and returns
     *  the appropriate vector expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T>        
    operator/ ( const VECTOR & v, 
                CONVERT s );

    /*
     *  The operator* takes a T and a vector expression, and returns
     *  the appropriate vector expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T>        
    operator* ( CONVERT        s, 
                const VECTOR & v );

    /* 
     *  The operator* takes a vector expression and a , and returns
     *  the appropriate vector expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T>        
    operator* ( const VECTOR & v, 
                CONVERT  s );

    /*
     *  The unary operator- takes a vector expression, and returns an
     *  appropriate vector expression.
     */
    EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,NegativeOp>       
    operator- ( const VECTOR & v );

    /*
     * The dot product is implemented as the operator| which takes two
     * vector expressions and returns a T.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline T                               
    operator| ( const VECTOR1 & v1, 
                const VECTOR2 & v2 );

    /*
     *  The dot product is also implemented as the operator* as takes
     *  two vector expressions and returns a T.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline T                               
    operator* ( const VECTOR1 & v1, 
                const VECTOR2 & v2 );

    /*
     *  The operator+ takes two matrix expressions and returns an
     *  appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,MATRIX1,PlusOp,MATRIX2>  
    operator+ ( const MATRIX1 & m1, 
                const MATRIX2 & m2 );

    /*
     *  The operator- takes two matrix expressions and returns an
     *  appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,MATRIX1,MinusOp,MATRIX2> 
    operator- ( const MATRIX1 & m1, 
                const MATRIX2 & m2 );

    /*
     *  The operator* takes two matrix expressions and returns an
     *  appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline Matrix<T,MATRIX1,TimesOp,MATRIX2> 
    operator* ( const MATRIX1 & m1, 
                const MATRIX2 & m2 );

    /*
     *  The operator* takes a matrix expressions and a vector
     *  expression, and returns the appropriate vector expression.
     */
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,MATRIX1,TimesOp,VECTOR2> 
    operator* ( const MATRIX1 & m,  
                const VECTOR2 & v );

    /*
     *  The operator* takes a T and a matrix expression, and returns
     *  the appropriate matrix expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Matrix<T,MATRIX,TimesOp,T>        
    operator* ( CONVERT s, 
                const MATRIX & m );

    /*
     *  The operator* takes a matrix expression and a T, and returns
     *  the appropriate matrix expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T>        
    operator* ( const MATRIX & m, 
                CONVERT        s );

    /*
     *  The operator/ takes a matrix expression and a T, and returns
     *  the appropriate matrix expression.
     */
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Matrix<T,MATRIX,TimesOp,T>        
    operator/ ( const MATRIX & m, 
                CONVERT s );

    /*
     *  The unary operator- takes a matrix expression and returns 
     *  an appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE 
    inline Matrix<T,MATRIX,NegativeOp>       
    operator- ( const MATRIX & m );

    /*
     *  The Transpose function, takes a matrix expression and returns
     *  an appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TransposeOp> 
    Transpose ( const MATRIX & m );

    /*
     *  The Dyadic function, takes two vector expressions and returns
     *  an appropriate matrix expression.
     */
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,VECTOR1,DyadicOp,VECTOR2>  
    Dyadic ( const VECTOR1 & v1, 
             const VECTOR2 & v2 );

    /*
     *  Return the distance between two vectors.
     */
    EXPRESSION_TEMPLATE_PAIR
    inline T dist ( const VECTOR1 & v1,
                    const VECTOR2 & v2 ); 

    /*
     *  Return the distance squared between two vectors.
     */
    EXPRESSION_TEMPLATE_PAIR
    inline T dist2 ( const VECTOR1 & v1, 
                     const VECTOR2 & v2 );

    /*
     *  Return the distance between two vectors with the first one shifted by s.
     */
    EXPRESSION_TEMPLATE_TRIPLET
    inline T distwithshift ( const VECTOR1 & v1, 
                             const VECTOR2 & v2, 
                             const VECTOR3 & v3 );

    /*
     *  Computes the inverse of a matrix expression.
     */
    EXPRESSION_TEMPLATE
    inline Matrix<T> 
    Inverse ( const MATRIX & m );

    /*
     *  Returns the rotation matrix around vector V by an angle vector.nrm().
     */
    EXPRESSION_TEMPLATE
    inline Matrix<T> 
    Rodrigues ( const VECTOR & v );

    /*
     *  Output to stream for vector expressions.
     */
    EXPRESSION_TEMPLATE 
    inline OSTREAM & 
    operator<< ( OSTREAM & o, 
                 const VECTOR & v );

    /*
     *  Output to stream for matrix expressions.
     */
    EXPRESSION_TEMPLATE 
    inline OSTREAM & 
    operator<< ( OSTREAM & o, 
                 const MATRIX & m );

/*
 * Convenient macros 
 */
#define crossProduct(a,b) ((a)^(b))
#define dotProduct(a,b)   ((a)|(b))
#define MTVmult(M,v)      (Transpose(M)*(v))

/*
 * Define the following for backward compatibility.
 */
#define returnVector(x,y,z) return Vector(x,y,z)
#define return_Vector(x,y,z) return Vector(x,y,z)
#define ZERO(v) v(0)
#define Vector_ZERO(v) Vector v(0)
#define Vector_INIT(v,x,y,z) Vector v(x,y,z)
#define DIFF(one,two) one-two
#define DIFFWITHSHIFT(one,two,shift) one-two+shift
#define Matrix_ZERO(m) Matrix m(0)
#define Matrix_INIT(m,xx,xy,xz,yx,yy,yz,zx,zy,zz) Matrix m(xx,xy,xz,yx,yy,yz,zx,zy,zz)
#define return_Matrix(xx,xy,xz,yx,yy,yz,zx,zy,zz) return Matrix(xx,xy,xz,yx,yy,yz,zx,zy,zz)
#define returnMatrix(xx,xy,xz,yx,yy,yz,zx,zy,zz) return Matrix(xx,xy,xz,yx,yy,yz,zx,zy,zz)

    /****************************/
    /*  INLINE IMPLEMENTATION  **/
    /****************************/

    // make a function for taking a square exists.
    template <typename T>
    inline T sqr(T x);

    template <typename T>
    inline T sqr(T x) 
    { 
        return x*x; 
    }

    //
    // member functions of the basic Vector type
    //
              
    // default constructor
    template <typename T>
    inline Vector<T>::Vector ( T x_, T y_, T z_ ) :
      x(x_), 
      y(y_), 
      z(z_) 
    {}

    // constuctor from a vector expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER 
    inline Vector<T>::Vector ( const VECTOR & v) :
      x(v.eval<0>()),
      y(v.eval<1>()),
      z(v.eval<2>())
    {}

    // assign zero to all elements
    template <typename T>
    inline void Vector<T>::zero() 
    {
        x = y = z = 0;
    }

    // access to elements through parenthesis
    // active: for assignment to ()
    template <typename T>
    inline T & Vector<T>::operator() ( const int i ) 
    { 
        return *(&x+i);
    }

    // passive
    template <typename T>
    inline const T & 
    Vector<T>::operator() ( const int i ) const 
    { 
        return *(&x+i);
    }

    // passive access to the elements through eval
    // (the basis of the template evaluation technique, really)
    template <typename T>
    template <int I> 
    inline T Vector<T>::eval() const 
    { 
        switch(I) {
        case 0: return x; 
        case 1: return y; 
        case 2: return z; 
        default: return 0;
        }
    }

    // norm squared
    template <typename T>
    inline T Vector<T>::nrm2() const 
    {
        return x*x+y*y+z*z;
    }

    // norm
    template <typename T>
    inline T Vector<T>::nrm() const 
    {
        T biggest = x<0?-x:x;
        T absval  = y<0?-y:y;
        if (absval>biggest) 
            biggest = absval;
        absval = z<0?-z:z;
        if (absval>biggest) 
            biggest = absval;
        if (biggest != 0)
            return biggest*(T)(sqrt(sqr(x/biggest)
                                    +sqr(y/biggest)
                                    +sqr(z/biggest)));
        else
            return 0;
    }

    // assignment from an expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER
    inline const Vector<T>& 
    Vector<T>::operator= ( const VECTOR & v ) 
    {
        T v_eval_1 = v.eval<1>();
        T v_eval_2 = v.eval<2>();
        x = v.eval<0>();
        y = v_eval_1;
        z = v_eval_2;
        return *this;
    }

    // assignment from a basic vector
    template <typename T>
    inline const Vector<T>& 
    Vector<T>::operator= ( const Vector<T> & v ) 
    {
        if (this != &v) {
            x = v.x;
            y = v.y;
            z = v.z;
        }
        return *this;
    }

    // addto-assignment from an expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER 
    inline Vector<T>& 
    Vector<T>::operator+= ( const VECTOR & v ) 
    {
        T v_eval_1 = v.eval<1>();
        T v_eval_2 = v.eval<2>();
        x += v.eval<0>();
        y += v_eval_1;
        z += v_eval_2;
        return *this;
    }

    // subtract-from assignment from an expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER 
    inline Vector<T>& 
    Vector<T>::operator-= ( const VECTOR & v ) 
    {
        T v_eval_1 = v.eval<1>();
        T v_eval_2 = v.eval<2>();
        x -= v.eval<0>();
        y -= v_eval_1;
        z -= v_eval_2;
        return *this;
    }

    // multiply-by assignment from a T
    template <typename T>
    CONVERTIBLE_TEMPLATE
    inline Vector<T>& 
    Vector<T>::operator*= ( const CONVERT s ) 
    {
        x *= s;
        y *= s;
        z *= s;
        return *this;
    }

    // divide-by assignment from a T
    template <typename T>
    CONVERTIBLE_TEMPLATE
    inline Vector<T>& 
    Vector<T>::operator/= ( const CONVERT s ) 
    {
        x /= s;
        y /= s;
        z /= s;
        return *this;
    }

    //                                           
    // member functions of the basic Matrix type 
    //                                           

    // constructor
    template <typename T>
    inline Matrix<T>::Matrix( T a, T b, T c,
                              T d, T e, T f,
                              T g, T h, T i ) :
      xx(a), xy(b), xz(c),
      yx(d), yy(e), yz(f),
      zx(g), zy(h), zz(i)
    {}

    // constructor from a Matrix expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER
    inline Matrix<T>::Matrix ( const MATRIX & m ) :
      xx(m.template eval<0,0>()), xy(m.template eval<0,1>()), xz(m.template eval<0,2>()),
      yx(m.template eval<1,0>()), yy(m.template eval<1,1>()), yz(m.template eval<1,2>()),
      zx(m.template eval<2,0>()), zy(m.template eval<2,1>()), zz(m.template eval<2,2>())
    {}

    // set all elements to zero
    template <typename T>
    inline void Matrix<T>::zero() 
    {
        xx = xy = xz = yx = yy = yz = zx = zy = zz = 0;
    }
 
    // turn this into an identity matrix
    template <typename T>
    inline void Matrix<T>::one() 
    {
        xx = yy = zz = 1;
        xy = xz = yx = yz = zx = zy = 0;
    }

    // set the elements of a row equal to those of a vector
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER 
    inline void Matrix<T>::setRow ( const int i, 
                                    const VECTOR & v ) 
    {
        switch(i) {
        case 0: xx = v.template eval<0>(); xy = v.template eval<1>(); xz = v.template eval<2>(); break;
        case 1: yx = v.template eval<0>(); yy = v.template eval<1>(); yz = v.template eval<2>(); break;
        case 2: zx = v.template eval<0>(); zy = v.template eval<1>(); zz = v.template eval<2>(); break;
        }
    }

    // set the elements of a column equal to those of a vector
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER 
    inline void Matrix<T>::setColumn ( const int j, 
                                       const VECTOR & v ) 
    {
        switch(j) {
        case 0: 
            xx = v.template eval<0>(); 
            yx = v.template eval<1>(); 
            zx = v.template eval<2>(); 
            break;
        case 1: 
            xy = v.template eval<0>(); 
            yy = v.template eval<1>(); 
            zy = v.template eval<2>(); 
            break;
        case 2: 
            xz = v.template eval<0>(); 
            yz = v.template eval<1>(); 
            zz = v.template eval<2>(); 
            break;
        }
    }

    // access to elements through parenthesis
    // active: for assignment to ()
    template <typename T>
    inline T & Matrix<T>::operator() ( const int i, 
                                       const int j ) 
    { 
        return *(&xx+3*i+j);
    }

    // passive:
    template <typename T>
    inline const T & Matrix<T>::operator() ( const int i, 
                                             const int j ) const 
    { 
        return *(&xx+3*i+j);
    }

    // passive access to the elements through eval function
    template <typename T>
    template <int I, int J>
    inline T Matrix<T>::eval() const 
    {
        switch (I) { 
        case 0: switch (J) { 
            case 0: return xx; 
            case 1: return xy; 
            case 2: return xz; 
            default: return 0; 
            } 
        case 1: switch (J) { 
            case 0: return yx; 
            case 1: return yy; 
            case 2: return yz; 
            default: return 0; 
            } 
        case 2: switch (J) { 
            case 0: return zx; 
            case 1: return zy; 
            case 2: return zz; 
            default: return 0; 
            } 
        } 
    }

    // trace
    template <typename T>
    inline T Matrix<T>::tr() const 
    { 
        return xx + yy + zz; 
    }

    // determinant
    template <typename T>
    inline T Matrix<T>::det() const 
    {
        return xx*(yy*zz-yz*zy)+xy*(yz*zx-yx*zz)+xz*(yx*zy-yy*zx);
    }

    // norm squared
    template <typename T>
    inline T Matrix<T>::nrm2() const 
    { 
        return sqr(xx) + sqr(xy) + sqr(xz) + sqr(yx) + sqr(yy) 
            + sqr(yz) + sqr(zx) + sqr(zy) + sqr(zz); 
    }

    // norm of the this matrix
    template <typename T>
    inline T Matrix<T>::nrm() const 
    {
        T absval;
        T biggest = xx<0?-xx:xx;
        absval = xy<0?-xy:xy;
        if (absval>biggest) 
            biggest = absval;
        absval = xz<0?-xz:xz;
        if (absval>biggest) 
            biggest = absval;
        absval = yx<0?-yx:yx;
        if (absval>biggest) 
            biggest = absval;
        absval = yy<0?-yy:yy;
        if (absval>biggest) 
            biggest = absval;
        absval = yz<0?-yz:yz;
        if (absval>biggest) 
            biggest = absval;
        absval = zx<0?-zx:zx;
        if (absval>biggest) 
            biggest = absval;
        absval = zy<0?-zy:zy;
        if (absval>biggest) 
            biggest = absval;
        absval = zz<0?-zz:zz;
        if (absval>biggest) 
            biggest = absval;
        if (biggest!= 0) 
            return biggest*(T)(sqrt(sqr(xx/biggest)
                                    +sqr(xy/biggest)
                                    +sqr(xz/biggest)
                                    +sqr(yx/biggest)
                                    +sqr(yy/biggest)
                                    +sqr(yz/biggest)
                                    +sqr(zx/biggest)
                                    +sqr(zy/biggest)
                                    +sqr(zz/biggest)));
        else 
            return 0;
    } 

    // assignment from a matrix expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER
    inline const Matrix<T>& 
    Matrix<T>::operator= ( const MATRIX & m ) 
    {
        T m_eval_xx, m_eval_xy, m_eval_xz, 
            m_eval_yx, m_eval_yy, m_eval_yz, 
            m_eval_zx, m_eval_zy;
        m_eval_xx = m.template eval<0,0>(); m_eval_xy = m.template eval<0,1>(); m_eval_xz = m.template eval<0,2>();
        m_eval_yx = m.template eval<1,0>(); m_eval_yy = m.template eval<1,1>(); m_eval_yz = m.template eval<1,2>();
        m_eval_zx = m.template eval<2,0>(); m_eval_zy = m.template eval<2,1>();
        zz = m.template eval<2,2>();
        xx = m_eval_xx; xy = m_eval_xy; xz = m_eval_xz;
        yx = m_eval_yx; yy = m_eval_yy; yz = m_eval_yz;
        zx = m_eval_zx; zy = m_eval_zy;
        return *this;
    }

    // assignment from a basic matrix
    template <typename T> 
    inline const Matrix<T>& 
    Matrix<T>::operator= ( const Matrix<T>& m ) 
    {
        if (this != &m) {
            xx = m.xx; yx = m.yx; zx = m.zx;
            xy = m.xy; yy = m.yy; zy = m.zy;
            xz = m.xz; yz = m.yz; zz = m.zz;
        }
        return *this;
    }

    // multiply-by assignment from a T
    template <typename T>
    CONVERTIBLE_TEMPLATE
    inline Matrix<T>& 
    Matrix<T>::operator*= ( const CONVERT s ) 
    {
        xx *= s; xy *= s; xz *= s;
        yx *= s; yy *= s; yz *= s;
        zx *= s; zy *= s; zz *= s;
        return *this;
    }

    // divide-by assignment from a T
    template <typename T>
    CONVERTIBLE_TEMPLATE
    inline Matrix<T>& 
    Matrix<T>::operator/= ( const CONVERT s ) 
    {
        xx /= s; xy /= s; xz /= s;
        yx /= s; yy /= s; yz /= s;
        zx /= s; zy /= s; zz /= s;
        return *this;
    }

    // add-to assignment from an expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER
    inline Matrix<T>& 
    Matrix<T>::operator+= ( const MATRIX & m ) 
    {
        T m_eval_xx, m_eval_xy, m_eval_xz, 
            m_eval_yx, m_eval_yy, m_eval_yz, 
            m_eval_zx, m_eval_zy;
        m_eval_xx = m.template eval<0,0>(); m_eval_xy = m.template eval<0,1>(); m_eval_xz = m.template eval<0,2>();
        m_eval_yx = m.template eval<1,0>(); m_eval_yy = m.template eval<1,1>(); m_eval_yz = m.template eval<1,2>();
        m_eval_zx = m.template eval<2,0>(); m_eval_zy = m.template eval<2,1>(); zz += m.template eval<2,2>();
        xx += m_eval_xx; xy += m_eval_xy; xz += m_eval_xz;
        yx += m_eval_yx; yy += m_eval_yy; yz += m_eval_yz;
        zx += m_eval_zx; zy += m_eval_zy;
        return *this;
    }

    // subtract-from assignment from an expression
    template <typename T> 
    EXPRESSION_TEMPLATE_MEMBER
    inline Matrix<T>& 
    Matrix<T>::operator-= ( const MATRIX & m ) 
    {
        T m_eval_xx, m_eval_xy, m_eval_xz, 
            m_eval_yx, m_eval_yy, m_eval_yz, 
            m_eval_zx, m_eval_zy;
        m_eval_xx = m.template eval<0,0>(); m_eval_xy = m.template eval<0,1>(); m_eval_xz = m.template eval<0,2>();
        m_eval_yx = m.template eval<1,0>(); m_eval_yy = m.template eval<1,1>(); m_eval_yz = m.template eval<1,2>();
        m_eval_zx = m.template eval<2,0>(); m_eval_zy = m.template eval<2,1>(); zz -= m.template eval<2,2>();
        xx -= m_eval_xx; xy -= m_eval_xy; xz -= m_eval_xz;
        yx -= m_eval_yx; yy -= m_eval_yy; yz -= m_eval_yz;
        zx -= m_eval_zx; zy -= m_eval_zy;
        return *this;
    }

    // multiply-by-matrix assignment from an expression
    template <typename T>
    EXPRESSION_TEMPLATE_MEMBER
    inline Matrix<T>& 
    Matrix<T>::operator*= ( const MATRIX & m ) 
    {
        Matrix<T> rhs(m);
        T x, y;
        x = xx; y = xy;
        xx *= rhs.xx; xx += y*rhs.yx+xz*rhs.zx;
        xy *= rhs.yy; xy += x*rhs.xy+xz*rhs.zy;
        xz *= rhs.zz; xz += x*rhs.xz+ y*rhs.yz;
        x = yx; y = yy;
        yx *= rhs.xx; yx += y*rhs.yx+yz*rhs.zx;
        yy *= rhs.yy; yy += x*rhs.xy+yz*rhs.zy;
        yz *= rhs.zz; yz += x*rhs.xz+ y*rhs.yz;
        x = zx; y = zy;
        zx *= rhs.xx; zx += y*rhs.yx+zz*rhs.zx;
        zy *= rhs.yy; zy += x*rhs.xy+zz*rhs.zy;
        zz *= rhs.zz; zz += x*rhs.xz+ y*rhs.yz;
        return *this;
    }


    //
    // member functions of Vector and Matrix
    //

    // define the body of const members function of Vector

    // to access the elements of a vector expression:
#define VECPARENTHESES operator() ( const int i ) const {\
        switch(i){                                          \
        case 0: return eval<0>();                           \
        case 1: return eval<1>();                           \
        case 2: return eval<2>();                           \
        default: return 0;                                  \
        }                                                   \
    }

    // to access the elements of a matrix expression
#define MATPARENTHESIS operator() ( const int i, \
                                    const int j ) const { \
        switch(i){                                           \
        case 0: switch(j){                                   \
            case 0: return eval<0,0>();                      \
            case 1: return eval<0,1>();                      \
            case 2: return eval<0,2>();                      \
            default: return 0;                               \
            }                                                \
        case 1: switch(j){                                   \
            case 0: return eval<1,0>();                      \
            case 1: return eval<1,1>();                      \
            case 2: return eval<1,2>();                      \
            default: return 0;                               \
            }                                                \
        case 2: switch(j){                                   \
            case 0: return eval<2,0>();                      \
            case 1: return eval<2,1>();                      \
            case 2: return eval<2,2>();                      \
            default: return 0;                               \
            }                                                \
        }                                                    \
    }

    // to get the norm of a vector
#define VECNRM nrm() const {                           \
        T x = eval<0>();                                  \
        T y = eval<1>();                                  \
        T z = eval<2>();                                  \
        T biggest = x<0?-x:x;                             \
        T absval = y<0?-y:y;                              \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = z<0?-z:z;                                \
        if (absval>biggest)                               \
            biggest = absval;                             \
        if (biggest!= 0)                                  \
            return biggest*(T)(sqrt(sqr(x/biggest)        \
                                    +sqr(y/biggest)       \
                                    +sqr(z/biggest)));    \
        else                                              \
            return 0;                                     \
    }

    // to get the norm of a matrix
#define MATNRM nrm() const {                           \
        T  xx = eval<0,0>(), xy = eval<0,1>(), xz = eval<0,2>(); \
        T  yx = eval<1,0>(), yy = eval<1,1>(), yz = eval<1,2>(); \
        T  zx = eval<2,0>(), zy = eval<2,1>(), zz = eval<2,2>(); \
        T  absval;                                        \
        T  biggest = xx<0?-xx:xx;                         \
        absval = xy<0?-xy:xy;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = xz<0?-xz:xz;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = yx<0?-yx:yx;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = yy<0?-yy:yy;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = yz<0?-yz:yz;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = zx<0?-zx:zx;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = zy<0?-zy:zy;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        absval = zz<0?-zz:zz;                             \
        if (absval>biggest)                               \
            biggest = absval;                             \
        if (biggest!= 0)                                  \
            return biggest*(T)(sqrt(sqr(xx/biggest)       \
                                    +sqr(xy/biggest)      \
                                    +sqr(xz/biggest)      \
                                    +sqr(yx/biggest)      \
                                    +sqr(yy/biggest)      \
                                    +sqr(yz/biggest)      \
                                    +sqr(zx/biggest)      \
                                    +sqr(zy/biggest)      \
                                    +sqr(zz/biggest)));   \
        else                                              \
            return 0;                                     \
    } 

    // to get the norm squared of a vector expression
#define VECNRM2 nrm2() const {                           \
        return sqr(eval<0>())+sqr(eval<1>())+sqr(eval<2>());\
    }

    // to get the norm squared of a matrix expression
#define MATNRM2 nrm2() const {                                 \
        return sqr(eval<0,0>())+sqr(eval<0,0>())+sqr(eval<0,0>()) \
            +sqr(eval<0,0>())+sqr(eval<0,0>())+sqr(eval<0,0>())   \
            +sqr(eval<0,0>())+sqr(eval<0,0>())+sqr(eval<0,0>());  \
    }

    // to get the trace of a matrix expression
#define MATTR tr() const {                             \
        return eval<0,0>() + eval<1,1>() + eval<2,2>();   \
    }

    // to get the determinant of a matrix expression
#define MATDET det() const {                                   \
        T  xx = eval<0,0>(), xy = eval<0,1>(), xz = eval<0,2>();  \
        T  yx = eval<1,0>(), yy = eval<1,1>(), yz = eval<1,2>();  \
        T  zx = eval<2,0>(), zy = eval<2,1>(), zz = eval<2,2>();  \
        return xx*(yy*zz-yz*zy)+xy*(yz*zx-yx*zz)+xz*(yx*zy-yy*zx);\
    }

    // to get the i-th row of a matrix expression
    // note: the class of the matrix expression has to be defined in CLASS
#define MATROW Vector<T,CLASS,RowOp> row ( const int i ) const {\
        return Vector<T,CLASS,RowOp>( *this, i );                  \
    }

    // to get the j-th column of a matrix expression
    // note: the class of the matrix expression has to be defined in CLASS
#define MATCOLUMN Vector<T,CLASS,ColOp> column ( const int j ) const {\
        return Vector<T,CLASS,ColOp>( *this, j );                        \
    }


    // Expression template class for '+'  between two Vector expressions
    EXPRESSION_TEMPLATE_PAIR
    class Vector<T,VECTOR1,PlusOp,VECTOR2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        { 
            return left->eval<I>() + right->eval<I>(); 
        }

        // constructor
        inline Vector ( const VECTOR1 & _left, 
                        const VECTOR2 & _right ) :
          left(&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const VECTOR1* left;
        const VECTOR2* right;
    };

  
    // Expression template class for '-'  between two Vector expressions
    EXPRESSION_TEMPLATE_PAIR
    class Vector<T,VECTOR1,MinusOp,VECTOR2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        { 
            return left->eval<I>() - right->eval<I>(); 
        }

        // constructor
        inline Vector( const VECTOR1 & _left, 
                       const VECTOR2 & _right ) :
          left(&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const VECTOR1* left;
        const VECTOR2* right;
    };


    // Expression template class for '^' between two Vector expressions
    EXPRESSION_TEMPLATE_PAIR
    class Vector<T,VECTOR1,TimesOp,VECTOR2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const  
        {
            switch (I) {
            case 0:
                return left->eval<1>() * right->eval<2>() - left->eval<2>() * right->eval<1>();
                break;
            case 1:
                return left->eval<2>() * right->eval<0>() - left->eval<0>() * right->eval<2>();
                break;
            case 2:
                return left->eval<0>() * right->eval<1>() - left->eval<1>() * right->eval<0>();
                break;
            default: 
                return 0;
                break;
            }
        }

        // constructor
        inline Vector( const VECTOR1 & _left, 
                       const VECTOR2 & _right ) :
          left(&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const VECTOR1* left;
        const VECTOR2* right;
    };


    // Expression template class for '*' between a Vector and a T
    EXPRESSION_TEMPLATE
    class Vector<T,VECTOR,TimesOp,T> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2  
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        { 
            return left->eval<I>() * right; 
        }

        // optimize a second multiplication with T
        CONVERTIBLE_TEMPLATE
        inline Vector & operator* ( CONVERT c )
        {
            right *= c;
            return *this;
        }

        // optimize a subsequent division by T
        CONVERTIBLE_TEMPLATE
        inline Vector & operator/ ( CONVERT c ) 
        {
            right /= c;
            return *this;
        }

        // constructor
        inline Vector( const VECTOR & _left, T _right ) :
          left(&_left), 
          right(_right) 
        {}

      private:

        // pointer to the subexpression and the T to multiply with
        const VECTOR* left;
        T right;

        // befriend multiplication operators that optimize further T
        // multiplication
        CONVERTIBLE_ANOTHER_EXPRESSION_TEMPLATE  
        friend Vector<T,ANOTHER_VECTOR,TimesOp,T>& 
        operator* (CONVERT b, 
                   Vector<T,ANOTHER_VECTOR,TimesOp,T> & a);
    };


    // Expression template class for unary '-' acting on a Vector expression
    EXPRESSION_TEMPLATE
    class Vector<T,VECTOR,NegativeOp> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        {
            return  - left->eval<I>();
        }

        // constructor
        inline Vector ( const VECTOR& _left ) :
          left(&_left) 
        {}

      private:

        // pointer to the subexpression
        const VECTOR* left;
    };


    // Expression template class for '+' between two Matrix expressions
    EXPRESSION_TEMPLATE_PAIR
    class Matrix<T,MATRIX1,PlusOp,MATRIX2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm(),(*this).te(), (*this).det(), element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS 

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX1,PlusOp,MATRIX2>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J>
        inline T eval() const 
        { 
            return left->eval<I,J>() + right->eval<I,J>(); 
        }

        // constructor
        inline Matrix( const MATRIX1 & _left, 
                       const MATRIX2 & _right ) :
          left(&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const MATRIX1* left;
        const MATRIX2* right;
    };


    // Expression template class for '-' between two Matrix expressions
    EXPRESSION_TEMPLATE_PAIR
    class Matrix<T,MATRIX1,MinusOp,MATRIX2> 
    {
      public:

        // define nrm2(), nrm(), tr(), det() and element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX1,MinusOp,MATRIX2>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J> 
        inline T eval() const 
        { 
            return left->eval<I,J>() - right->eval<I,J>(); 
        }

        // constructor
        inline Matrix( const MATRIX1 & _left,
                       const MATRIX2 & _right ) :
          left (&_left), 
          right(&_right) 
        {}
        
      private:

        // pointers to the sub-expressions
        const MATRIX1* left;
        const MATRIX2* right;
    };


    // Expression template class for '*' between a Matrix expression and a T
    EXPRESSION_TEMPLATE
    class Matrix<T,MATRIX,TimesOp,T> 
    {
      public:

        // define nrm2(), nrm(), tr(), det() and element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX,TimesOp,T>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // constructor
        inline  Matrix ( const MATRIX & _left, 
                         T _right ) :
          left(&_left), 
          right(_right) 
        {}

        // define what this operation evaluates to
        template <int I, int J> 
        inline T eval() const 
        { 
            return left->eval<I,J>() * right; 
        }

        CONVERTIBLE_TEMPLATE
        inline Matrix & operator* ( CONVERT c ) 
        {
            right *= c;
            return *this;
        }

        CONVERTIBLE_TEMPLATE
        inline Matrix & operator/ ( CONVERT c ) 
        {
            right /= c;
            return *this;
        }

      private:

        // pointer to the sub-expression and the T to multiply with
        const MATRIX* left;
        T right;

        // be-friend multiplication operators that optimize further T
        // multiplication
        CONVERTIBLE_ANOTHER_EXPRESSION_TEMPLATE 
        friend Matrix<T,ANOTHER_MATRIX,TimesOp,T>& 
        operator*(CONVERT b, 
                  Matrix<T,ANOTHER_MATRIX,TimesOp,T>& a);
    };

    // Expression template class for '*' between two Matrix expressions
    EXPRESSION_TEMPLATE_PAIR
    class Matrix<T,MATRIX1,TimesOp,MATRIX2> 
    {
      public:

        // define nrm2(), nrm(), tr(), det() and element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX1,TimesOp,MATRIX2>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J> inline T eval() const 
        {
            return  left->eval<I,0>() * right->eval<0,J>()
                  + left->eval<I,1>() * right->eval<1,J>()
                  + left->eval<I,2>() * right->eval<2,J>();
        }

        // constructor
        inline Matrix ( const MATRIX1 & _left, 
                        const MATRIX2 & _right) :
          left (&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const MATRIX1* left;
        const MATRIX2* right;
    };


    // Expression template class for unary '-' acting on a Matrix expression 
    EXPRESSION_TEMPLATE
    class Matrix<T,MATRIX,NegativeOp> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm(),(*this).te(), (*this).det(), element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX,NegativeOp>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J> 
        inline T eval() const 
        {
            return - right->eval<I,J>;
        }

        // constructor
        inline Matrix ( const MATRIX & _right ) :
          right(&_right) 
        {}

      private:

        // pointer to the sub-expression
        const MATRIX* right;
    };


    // Expression template class for transpose function 
    EXPRESSION_TEMPLATE
    class Matrix<T,MATRIX,TransposeOp> 
    {
      public:
        // define (*this).nrm2(), (*this).nrm(),(*this).te(), (*this).det(), element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,MATRIX,TransposeOp>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J> 
        inline T eval() const 
        {
            return  left->eval<J,I>();
        }

        // constructor
        inline Matrix ( const MATRIX & _left ) :
          left(&_left) 
        {}

      private:

        // pointer to the sub-expression
        const MATRIX* left;
    };


    // Expression template class for Dyadic operation between two Vector expressions
    EXPRESSION_TEMPLATE_PAIR
    class Matrix<T,VECTOR1,DyadicOp,VECTOR2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm(),(*this).te(), (*this).det(), element access (*this)(i,j)
        inline T MATNRM2
        inline T MATNRM
        inline T MATDET
        inline T MATTR
        inline T MATPARENTHESIS

        // define row and column access as if they were Vectors
        // CLASS is used in the MATROW and MATCOLUMN macros
#define CLASS Matrix<T,VECTOR1,DyadicOp,VECTOR2>
        inline MATROW
        inline MATCOLUMN
#undef CLASS

        // define what this operation evaluates to
        template <int I, int J> 
        inline T eval() const 
        {
            return left->eval<I>() * right->eval<J>();
        }

        // constructor
        inline Matrix( const VECTOR1 & _left, 
                       const VECTOR2 & _right) :
          left (&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expression
        const VECTOR1* left;
        const VECTOR2* right;
    };


    // Expression template class for row operation acting on a Matrix expression 
    EXPRESSION_TEMPLATE
    class Vector<T,MATRIX,RowOp> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int J> 
        inline T eval() const 
        {
            switch(i) {
            case 0: return left->eval<0,J>();
            case 1: return left->eval<1,J>();
            case 2: return left->eval<2,J>();
            default: return 0;
            }
        }

        // constructor
        inline Vector ( const MATRIX & _left, 
                        int _i) :
          left(&_left), 
          i(_i) 
        {}

      private:

        // pointer to the sub-expression and the index of the row
        const MATRIX* left;
        const int     i;
    };


    // Expression template class for the column operation acting on a Matrix expression
    EXPRESSION_TEMPLATE
    class Vector<T,MATRIX,ColOp> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        {
            switch(j) {
            case 0: return left->eval<I,0>();
            case 1: return left->eval<I,1>();
            case 2: return left->eval<I,2>();
            default: return 0;
            }
        }

        // constructor
        inline Vector ( const MATRIX & _left, int _j ) :
          left(&_left), 
          j(_j) 
        {}

      private:

        // pointer to the sub-expression and the index of the column
        const MATRIX* left;
        const int     j;
    };


    // Expression template class for '*' between a Matrix expression and a
    // Vector expression
    EXPRESSION_TEMPLATE_PAIR
    class Vector<T,MATRIX1,TimesOp,VECTOR2> 
    {
      public:

        // define (*this).nrm2(), (*this).nrm() and element access (*this)(i) 
        inline T VECNRM2
        inline T VECNRM
        inline T VECPARENTHESES

        // define what this operation evaluates to
        template <int I> 
        inline T eval() const 
        {
            return left->eval<I,0>() * right->eval<0>()
                + left->eval<I,1>() * right->eval<1>()
                + left->eval<I,2>() * right->eval<2>();
        }

        // constructor
        inline Vector ( const MATRIX1 & _left, 
                        const VECTOR2 & _right) :
          left (&_left), 
          right(&_right) 
        {}

      private:

        // pointers to the sub-expressions
        const MATRIX1* left;
        const VECTOR2* right;
    };


    // helper class to implemement a comma list assignment
    template <typename T>
    class CommaOp 
    {
      public:

        // the comma operator sets the first element and moves the
        // pointer to the next element
        inline CommaOp & operator, ( T e )
        { 
            // make sure the vector or matrix is not completely filled yet
            if (ptr <= end) 
                *ptr++ = e;  
            return *this; 
        }

        // destructor: fill whatever elements remain with zeros
        inline ~CommaOp()
        { 
            while (ptr <= end) 
                *ptr++ = 0;
        }

      private:

        // store pointers to the next element to be filled and the last
        // element fillable.
        T* ptr;
        T* const end; 

        // constructor: set pointer to the element to be filled next and
        // to the last element
        inline CommaOp ( T & e1, 
                         T & e2 ) :
          ptr(&e1), 
          end(&e2) 
        {}

        // this constructor is private, and so that the CommaOp class can
        // only be used by the friend member functions operator= of
        // Vector<T> and Matrix<T>
        friend CommaOp<T> Vector<T>::operator= (T e);
        friend CommaOp<T> Matrix<T>::operator= (T e);
    };

    // assignment operator of Vector<T> which triggers a comma separated
    // initializer list
    template <typename T>
    inline CommaOp<T> 
    Vector<T>::operator= (T e) 
    {
        x = e;
        return CommaOp<T>(y, z);
    }

    // assignment operatorof Matrix<T> that triggers a comma separated
    // initializer list
    template <typename T>
    inline CommaOp<T> 
    Matrix<T>::operator= (T e) 
    {
        xx = e;
        return CommaOp<T>(xy, zz);
    }

     //
    // definition of the operators
     //

    // Vector + Vector
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,PlusOp,VECTOR2> 
    operator+ ( const VECTOR1 & v1, 
                const VECTOR2 & v2) 
    { 
        return Vector<T,VECTOR1,PlusOp,VECTOR2>(v1, v2);
    }

    // Vector - Vector
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,MinusOp,VECTOR2> 
    operator- ( const VECTOR1 & v1, 
                const VECTOR2 & v2 ) 
    { 
        return Vector<T,VECTOR1,MinusOp,VECTOR2>(v1, v2);
    }

    // Vector ^ Vector
    EXPRESSION_TEMPLATE_PAIR 
    inline Vector<T,VECTOR1,TimesOp,VECTOR2> 
    operator^ ( const VECTOR1 & v1, 
                const VECTOR2 & v2 ) 
    { 
        return Vector<T,VECTOR1,TimesOp,VECTOR2>(v1, v2);
    }

    // T * Vector
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T> 
    operator* ( CONVERT s, 
                const VECTOR & v ) 
    { 
        return Vector<T,VECTOR,TimesOp,T>(v, s);
    }

    // Vector * T
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T> 
    operator* ( const VECTOR & v, 
                CONVERT s ) 
    { 
        return Vector<T,VECTOR,TimesOp,T>(v, s);
    }

    // Vector / T
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T> 
    operator/ ( const VECTOR & v, 
                CONVERT s ) 
    { 
        return Vector<T,VECTOR,TimesOp,T>(v, 1/s);
    }

    // Vector * T * T
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T>& 
    operator* ( CONVERT s, 
                Vector<T,VECTOR,TimesOp,T> & v ) 
    { 
        v.right *= s; 
        return v; 
    }

    // Vector * T / T
    CONVERTIBLE_EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,TimesOp,T>& 
    operator/ ( Vector<T,VECTOR,TimesOp, T> & v, 
                CONVERT s ) 
    { 
        v.right /= s; 
        return v; 
    }

    // - Vector
    EXPRESSION_TEMPLATE 
    inline Vector<T,VECTOR,NegativeOp> 
    operator- ( const VECTOR & v ) 
    { 
        return Vector<T,VECTOR,NegativeOp>(v);
    }

    // Vector | Vector
    EXPRESSION_TEMPLATE_PAIR
    inline T 
    operator| ( const VECTOR1 & v1, 
                const VECTOR2 & v2 ) 
    {
        return v1.template eval<0>()*v2.template eval<0>() 
            + v1.template eval<1>()*v2.template eval<1>() 
            + v1.template eval<2>()*v2.template eval<2>();
    }

    // Vector * Vector
    EXPRESSION_TEMPLATE_PAIR
    inline T 
    operator* ( const VECTOR1 & v1, 
                const VECTOR2 & v2 ) 
    {
        return v1.template eval<0>()*v2.template eval<0>() 
            + v1.template eval<1>()*v2.template eval<1>() 
            + v1.template eval<2>()*v2.template eval<2>();
    }

    // Matrix + Matrix
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,MATRIX1,PlusOp,MATRIX2> 
    operator+ ( const MATRIX1 & v1, 
                const MATRIX2 & v2 ) 
    { 
        return Matrix<T,MATRIX1,PlusOp,MATRIX2>(v1, v2);
    }

    // Matrix - Matrix
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,MATRIX1,MinusOp,MATRIX2>
    operator- ( const MATRIX1 & v1, 
                const MATRIX2 & v2 ) 
    { 
        return Matrix<T,MATRIX1,MinusOp,MATRIX2>(v1, v2);
    }

    // T * Matrix
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T> 
    operator* ( CONVERT s, 
                const MATRIX & m ) 
    { 
        return Matrix<T,MATRIX,TimesOp,T>(m, s);
    }

    // Matrix * T
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T> 
    operator* ( const MATRIX & m, 
                CONVERT s ) 
    { 
        return Matrix<T,MATRIX,TimesOp,T>(m, s);
    }

    // Matrix / T
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T> 
    operator/ ( const MATRIX & m, 
                CONVERT s ) 
    { 
        return Matrix<T,MATRIX,TimesOp,T>(m, 1.0/s);
    }

    // Matrix * T * T
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T>& 
    operator* ( CONVERT s, 
                Matrix<T,MATRIX,TimesOp,T> & m ) 
    { 
        m.right *= s; 
        return m; 
    }

    // Matrix * T / T
    CONVERTIBLE_EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TimesOp,T>& 
    operator/ ( Matrix<T,MATRIX,TimesOp,T> & m, 
                CONVERT s ) 
    { 
        m.right /= s; 
        return m; 
    }

    // Matrix * Matrix
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,MATRIX1,TimesOp,MATRIX2> 
    operator* ( const MATRIX1 & a, 
                const MATRIX2 & b ) 
    { 
        return Matrix<T,MATRIX1,TimesOp,MATRIX2>(a, b);
    }

    // -Matrix 
    EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,NegativeOp> 
    operator- ( const MATRIX & m ) 
    { 
        return Matrix<T,MATRIX,NegativeOp>(m);
    }

    // Matrix * Vector
    EXPRESSION_TEMPLATE_PAIR
    inline Vector<T,MATRIX1,TimesOp,VECTOR2> 
    operator* ( const MATRIX1 & m, 
                const VECTOR2 & v ) 
    { 
        return Vector<T,MATRIX1,TimesOp,VECTOR2> (m, v);
    }

     //
    // Implementation of the row and column member functions of a Matrix
     //

    template <typename T>
    inline Vector<T,Matrix<T>,RowOp> 
    Matrix<T>::row ( const int i ) const 
    {
        return Vector<T,Matrix<T>,RowOp>(*this,i);
    }

    template <typename T>
    inline Vector<T,Matrix<T>,ColOp> 
    Matrix<T>::column ( const int j ) const 
    {
        return Vector<T,Matrix<T>,ColOp>(*this,j);
    }

    // Can now implement Gramm-Schmidt orthogonalization of the rows of the matrix
    template <typename T>
    inline void Matrix<T>::reorthogonalize() 
    {
        T z;
        int num = 10;
        while ( fabs(det()-1)> 1E-16 and --num ) {
            z = 1/row(0).nrm();   
            xx *= z;    
            xy *= z;    
            xz *= z;    
            z = xx*yx + xy*yy + xz*yz;          
            yx -= z*xx; 
            yy -= z*xy; 
            yz -= z*xz; 
            z = 1/row(1).nrm();                
            yx *= z;  
            yy *= z;  
            yz *= z;        
            zx = xy*yz - xz*yy;                 
            zy = xz*yx - xx*yz;
            zz = xx*yy - xy*yx;
            z = 1/row(2).nrm();                
            zx *= z;  
            zy *= z;  
            zz *= z;        
        }
    }

     //
    // Non-member functions
     //

    // return |a-b|
    EXPRESSION_TEMPLATE_PAIR
    inline T dist ( const VECTOR1 & v1,
                    const VECTOR2 & v2 ) 
    {
        T x = fabs(v1.template eval<0>() - v2.template eval<0>());
        T y = fabs(v1.template eval<1>() - v2.template eval<1>());
        T z = fabs(v1.template eval<2>() - v2.template eval<2>());
        T biggest = x;  
        if (y > biggest) 
            biggest = y;
        if (z > biggest) 
            biggest = z;   
        if (biggest != 0) {
            x /= biggest;
            y /= biggest;
            z /= biggest;
            return biggest*sqrt(x*x+y*y+z*z);
        }
        else 
            return 0;
    }

    // return (a-b)|(a-b)
    EXPRESSION_TEMPLATE_PAIR
    inline T dist2 ( const VECTOR1 & v1, 
                     const VECTOR2 & v2 ) 
    {
        T  d = v1.template eval<0>() - v2.template eval<0>();
        T  c = d * d;
        d = v1.template eval<1>() - v2.template eval<1>();
        c += d * d;
        d = v1.template eval<2>() - v2.template eval<2>();
        return c + d * d;
    }

    // return |a+s-b|
    EXPRESSION_TEMPLATE_TRIPLET
    inline T distwithshift ( const VECTOR1 & v1, 
                             const VECTOR2 & v2, 
                             const VECTOR3 & v3 ) 
    {
        T  x = fabs(v3.template eval<0>()+v1.template eval<0>()-v2.template eval<0>());
        T  y = fabs(v3.template eval<1>()+v1.template eval<1>()-v2.template eval<1>());
        T  z = fabs(v3.template eval<2>()+v1.template eval<2>()-v2.template eval<2>());
        T  biggest = x;
        if (y>biggest) biggest = y;
        if (z>biggest) biggest = z;   
        if (biggest != 0) {
            x /= biggest;
            y /= biggest;
            z /= biggest;
            return biggest*sqrt(x*x+y*y+z*z);
        }
        else 
            return 0;
    }

    // the inverse of a matrix
    EXPRESSION_TEMPLATE
    inline Matrix<T> 
    Inverse ( const MATRIX & me ) 
    {
        Matrix<T> m = me;
        T s = 1/m.det();
        return Matrix<T>(
            (m.yy*m.zz-m.yz*m.zy)*s, -(m.xy*m.zz-m.xz*m.zy)*s,  (m.xy*m.yz-m.xz*m.yy)*s,
            -(m.yx*m.zz-m.yz*m.zx)*s,  (m.xx*m.zz-m.xz*m.zx)*s, -(m.xx*m.yz-m.xz*m.yx)*s,
            (m.yx*m.zy-m.yy*m.zx)*s, -(m.xx*m.zy-m.xy*m.zx)*s,  (m.xx*m.yy-m.xy*m.yx)*s);
    }

    // rotation matrix built using the Rodrigues formula
    EXPRESSION_TEMPLATE
    inline Matrix<T> 
    Rodrigues ( const VECTOR & ve ) 
    {
        Vector<T> v = ve;
        T theta = v.nrm();
        if (theta != 0) {
            T s, c;
            #ifdef SINCOS
            SINCOS(theta,&s,&c);
            #else
            s = sin(theta);
            c = cos(theta);
            #endif
            T inrm = 1/theta;
            T wx = v.x*inrm;
            T wy = v.y*inrm;
            T wz = v.z*inrm;
            T oneminusc = 1-c;
            T wxwy1mc = wx*wy*oneminusc;
            T wxwz1mc = wx*wz*oneminusc;
            T wywz1mc = wy*wz*oneminusc;
            T wxs = wx*s;
            T wys = wy*s;
            T wzs = wz*s;
            return Matrix<T>(c+wx*wx*oneminusc, wxwy1mc-wzs,       wxwz1mc+wys,
                             wxwy1mc+wzs,       c+wy*wy*oneminusc, wywz1mc-wxs,
                             wxwz1mc-wys,       wywz1mc+wxs,       c+wz*wz*oneminusc);
        } else {
            return Matrix<T>(1,0,0,0,1,0,0,0,1);
        }
    }

    // transpose of a matrix
    EXPRESSION_TEMPLATE
    inline Matrix<T,MATRIX,TransposeOp> 
    Transpose ( const MATRIX & m ) 
    { 
        return Matrix<T,MATRIX,TransposeOp>(m);
    }

    // dyadic product of two vectors
    EXPRESSION_TEMPLATE_PAIR
    inline Matrix<T,VECTOR1,DyadicOp,VECTOR2> 
    Dyadic ( const VECTOR1 & a, 
             const VECTOR2 & b) 
    {
        return Matrix<T,VECTOR1,DyadicOp,VECTOR2> (a,b);
    }
    
    //
    // Output
    //

    // vectors
    EXPRESSION_TEMPLATE
    inline OSTREAM & operator<< ( OSTREAM & o, 
                                  const VECTOR & v ) 
    {
        return o << v.template eval<0>() << " " 
                 << v.template eval<1>() << " " 
                 << v.template eval<2>();
    }

    // matrices
    EXPRESSION_TEMPLATE
    OSTREAM & operator<< ( OSTREAM & o, 
                           const MATRIX & m ) 
    {
        return o << "\n"<< m.template eval<0,0>() 
                 << " " << m.template eval<0,1>() 
                 << " " << m.template eval<0,2>()
                 << "\n"<< m.template eval<1,0>() 
                 << " " << m.template eval<1,1>() 
                 << " " << m.template eval<1,2>()
                 << "\n"<< m.template eval<2,0>() 
                 << " " << m.template eval<2,1>() 
                 << " " << m.template eval<2,2>() 
                 << "\n";
    }

#undef OSTREAM
#undef VECNRM
#undef VECNRM2
#undef VECPARENTHESIS
#undef MATNRM
#undef MATNRM2
#undef MATPARENTHESIS
#undef MATTR
#undef MATDET
#undef MATROW
#undef MATCOLUMN
#undef EXPRESSION_TEMPLATE     
#undef ANOTHER_EXPRESSION_TEMPLATE 
#undef EXPRESSION_TEMPLATE_PAIR   
#undef EXPRESSION_TEMPLATE_TRIPLET
#undef EXPRESSION_TEMPLATE_MEMBER
#undef CONVERTIBLE_EXPRESSION_TEMPLATE     
#undef CONVERTIBLE_ANOTHER_EXPRESSION_TEMPLATE 
#undef VECTOR 
#undef VECTOR1
#undef VECTOR2
#undef VECTOR3
#undef MATRIX 
#undef MATRIX1
#undef MATRIX2
#undef ENODE
}  // end namespace vecmat3

#ifndef NOVECMAT3DEF

// Backwards compatible type definitions of Vector and Matrix
// To not use these, state:
// #define NOVECMAT3DEF
// #include "vecmat3.h"

#ifndef DOUBLE
// default type of elements is double
#define DOUBLE double
#endif

typedef vecmat3::Vector<DOUBLE> Vector; 
typedef vecmat3::Matrix<DOUBLE> Matrix;

#endif

#endif

//
//  On the implementations of the template expressions:
// 
//  To have e.g. 'c = a+b' be unrolled to 
//  'c.x = a.x+b.x; c.y = a.y+b.y; c.z = a.z+b.z;',
//  an expression template class is defined with the following properties:
// 
//  - It contains pointers to the vectors 'a' and 'b'.
// 
//  - The constructor of this class sets up these pointers
// 
//  - It has an int template parameter indicating the type of operation.
//    The possible operations are listed in the unnamed 'enum' below.
// 
//  - It contains a function 'eval<I>()', where I is a integer
//    parameter, 0, 1 or 2, which returns the value of the Ith element
//    of the vector expression. [For matrix expressions this is replaced
//    by 'eval<I,J>()'.]
// 
//  - The 'operator+' is overloaded to return an instance of this class
// 
//  - The 'operator= ' is overloaded to call the 'eval<I>()' function.
// 
//  - Two class template-parameters are needed to specify the operant types.
//    This way one can distinguish different operations, e.g. vector-vector 
//    multiplication from vector-double multiplication. 
// 
//  - There are separate expression templates Vector-valued and Matrix-valued 
//    expressions, such that one can distinguish e.g. Matrix-Vector from 
//    Matrix-Matrix multiplication.
// 
//  The general expression templates are defined as follows:
//    template <typename T, class A, operation B, class C>
//    Vector<T,A,B,C>; 
//  and
//    template <typename T, class A, operation B, class C>
//    Matrix<T,A,B,C>;
//  
//  The general template classes are not defined, only special instances for 
//  allowed operations B. 
// 
//  The default values for the templates arguments are 'Base' for 'A'
//  and 'C', and 'NO_OPERATION' for 'B'. The templates with these
//  default values, denoted as 'Vector<T>' and 'Matrix<T>' serve as
//  the actual Vector and Matrix class. They are partially specialized
//  template classes in the vecmat3 namespace, with the template
//  parameter T signifying the data type.  For convenience and
//  backward compatibility, Vector<DOUBLE> and Matrix<DOUBLE> and are
//  'typedef'ed as the default Vector and Matrix classes.  The
//  specializations of these default templates therefore contain the
//  actual elements and much of the basic functionality of vectors and
//  matrices.
// 
//  Note that the classes A and C can themselved be expression classes,
//  and so nested expressions such as (a+b)*(d+2*c) are perfectly
//  possible.  This recursiveness does make the notation somewhat
//  involved.
// 
//  For further information on the technique of expression templates
//  in c++, see:
//
// -ubiety.uwaterloo.ca/~tveldhui/papers/Expression-Templates/expre_eval_l.html
//
// -www.oonumerics.org/blitz
//
// -tvmet.sourceforge.net
//