gee_lwo.hpp

Go to the documentation of this file.

#ifndef __VECTORMATRIX_HPP
#define __VECTORMATRIX_HPP

//changes history
//[021027 jussi] added vector comparisons with epsilon (uses radius, not component-wise threshold)
//[021004 jussi] added const return values to all functions returning objects
//[030128 jussi] added Matrix2x3 and Matrix3x2 classes
//[030218 jussi] fixed a bug in Matrix2 invert

// Review History:
// 020726 jan
// 020729 jan
// 020730 jan
// 020731 jan
// \todo [020726 jan] doxygen comments
// \todo [020726 jan] separate .inl file for inline functions
//       -jussi: no, because we want to have as few files as possible for easy ripping

// \todo [020726 jan] in general, write operator@(a,b) as { return T(a) @= b; }
//                                              - preferred by Stroustrup and Meyers
// jussi: NOT DONE: produces identical code at least on msvc (with optimizations),
//                  + proposed way is unintuitive + forces all member functions to return reference
//                    to *this (for example invert couldn't return a bool anymore)
// jan:                         I will check this out. It's most likely MSVC inliner doing the optimization
//
// \todo [020726 jan] consider private with public accessors for Vectors. Annoying, yes, but you can
//                                              always use operator[]
//                                              - allows data validity checks in debug build
// jussi: NOT DONE: VERY annoying to use, access to Vector.x/y/z is very intuitive and used all the time.
//                  what would be the data validity checks to be performed? check for NAN/INF?
//                                      what is the use of these tests? aren't they dependent on application and thus
//                                      couldn't be generic?
//      jan:                    NAN/DENORMAL/INF should be checked in (paranoid) debug build. Can be asserted in math functions as well.

//TODO check the code produced by m[][] access (they are used heavily internally)

//TODO is there a better name for vector::project that indicates which vector is projected on which?

//TODO should we have defaultMultiplicationOrder and defaultHandedness global variables to cut down
//     the number of parameters the user needs to give?
//  jan: no. #define DEFAULT_HANDEDNESS etc. could be defined prior header #include directive,
//                      but may cause grey hair. if not defined, then define it as you like


//TODO add Matrix::skewX,skewY and skewZ?

//TODO initialize Vectors in debug build to NaN






//MatrixRxC, R = number of rows, C = number of columns
//indexing: matrix[row][column]

// left multiply : M * columnVector
// right multiply: rowVector * N
// the matrices are related by M = N^T <=> M^T = N              (^T = transpose)
// for left multiply, the columns of the matrix are the basis vectors of the coordinate system
// for right multiply, the rows of the matrix are the basis vectors of the coordinate system

//matrices formed by RIGHT_MULTIPLY correspond to OpenGL input format (column major order)
//matrices formed by LEFT_MULTIPLY must be transposed before passing them to OpenGL
//in addition to this, OpenGL uses right-handed matrices (z-axis points towards the viewer)
// and our matrices are left-handed (z-axis points away from the viewer). this can be confusing when
// mixing OpenGL matrix operations with ours.
// calculating everything using our matrices (including the projection matrix)
// and passing the result (which is a left handed matrix) directly to OpenGL works.

//to convert between left and right handed coordinate systems, negate the z-axis (for left multiply
// it is the third column, for right multiply the third row)

#include <math.h>               //for sqrt, sin and cos
#include <stdlib.h>             //for rand
#include <assert.h>

namespace VECTORMATRIX  //TODO is this a proper name?
{

typedef float T;        //TODO is this enough, or do we need to use templates?

#ifndef ASSERT
#define ASSERT assert
#endif

inline T                                                RAND()                          { return T(rand()) * T(2) / T(RAND_MAX) - T(1); }       //random number [-1,1], used by VectorFactories
template<class TT> inline TT    OOSQRT( TT a )          { return TT(1) / TT(sqrt(a)); } //used by Vector::normalize
template<class TT> inline TT    SQRT( TT a )            { return TT(sqrt(a)); }                 //used by Vector::length
template<class TT> inline TT    SQR( TT a )                     { return a*a; }                                 //used by vector isEqual
template<class TT> inline TT    ABS( TT a )                     { return (a < TT(0)) ? -a : a; }//used by VectorFactory::perpendicular
template<class TT> inline TT    SIN( TT a )                     { return TT(sin(double(a))); }  //used by Matrix::rotate
template<class TT> inline TT    COS( TT a )                     { return TT(cos(double(a))); }  //used by Matrix::rotate
template<class TT> inline TT    TAN( TT a )                     { return TT(tan(double(a))); }  //used by MatrixFactory::frustum, frustum01
template<class TT> inline void  SWAP( TT &a, TT &b ){ TT tmp = a; a = b; b = tmp; }     //used by Matrix::transpose

enum MultiplicationOrder
{
        LEFT_MULTIPLY,          //for column vectors to be multiplied from the left
        RIGHT_MULTIPLY          //for row vectors to be multiplied from the right
};

class Matrix2;
class Matrix3;
class Matrix4;
class Matrix2x3;        //2 rows, 3 columns, for column vector transformations with implied 3rd row [0,0,1]
class Matrix3x2;        //3 rows, 2 columns, for row vector transformations with implied 3rd column [0,0,1]
class Matrix3x4;        //3 rows, 4 columns, for column vector transformations with implied 4th row [0,0,0,1]
class Matrix4x3;        //4 rows, 3 columns, for row vector transformations with implied 4th column [0,0,0,1]

class Vector2;
class Vector3;
class Vector4;

/*
design principles:
--------------------------------------------------------------------------
-orthogonal set of functions
 *except not all combinations of vector * matrix and matrix * matrix are present since these are mostly
   useless
 *except Matrix3x4 does not have right multiply since in practice it is used for left multiply matrices
   only. the same applies for Matrix4x3
 *different types of matrices can be multiplied through the use of constructors and explicit conversion
   functions

-follows mathematical notation
 *the order of multiplications was considered important. for example Matrix3x4, which has only left
   multiply, does not have operator*=

-allows row vectors with matrix multiplication from right, and column vectors with multiplication
  from left

-not dependent on proprietary code
 *only a few operations need to be defined
 *all of these are defined above, so replacing them with faster variants is easy

-usability was considered more important than speed
 *with less usable functions some temporaries could be avoided
   for example void cross(Vector& dest, const Vector& v1, const Vector& v2) instead of
             Vector cross( const Vector& v1, const Vector& v2 )
   but then it wouldn't be possible to write Vector a( normalize( cross( b, c ) ) );

-avoids redundancy
 *except some functions are both globals and members. there are no strict rules for this though.

-does not return references because Vector3 a; a.normalize().negate(); is bad style compared to
 Vector3 a; a = -normalize( a ); even though the former avoids the use of temporaries.
 note that it is still possible to write Vector3 a; a.normalize(); a.negate(); to avoid temporaries
 it would also be impossible to return a boolean telling if normalize and invert succeeded

-single header + cpp
 *easy to copy to other projects
   TODO for readability and ease of checking the available functions a separate .inl file could be handy.
        functions from the .cpp could be moved there?

-handling of mathematical singularities
 *no assertions, the user is expected to do the checks when necessary
 *different behaviour for members and globals (members return a boolean indicating success or failure)
 *if an operation fails, the operand is not changed
 *zero vector normalization returns a zero vector
   TODO should there be two versions of normalize, one that quietly outputs a zero vector if the input
        is a zero vector, and another that asserts that the input is non-zero?
 *projection on a zero vector returns a zero vector
   TODO vector::project: what the behaviour should be if the second vector is zero? (assertion or return 0?)
 *global inversion of a singular matrix returns the original matrix
   TODO should we assert that the input is not singular?

example of multiplication from right:
Matrix4x3 objectToWorld,cameraToWorld;
Matrix4x3 worldToCamera = invert(cameraToWorld);
Matrix4x3 objectToCamera = objectToWorld * worldToCamera;
Vector3 vertexPosition(1,2,3);
Vector4 result = Vector4(vertexPosition,1) * objectToCamera; //implicit conversion to Matrix4

example of multiplication from left:
Matrix3x4 objectToWorld,cameraToWorld;
Matrix3x4 worldToCamera = invert(cameraToWorld);
Matrix3x4 objectToCamera = worldToCamera * objectToWorld;
Vector3 vertexPosition(1,2,3);
Vector4 result = objectToCamera * Vector4(vertexPosition,1); //implicit conversion to Matrix4

transformations:
point p (Vector3) by Matrix4 m:                Vector4 r = Vector4(p,1) * m;
normal n (Vector3) by Matrix4 m:               Vector3 r = n * transpose(invert(m.getRotation()));
tangent or direction t (Vector3) by Matrix4 m: Vector3 r = t * m.getRotation();
*/

//==============================================================================================
//==============================================================================================

class Matrix2
{
public:
        inline                                  Matrix2                 ();                                             //initialize to identity
        inline                                  Matrix2                 ( const Matrix2& m );   //copy
        //construction from other types is explicit for safer use (see makeMatrix2 functions below)
        inline                                  Matrix2                 ( T m00, T m01, T m10, T m11 );
        inline                                  ~Matrix2                ();
        inline Matrix2&                 operator=               ( const Matrix2& m );
        inline Vector2&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector2&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m10, T m11 );
        inline const Vector2    getRow                  ( int i ) const;
        inline const Vector2    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector2& v );
        inline void                             setColumn               ( int i, const Vector2& v );
        inline void                             operator*=              ( const Matrix2& m );
        inline void                             multiply                ( const Matrix2& m, MultiplicationOrder o );
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix2& m );
        inline void                             operator-=              ( const Matrix2& m );
        inline const Matrix2    operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        inline T                                det                             () const;
        void                                    rotate                  ( T radians, MultiplicationOrder o );   //counterclockwise
        inline void                             scale                   ( const Vector2& v, MultiplicationOrder o );

private:
        T                                               matrix[2][2];
};

//==============================================================================================
//==============================================================================================

class Matrix3
{
public:
        inline                                  Matrix3                 ();                                             //initialize to identity
        inline                                  Matrix3                 ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3                 ( const Matrix2x3& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3                 ( const Matrix3x2& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3                 ( const Matrix3& m );   //copy
        //construction from other types is explicit for safer use (see makeMatrix3 functions below)
        inline                                  Matrix3                 ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22 );
        inline                                  ~Matrix3                ();
        inline Matrix3&                 operator=               ( const Matrix3& m );
        inline Vector3&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector3&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22 );
        inline const Vector3    getRow                  ( int i ) const;
        inline const Vector3    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector3& v );
        inline void                             setColumn               ( int i, const Vector3& v );
        inline void                             operator*=              ( const Matrix3& m );
        inline void                             multiply                ( const Matrix3& m, MultiplicationOrder o );
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix3& m );
        inline void                             operator-=              ( const Matrix3& m );
        inline const Matrix3    operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        inline T                                det                             () const;
        void                                    rotate                  ( T radians, const Vector3& aboutThis, MultiplicationOrder o );
        inline void                             scale                   ( const Vector3& v, MultiplicationOrder o );

private:
        T                                               matrix[3][3];
};

//==============================================================================================
//==============================================================================================

class Matrix4
{
public:
        inline                                  Matrix4                 ();                                             //initialize to identity
        inline                                  Matrix4                 ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix2x3& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix3x2& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix3& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix3x4& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix4x3& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4                 ( const Matrix4& m );   //copy
        inline                                  Matrix4                 ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33 );
        inline                                  ~Matrix4                ();
        inline Matrix4&                 operator=               ( const Matrix4& m );
        inline Vector4&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector4&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33 );
        inline const Vector4    getRow                  ( int i ) const;
        inline const Vector4    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector4& v );
        inline void                             setColumn               ( int i, const Vector4& v );
        inline const Matrix3    getRotation             () const;       //upper-left submatrix
        inline const Vector3    getTranslation  ( MultiplicationOrder o ) const;        //third row or column depending on the multiplication order
        inline void                             setRotation             ( const Matrix3& m );   //upper-left submatrix
        inline void                             setTranslation  ( const Vector3& v, MultiplicationOrder o );    //third row or column depending on the multiplication order
        inline void                             operator*=              ( const Matrix4& m );
        inline void                             multiply                ( const Matrix4& m, MultiplicationOrder o );
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix4& m );
        inline void                             operator-=              ( const Matrix4& m );
        inline const Matrix4    operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        T                                               det                             () const;
        void                                    rotate                  ( T radians, const Vector3& aboutThis, MultiplicationOrder o );
        inline void                             translate               ( const Vector3& v, MultiplicationOrder o );
        inline void                             scale                   ( const Vector3& v, MultiplicationOrder o );

private:
        T                                               matrix[4][4];
};

//==============================================================================================
//==============================================================================================

//2 rows, 3 columns, for column vector transformations with implied 3rd row [0,0,1]
//uses always left multiply
class Matrix2x3
{
public:
        inline                                  Matrix2x3               ();                                             //initialize to identity
        inline                                  Matrix2x3               ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix2x3               ( const Matrix2x3& m ); //copy
        //construction from other types is explicit for safer use (see makeMatrix2x3 functions below)
        inline                                  Matrix2x3               ( T m00, T m01, T m02, T m10, T m11, T m12 );
        inline                                  ~Matrix2x3              ();
        inline Matrix2x3&               operator=               ( const Matrix2x3& m );
        inline Vector3&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector3&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m02, T m10, T m11, T m12 );
        inline const Vector3    getRow                  ( int i ) const;
        inline const Vector2    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector3& v );
        inline void                             setColumn               ( int i, const Vector2& v );
        inline const Matrix2    getRotation             () const;
        inline const Vector2    getTranslation  () const;
        inline void                             setRotation             ( const Matrix2& m );
        inline void                             setTranslation  ( const Vector2& v );
        inline void                             multiply                ( const Matrix2x3& m ); // this = m * this (left multiply)
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix2x3& m );
        inline void                             operator-=              ( const Matrix2x3& m );
        inline const Matrix2x3  operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        T                                               det                             () const;
        void                                    rotate                  ( T radians );
        inline void                             translate               ( const Vector2& v );
        inline void                             scale                   ( const Vector2& v );

private:
        T                                               matrix[2][3];
};

//==============================================================================================
//==============================================================================================

//3 rows, 2 columns, for row vector transformations with implied 3rd column [0,0,1]
//uses always right multiply
class Matrix3x2
{
public:
        inline                                  Matrix3x2               ();                                             //initialize to identity
        inline                                  Matrix3x2               ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3x2               ( const Matrix3x2& m ); //copy
        //construction from other types is explicit for safer use (see makeMatrix3x2 functions below)
        inline                                  Matrix3x2               ( T m00, T m01, T m10, T m11, T m20, T m21 );
        inline                                  ~Matrix3x2              ();
        inline Matrix3x2&               operator=               ( const Matrix3x2& m );
        inline Vector2&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector2&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m10, T m11, T m20, T m21 );
        inline const Vector2    getRow                  ( int i ) const;
        inline const Vector3    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector2& v );
        inline void                             setColumn               ( int i, const Vector3& v );
        inline const Matrix2    getRotation             () const;
        inline const Vector2    getTranslation  () const;
        inline void                             setRotation             ( const Matrix2& m );
        inline void                             setTranslation  ( const Vector2& v );
        inline void                             operator*=              ( const Matrix3x2& m ); // this = this * m (right multiply)
        inline void                             multiply                ( const Matrix3x2& m ); // this = this * m (right multiply)
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix3x2& m );
        inline void                             operator-=              ( const Matrix3x2& m );
        inline const Matrix3x2  operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        inline T                                det                             () const;
        void                                    rotate                  ( T radians );
        inline void                             translate               ( const Vector2& v );
        inline void                             scale                   ( const Vector2& v );

private:
        T                                               matrix[3][2];
};

//==============================================================================================
//==============================================================================================

//3 rows, 4 columns, for column vector transformations with implied 4th row [0,0,0,1]
//uses always left multiply
class Matrix3x4
{
public:
        inline                                  Matrix3x4               ();                                             //initialize to identity
        inline                                  Matrix3x4               ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3x4               ( const Matrix2x3& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3x4               ( const Matrix3x2& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3x4               ( const Matrix3& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix3x4               ( const Matrix3x4& m ); //copy
        //construction from other types is explicit for safer use (see makeMatrix3x4 functions below)
        inline                                  Matrix3x4               ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23 );
        inline                                  ~Matrix3x4              ();
        inline Matrix3x4&               operator=               ( const Matrix3x4& m );
        inline Vector4&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector4&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23 );
        inline const Vector4    getRow                  ( int i ) const;
        inline const Vector3    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector4& v );
        inline void                             setColumn               ( int i, const Vector3& v );
        inline const Matrix3    getRotation             () const;
        inline const Vector3    getTranslation  () const;
        inline void                             setRotation             ( const Matrix3& m );
        inline void                             setTranslation  ( const Vector3& v );
        inline void                             multiply                ( const Matrix3x4& m ); // this = m * this (left multiply)
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix3x4& m );
        inline void                             operator-=              ( const Matrix3x4& m );
        inline const Matrix3x4  operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        T                                               det                             () const;
        void                                    rotate                  ( T radians, const Vector3& aboutThis );
        inline void                             translate               ( const Vector3& v );
        inline void                             scale                   ( const Vector3& v );

private:
        T                                               matrix[3][4];
};

//==============================================================================================
//==============================================================================================

//4 rows, 3 columns, for row vector transformations with implied 4th column [0,0,0,1]
//uses always right multiply
class Matrix4x3
{
public:
        inline                                  Matrix4x3               ();                                             //initialize to identity
        inline                                  Matrix4x3               ( const Matrix2& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4x3               ( const Matrix2x3& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4x3               ( const Matrix3x2& m ); //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4x3               ( const Matrix3& m );   //extend: 1 to the diagonal and 0 to off-diagonal
        inline                                  Matrix4x3               ( const Matrix4x3& m ); //copy
        //construction from other types is explicit for safer use (see makeMatrix4x3 functions below)
        inline                                  Matrix4x3               ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22, T m30, T m31, T m32 );
        inline                                  ~Matrix4x3              ();
        inline Matrix4x3&               operator=               ( const Matrix4x3& m );
        inline Vector3&                 operator[]              ( int i );                              //returns a row vector
        inline const Vector3&   operator[]              ( int i ) const;
        inline void                             set                             ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22, T m30, T m31, T m32 );
        inline const Vector3    getRow                  ( int i ) const;
        inline const Vector4    getColumn               ( int i ) const;
        inline void                             setRow                  ( int i, const Vector3& v );
        inline void                             setColumn               ( int i, const Vector4& v );
        inline const Matrix3    getRotation             () const;
        inline const Vector3    getTranslation  () const;
        inline void                             setRotation             ( const Matrix3& m );
        inline void                             setTranslation  ( const Vector3& v );
        inline void                             operator*=              ( const Matrix4x3& m ); // this = this * m (right multiply)
        inline void                             multiply                ( const Matrix4x3& m ); // this = this * m (right multiply)
        inline void                             operator*=              ( T f );
        inline void                             operator+=              ( const Matrix4x3& m );
        inline void                             operator-=              ( const Matrix4x3& m );
        inline const Matrix4x3  operator-               () const;
        inline void                             identity                ();
        inline void                             transpose               ();
        bool                                    invert                  ();     //if the matrix is singular, returns false and leaves it unmodified
        inline T                                det                             () const;
        void                                    rotate                  ( T radians, const Vector3& aboutThis );
        inline void                             translate               ( const Vector3& v );
        inline void                             scale                   ( const Vector3& v );

private:
        T                                               matrix[4][3];
};

//==============================================================================================
//==============================================================================================

class Vector2
{
public:
        //not initialized (avoids some of the cost of temporary vectors)
        inline                                  Vector2                 ()                                                              {}
        inline                                  Vector2                 ( const Vector2& v )                    { x = v.x; y = v.y; }
        inline                                  Vector2                 ( T fx, T fy )                                  { x = fx; y = fy; }
        inline                                  ~Vector2                ()                                                              {}
        inline Vector2&                 operator=               ( const Vector2& v )                    { x = v.x; y = v.y; return *this; }
        inline T&                               operator[]              ( int i )                                               { ASSERT(i>=0&&i<2); return (&x)[i]; }
        inline const T&                 operator[]              ( int i ) const                                 { ASSERT(i>=0&&i<2); return (&x)[i]; }
        inline void                             set                             ( T fx, T fy )                                  { x = fx; y = fy; }
        inline void                             operator*=              ( T f )                                                 { x *= f; y *= f; }
        //row vector * matrix (= RIGHT_MULTIPLY)
        inline void                             operator*=              ( const Matrix2& m )                    { T tx = x,ty = y; x = tx*m[0][0]+ty*m[1][0]; y = tx*m[0][1]+ty*m[1][1]; }
        //right multiply with implicit 1 in the third component
        inline void                             operator*=              ( const Matrix3x2& m )                  { T tx = x,ty = y; x = tx*m[0][0]+ty*m[1][0]+m[2][0]; y = tx*m[0][1]+ty*m[1][1]+m[2][1]; }
        //no default parameter for safer use
        inline void                             multiply                ( const Matrix2& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) { T tx = x,ty = y; x = tx*m[0][0]+ty*m[0][1]; y = tx*m[1][0]+ty*m[1][1]; } else { *this *= m; } }
        //right multiply with implicit 1 in the third component
        inline void                             multiply                ( const Matrix3x2& m )                  { T tx = x,ty = y; x = tx*m[0][0]+ty*m[1][0]+m[2][0]; y = tx*m[0][1]+ty*m[1][1]+m[2][1]; }
        //left multiply with implicit 1 in the third component
        inline void                             multiply                ( const Matrix2x3& m );
        inline void                             operator+=              ( const Vector2& v )                    { x += v.x; y += v.y; }
        inline void                             operator-=              ( const Vector2& v )                    { x -= v.x; y -= v.y; }
        inline const Vector2    operator-               () const                                                { return Vector2(-x,-y); }
        //if the vector is zero, returns false and leaves it unmodified
        inline bool                             normalize               ()                                                              { T l = x*x+y*y; if( l == T(0) ) return false; l = OOSQRT(l); x *= l; y *= l; return true; }
        inline T                                length                  () const                                                { return SQRT(x*x+y*y); }
        inline void                             scale                   ( const Vector2& v )                    { x *= v.x; y *= v.y; } //component-wise scale
        inline void                             negate                  ()                                                              { x = -x; y = -y; }

        T                                               x,y;
};

//==============================================================================================
//==============================================================================================

class Vector3
{
public:
        //not initialized (avoids some of the cost of temporary vectors)
        inline                                  Vector3                 ()                                                              {}
        inline                                  Vector3                 ( const Vector2& v, T fz )              { x = v.x; y = v.y; z = fz; }
        inline                                  Vector3                 ( const Vector3& v )                    { x = v.x; y = v.y; z = v.z; }
        inline                                  Vector3                 ( T fx, T fy, T fz )                    { x = fx; y = fy; z = fz; }
        inline                                  ~Vector3                ()                                                              {}
        inline Vector3&                 operator=               ( const Vector3& v )                    { x = v.x; y = v.y; z = v.z; return *this; }
        inline T&                               operator[]              ( int i )                                               { ASSERT(i>=0&&i<3); return (&x)[i]; }
        inline const T&                 operator[]              ( int i ) const                                 { ASSERT(i>=0&&i<3); return (&x)[i]; }
        inline void                             set                             ( T fx, T fy, T fz )                    { x = fx; y = fy; z = fz; }
        inline void                             operator*=              ( T f )                                                 { x *= f; y *= f; z *= f; }
        //row vector * matrix (= RIGHT_MULTIPLY)
        inline void                             operator*=              ( const Matrix3& m )                    { T tx = x,ty = y,tz = z; x = tx*m[0][0]+ty*m[1][0]+tz*m[2][0]; y = tx*m[0][1]+ty*m[1][1]+tz*m[2][1]; z = tx*m[0][2]+ty*m[1][2]+tz*m[2][2]; }
        //right multiply with implicit 1 in the fourth component
        inline void                             operator*=              ( const Matrix4x3& m )                  { T tx = x,ty = y,tz = z; x = tx*m[0][0]+ty*m[1][0]+tz*m[2][0]+m[3][0]; y = tx*m[0][1]+ty*m[1][1]+tz*m[2][1]+m[3][1]; z = tx*m[0][2]+ty*m[1][2]+tz*m[2][2]+m[3][2]; }
        //no default parameter for safer use
        inline void                             multiply                ( const Matrix3& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) { T tx = x,ty = y,tz = z; x = tx*m[0][0]+ty*m[0][1]+tz*m[0][2]; y = tx*m[1][0]+ty*m[1][1]+tz*m[1][2]; z = tx*m[2][0]+ty*m[2][1]+tz*m[2][2]; } else { *this *= m; } }
        //right multiply with implicit 1 in the fourth component
        inline void                             multiply                ( const Matrix4x3& m )                  { T tx = x,ty = y,tz = z; x = tx*m[0][0]+ty*m[1][0]+tz*m[2][0]+m[3][0]; y = tx*m[0][1]+ty*m[1][1]+tz*m[2][1]+m[3][1]; z = tx*m[0][2]+ty*m[1][2]+tz*m[2][2]+m[3][2]; }
        //left multiply with implicit 1 in the fourth component
        inline void                             multiply                ( const Matrix3x4& m );
        inline void                             operator+=              ( const Vector3& v )                    { x += v.x; y += v.y; z += v.z; }
        inline void                             operator-=              ( const Vector3& v )                    { x -= v.x; y -= v.y; z -= v.z; }
        inline const Vector3    operator-               () const                                                { return Vector3(-x,-y,-z); }
        //if the vector is zero, returns false and leaves it unmodified
        inline bool                             normalize               ()                                                              { T l = x*x+y*y+z*z; if( l == T(0) ) return false; l = OOSQRT(l); x *= l; y *= l; z *= l; return true; }
        inline T                                length                  () const                                                { return SQRT(x*x+y*y+z*z); }
        inline void                             scale                   ( const Vector3& v )                    { x *= v.x; y *= v.y; z *= v.z; }       //component-wise scale
        inline void                             negate                  ()                                                              { x = -x; y = -y; z = -z; }

        T                                               x,y,z;
};

//==============================================================================================
//==============================================================================================

class Vector4
{
public:
        //not initialized (avoids some of the cost of temporary vectors)
        inline                                  Vector4                 ()                                                              {}
        inline                                  Vector4                 ( const Vector2& v, T fz, T fw ){ x = v.x; y = v.y; z = fz; w = fw; }
        inline                                  Vector4                 ( const Vector3& v, T fw )              { x = v.x; y = v.y; z = v.z; w = fw; }
        inline                                  Vector4                 ( const Vector4& v )                    { x = v.x; y = v.y; z = v.z; w = v.w; }
        inline                                  Vector4                 ( T fx, T fy, T fz, T fw )              { x = fx; y = fy; z = fz; w = fw; }
        inline                                  ~Vector4                ()                                                              {}
        inline Vector4&                 operator=               ( const Vector4& v )                    { x = v.x; y = v.y; z = v.z; w = v.w; return *this; }
        inline T&                               operator[]              ( int i )                                               { ASSERT(i>=0&&i<4); return (&x)[i]; }
        inline const T&                 operator[]              ( int i ) const                                 { ASSERT(i>=0&&i<4); return (&x)[i]; }
        inline void                             set                             ( T fx, T fy, T fz, T fw )              { x = fx; y = fy; z = fz; w = fw; }
        inline void                             operator*=              ( T f )                                                 { x *= f; y *= f; z *= f; w *= f; }
        //row vector * matrix (= RIGHT_MULTIPLY)
        inline void                             operator*=              ( const Matrix4& m )                    { T tx = x,ty = y,tz = z,tw = w; x = tx*m[0][0]+ty*m[1][0]+tz*m[2][0]+tw*m[3][0]; y = tx*m[0][1]+ty*m[1][1]+tz*m[2][1]+tw*m[3][1]; z = tx*m[0][2]+ty*m[1][2]+tz*m[2][2]+tw*m[3][2]; w = tx*m[0][3]+ty*m[1][3]+tz*m[2][3]+tw*m[3][3]; }
        //no default parameter for safer use
        inline void                             multiply                ( const Matrix4& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) { T tx = x,ty = y,tz = z,tw = w; x = tx*m[0][0]+ty*m[0][1]+tz*m[0][2]+tw*m[0][3]; y = tx*m[1][0]+ty*m[1][1]+tz*m[1][2]+tw*m[1][3]; z = tx*m[2][0]+ty*m[2][1]+tz*m[2][2]+tw*m[2][3]; w = tx*m[3][0]+ty*m[3][1]+tz*m[3][2]+tw*m[3][3]; } else { *this *= m; } }
        inline void                             operator+=              ( const Vector4& v )                    { x += v.x; y += v.y; z += v.z; w += v.w; }
        inline void                             operator-=              ( const Vector4& v )                    { x -= v.x; y -= v.y; z -= v.z; w -= v.w; }
        inline const Vector4    operator-               () const                                                { return Vector4(-x,-y,-z,-w); }
        //if the vector is zero, returns false and leaves it unmodified
        inline bool                             normalize               ()                                                              { T l = x*x+y*y+z*z+w*w; if( l == T(0) ) return false; l = OOSQRT(l); x *= l; y *= l; z *= l; w *= l; return true; }
        inline T                                length                  () const                                                { return SQRT(x*x+y*y+z*z+w*w); }
        inline void                             scale                   ( const Vector4& v )                    { x *= v.x; y *= v.y; z *= v.z; w *= v.w; }     //component-wise scale
        inline void                             negate                  ()                                                              { x = -x; y = -y; z = -z; w = -w; }

        T                                               x,y,z,w;
};

//==============================================================================================
//==============================================================================================

//Vector2 global functions
inline bool                             operator==      ( const Vector2& v1, const Vector2& v2 )        { return (v1.x == v2.x) && (v1.y == v2.y); }
inline bool                             operator!=      ( const Vector2& v1, const Vector2& v2 )        { return (v1.x != v2.x) || (v1.y != v2.y); }
inline bool                             isEqual         ( const Vector2& v1, const Vector2& v2, T epsilon )     { return SQR(v2.x-v1.x) + SQR(v2.y-v1.y) <= epsilon*epsilon; }
inline const Vector2    operator*       ( T f, const Vector2& v )                                       { return Vector2(v.x*f,v.y*f); }
inline const Vector2    operator*       ( const Vector2& v, T f )                                       { return Vector2(v.x*f,v.y*f); }
inline const Vector2    operator+       ( const Vector2& v1, const Vector2& v2 )        { return Vector2(v1.x+v2.x, v1.y+v2.y); }
inline const Vector2    operator-       ( const Vector2& v1, const Vector2& v2 )        { return Vector2(v1.x-v2.x, v1.y-v2.y); }
inline T                                dot                     ( const Vector2& v1, const Vector2& v2 )        { return v1.x*v2.x+v1.y*v2.y; }
//if v is a zero vector, returns a zero vector
inline const Vector2    normalize       ( const Vector2& v )                                            { T l = dot(v,v); if( l != T(0) ) l = OOSQRT(l); return v * l; }
//if onThis is a zero vector, returns a zero vector
inline const Vector2    project         ( const Vector2& v, const Vector2& onThis ) { T l = dot(onThis,onThis); if( l != T(0) ) l = dot(v, onThis)/l; return onThis * l; }
inline const Vector2    lerp            ( const Vector2& v1, const Vector2& v2, T ratio )       { return v1 + ratio * (v2 - v1); }
inline const Vector2    scale           ( const Vector2& v1, const Vector2& v2 )        { return Vector2(v1.x*v2.x, v1.y*v2.y); }

//==============================================================================================
//==============================================================================================

//Vector3 global functions
inline bool                             operator==      ( const Vector3& v1, const Vector3& v2 )        { return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z); }
inline bool                             operator!=      ( const Vector3& v1, const Vector3& v2 )        { return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z); }
inline bool                             isEqual         ( const Vector3& v1, const Vector3& v2, T epsilon )     { return SQR(v2.x-v1.x) + SQR(v2.y-v1.y) + SQR(v2.z-v1.z) <= epsilon*epsilon; }
inline const Vector3    operator*       ( T f, const Vector3& v )                                       { return Vector3(v.x*f,v.y*f,v.z*f); }
inline const Vector3    operator*       ( const Vector3& v, T f )                                       { return Vector3(v.x*f,v.y*f,v.z*f); }
inline const Vector3    operator+       ( const Vector3& v1, const Vector3& v2 )        { return Vector3(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z); }
inline const Vector3    operator-       ( const Vector3& v1, const Vector3& v2 )        { return Vector3(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z); }
inline T                                dot                     ( const Vector3& v1, const Vector3& v2 )        { return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z; }
inline const Vector3    cross           ( const Vector3& v1, const Vector3& v2 )        { return Vector3( v1.y*v2.z-v1.z*v2.y, v1.z*v2.x-v1.x*v2.z, v1.x*v2.y-v1.y*v2.x ); }
//if v is a zero vector, returns a zero vector
inline const Vector3    normalize       ( const Vector3& v )                                            { T l = dot(v,v); if( l != T(0) ) l = OOSQRT(l); return v * l; }
//if onThis is a zero vector, returns a zero vector
inline const Vector3    project         ( const Vector3& v, const Vector3& onThis ) { T l = dot(onThis,onThis); if( l != T(0) ) l = dot(v, onThis)/l; return onThis * l; }
inline const Vector3    lerp            ( const Vector3& v1, const Vector3& v2, T ratio )       { return v1 + ratio * (v2 - v1); }
inline const Vector3    scale           ( const Vector3& v1, const Vector3& v2 )        { return Vector3(v1.x*v2.x, v1.y*v2.y, v1.z*v2.z); }

//==============================================================================================
//==============================================================================================

//Vector4 global functions
inline bool                             operator==      ( const Vector4& v1, const Vector4& v2 )        { return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z) && (v1.w == v2.w); }
inline bool                             operator!=      ( const Vector4& v1, const Vector4& v2 )        { return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z) || (v1.w != v2.w); }
inline bool                             isEqual         ( const Vector4& v1, const Vector4& v2, T epsilon )     { return SQR(v2.x-v1.x) + SQR(v2.y-v1.y) + SQR(v2.z-v1.z) + SQR(v2.w-v1.w) <= epsilon*epsilon; }
inline const Vector4    operator*       ( T f, const Vector4& v )                                       { return Vector4(v.x*f,v.y*f,v.z*f,v.w*f); }
inline const Vector4    operator*       ( const Vector4& v, T f )                                       { return Vector4(v.x*f,v.y*f,v.z*f,v.w*f); }
inline const Vector4    operator+       ( const Vector4& v1, const Vector4& v2 )        { return Vector4(v1.x+v2.x, v1.y+v2.y, v1.z+v2.z, v1.w+v2.w); }
inline const Vector4    operator-       ( const Vector4& v1, const Vector4& v2 )        { return Vector4(v1.x-v2.x, v1.y-v2.y, v1.z-v2.z, v1.w-v2.w); }
inline T                                dot                     ( const Vector4& v1, const Vector4& v2 )        { return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z+v1.w*v2.w; }
//if v is a zero vector, returns a zero vector
inline const Vector4    normalize       ( const Vector4& v )                                            { T l = dot(v,v); if( l != T(0) ) l = OOSQRT(l); return v * l; }
//if onThis is a zero vector, returns a zero vector
inline const Vector4    project         ( const Vector4& v, const Vector4& onThis ) { T l = dot(onThis,onThis); if( l != T(0) ) l = dot(v, onThis)/l; return onThis * l; }
inline const Vector4    lerp            ( const Vector4& v1, const Vector4& v2, T ratio )       { return v1 + ratio * (v2 - v1); }
inline const Vector4    scale           ( const Vector4& v1, const Vector4& v2 )        { return Vector4(v1.x*v2.x, v1.y*v2.y, v1.z*v2.z, v1.w*v2.w); }

//==============================================================================================
//==============================================================================================

//Vector2 inline functions (cannot be defined in the class because of dependency issues)
inline void                             Vector2::multiply               ( const Matrix2x3& m )                  { T tx = x,ty = y; x = tx*m[0][0]+ty*m[0][1]+m[0][2]; y = tx*m[1][0]+ty*m[1][1]+m[1][2]; }

//global Vector2-Matrix2 multiplications
//row vector * matrix
inline const Vector2    operator*       ( const Vector2& v, const Matrix2& m)           { return Vector2( v.x*m[0][0]+v.y*m[1][0], v.x*m[0][1]+v.y*m[1][1] ); }
//row vector * matrix with implied 1 in the third component of the vector
inline const Vector2    operator*       ( const Vector2& v, const Matrix3x2& m)         { return Vector2( v.x*m[0][0]+v.y*m[1][0]+m[2][0], v.x*m[0][1]+v.y*m[1][1]+m[2][1] ); }
//matrix * column vector
inline const Vector2    operator*       ( const Matrix2& m, const Vector2& v)           { return Vector2( v.x*m[0][0]+v.y*m[0][1], v.x*m[1][0]+v.y*m[1][1] ); }
//matrix * column vector with implied 1 in the third component of the vector
inline const Vector2    operator*       ( const Matrix2x3& m, const Vector2& v)         { return Vector2( v.x*m[0][0]+v.y*m[0][1]+m[0][2], v.x*m[1][0]+v.y*m[1][1]+m[1][2] ); }

//==============================================================================================
//==============================================================================================

//Vector3 inline functions (cannot be defined in the class because of dependency issues)
inline void                             Vector3::multiply               ( const Matrix3x4& m )                  { T tx = x,ty = y,tz = z; x = tx*m[0][0]+ty*m[0][1]+tz*m[0][2]+m[0][3]; y = tx*m[1][0]+ty*m[1][1]+tz*m[1][2]+m[1][3]; z = tx*m[2][0]+ty*m[2][1]+tz*m[2][2]+m[2][3]; }

//global Vector3-Matrix3 multiplications
//row vector * matrix
inline const Vector3    operator*       ( const Vector3& v, const Matrix3& m)           { return Vector3( v.x*m[0][0]+v.y*m[1][0]+v.z*m[2][0], v.x*m[0][1]+v.y*m[1][1]+v.z*m[2][1], v.x*m[0][2]+v.y*m[1][2]+v.z*m[2][2] ); }
//row vector * matrix with implied 1 in the fourth component of the vector
inline const Vector3    operator*       ( const Vector3& v, const Matrix4x3& m)         { return Vector3( v.x*m[0][0]+v.y*m[1][0]+v.z*m[2][0]+m[3][0], v.x*m[0][1]+v.y*m[1][1]+v.z*m[2][1]+m[3][1], v.x*m[0][2]+v.y*m[1][2]+v.z*m[2][2]+m[3][2] ); }
//matrix * column vector
inline const Vector3    operator*       ( const Matrix3& m, const Vector3& v)           { return Vector3( v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2], v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2], v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2] ); }
//matrix * column vector with implied 1 in the fourth component of the vector
inline const Vector3    operator*       ( const Matrix3x4& m, const Vector3& v)         { return Vector3( v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2]+m[0][3], v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2]+m[1][3], v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2]+m[2][3] ); }

//==============================================================================================
//==============================================================================================

//global Vector4-Matrix4 multiplications
//row vector * matrix
inline const Vector4    operator*       ( const Vector4& v, const Matrix4& m)           { return Vector4( v.x*m[0][0]+v.y*m[1][0]+v.z*m[2][0]+v.w*m[3][0], v.x*m[0][1]+v.y*m[1][1]+v.z*m[2][1]+v.w*m[3][1], v.x*m[0][2]+v.y*m[1][2]+v.z*m[2][2]+v.w*m[3][2], v.x*m[0][3]+v.y*m[1][3]+v.z*m[2][3]+v.w*m[3][3] ); }
//matrix * column vector
inline const Vector4    operator*       ( const Matrix4& m, const Vector4& v)           { return Vector4( v.x*m[0][0]+v.y*m[0][1]+v.z*m[0][2]+v.w*m[0][3], v.x*m[1][0]+v.y*m[1][1]+v.z*m[1][2]+v.w*m[1][3], v.x*m[2][0]+v.y*m[2][1]+v.z*m[2][2]+v.w*m[2][3], v.x*m[3][0]+v.y*m[3][1]+v.z*m[3][2]+v.w*m[3][3] ); }

//==============================================================================================
//==============================================================================================

//Matrix2 global functions
inline bool                             operator==      ( const Matrix2& m1, const Matrix2& m2 )        { for(int i=0;i<2;i++) for(int j=0;j<2;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix2& m1, const Matrix2& m2 )        { return !(m1 == m2); }
inline const Matrix2    operator*       ( const Matrix2& m1, const Matrix2& m2 )        { Matrix2 t; for(int i=0;i<2;i++) for(int j=0;j<2;j++) t[i][j] = m1[i][0] * m2[0][j] + m1[i][1] * m2[1][j]; return t; }
inline const Matrix2    operator*       ( T f, const Matrix2& m )                                       { Matrix2 t(m); t *= f; return t; }
inline const Matrix2    operator*       ( const Matrix2& m, T f )                                       { Matrix2 t(m); t *= f; return t; }
inline const Matrix2    operator+       ( const Matrix2& m1, const Matrix2& m2 )        { Matrix2 t(m1); t += m2; return t; }
inline const Matrix2    operator-       ( const Matrix2& m1, const Matrix2& m2 )        { Matrix2 t(m1); t -= m2; return t; }
inline const Matrix2    transpose       ( const Matrix2& m )                                            { Matrix2 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix2    invert          ( const Matrix2& m )                                            { Matrix2 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix2 inline functions (cannot be inside the class because Vector2 is not defined yet when Matrix2 is defined)
inline                                  Matrix2::Matrix2        ()                                                                      { identity(); }
inline                                  Matrix2::Matrix2        ( const Matrix2& m )                            { *this = m; }
inline                                  Matrix2::Matrix2        ( T m00, T m01, T m10, T m11 )          { set(m00,m01,m10,m11); }
inline                                  Matrix2::~Matrix2       ()                                                                      {}
inline Matrix2&                 Matrix2::operator=      ( const Matrix2& m )                            { for(int i=0;i<2;i++) for(int j=0;j<2;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector2&                 Matrix2::operator[]     ( int i )                                                       { ASSERT(i>=0&&i<2); return (Vector2&)matrix[i][0]; }
inline const Vector2&   Matrix2::operator[]     ( int i ) const                                         { ASSERT(i>=0&&i<2); return (const Vector2&)matrix[i][0]; }
inline void                             Matrix2::set            ( T m00, T m01, T m10, T m11 )          { matrix[0][0] = m00; matrix[0][1] = m01; matrix[1][0] = m10; matrix[1][1] = m11; }
inline const Vector2    Matrix2::getRow         ( int i ) const                                         { ASSERT(i>=0&&i<2); return Vector2(matrix[i][0],matrix[i][1]); }
inline const Vector2    Matrix2::getColumn      ( int i ) const                                         { ASSERT(i>=0&&i<2); return Vector2(matrix[0][i],matrix[1][i]); }
inline void                             Matrix2::setRow         ( int i, const Vector2& v )                     { ASSERT(i>=0&&i<2); matrix[i][0] = v.x; matrix[i][1] = v.y; }
inline void                             Matrix2::setColumn      ( int i, const Vector2& v )                     { ASSERT(i>=0&&i<2); matrix[0][i] = v.x; matrix[1][i] = v.y; }
inline void                             Matrix2::operator*=     ( const Matrix2& m )                            { *this = *this * m; }
inline void                             Matrix2::multiply       ( const Matrix2& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) *this = m * (*this); else *this *= m; }
inline void                             Matrix2::operator*=     ( T f )                                                         { for(int i=0;i<2;i++) for(int j=0;j<2;j++) matrix[i][j] *= f; }
inline void                             Matrix2::operator+=     ( const Matrix2& m )                            { for(int i=0;i<2;i++) for(int j=0;j<2;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix2::operator-=     ( const Matrix2& m )                            { for(int i=0;i<2;i++) for(int j=0;j<2;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix2    Matrix2::operator-      () const                                                        { return Matrix2( -matrix[0][0],-matrix[0][1],-matrix[1][0],-matrix[1][1]); }
inline void                             Matrix2::identity       ()                                                                      { for(int i=0;i<2;i++) for(int j=0;j<2;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix2::transpose      ()                                                                      { SWAP(matrix[1][0], matrix[0][1]); }
inline T                                Matrix2::det            () const                                                        { return matrix[0][0] * matrix[1][1] - matrix[1][0]*matrix[0][1]; }
inline void                             Matrix2::scale          ( const Vector2& v, MultiplicationOrder o ) { multiply( Matrix2( v.x, T(0), T(0),  v.y ), o ); }

//==============================================================================================
//==============================================================================================

//Matrix3 global functions
inline bool                             operator==      ( const Matrix3& m1, const Matrix3& m2 )        { for(int i=0;i<3;i++) for(int j=0;j<3;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix3& m1, const Matrix3& m2 )        { return !(m1 == m2); }
inline const Matrix3    operator*       ( const Matrix3& m1, const Matrix3& m2 )        { Matrix3 t; for(int i=0;i<3;i++) for(int j=0;j<3;j++) t[i][j] = m1[i][0] * m2[0][j] + m1[i][1] * m2[1][j] + m1[i][2] * m2[2][j]; return t; }
inline const Matrix3    operator*       ( T f, const Matrix3& m )                                       { Matrix3 t(m); t *= f; return t; }
inline const Matrix3    operator*       ( const Matrix3& m, T f )                                       { Matrix3 t(m); t *= f; return t; }
inline const Matrix3    operator+       ( const Matrix3& m1, const Matrix3& m2 )        { Matrix3 t(m1); t += m2; return t; }
inline const Matrix3    operator-       ( const Matrix3& m1, const Matrix3& m2 )        { Matrix3 t(m1); t -= m2; return t; }
inline const Matrix3    transpose       ( const Matrix3& m )                                            { Matrix3 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix3    invert          ( const Matrix3& m )                                            { Matrix3 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix3 inline functions (cannot be inside the class because Vector3 is not defined yet when Matrix3 is defined)
inline                                  Matrix3::Matrix3        ()                                                                      { identity(); }
inline                                  Matrix3::Matrix3        ( const Matrix2& m )                            { set( m[0][0], m[0][1], T(0), m[1][0], m[1][1], T(0), T(0), T(0), T(1) ); }
inline                                  Matrix3::Matrix3        ( const Matrix2x3& m )                          { set( m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], T(0), T(0), T(1) ); }
inline                                  Matrix3::Matrix3        ( const Matrix3x2& m )                          { set( m[0][0], m[0][1], T(0), m[1][0], m[1][1], T(0), m[2][0], m[2][1], T(1) ); }
inline                                  Matrix3::Matrix3        ( const Matrix3& m )                            { *this = m; }
inline                                  Matrix3::Matrix3        ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22 )       { set(m00,m01,m02,m10,m11,m12,m20,m21,m22); }
inline                                  Matrix3::~Matrix3       ()                                                                      {}
inline Matrix3&                 Matrix3::operator=      ( const Matrix3& m )                            { for(int i=0;i<3;i++) for(int j=0;j<3;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector3&                 Matrix3::operator[]     ( int i )                                                       { ASSERT(i>=0&&i<3); return (Vector3&)matrix[i][0]; }
inline const Vector3&   Matrix3::operator[]     ( int i ) const                                         { ASSERT(i>=0&&i<3); return (const Vector3&)matrix[i][0]; }
inline void                             Matrix3::set            ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22 ) {     matrix[0][0] = m00; matrix[0][1] = m01; matrix[0][2] = m02; matrix[1][0] = m10; matrix[1][1] = m11; matrix[1][2] = m12; matrix[2][0] = m20; matrix[2][1] = m21; matrix[2][2] = m22; }
inline const Vector3    Matrix3::getRow         ( int i ) const                                         { ASSERT(i>=0&&i<3); return Vector3(matrix[i][0],matrix[i][1],matrix[i][2]); }
inline const Vector3    Matrix3::getColumn      ( int i ) const                                         { ASSERT(i>=0&&i<3); return Vector3(matrix[0][i],matrix[1][i],matrix[2][i]); }
inline void                             Matrix3::setRow         ( int i, const Vector3& v )                     { ASSERT(i>=0&&i<3); matrix[i][0] = v.x; matrix[i][1] = v.y; matrix[i][2] = v.z; }
inline void                             Matrix3::setColumn      ( int i, const Vector3& v )                     { ASSERT(i>=0&&i<3); matrix[0][i] = v.x; matrix[1][i] = v.y; matrix[2][i] = v.z; }
inline void                             Matrix3::operator*=     ( const Matrix3& m )                            { *this = *this * m; }
inline void                             Matrix3::multiply       ( const Matrix3& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) *this = m * (*this); else *this *= m; }
inline void                             Matrix3::operator*=     ( T f )                                                         { for(int i=0;i<3;i++) for(int j=0;j<3;j++) matrix[i][j] *= f; }
inline void                             Matrix3::operator+=     ( const Matrix3& m )                            { for(int i=0;i<3;i++) for(int j=0;j<3;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix3::operator-=     ( const Matrix3& m )                            { for(int i=0;i<3;i++) for(int j=0;j<3;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix3    Matrix3::operator-      () const                                                        { return Matrix3( -matrix[0][0],-matrix[0][1],-matrix[0][2], -matrix[1][0],-matrix[1][1],-matrix[1][2], -matrix[2][0],-matrix[2][1],-matrix[2][2]); }
inline void                             Matrix3::identity       ()                                                                      { for(int i=0;i<3;i++) for(int j=0;j<3;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix3::transpose      ()                                                                      { SWAP(matrix[1][0], matrix[0][1]); SWAP(matrix[2][0], matrix[0][2]); SWAP(matrix[2][1], matrix[1][2]); }
inline T                                Matrix3::det            () const                                                        { return matrix[0][0] * (matrix[1][1]*matrix[2][2] - matrix[2][1]*matrix[1][2]) + matrix[0][1] * (matrix[2][0]*matrix[1][2] - matrix[1][0]*matrix[2][2]) + matrix[0][2] * (matrix[1][0]*matrix[2][1] - matrix[2][0]*matrix[1][1]); }
inline void                             Matrix3::scale          ( const Vector3& v, MultiplicationOrder o )     { multiply( Matrix3( v.x, T(0), T(0), T(0), v.y, T(0), T(0), T(0), v.z ), o ); }

//==============================================================================================
//==============================================================================================

//Matrix4 global functions
inline bool                             operator==      ( const Matrix4& m1, const Matrix4& m2 )        { for(int i=0;i<4;i++) for(int j=0;j<4;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix4& m1, const Matrix4& m2 )        { return !(m1 == m2); }
inline const Matrix4    operator*       ( const Matrix4& m1, const Matrix4& m2 )        { Matrix4 t; for(int i=0;i<4;i++) for(int j=0;j<4;j++) t[i][j] = m1[i][0] * m2[0][j] + m1[i][1] * m2[1][j] + m1[i][2] * m2[2][j] + m1[i][3] * m2[3][j]; return t; }
inline const Matrix4    operator*       ( T f, const Matrix4& m )                                       { Matrix4 t(m); t *= f; return t; }
inline const Matrix4    operator*       ( const Matrix4& m, T f )                                       { Matrix4 t(m); t *= f; return t; }
inline const Matrix4    operator+       ( const Matrix4& m1, const Matrix4& m2 )        { Matrix4 t(m1); t += m2; return t; }
inline const Matrix4    operator-       ( const Matrix4& m1, const Matrix4& m2 )        { Matrix4 t(m1); t -= m2; return t; }
inline const Matrix4    transpose       ( const Matrix4& m )                                            { Matrix4 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix4    invert          ( const Matrix4& m )                                            { Matrix4 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix4 inline functions (cannot be inside the class because Vector4 is not defined yet when Matrix4 is defined)
inline                                  Matrix4::Matrix4                ()                                                                      { identity(); }
inline                                  Matrix4::Matrix4                ( const Matrix2& m )                            { set(m[0][0],m[0][1],T(0),T(0),m[1][0],m[1][1],T(0),T(0),T(0),T(0),T(1),T(0),T(0),T(0),T(0),T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix2x3& m )                          { set(m[0][0],m[0][1],m[0][2],T(0),m[1][0],m[1][1],m[1][2],T(0),T(0),T(0),T(1),T(0),T(0),T(0),T(0),T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix3x2& m )                          { set(m[0][0],m[0][1],T(0),T(0),m[1][0],m[1][1],T(0),T(0),m[2][0],m[2][1],T(1),T(0),T(0),T(0),T(0),T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix3& m )                            { set(m[0][0],m[0][1],m[0][2],T(0),m[1][0],m[1][1],m[1][2],T(0),m[2][0],m[2][1],m[2][2],T(0),T(0),T(0),T(0),T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix3x4& m )                          { set(m[0][0],m[0][1],m[0][2],m[0][3],m[1][0],m[1][1],m[1][2],m[1][3],m[2][0],m[2][1],m[2][2],m[2][3],T(0),T(0),T(0),T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix4x3& m )                          { set(m[0][0],m[0][1],m[0][2],T(0),m[1][0],m[1][1],m[1][2],T(0),m[2][0],m[2][1],m[2][2],T(0),m[3][0],m[3][1],m[3][2],T(1)); }
inline                                  Matrix4::Matrix4                ( const Matrix4& m )                            { *this = m; }
inline                                  Matrix4::Matrix4                ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33 )      { set(m00,m01,m02,m03,m10,m11,m12,m13,m20,m21,m22,m23,m30,m31,m32,m33); }
inline                                  Matrix4::~Matrix4               ()                                                                      {}
inline Matrix4&                 Matrix4::operator=              ( const Matrix4& m )                            { for(int i=0;i<4;i++) for(int j=0;j<4;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector4&                 Matrix4::operator[]             ( int i )                                                       { ASSERT(i>=0&&i<4); return (Vector4&)matrix[i][0]; }
inline const Vector4&   Matrix4::operator[]             ( int i ) const                                         { ASSERT(i>=0&&i<4); return (const Vector4&)matrix[i][0]; }
inline void                             Matrix4::set                    ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33  ) { matrix[0][0] = m00; matrix[0][1] = m01; matrix[0][2] = m02; matrix[0][3] = m03; matrix[1][0] = m10; matrix[1][1] = m11; matrix[1][2] = m12;  matrix[1][3] = m13; matrix[2][0] = m20; matrix[2][1] = m21; matrix[2][2] = m22; matrix[2][3] = m23; matrix[3][0] = m30; matrix[3][1] = m31; matrix[3][2] = m32; matrix[3][3] = m33; }
inline const Vector4    Matrix4::getRow                 ( int i ) const                                         { ASSERT(i>=0&&i<4); return Vector4(matrix[i][0],matrix[i][1],matrix[i][2],matrix[i][3]); }
inline const Vector4    Matrix4::getColumn              ( int i ) const                                         { ASSERT(i>=0&&i<4); return Vector4(matrix[0][i],matrix[1][i],matrix[2][i],matrix[3][i]); }
inline void                             Matrix4::setRow                 ( int i, const Vector4& v )                     { ASSERT(i>=0&&i<4); matrix[i][0] = v.x; matrix[i][1] = v.y; matrix[i][2] = v.z; matrix[i][3] = v.w; }
inline void                             Matrix4::setColumn              ( int i, const Vector4& v )                     { ASSERT(i>=0&&i<4); matrix[0][i] = v.x; matrix[1][i] = v.y; matrix[2][i] = v.z; matrix[3][i] = v.w; }
inline const Matrix3    Matrix4::getRotation    () const                                                        { return Matrix3(matrix[0][0],matrix[0][1],matrix[0][2],matrix[1][0],matrix[1][1],matrix[1][2],matrix[2][0],matrix[2][1],matrix[2][2]); }
inline const Vector3    Matrix4::getTranslation ( MultiplicationOrder o ) const         { if( o == LEFT_MULTIPLY ) return Vector3(matrix[0][3],matrix[1][3],matrix[2][3]); else return Vector3(matrix[3][0],matrix[3][1],matrix[3][2]); }
inline void                             Matrix4::setRotation    ( const Matrix3& m )                            { matrix[0][0] = m[0][0]; matrix[0][1] = m[0][1]; matrix[0][2] = m[0][2]; matrix[1][0] = m[1][0]; matrix[1][1] = m[1][1]; matrix[1][2] = m[1][2]; matrix[2][0] = m[2][0]; matrix[2][1] = m[2][1]; matrix[2][2] = m[2][2]; }
inline void                             Matrix4::setTranslation ( const Vector3& v, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) { matrix[0][3] = v.x; matrix[1][3] = v.y; matrix[2][3] = v.z; } else { matrix[3][0] = v.x; matrix[3][1] = v.y; matrix[3][2] = v.z; } }
inline void                             Matrix4::operator*=             ( const Matrix4& m )                            { *this = *this * m; }
inline void                             Matrix4::multiply               ( const Matrix4& m, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) *this = m * (*this); else *this *= m; }
inline void                             Matrix4::operator*=             ( T f )                                                         { for(int i=0;i<4;i++) for(int j=0;j<4;j++) matrix[i][j] *= f; }
inline void                             Matrix4::operator+=             ( const Matrix4& m )                            { for(int i=0;i<4;i++) for(int j=0;j<4;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix4::operator-=             ( const Matrix4& m )                            { for(int i=0;i<4;i++) for(int j=0;j<4;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix4    Matrix4::operator-              () const                                                        { return Matrix4( -matrix[0][0],-matrix[0][1],-matrix[0][2],-matrix[0][3], -matrix[1][0],-matrix[1][1],-matrix[1][2],-matrix[1][3], -matrix[2][0],-matrix[2][1],-matrix[2][2],-matrix[2][3],-matrix[3][0],-matrix[3][1],-matrix[3][2],-matrix[3][3]); }
inline void                             Matrix4::identity               ()                                                                      { for(int i=0;i<4;i++) for(int j=0;j<4;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix4::transpose              ()                                                                      { SWAP(matrix[1][0], matrix[0][1]); SWAP(matrix[2][0], matrix[0][2]); SWAP(matrix[3][0], matrix[0][3]); SWAP(matrix[2][1], matrix[1][2]); SWAP(matrix[3][1], matrix[1][3]); SWAP(matrix[3][2], matrix[2][3]); }
inline void                             Matrix4::translate              ( const Vector3& v, MultiplicationOrder o )     { if( o == LEFT_MULTIPLY ) { matrix[0][3] += v.x; matrix[1][3] += v.y; matrix[2][3] += v.z; } else { matrix[3][0] += v.x; matrix[3][1] += v.y; matrix[3][2] += v.z; } }
inline void                             Matrix4::scale                  ( const Vector3& v, MultiplicationOrder o )     { multiply( Matrix4( v.x, T(0), T(0), T(0), T(0), v.y, T(0), T(0), T(0), T(0), v.z, T(0), T(0), T(0), T(0), T(1) ), o ); }

//==============================================================================================
//==============================================================================================

//Matrix2x3 global functions
inline bool                             operator==      ( const Matrix2x3& m1, const Matrix2x3& m2 )    { for(int i=0;i<2;i++) for(int j=0;j<3;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix2x3& m1, const Matrix2x3& m2 )    { return !(m1 == m2); }
const Matrix2x3                 operator*       ( const Matrix2x3& m1, const Matrix2x3& m2 );   //in .cpp
inline const Matrix2x3  operator*       ( T f, const Matrix2x3& m )                                             { Matrix2x3 t(m); t *= f; return t; }
inline const Matrix2x3  operator*       ( const Matrix2x3& m, T f )                                             { Matrix2x3 t(m); t *= f; return t; }
inline const Matrix2x3  operator+       ( const Matrix2x3& m1, const Matrix2x3& m2 )    { Matrix2x3 t(m1); t += m2; return t; }
inline const Matrix2x3  operator-       ( const Matrix2x3& m1, const Matrix2x3& m2 )    { Matrix2x3 t(m1); t -= m2; return t; }
inline const Matrix2x3  transpose       ( const Matrix2x3& m )                                                  { Matrix2x3 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix2x3  invert          ( const Matrix2x3& m )                                                  { Matrix2x3 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix2x3 inline functions
inline                                  Matrix2x3::Matrix2x3            ()                                                              { identity(); }
inline                                  Matrix2x3::Matrix2x3            ( const Matrix2& m )                    { set(m[0][0],m[0][1],T(0),m[1][0],m[1][1],T(0)); }
inline                                  Matrix2x3::Matrix2x3            ( const Matrix2x3& m )                  { *this = m; }
inline                                  Matrix2x3::Matrix2x3            ( T m00, T m01, T m02, T m10, T m11, T m12 )    { set(m00,m01,m02,m10,m11,m12); }
inline                                  Matrix2x3::~Matrix2x3           ()                                                              {}
inline Matrix2x3&               Matrix2x3::operator=            ( const Matrix2x3& m )                  { for(int i=0;i<2;i++) for(int j=0;j<3;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector3&                 Matrix2x3::operator[]           ( int i )                                               { ASSERT(i>=0&&i<2); return (Vector3&)matrix[i][0]; }
inline const Vector3&   Matrix2x3::operator[]           ( int i ) const                                 { ASSERT(i>=0&&i<2); return (const Vector3&)matrix[i][0]; }
inline void                             Matrix2x3::set                          ( T m00, T m01, T m02, T m10, T m11, T m12 ) { matrix[0][0] = m00; matrix[0][1] = m01; matrix[0][2] = m02; matrix[1][0] = m10; matrix[1][1] = m11; matrix[1][2] = m12;  }
inline const Vector3    Matrix2x3::getRow                       ( int i ) const                                 { ASSERT(i>=0&&i<2); return Vector3(matrix[i][0],matrix[i][1],matrix[i][2]); }
inline const Vector2    Matrix2x3::getColumn            ( int i ) const                                 { ASSERT(i>=0&&i<3); return Vector2(matrix[0][i],matrix[1][i]); }
inline void                             Matrix2x3::setRow                       ( int i, const Vector3& v )             { ASSERT(i>=0&&i<2); matrix[i][0] = v.x; matrix[i][1] = v.y; matrix[i][2] = v.z; }
inline void                             Matrix2x3::setColumn            ( int i, const Vector2& v )             { ASSERT(i>=0&&i<3); matrix[0][i] = v.x; matrix[1][i] = v.y; }
inline const Matrix2    Matrix2x3::getRotation          () const                                                { return Matrix2(matrix[0][0],matrix[0][1],matrix[1][0],matrix[1][1]); }
inline const Vector2    Matrix2x3::getTranslation       () const                                                { return Vector2(matrix[0][2],matrix[1][2]); }
inline void                             Matrix2x3::setRotation          ( const Matrix2& m )                    { matrix[0][0] = m[0][0]; matrix[0][1] = m[0][1]; matrix[1][0] = m[1][0]; matrix[1][1] = m[1][1]; }
inline void                             Matrix2x3::setTranslation       ( const Vector2& v )                    { matrix[0][2] = v[0]; matrix[1][2] = v[1]; }
inline void                             Matrix2x3::multiply                     ( const Matrix2x3& m )                  { *this = m * (*this); }
inline void                             Matrix2x3::operator*=           ( T f )                                                 { for(int i=0;i<2;i++) for(int j=0;j<3;j++) matrix[i][j] *= f; }
inline void                             Matrix2x3::operator+=           ( const Matrix2x3& m )                  { for(int i=0;i<2;i++) for(int j=0;j<3;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix2x3::operator-=           ( const Matrix2x3& m )                  { for(int i=0;i<2;i++) for(int j=0;j<3;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix2x3  Matrix2x3::operator-            () const                                                { return Matrix2x3( -matrix[0][0],-matrix[0][1],-matrix[0][2], -matrix[1][0],-matrix[1][1],-matrix[1][2]); }
inline void                             Matrix2x3::identity                     ()                                                              { for(int i=0;i<2;i++) for(int j=0;j<3;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix2x3::transpose            ()                                                              { SWAP(matrix[1][0], matrix[0][1]); }
inline T                                Matrix2x3::det                          () const                                                { return getRotation().det( ); }
inline void                             Matrix2x3::translate            ( const Vector2& v )                    { matrix[0][2] += v.x; matrix[1][2] += v.y; }
inline void                             Matrix2x3::scale                        ( const Vector2& v )                    { multiply( Matrix2x3( v.x, T(0), T(0), T(0), v.y, T(0) ) ); }

//==============================================================================================
//==============================================================================================

//Matrix3x2 global functions
inline bool                             operator==      ( const Matrix3x2& m1, const Matrix3x2& m2 )    { for(int i=0;i<3;i++) for(int j=0;j<2;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix3x2& m1, const Matrix3x2& m2 )    { return !(m1 == m2); }
const Matrix3x2                 operator*       ( const Matrix3x2& m1, const Matrix3x2& m2 );   //in .cpp
inline const Matrix3x2  operator*       ( T f, const Matrix3x2& m )                                             { Matrix3x2 t(m); t *= f; return t; }
inline const Matrix3x2  operator*       ( const Matrix3x2& m, T f )                                             { Matrix3x2 t(m); t *= f; return t; }
inline const Matrix3x2  operator+       ( const Matrix3x2& m1, const Matrix3x2& m2 )    { Matrix3x2 t(m1); t += m2; return t; }
inline const Matrix3x2  operator-       ( const Matrix3x2& m1, const Matrix3x2& m2 )    { Matrix3x2 t(m1); t -= m2; return t; }
inline const Matrix3x2  transpose       ( const Matrix3x2& m )                                                  { Matrix3x2 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix3x2  invert          ( const Matrix3x2& m )                                                  { Matrix3x2 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix3x2 inline functions
inline                                  Matrix3x2::Matrix3x2            ()                                                              { identity(); }
inline                                  Matrix3x2::Matrix3x2            ( const Matrix2& m )                    { set(m[0][0],m[0][1], m[1][0],m[1][1], T(0),T(0)); }
inline                                  Matrix3x2::Matrix3x2            ( const Matrix3x2& m )                  { *this = m; }
inline                                  Matrix3x2::Matrix3x2            ( T m00, T m01, T m10, T m11, T m20, T m21 )    { set(m00,m01,m10,m11,m20,m21); }
inline                                  Matrix3x2::~Matrix3x2           ()                                                              {}
inline Matrix3x2&               Matrix3x2::operator=            ( const Matrix3x2& m )                  { for(int i=0;i<3;i++) for(int j=0;j<2;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector2&                 Matrix3x2::operator[]           ( int i )                                               { ASSERT(i>=0&&i<3); return (Vector2&)matrix[i][0]; }
inline const Vector2&   Matrix3x2::operator[]           ( int i ) const                                 { ASSERT(i>=0&&i<3); return (const Vector2&)matrix[i][0]; }
inline void                             Matrix3x2::set                          ( T m00, T m01, T m10, T m11, T m20, T m21 ) { matrix[0][0] = m00; matrix[0][1] = m01; matrix[1][0] = m10; matrix[1][1] = m11; matrix[2][0] = m20; matrix[2][1] = m21; }
inline const Vector2    Matrix3x2::getRow                       ( int i ) const                                 { ASSERT(i>=0&&i<3); return Vector2(matrix[i][0],matrix[i][1]); }
inline const Vector3    Matrix3x2::getColumn            ( int i ) const                                 { ASSERT(i>=0&&i<2); return Vector3(matrix[0][i],matrix[1][i],matrix[2][i]); }
inline void                             Matrix3x2::setRow                       ( int i, const Vector2& v )             { ASSERT(i>=0&&i<3); matrix[i][0] = v.x; matrix[i][1] = v.y; }
inline void                             Matrix3x2::setColumn            ( int i, const Vector3& v )             { ASSERT(i>=0&&i<2); matrix[0][i] = v.x; matrix[1][i] = v.y; matrix[2][i] = v.z; }
inline const Matrix2    Matrix3x2::getRotation          () const                                                { return Matrix2(matrix[0][0],matrix[0][1],matrix[1][0],matrix[1][1]); }
inline const Vector2    Matrix3x2::getTranslation       () const                                                { return (*this)[2]; }
inline void                             Matrix3x2::setRotation          ( const Matrix2& m )                    { matrix[0][0] = m[0][0]; matrix[0][1] = m[0][1]; matrix[1][0] = m[1][0]; matrix[1][1] = m[1][1]; }
inline void                             Matrix3x2::setTranslation       ( const Vector2& v )                    { (*this)[2] = v; }
inline void                             Matrix3x2::operator*=           ( const Matrix3x2& m )                  { *this = *this * m; }
inline void                             Matrix3x2::multiply                     ( const Matrix3x2& m )                  { *this *= m; }
inline void                             Matrix3x2::operator*=           ( T f )                                                 { for(int i=0;i<3;i++) for(int j=0;j<2;j++) matrix[i][j] *= f; }
inline void                             Matrix3x2::operator+=           ( const Matrix3x2& m )                  { for(int i=0;i<3;i++) for(int j=0;j<2;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix3x2::operator-=           ( const Matrix3x2& m )                  { for(int i=0;i<3;i++) for(int j=0;j<2;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix3x2  Matrix3x2::operator-            () const                                                { return Matrix3x2( -matrix[0][0],-matrix[0][1], -matrix[1][0],-matrix[1][1], -matrix[2][0],-matrix[2][1]); }
inline void                             Matrix3x2::identity                     ()                                                              { for(int i=0;i<3;i++) for(int j=0;j<2;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix3x2::transpose            ()                                                              { SWAP(matrix[1][0], matrix[0][1]); }
inline T                                Matrix3x2::det                          () const                                                { return ((Matrix2&)(*this)).det( ); }  //works because upperleft corner of Matrix3x2 is Matrix2x2
inline void                             Matrix3x2::translate            ( const Vector2& v )                    { matrix[2][0] += v.x; matrix[2][1] += v.y; }
inline void                             Matrix3x2::scale                        ( const Vector2& v )                    { multiply( Matrix3x2( v.x, T(0), T(0), v.y, T(0), T(0) ) ); }

//==============================================================================================
//==============================================================================================

//Matrix3x4 global functions
inline bool                             operator==      ( const Matrix3x4& m1, const Matrix3x4& m2 )    { for(int i=0;i<3;i++) for(int j=0;j<4;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix3x4& m1, const Matrix3x4& m2 )    { return !(m1 == m2); }
const Matrix3x4                 operator*       ( const Matrix3x4& m1, const Matrix3x4& m2 );   //in .cpp
inline const Matrix3x4  operator*       ( T f, const Matrix3x4& m )                                             { Matrix3x4 t(m); t *= f; return t; }
inline const Matrix3x4  operator*       ( const Matrix3x4& m, T f )                                             { Matrix3x4 t(m); t *= f; return t; }
inline const Matrix3x4  operator+       ( const Matrix3x4& m1, const Matrix3x4& m2 )    { Matrix3x4 t(m1); t += m2; return t; }
inline const Matrix3x4  operator-       ( const Matrix3x4& m1, const Matrix3x4& m2 )    { Matrix3x4 t(m1); t -= m2; return t; }
inline const Matrix3x4  transpose       ( const Matrix3x4& m )                                                  { Matrix3x4 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix3x4  invert          ( const Matrix3x4& m )                                                  { Matrix3x4 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix3x4 inline functions
inline                                  Matrix3x4::Matrix3x4            ()                                                              { identity(); }
inline                                  Matrix3x4::Matrix3x4            ( const Matrix2& m )                    { set(m[0][0],m[0][1],T(0),T(0),m[1][0],m[1][1],T(0),T(0),T(0),T(0),T(1),T(0)); }
inline                                  Matrix3x4::Matrix3x4            ( const Matrix2x3& m )                  { set(m[0][0],m[0][1],m[0][2],T(0),m[1][0],m[1][1],m[1][2],T(0),T(0),T(0),T(1),T(0)); }
inline                                  Matrix3x4::Matrix3x4            ( const Matrix3x2& m )                  { set(m[0][0],m[0][1],T(0),T(0),m[1][0],m[1][1],T(0),T(0),m[2][0],m[2][1],T(1),T(0)); }
inline                                  Matrix3x4::Matrix3x4            ( const Matrix3& m )                    { set(m[0][0],m[0][1],m[0][2],T(0),m[1][0],m[1][1],m[1][2],T(0),m[2][0],m[2][1],m[2][2],T(0)); }
inline                                  Matrix3x4::Matrix3x4            ( const Matrix3x4& m )                  { *this = m; }
inline                                  Matrix3x4::Matrix3x4            ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23 )  { set(m00,m01,m02,m03,m10,m11,m12,m13,m20,m21,m22,m23); }
inline                                  Matrix3x4::~Matrix3x4           ()                                                              {}
inline Matrix3x4&               Matrix3x4::operator=            ( const Matrix3x4& m )                  { for(int i=0;i<3;i++) for(int j=0;j<4;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector4&                 Matrix3x4::operator[]           ( int i )                                               { ASSERT(i>=0&&i<3); return (Vector4&)matrix[i][0]; }
inline const Vector4&   Matrix3x4::operator[]           ( int i ) const                                 { ASSERT(i>=0&&i<3); return (const Vector4&)matrix[i][0]; }
inline void                             Matrix3x4::set                          ( T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23 ) { matrix[0][0] = m00; matrix[0][1] = m01; matrix[0][2] = m02; matrix[0][3] = m03; matrix[1][0] = m10; matrix[1][1] = m11; matrix[1][2] = m12;  matrix[1][3] = m13; matrix[2][0] = m20; matrix[2][1] = m21; matrix[2][2] = m22; matrix[2][3] = m23; }
inline const Vector4    Matrix3x4::getRow                       ( int i ) const                                 { ASSERT(i>=0&&i<3); return Vector4(matrix[i][0],matrix[i][1],matrix[i][2],matrix[i][3]); }
inline const Vector3    Matrix3x4::getColumn            ( int i ) const                                 { ASSERT(i>=0&&i<4); return Vector3(matrix[0][i],matrix[1][i],matrix[2][i]); }
inline void                             Matrix3x4::setRow                       ( int i, const Vector4& v )             { ASSERT(i>=0&&i<3); matrix[i][0] = v.x; matrix[i][1] = v.y; matrix[i][2] = v.z; matrix[i][3] = v.w; }
inline void                             Matrix3x4::setColumn            ( int i, const Vector3& v )             { ASSERT(i>=0&&i<4); matrix[0][i] = v.x; matrix[1][i] = v.y; matrix[2][i] = v.z; }
inline const Matrix3    Matrix3x4::getRotation          () const                                                { return Matrix3(matrix[0][0],matrix[0][1],matrix[0][2],matrix[1][0],matrix[1][1],matrix[1][2],matrix[2][0],matrix[2][1],matrix[2][2]); }
inline const Vector3    Matrix3x4::getTranslation       () const                                                { return Vector3(matrix[0][3],matrix[1][3],matrix[2][3]); }
inline void                             Matrix3x4::setRotation          ( const Matrix3& m )                    { matrix[0][0] = m[0][0]; matrix[0][1] = m[0][1]; matrix[0][2] = m[0][2]; matrix[1][0] = m[1][0]; matrix[1][1] = m[1][1]; matrix[1][2] = m[1][2]; matrix[2][0] = m[2][0]; matrix[2][1] = m[2][1]; matrix[2][2] = m[2][2]; }
inline void                             Matrix3x4::setTranslation       ( const Vector3& v )                    { matrix[0][3] = v[0]; matrix[1][3] = v[1]; matrix[2][3] = v[2]; }
inline void                             Matrix3x4::multiply                     ( const Matrix3x4& m )                  { *this = m * (*this); }
inline void                             Matrix3x4::operator*=           ( T f )                                                 { for(int i=0;i<3;i++) for(int j=0;j<4;j++) matrix[i][j] *= f; }
inline void                             Matrix3x4::operator+=           ( const Matrix3x4& m )                  { for(int i=0;i<3;i++) for(int j=0;j<4;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix3x4::operator-=           ( const Matrix3x4& m )                  { for(int i=0;i<3;i++) for(int j=0;j<4;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix3x4  Matrix3x4::operator-            () const                                                { return Matrix3x4( -matrix[0][0],-matrix[0][1],-matrix[0][2],-matrix[0][3], -matrix[1][0],-matrix[1][1],-matrix[1][2],-matrix[1][3], -matrix[2][0],-matrix[2][1],-matrix[2][2],-matrix[2][3]); }
inline void                             Matrix3x4::identity                     ()                                                              { for(int i=0;i<3;i++) for(int j=0;j<4;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix3x4::transpose            ()                                                              { SWAP(matrix[1][0], matrix[0][1]); SWAP(matrix[2][0], matrix[0][2]); SWAP(matrix[2][1], matrix[1][2]); }
inline T                                Matrix3x4::det                          () const                                                { return getRotation().det( ); }
inline void                             Matrix3x4::translate            ( const Vector3& v )                    { matrix[0][3] += v.x; matrix[1][3] += v.y; matrix[2][3] += v.z; }
inline void                             Matrix3x4::scale                        ( const Vector3& v )                    { multiply( Matrix3x4( v.x, T(0), T(0), T(0), T(0), v.y, T(0), T(0), T(0), T(0), v.z, T(0)) ); }

//==============================================================================================
//==============================================================================================

//Matrix4x3 global functions
inline bool                             operator==      ( const Matrix4x3& m1, const Matrix4x3& m2 )    { for(int i=0;i<4;i++) for(int j=0;j<3;j++) if( m1[i][j] != m2[i][j] ) return false; return true; }
inline bool                             operator!=      ( const Matrix4x3& m1, const Matrix4x3& m2 )    { return !(m1 == m2); }
const Matrix4x3                 operator*       ( const Matrix4x3& m1, const Matrix4x3& m2 );   //in .cpp
inline const Matrix4x3  operator*       ( T f, const Matrix4x3& m )                                             { Matrix4x3 t(m); t *= f; return t; }
inline const Matrix4x3  operator*       ( const Matrix4x3& m, T f )                                             { Matrix4x3 t(m); t *= f; return t; }
inline const Matrix4x3  operator+       ( const Matrix4x3& m1, const Matrix4x3& m2 )    { Matrix4x3 t(m1); t += m2; return t; }
inline const Matrix4x3  operator-       ( const Matrix4x3& m1, const Matrix4x3& m2 )    { Matrix4x3 t(m1); t -= m2; return t; }
inline const Matrix4x3  transpose       ( const Matrix4x3& m )                                                  { Matrix4x3 t(m); t.transpose(); return t; }
// if the matrix is singular, returns it unmodified
inline const Matrix4x3  invert          ( const Matrix4x3& m )                                                  { Matrix4x3 t(m); t.invert(); return t; }

//==============================================================================================
//==============================================================================================

//Matrix4x3 inline functions
inline                                  Matrix4x3::Matrix4x3            ()                                                              { identity(); }
inline                                  Matrix4x3::Matrix4x3            ( const Matrix2& m )                    { set(m[0][0],m[0][1],T(0),m[1][0],m[1][1],T(0),T(0),T(0),T(1),T(0),T(0),T(0)); }
inline                                  Matrix4x3::Matrix4x3            ( const Matrix2x3& m )                  { set(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],T(0),T(0),T(1),T(0),T(0),T(0)); }
inline                                  Matrix4x3::Matrix4x3            ( const Matrix3x2& m )                  { set(m[0][0],m[0][1],T(0),m[1][0],m[1][1],T(0),m[2][0],m[2][1],T(1),T(0),T(0),T(0)); }
inline                                  Matrix4x3::Matrix4x3            ( const Matrix3& m )                    { set(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2],T(0),T(0),T(0)); }
inline                                  Matrix4x3::Matrix4x3            ( const Matrix4x3& m )                  { *this = m; }
inline                                  Matrix4x3::Matrix4x3            ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22, T m30, T m31, T m32 )  { set(m00,m01,m02,m10,m11,m12,m20,m21,m22,m30,m31,m32); }
inline                                  Matrix4x3::~Matrix4x3           ()                                                              {}
inline Matrix4x3&               Matrix4x3::operator=            ( const Matrix4x3& m )                  { for(int i=0;i<4;i++) for(int j=0;j<3;j++) matrix[i][j] = m.matrix[i][j]; return *this; }
inline Vector3&                 Matrix4x3::operator[]           ( int i )                                               { ASSERT(i>=0&&i<4); return (Vector3&)matrix[i][0]; }
inline const Vector3&   Matrix4x3::operator[]           ( int i ) const                                 { ASSERT(i>=0&&i<4); return (const Vector3&)matrix[i][0]; }
inline void                             Matrix4x3::set                          ( T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22, T m30, T m31, T m32 ) { matrix[0][0] = m00; matrix[0][1] = m01; matrix[0][2] = m02; matrix[1][0] = m10; matrix[1][1] = m11; matrix[1][2] = m12; matrix[2][0] = m20; matrix[2][1] = m21; matrix[2][2] = m22; matrix[3][0] = m30; matrix[3][1] = m31; matrix[3][2] = m32; }
inline const Vector3    Matrix4x3::getRow                       ( int i ) const                                 { ASSERT(i>=0&&i<4); return Vector3(matrix[i][0],matrix[i][1],matrix[i][2]); }
inline const Vector4    Matrix4x3::getColumn            ( int i ) const                                 { ASSERT(i>=0&&i<3); return Vector4(matrix[0][i],matrix[1][i],matrix[2][i],matrix[3][i]); }
inline void                             Matrix4x3::setRow                       ( int i, const Vector3& v )             { ASSERT(i>=0&&i<4); matrix[i][0] = v.x; matrix[i][1] = v.y; matrix[i][2] = v.z; }
inline void                             Matrix4x3::setColumn            ( int i, const Vector4& v )             { ASSERT(i>=0&&i<3); matrix[0][i] = v.x; matrix[1][i] = v.y; matrix[2][i] = v.z; matrix[3][i] = v.w; }
inline const Matrix3    Matrix4x3::getRotation          () const                                                { return Matrix3(matrix[0][0],matrix[0][1],matrix[0][2],matrix[1][0],matrix[1][1],matrix[1][2],matrix[2][0],matrix[2][1],matrix[2][2]); }
inline const Vector3    Matrix4x3::getTranslation       () const                                                { return (*this)[3]; }
inline void                             Matrix4x3::setRotation          ( const Matrix3& m )                    { matrix[0][0] = m[0][0]; matrix[0][1] = m[0][1]; matrix[0][2] = m[0][2]; matrix[1][0] = m[1][0]; matrix[1][1] = m[1][1]; matrix[1][2] = m[1][2]; matrix[2][0] = m[2][0]; matrix[2][1] = m[2][1]; matrix[2][2] = m[2][2]; }
inline void                             Matrix4x3::setTranslation       ( const Vector3& v )                    { (*this)[3] = v; }
inline void                             Matrix4x3::operator*=           ( const Matrix4x3& m )                  { *this = *this * m; }
inline void                             Matrix4x3::multiply                     ( const Matrix4x3& m )                  { *this *= m; }
inline void                             Matrix4x3::operator*=           ( T f )                                                 { for(int i=0;i<4;i++) for(int j=0;j<3;j++) matrix[i][j] *= f; }
inline void                             Matrix4x3::operator+=           ( const Matrix4x3& m )                  { for(int i=0;i<4;i++) for(int j=0;j<3;j++) matrix[i][j] += m.matrix[i][j]; }
inline void                             Matrix4x3::operator-=           ( const Matrix4x3& m )                  { for(int i=0;i<4;i++) for(int j=0;j<3;j++) matrix[i][j] -= m.matrix[i][j]; }
inline const Matrix4x3  Matrix4x3::operator-            () const                                                { return Matrix4x3( -matrix[0][0],-matrix[0][1],-matrix[0][2], -matrix[1][0],-matrix[1][1],-matrix[1][2], -matrix[2][0],-matrix[2][1],-matrix[2][2],-matrix[3][0],-matrix[3][1],-matrix[3][2]); }
inline void                             Matrix4x3::identity                     ()                                                              { for(int i=0;i<4;i++) for(int j=0;j<3;j++) matrix[i][j] = (i == j) ? T(1) : T(0); }
inline void                             Matrix4x3::transpose            ()                                                              { SWAP(matrix[1][0], matrix[0][1]); SWAP(matrix[2][0], matrix[0][2]); SWAP(matrix[2][1], matrix[1][2]); }
inline T                                Matrix4x3::det                          () const                                                { return ((Matrix3&)(*this)).det( ); }  //works because upperleft corner of Matrix4x3 is Matrix3x3
inline void                             Matrix4x3::translate            ( const Vector3& v )                    { matrix[3][0] += v.x; matrix[3][1] += v.y; matrix[3][2] += v.z; }
inline void                             Matrix4x3::scale                        ( const Vector3& v )                    { multiply( Matrix4x3( v.x, T(0), T(0), T(0), v.y, T(0), T(0), T(0), v.z, T(0), T(0), T(0)) ); }

//==============================================================================================
//==============================================================================================

//global matrix conversions
//NOTE: these could be constructors of respective classes, but as the compiler would use constructors
// in unexpected places, it was considered safer to make the conversion explicit.
// the safe conversions (from smaller to larger) are written as constructors

//there are no conversion for vectors, as these are easy enough to write with constructors:
// Vector3 a(1,2,3); Vector2 b(a.x, a.y);

//these transpose the matrix
inline const Matrix2x3  makeMatrix2x3   ( const Matrix3x2& m )                                          { Matrix2x3 t; for(int i=0;i<3;i++) for(int j=0;j<2;j++) t[j][i] = m[i][j]; return t; }
inline const Matrix3x2  makeMatrix3x2   ( const Matrix2x3& m )                                          { Matrix3x2 t; for(int i=0;i<3;i++) for(int j=0;j<2;j++) t[i][j] = m[j][i]; return t; }

inline const Matrix3x4  makeMatrix3x4   ( const Matrix4x3& m )                                          { Matrix3x4 t; for(int i=0;i<4;i++) for(int j=0;j<3;j++) t[j][i] = m[i][j]; return t; }
inline const Matrix4x3  makeMatrix4x3   ( const Matrix3x4& m )                                          { Matrix4x3 t; for(int i=0;i<4;i++) for(int j=0;j<3;j++) t[i][j] = m[j][i]; return t; }

//these take the upper-left submatrix of the larger matrix
inline const Matrix2    makeMatrix2             ( const Matrix2x3& m )                                          { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }
inline const Matrix2    makeMatrix2             ( const Matrix3x2& m )                                          { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }
inline const Matrix2    makeMatrix2             ( const Matrix3& m )                                            { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }
inline const Matrix2    makeMatrix2             ( const Matrix3x4& m )                                          { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }
inline const Matrix2    makeMatrix2             ( const Matrix4x3& m )                                          { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }
inline const Matrix2    makeMatrix2             ( const Matrix4& m )                                            { return Matrix2(m[0][0],m[0][1],m[1][0],m[1][1]); }

inline const Matrix3x2  makeMatrix3x2   ( const Matrix3& m )                                            { return Matrix3x2(m[0][0],m[0][1],m[1][0],m[1][1],m[2][0],m[2][1]); }
inline const Matrix2x3  makeMatrix2x3   ( const Matrix3& m )                                            { return Matrix2x3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2]); }

inline const Matrix3    makeMatrix3             ( const Matrix3x4& m )                                          { return Matrix3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2]); }
inline const Matrix3    makeMatrix3             ( const Matrix4x3& m )                                          { return Matrix3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2]); }
inline const Matrix3    makeMatrix3             ( const Matrix4& m )                                            { return Matrix3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2]); }

inline const Matrix3x2  makeMatrix3x2   ( const Matrix4& m )                                            { return Matrix3x2(m[0][0],m[0][1],m[1][0],m[1][1],m[2][0],m[2][1]); }
inline const Matrix2x3  makeMatrix2x3   ( const Matrix4& m )                                            { return Matrix2x3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2]); }

inline const Matrix4x3  makeMatrix4x3   ( const Matrix4& m )                                            { return Matrix4x3(m[0][0],m[0][1],m[0][2],m[1][0],m[1][1],m[1][2],m[2][0],m[2][1],m[2][2],m[3][0],m[3][1],m[3][2]); }
inline const Matrix3x4  makeMatrix3x4   ( const Matrix4& m )                                            { return Matrix3x4(m[0][0],m[0][1],m[0][2],m[0][3],m[1][0],m[1][1],m[1][2],m[1][3],m[2][0],m[2][1],m[2][2],m[2][3]); }

//==============================================================================================
//==============================================================================================





//==============================================================================================
//==============================================================================================

class Vector2Factory
{
public:
        static const Vector2    randInSquare    ( T edgeLength = T(2) );
        static const Vector2    randInCircle    ( T radius = T(1) );
        static const Vector2    randOnCircle    ( T radius = T(1) );
        //returns a vector perpendicular to v in the counterclockwise direction
        static const Vector2    perpendicular   ( const Vector2& v );
};

//==============================================================================================
//==============================================================================================

class Vector3Factory
{
public:
        static const Vector3    randInCube              ( T edgeLength = T(2) );
        static const Vector3    randInSphere    ( T radius = T(1) );
        static const Vector3    randOnSphere    ( T radius = T(1) );
        //returns a "random" vector perpendicular to v
        //for more randomness, rotate the result by a random amount about v
        static const Vector3    perpendicular   ( const Vector3& v );
};

//==============================================================================================
//==============================================================================================

class Vector4Factory
{
public:
        static const Vector4    randInCube              ( T edgeLength = T(2) );
        static const Vector4    randInSphere    ( T radius = T(1) );
        static const Vector4    randOnSphere    ( T radius = T(1) );
};

//==============================================================================================
//==============================================================================================

//constructs left-handed coordinate systems
//TODO add support for right-handed coordinate systems
//TODO add random rotation matrix?

// \todo [020731 jan] I'd prefer complete set of functionality for both D3D and OpenGL.
//                                              - enum ProjectionModel { OPENGL, DIRECT3D, WHATNOT };
//                                              - need frustum, ortho, frustumFOV

class MatrixFactory
{
public:
        //given the position of the viewer and the target, and the up-vector, makes an orthonormal rotation
        //matrix where z-axis points towards the target
        static const Matrix3    target          ( const Vector3& pos, const Vector3& tgt, const Vector3& up, MultiplicationOrder o );

        //The (left, bottom, znear) and (right, top, znear) parameters specify the points on the
        //near clipping plane that are mapped to the lower-left and upper-right corners of the window,
        //respectively, assuming that the eye is located at (0, 0, 0)
        //this maps the view frustum to a box with dimensions [-1,1]
        //this corresponds to OpenGL projection matrix format
        static const Matrix4    frustum         ( T left, T right, T bottom, T top, T zNear, T zFar, MultiplicationOrder o );

        //this maps the view frustum to a box with dimensions [-1,1]
        //this corresponds to OpenGL projection matrix format
        //FOV is given in radians
        static const Matrix4    frustum         ( T horizontalFOV, T verticalFOV, T zNear, T zFar, MultiplicationOrder o );

        //The (left, bottom, near) and (right, top, near) parameters specify the points on the
        //near clipping plane that are mapped to the lower-left and upper-right corners of the window,
        //respectively, assuming that the eye is located at (0, 0, 0)
        //this maps the view frustum to a box with dimensions [-1,1]
        //this corresponds to OpenGL projection matrix format
        static const Matrix4    ortho           ( T left, T right, T bottom, T top, T zNear, T zFar, MultiplicationOrder o );

        //this maps the view frustum to a box with dimensions [-1,1] for x- and y-axis, and [0,1] for z-axis
        //this corresponds to Direct3D projection matrix format
        //FOV is given in radians
        static const Matrix4    frustum01       ( T horizontalFOV, T verticalFOV, T zNear, T zFar, MultiplicationOrder o );

};

//==============================================================================================
//==============================================================================================

}      //namespace VECTORMATRIX

#endif //__VECTORMATRIX_HPP


// NEW FILE

#ifndef SBDEF_H
#define SBDEF_H

//=======================================================================
//=======================================================================

//#include <float.h>            //for FLT_MAX
#include <math.h>               //for sqrt, sin and cos
#include <stdlib.h>             //for rand
#include <assert.h>
#include <string.h>

//TODO add Timer class?

//=======================================================================
//=======================================================================

#define PI                                              3.141592654f
template<class TT> inline TT    RAD( TT x )                                                             { return x * TT(PI) / TT(180); }
template<class TT> inline TT    DEG( TT x )                                                             { return x * TT(180) / TT(PI); }
template<class TT> inline TT    MAX( TT a, TT b )                                               { return ( a > b ) ? a : b; }
template<class TT> inline TT    MIN( TT a, TT b )                                               { return ( a < b ) ? a : b; }
template<class TT> inline TT    CLAMP( TT x, TT min, TT max )                   { ASSERT( min <= max ); return (x >= max) ? max : (x <= min) ? min : x; }
template<class TT> inline TT    STEP( TT x, TT edge )                                   { if( x < edge ) return TT(0); return TT(1); }
template<class TT> inline TT    SMOOTHSTEP( TT x, TT edge0, TT edge1 )  { ASSERT( edge0 <= edge1 ); if(x <= edge0) return TT(0); if(x >= edge1) return TT(1); TT t = (x - edge0) / (edge1 - edge0); return t * t * (TT(3) - TT(2) * t); }
template<class TT> inline int   SIGN( TT x )                                                    { unsigned int a = *reinterpret_cast<const unsigned int*>(&x); if( !(a & 0x7fffffff) ) return 0; if( a & 0x80000000 ) return -1; return 1; }
template<class TT> inline TT    LERP( TT x, TT y, TT r )                                { return x + r * (y - x); }
inline float                                    RAND01()                                                                { return float(rand()) / float(RAND_MAX); }     //random number [0,1]
inline int                                              FLOOR( const float& f )
{ 
        int a                   = *reinterpret_cast<const int*>(&f);                                                                    // take bit pattern of float into a register
        int sign                = (a>>31);                                                                                                                              // sign = 0xFFFFFFFF if original value is negative, 0 if positive
        a&=0x7fffffff;                                                                                                                                                  // we don't need the sign any more
        int exponent    = (a>>23)-127;                                                                                                                  // extract the exponent
        int expsign             = ~(exponent>>31);                                                                                                              // 0xFFFFFFFF if exponent is positive, 0 otherwise
        int imask               = ( (1<<(31-(exponent))))-1;                                                                                    // mask for true integer values
        int mantissa    = (a&((1<<23)-1));                                                                                                              // extract mantissa (without the hidden bit)
        int r                   = ((unsigned int)(mantissa|(1<<23))<<8)>>(31-exponent);                                 // ((1<<exponent)*(mantissa|hidden bit))>>24 -- (we know that mantissa > (1<<24))
        r = ((r & expsign) ^ (sign)) + ((!((mantissa<<8)&imask)&(expsign^((a-1)>>31)))&sign);   // if (fabs(value)<1.0) value = 0; copy sign; if (value < 0 && value==(int)(value)) value++; 
        return r;
}
inline int                                              CEIL( const float& f )
{ 
        int a                   = *reinterpret_cast<const int*>(&f) ^ 0x80000000;                                               // take bit pattern of float into a register
        int sign                = (a>>31);                                                                                                                              // sign = 0xFFFFFFFF if original value is negative, 0 if positive
        a&=0x7fffffff;                                                                                                                                                  // we don't need the sign any more
        int exponent    = (a>>23)-127;                                                                                                                  // extract the exponent
        int expsign             = ~(exponent>>31);                                                                                                              // 0xFFFFFFFF if exponent is positive, 0 otherwise
        int imask               = ( (1<<(31-(exponent))))-1;                                                                                    // mask for true integer values
        int mantissa    = (a&((1<<23)-1));                                                                                                              // extract mantissa (without the hidden bit)
        int r                   = ((unsigned int)(mantissa|(1<<23))<<8)>>(31-exponent);                                 // ((1<<exponent)*(mantissa|hidden bit))>>24 -- (we know that mantissa > (1<<24))
        r = ((r & expsign) ^ (sign)) + ((!((mantissa<<8)&imask)&(expsign^((a-1)>>31)))&sign);   // if (fabs(value)<1.0) value = 0; copy sign; if (value < 0 && value==(int)(value)) value++; 
        return -r;
}
inline float                                    FRACT( float x )                                                { return x - float(FLOOR(x)); }
inline unsigned int                             CEILPOW2( unsigned int x )
{
        x--;
        x |= (x >>  1);
        x |= (x >>  2);
        x |= (x >>  4);
        x |= (x >>  8);
        x |= (x >> 16);
        return x + 1;
}

//=======================================================================
//=======================================================================

class Surface
{
public:
        Surface( int w, int h )
        {
                m_data = NULL;
                m_width = 0;
                m_height = 0;
                resize(w,h);
        }
        ~Surface()
        {
                delete[] m_data;
        }

        void            resize( int w, int h )
        {
                assert(w>0 && h>0);
                if( w == m_width && h == m_height )
                        return;
                delete[] m_data;
                m_width = w;
                m_height = h;
                m_data = new unsigned int[w*h];
        }

        void            clear()
        {
                memset(m_data,0,sizeof(unsigned int)*m_width*m_height);
        }

        int                             getWidth() const                { return m_width; }
        int                             getHeight() const               { return m_height; }
        unsigned int*   getSurfaceData() const  { return m_data; }
private:
        unsigned int*   m_data;
        int                             m_width;
        int                             m_height;
};

//=======================================================================
//=======================================================================

//#ifdef DETECT_MEMORY_LEAKS ???
// memory leak detection
/*void* operator new            (size_t s);
void*   operator new[]          (size_t s);
void    operator delete         (void* p);
void    operator delete[]       (void* p);
*/
//#endif
/*      //put these in main.cpp
static int s_allocations = 0, s_allocationsDone = 0;
void*   operator new            (size_t s)      { ++s_allocations; ++s_allocationsDone; return malloc(s);       }
void*   operator new[]          (size_t s)      { ++s_allocations; ++s_allocationsDone; return malloc(s);       }
void    operator delete         (void* p)       { if (p) { free(p); --s_allocations;}   }
void    operator delete[]       (void* p)       { if (p) { free(p); --s_allocations;}   }
*/

//=======================================================================
//=======================================================================

// #include "VectorMatrix.hpp"
using namespace VECTORMATRIX;

//=======================================================================
//=======================================================================

//=======================================================================
//=======================================================================

#endif //SBDEF_H


#ifndef __FILE_H
#define __FILE_H

bool                    fileExists( const char *name );
char*                   fileLoad( const char *name );
void                    fileSave( const char *name, const void *s, unsigned int l );
unsigned int    fileSize( const char *name );
void                    fileFree( char *data );

#endif //__FILE_H
#ifndef MESH_H
#define MESH_H

// #include "def.h"

#ifdef _MSC_VER         //ms visual studio
#pragma warning(disable:4786)   //debug info was truncated to 255 characters
#endif

#include <string>
#include <vector>

//TODO test discontinuos vertex maps (VMAD) with more than one point chunk

//=============================================================================
//=============================================================================

//fast next and prev functions with mod 3
//the input must be 0,1 or 2
inline int nextMod3( int i )    { assert(i>=0&&i<3); return (1 << i) & 3; }     // (i+1) % 3
inline int prevMod3( int i )    { assert(i>=0&&i<3); return (0x102 >> (i<<2)) & 0xf; }  // (i-1+3) % 3

//=============================================================================
//=============================================================================

class Mesh
{
public:
#pragma pack( push, 4 )
        struct Vertex
        {
                void                    null()
                {
                        m_position.set(0,0,0);
                        m_normal.set(0,0,0);
                        m_t1.set(0,0,0);
                        m_t2.set(0,0,0);
                        m_uv.set(0,0);
                        m_diffuse.set(0,0,0,0);
                }
                void                    add( const Vertex& v )
                {
                        m_position += v.m_position;
                        m_normal += v.m_normal;
                        m_t1 += v.m_t1;
                        m_t2 += v.m_t2;
                        m_uv += v.m_uv;
                        m_diffuse += v.m_diffuse;
                }
                void                    mad( const Vertex& v, float c )
                {
                        m_position += v.m_position * c;
                        m_normal += v.m_normal * c;
                        m_t1 += v.m_t1 * c;
                        m_t2 += v.m_t2 * c;
                        m_uv += v.m_uv * c;
                        m_diffuse += v.m_diffuse * c;
                }
                void                    normalize()
                {
                        m_normal.normalize();
                        m_t1.normalize();
                        m_t2.normalize();
                }
                Vector3                 m_position;
                Vector3                 m_normal;
                Vector3                 m_t1;                           //tangent 1
                Vector3                 m_t2;                           //tangent 2
                int                             m_triangleUsing;        //30msb polyindex 2lsb edgeindex
                Vector2                 m_uv;
                Vector4                 m_diffuse;
                bool                    m_reliable;                     //true if it was possible to calculate the tangent field reliably using curvature
        };

        struct Triangle
        {
                Vector3                 m_normal;
                float                   m_area;
                int                             m_surface;                      //index into the surface array
                int                             m_v[3];                         //vertex indices into the vertex array
                int                             m_n[3];                         //neighbour indices into the triangle array
                char                    m_nedge[3];                     //corresponding neighbours' edge numbers
                char                    m_dummy;                        //for alignment
                Vector2                 m_uv[3];
                int                             m_baseTriangle;         //original triangle index before subdivision
        };

        struct Edge
        {
                bool operator<( const Edge& e ) const   { if( m_v[0] < e.m_v[0] ) return true; else if( m_v[0] > e.m_v[0] ) return false; else return (m_v[1] < e.m_v[1]) ? true : false; }
                bool operator==( const Edge& e ) const  { if( m_v[0] == e.m_v[0] && m_v[1] == e.m_v[1] ) return true; return false; }
                int                             m_v[2];
                int                             m_t[2];
        };

        struct Surface
        {
                Vector4                 m_color;
                int                             m_polygons;                     //polygons using this surface material
                std::string             m_name;
        };
#pragma pack( pop )

        Mesh();
        Mesh( int vertices, int triangles );
        Mesh( const char* fileName );
        Mesh( const Mesh& m );
        ~Mesh();

        void                                    load( const char* fileName );
        void                                    save( const char* fileName );
        void                                    copy( const Mesh& m );
        void                                    clear( int newVertices, int newTriangles );

        Vertex&                                 getVertex( int i )                              { assert(i>=0&&i<(int)m_vertices.size()); return m_vertices[i]; }
        Triangle&                               getTriangle( int i )                    { assert(i>=0&&i<(int)m_triangles.size()); return m_triangles[i]; }
        Edge&                                   getEdge( int i )                                { assert(i>=0&&i<(int)m_edges.size()); return m_edges[i]; }
        Surface&                                getSurface( int i )                             { assert(i>=0&&i<(int)m_surfaces.size()); return m_surfaces[i]; }

        int                                             getVertexCount() const                  { return m_vertices.size(); }
        int                                             getTriangleCount() const                { return m_triangles.size(); }
        int                                             getEdgeCount() const                    { return m_edges.size(); }
        int                                             getSurfaceCount() const                 { return m_surfaces.size(); }

        std::vector<Vertex>&    getVertices()                                   { return m_vertices; }
        std::vector<Triangle>&  getTriangles()                                  { return m_triangles; }
        std::vector<Edge>&              getEdges()                                              { return m_edges; }
        std::vector<Surface>&   getSurfaces()                                   { return m_surfaces; }

        const Vector3&                  getBoundingBoxMin() const               { return m_boundingBoxMin; }
        const Vector3&                  getBoundingBoxMax() const               { return m_boundingBoxMax; }
        float                                   getBoundingRadius() const               { return m_boundingRadius; }

        void                                    calculateAdjacency();   //calculates also edges and sets triangleUsing
        void                                    recalculateAdjacency(); //forces recalculation of adjacency
        void                                    calculateCurvature();
        void                                    calculateNormals();
        void                                    setParametrization();
        void                                    interpolateTangents( float stddev );
        void                                    singleSurfaceVertices();        //duplicates vertices which have more than 1 surface
        void                                    silhouette( std::vector<int>& lines, const Matrix4x3& objectToWorldMatrix, const Vector3& cameraInWorld );
        void                                    calculateBoundingBox();
        void                                    duplicateVertices();    //makes each triangle to have its own vertices

        //vertices, triangles and edges at the current subdivision level will be passed in the
        // displacement function. if original vertices etc. are needed, they can be requested
        // as normally through mesh->getVertices() etc.
        typedef void DisplacementFunc( Mesh* mesh, std::vector<Vertex>& vertices, const std::vector<Triangle>& triangles, const std::vector<Edge>& edges, int level );
        void                                    subdivide( int n, DisplacementFunc* displacement = NULL );

        struct ventry
        {
                int             m_vertex;
                int             m_i;
        };
        const std::vector<ventry>&      get1Neighbourhood( int ver, int& valence, bool& boundary, int& left, int& right ) const;
        const std::vector<int>&         get1NeighbourhoodPolygons( int ver ) const;
protected:
        std::vector<Vertex>             m_vertices;
        std::vector<Triangle>   m_triangles;
        std::vector<Edge>               m_edges;
        std::vector<Surface>    m_surfaces;
        std::string                             m_filename;

        Vector3                                 m_boundingBoxMin;
        Vector3                                 m_boundingBoxMax;
        float                                   m_boundingRadius;
        int                                             m_ambiguousNeighbours;
        bool                                    m_adjacencyCalculated;

        static unsigned char*   m_fileptr;
        static unsigned char*   m_nextchunkptr;
        static unsigned char*   m_nextsubchunkptr;
        static unsigned char*   m_fileend;
        static unsigned char*   m_readptr;
        static unsigned char*   m_writeptr;

        unsigned int                    readID4();
        int                                             readI1();
        int                                             readI2();
        int                                             readI4();
        unsigned char                   readU1();
        unsigned short                  readU2();
        unsigned int                    readU4();
        float                                   readF4();
        const char*                             readS0();
        unsigned int                    readVX();
        Vector3                                 readCOL12();
        Vector3                                 readVEC12();

        void                                    writeID4( unsigned int i );
        void                                    writeI1( signed char i );
        void                                    writeI2( signed short i );
        void                                    writeI4( int i );
        void                                    writeU1( unsigned char i );
        void                                    writeU2( unsigned short i );
        void                                    writeU4( unsigned int i );
        void                                    writeF4( float f );
        void                                    writeS0( const char* s );
        void                                    writeVX( unsigned int v );
        void                                    writeCOL12( const Vector3& c );
        void                                    writeVEC12( const Vector3& v );

        void                                    loadLW5();
        void                                    loadLW6();
        void                                    saveLW6();
};

//=============================================================================
//=============================================================================

#endif

---