#ifndef SOMVECTOR_INCLUDED
#define SOMVECTOR_INCLUDED

#ifndef SOMVECTOR_MATH
#define SOMVECTOR_MATH std::
#endif
#ifndef SOMVECTOR_SIN
#define SOMVECTOR_SIN SOMVECTOR_MATH sin
#endif
#ifndef SOMVECTOR_TAN
#define SOMVECTOR_TAN SOMVECTOR_MATH tan
#endif
#ifndef SOMVECTOR_COS
#define SOMVECTOR_COS SOMVECTOR_MATH cos
#endif
#ifndef SOMVECTOR_ACOS
#define SOMVECTOR_ACOS SOMVECTOR_MATH acos
#endif
#ifndef SOMVECTOR_SQRT
#define SOMVECTOR_SQRT SOMVECTOR_MATH sqrt
#endif
#ifndef SOMVECTOR_PI
#define SOMVECTOR_PI 3.14159265358979323846
#endif

  /*///////////////// WARNING //////////////////

	Deriving from these classes requires Empty
	Base Optimization. Double check alignments
	We cannot use compilers that do not do EBO

  /*////////////////////////////////////////////

class Somvector //namespace
{	
public:	//matrices are column-major order
	   //So the rules you'll find in most
	  //books are going to be upside down

	template<class,size_t,size_t,size_t> class matrix;

	////WARNING for MSVC2005 developers //////////////
	//
	// See vector::premultiply for an example of a bug
	// that may lie dormant among the unused templates

	//vector (zero sized container)
	template<class S, size_t N> class vector
	{		
	private: friend class Somvector; friend class vector;

		typedef S scalarN[N]; typedef vector<S,N> vectorN; 

		//these do not enforce const-ness and will conflict
		//with any new operators--which we don't want anyway
		//[enforcing const-ness will break builtin operators]
		inline operator scalarN&(){ return *(scalarN*)this; }
		inline operator scalarN&()const{ return *(scalarN*)this; }

	protected: vector<S,N>(){} vector<S,N>(const vector<S,N>&){}

	public:	typedef S scalar; static const size_t N = N; 
		
		//sub element (single) selector
		template<size_t a> inline scalar &se() //const 
		{
			int compile[a<N]; return (*this)[a];
		}
		template<size_t a> inline const scalar &se()const
		{
			int compile[a<N]; return (*this)[a];
		}	

		//sub element range (serial) selector
		template<size_t a, size_t z> inline vector<S,z-a+1> &ser(size_t len) //const
		{
			int compile[z-a<N]; assert(len==z-a+1); return *(vector<S,z-a+1>*)(*this+a);
		}
		template<size_t a, size_t z> inline const vector<S,z-a+1> &ser(size_t len)const
		{
			int compile[z-a<N]; assert(len==z-a+1); return *(vector<S,z-a+1>*)(*this+a);
		}	
		
		//static method (map a vector to memory)
		template<size_t Z> static inline vectorN &map(scalar (&mem)[Z]) //const
		{
			int compile[Z>=N]; return *(vectorN*)mem; 
		}
		template<size_t Z> static inline const vectorN &map(const scalar (&mem)[Z])
		{
			int compile[Z>=N]; return *(vectorN*)mem; 
		}
		template<size_t Z> static inline vectorN &map(scalar* &mem) //const
		{
			int compile[Z>=N]; return *(vectorN*)mem; 
		}
		template<size_t Z> static inline const vectorN &map(const scalar* &mem)
		{
			int compile[Z>=N]; return *(vectorN*)mem; 
		} 	

		//map acquisition (non-static)
		template<size_t Z> inline scalarN &map() //const
		{
			int compile[Z==N]; return *this;  
		}
		template<size_t Z> inline const scalarN &map()const
		{
			int compile[Z==N]; return *this;			
		}

		//WARNING: traditionally set is used to setup the values of the vector
		//That is NOT what is going on here. Unfortunately "objectifying" data 
		//introduces a question of semantics. Is the object copying the input?
		//Or is the input copying the object? So we're trying to be consistent
		//here (an entire afternoon was burned up in search of a better name!)

		//write to memory 
		template<size_t n, size_t Z, class T> inline void set(T (&v)[Z])const
		{
			int compile[n<=N&&z<=Z]; for(size_t i=0;i<n;i++) v[i] = (*this)[i]; 
		}	
		template<size_t n, size_t Z, class T> inline void set(T *v)const
		{
			int compile[n<=N&&z<=Z]; for(size_t i=0;i<n;i++) v[i] = (*this)[i]; 
		}	

		//extra friendly copy, move, scale
		#define SONVECTOR_LOOP_OP_(func,op) \
		template<size_t n, class T> inline vectorN &func(const T &v)\
		{\
			int compile[n<=N&&n<=T::N]; for(size_t i=0;i<n;i++) (*this)[i] op v[i]; return *this;\
		}\
		template<size_t n, size_t Z, class T> inline vectorN &func(const T (&v)[Z])\
		{\
			int compile[n<=N&&n<=Z]; for(size_t i=0;i<n;i++) (*this)[i] op v[i]; return *this;\
		}\
		template<size_t n, size_t Z, class T> inline vectorN &func(const T *v)\
		{\
			int compile[n<=N&&n<=Z]; for(size_t i=0;i<n;i++) (*this)[i] op v[i]; return *this;\
		}\
		/*scalar extended convention (the naming is deliberately misleading)*/\
		template<size_t n, class T> inline vectorN &se_##func(const T &s)\
		{\
			int compile[n<=N]; for(size_t i=0;i<n;i++) (*this)[i] op s; return *this;\
		}	
		SONVECTOR_LOOP_OP_(copy,=)
		SONVECTOR_LOOP_OP_(move,+=)
		SONVECTOR_LOOP_OP_(scale,*=)
		#undef SONVECTOR_LOOP_OP_

		//negation
		template<size_t n> inline vectorN &flip()
		{
			int compile[n<=N]; for(size_t i=0;i<n;i++) (*this)[i] = -(*this)[i]; return *this;
		}

		//limitation
		template<size_t n, class T> inline vectorN &limit(const T &min, const T &max)
		{
			int compile[n<=N]; for(size_t i=0;i<n;i++) 				
			if((*this)[i]<min) (*this)[i] = min; else if((*this)[i]>max) (*this)[i] = max;			
			return *this;
		}
		template<size_t n, class T> inline vectorN &limit(const T &symmetric_min_max)
		{
			return limit<n>(-symmetric_min_max,+symmetric_min_max);
		}

		//cross product
		template<size_t n, class T> inline vectorN &cross(const T &v)
		{
			//// unsure if this generalizes to other dimensions ////
		
			scalar compile[n==3&&n<=N&&n<=v.N], w[3]; map(w).copy(*this); 
			
			(*this)[0] = w[1]*v[2]-w[2]*v[1]; 
			(*this)[1] = w[2]*v[0]-w[0]*v[2]; 
			(*this)[2] = w[0]*v[1]-w[1]*v[0]; return *this;
		}		
		
		//dot product
		template<size_t n, class T> inline scalar dot(const T &v)const
		{
			scalar compile[n<=N&&n<=v.N], o = 0;
			
			for(size_t i=0;i<n;i++) o+=(*this)[i]*v[i]; return o;
		}	

		//normalization	
		template<size_t n> inline vectorN &unit()
		{		
			scalar rcp = length<n>(); if(rcp) rcp = scalar(1)/rcp;

			for(size_t i=0;i<n;i++) (*this)[i]*=rcp; return *this;
		}	
		template<size_t n> inline scalar length()const
		{
			return SOMVECTOR_SQRT(dot<n>(*this));
		}
								
		//transformation
		template<size_t n, class T> inline vectorN &premultiply(const T &t)
		{
			scalar cp[N]; Somvector::multiply<n>(map(cp).copy<N>(*this),t,*this);
			return *this;
		}
		template<size_t n, class T> inline vectorN &postmultiply(const T &t)
		{
			int compile[n==N]; //TODO: implement homogeneous extension

			scalar cp[t.N][1]; matrix<S,t.N,1> &col = transpose<t.N>(); 

			 //MSVC2005 (bug)
			//Somvector::multiply<n,1>(t,col.map(cp).copy<t.N,1>(col),col); 
			Somvector::multiply<n,1>(t,col.map(cp).copy<T::N,1>(col),col); 
			return *this;
		}	
		
		//transposition
		template<size_t m> inline matrix<S,m,1,1> &transpose() //const
		{
			int compile[m<=N]; return matrix<S,m,1,1>::map<m,1>(*this);
		}
		template<size_t m> inline const matrix<S,m,1,1> &transpose()const
		{
			int compile[m<=N]; return matrix<S,m,1,1>::map<m,1>(*this);
		}

		/*//// quaternions (taking "quat" to imply 4D) ////////*/

		//TODO: add generalized solutions from Somimp's Euler.hpp 

		//standard construction methods
		template<class S, size_t Z> inline vectorN &quaternion(const vector<S,Z> &Euler)
		{
			int compile[N==4&&Z>=3]; const scalar _5 = .5, //Euler = { yaw, pitch, roll }

			s0(SOMVECTOR_SIN(_5*Euler[0])), s1(SOMVECTOR_SIN(_5*Euler[1])),	//in radians
			s2(SOMVECTOR_SIN(_5*Euler[2])), c0(SOMVECTOR_COS(_5*Euler[0])), s0_s2(s0*s2),			
			c1(SOMVECTOR_COS(_5*Euler[1])), c2(SOMVECTOR_COS(_5*Euler[2])), c0_c2(c0*c2);

			(*this)[0] = s1*c0_c2-c1*s0_s2;	(*this)[1] = c1*s0*c2+s1*c0*s2;
			(*this)[2] = c1*c0*s2-s1*s0*c2;	(*this)[3] = c1*c0_c2+s1*s0_s2; return *this;
		}
		template<class S, size_t Z> inline vectorN &quaternion(const S (&Euler)[Z])
		{
			return quaternion(vector<S,Z>::map(Euler)); //treat scalar arrays like Euler angles
		}
		template<class S, class T, size_t Z> inline vectorN &quaternion(const vector<S,Z> &unit_axis, const T angle)
		{
			int compile[N==4&&Z>=3]; const scalar rad = scalar(0.5)*angle; //in radians

			(*this)[3] = SOMVECTOR_COS(rad); return copy<3>(unit_axis).se_scale<3>(SOMVECTOR_SIN(rad));
		}
		template<class S, class T, size_t W, size_t Z> inline vectorN &quaternion(const vector<S,W> &s, const vector<T,Z> &t)
		{			
			scalar cp[4]; return quaternion(map(cp).copy<3>(s).cross(t).unit<3>(),SOMVECTOR_ACOS(s.dot(t)));
		}		
		template<class T, size_t W, size_t Z, size_t ZZ> inline vectorN &quaternion(const matrix<T,W,Z,ZZ> &t) 
		{			
			int compile[N==4&&t.M>=3&&t.N>=3]; //TODO: non-trivial test

			scalar trace = t.trace<3>(), square = SOMVECTOR_SQRT(trace)*2;

			if(trace>0) //non-special case
			{
				(*this)[3] = square*scalar(0.25); square = scalar(1)/square;

				(*this)[0] = square*(t[1][2]-t[2][1]); (*this)[1] = square*(t[2][0]-t[0][2]);
				(*this)[2] = square*(t[0][1]-t[1][0]); 
			}
			else if(t[0][0]>t[1][1]&&t[0][0]>t[2][2])
			{	
				square = SOMVECTOR_SQRT(scalar(1)+t[0][0]-t[1][1]-t[2][2])*2;

				(*this)[0] = square*scalar(0.25); square = scalar(1)/square;

				(*this)[1] = square*(t[0][1]+t[1][0]); (*this)[2] = square*(t[2][0]+t[0][2]);
				(*this)[3] = square*(t[1][2]-t[2][1]);
			} 
			else if(t[1][1]>t[2][2]) 
			{ 
				square = SOMVECTOR_SQRT(scalar(1)+t[1][1]-t[0][0]-t[2][2])*2;

				(*this)[1] = square*scalar(0.25); square = scalar(1)/square;

				(*this)[0] = square*(t[0][1]+t[1][0]); (*this)[2] = square*(t[1][2]+t[2][1]);
				(*this)[3] = square*(t[2][0]-t[0][2]);
			} 
			else 
			{ 
				square = SOMVECTOR_SQRT(scalar(1)+t[2][2]-t[0][0]-t[1][1])*2;

				(*this)[2] = square*scalar(0.25); square = scalar(1)/square;

				(*this)[0] = square*(t[2][0]+t[0][2]); (*this)[1] = square*(t[1][2]+t[2][1]);
				(*this)[3] = square*(t[0][1]-t[1][0]);
			}
			return *this;
		}		

		//rotation about (a quaternion)
		template<size_t n, class T> inline vectorN &rotate(const vector<T,4> &q)
		{
			scalar cp[4], cq[4] = {-q[0],-q[1],-q[2], q[3]}; //conjugate

			int compile[n==3]; vector<S,3> &_this = ser<0,2>(n); vector<S,4> s4; //tricky

			Somvector::multiply(Somvector::multiply(q,_this,s4.map(cp)),s4.map(cq),_this);

			return *this;
		}

		//quaternion quaternion product
		template<class T> inline vectorN &postmultiply(const vector<T,4> &t)
		{				
			scalar cp[4]; return Somvector::multiply(t,map(cp).copy<4>(*this),*this);
		}
		template<class T> inline vectorN &premultiply(const vector<T,4> &t)
		{
			scalar cp[4]; return Somvector::multiply(map(cp).copy<4>(*this),t,*this);
		}
	};

	//matrix (zero sized container)
	template<class S, size_t M, size_t N, size_t NN=N> class matrix
	{
	private: friend class Somvector; friend class matrix;

		typedef S scalarMxN[M][NN]; typedef matrix<S,M,N,NN> matrixMxN; 

		//these do not enforce const-ness and will conflict
		//with any new operators--which we don't want anyway
		//[enforcing const-ness will break builtin operators]
		inline operator scalarMxN&(){ return *(scalarMxN*)this; }
		inline operator scalarMxN&()const{ return *(scalarMxN*)this; }

	protected: matrix<S,M,N,NN>(){ int compile[N<=NN]; }; 
		
		matrix<S,M,N,NN>(const matrix<S,M,N,NN>&){ int compile[N<=NN]; }; 
		
	public:	typedef S scalar; static const size_t M = M, N = N, NN = NN;

		//sub element (single) selector
		template<size_t a, size_t aa> inline scalar &se() //const 
		{
			int compile[a<M&&aa<N]; return (*this)[a][aa];
		}
		template<size_t a, size_t aa> inline const scalar &se()const
		{
			int compile[a<M&&aa<N]; return (*this)[a][aa];
		}	

		//sub element range (serial) selector
		template<size_t a, size_t z, size_t aa, size_t zz> 
		inline  matrix<S,z-a+1,zz-aa+1,NN> &ser(size_t mm, size_t nn) //const
		{
			int compile[z-a<M&&zz-aa<N]; assert(mm==z-a+1&&nn==zz-aa+1); 
			
			return *(matrix<S,z-a+1,zz-aa+1,NN>*)&(*this)[a][aa];
		}
		template<size_t a, size_t z, size_t aa, size_t zz> 
		inline const matrix<S,z-a+1,zz-aa+1,NN> &ser(size_t mm, size_t nn)const 
		{
			int compile[z-a<M&&zz-aa<N]; assert(mm==z-a+1&&nn==zz-aa+1); 
			
			return *(matrix<S,z-a+1,zz-aa+1,NN>*)&(*this)[a][aa];
		}

		//map a matrix to memory (static)
		#define SOMVECTOR_MATRIX_MAP_(T,Test,...) \
		template<size_t W, size_t Z> static inline __VA_ARGS__ &map(T)\
		{\
			int compile[W>=M&&Z>=N&&Test]; return *(__VA_ARGS__*)mem;\
		}\
		template<size_t W, size_t Z> static inline const __VA_ARGS__ &map(const T)\
		{\
			int compile[W>=M&&Z>=N&&Test]; return *(__VA_ARGS__*)mem;\
		}
		//(if you derive a subclass you may need to copy these)
		SOMVECTOR_MATRIX_MAP_(scalar (&mem)[W][Z],1,matrix<S,M,N,Z>)
		SOMVECTOR_MATRIX_MAP_(scalar (&mem)[W],N==1,matrix<S,M,N,Z>)
		SOMVECTOR_MATRIX_MAP_(scalar* &mem,     W*Z,matrix<S,M,N,Z>)
			
		//map acquisition (non-static)
		template<size_t W, size_t Z> inline scalarMxN &map() //const
		{
			int compile[W==M&&Z==NN]; return *this;  
		}
		template<size_t W, size_t Z> inline const scalarMxN &map()const
		{
			int compile[W==M&&Z==NN]; return *this; 
		}

		#define W_Z_T size_t W,size_t Z,class T //lazy: for legibility sake

		//WARNING: traditionally set is used to setup the values of the matrix
		//That is NOT what is going on here. Unfortunately "objectifying" data 
		//introduces a question of semantics. Is the object copying the input?
		//Or is the input copying the object? So we're trying to be consistent
		//here (an entire afternoon was burned up in search of a better name!)

		//write matrix to memory 
		template<size_t m, size_t n, W_Z_T> inline void set(T (&t)[W][Z])const
		{
			int compile[m<=M&&m<=W&&n<=N&&n<=Z];		
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) t[i][j] = (*this)[i][j]; 
		}	
		template<size_t m, size_t n, W_Z_T> inline void set(T t[W*Z])const
		{
			int compile[m<=M&&m<=W&&n<=N&&n<=Z];		
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) t[Z*i+j] = (*this)[i][j]; 
		}	
		
		//matrix matrix copy assignment
		template<size_t m, size_t n, class T> inline matrixMxN &copy(const T &t)
		{	
			int compile[m<=M&&m<=t.M&&n<=N&&n<=t.N];
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) (*this)[i][j] = t[i][j]; 		
			return *this;
		}	
		template<size_t m, size_t n, W_Z_T> inline matrixMxN &copy(const T (&t)[W][Z])
		{	
			int compile[m<=M&&m<=W&&n<=N&&n<=Z];
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) (*this)[i][j] = t[i][j];
			return *this;
		}	
		template<size_t m, size_t n, W_Z_T> inline matrixMxN &copy(const T t[W*Z])
		{	
			int compile[m<=M&&m<=W&&n<=N&&n<=Z];
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) (*this)[i][j] = t[W*i+j];
			return *this;
		}	

		//matrix quaternion copy conversion
		template<size_t m, size_t n, class T> inline matrixMxN &copy_quaternion(const vector<T,4> &q)
		{	
			int compile[m<=M&&m<=4&&n<=N&&n<=4]; 
			
			if(m>=1){
			if(n>=1) (*this)[0][0] = scalar(1)-scalar(2)*(q[1]*q[1]+q[2]*q[2]);
			if(n>=2) (*this)[0][1] = scalar(2)          *(q[0]*q[1]+q[2]*q[3]);
			if(n>=3) (*this)[0][2] = scalar(2)          *(q[0]*q[2]-q[1]*q[3]);
			if(n>=4) (*this)[0][3] = 0;}
			if(m>=2){
			if(n>=1) (*this)[1][0] = scalar(2)          *(q[0]*q[1]-q[2]*q[3]);
			if(n>=2) (*this)[1][1] = scalar(1)-scalar(2)*(q[0]*q[0]+q[2]*q[2]);
			if(n>=3) (*this)[1][2] = scalar(2)          *(q[1]*q[2]+q[0]*q[3]);
			if(n>=4) (*this)[1][3] = 0;}
			if(m>=3){
			if(n>=1) (*this)[2][0] = scalar(2)          *(q[0]*q[2]+q[1]*q[3]);
			if(n>=2) (*this)[2][1] = scalar(2)          *(q[1]*q[2]-q[0]*q[3]);			
			if(n>=3) (*this)[2][2] = scalar(1)-scalar(2)*(q[0]*q[0]+q[1]*q[1]);
			if(n>=4) (*this)[2][3] = 0;}
			if(m>=4){
			if(n>=1) (*this)[3][0] = 0;	if(n>=2) (*this)[3][1] = 0;
			if(n>=3) (*this)[3][2] = 0;	if(n>=4) (*this)[3][3] = 1;}
			
			return *this; 
		}	

		//identity matrix
		template<size_t m, size_t n> inline matrixMxN &identity() 
		{				
			int compile[m<=M&&n<=N]; 
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) (*this)[i][j] = i==j?1:0; 
			return *this;
		}
		template<size_t m> inline scalar trace()const
		{
			scalar sum = 0; for(size_t i=0;i<m;i++) sum+=(*this)[i][i]; return sum;
		}

		//matrix matrix product
		template<size_t m, size_t n, class T> inline matrixMxN &premultiply(const T &t)
		{
			scalar cp[m][N]; Somvector::multiply<m,n>(matrix<S,m,N>::map(cp).copy<m,N>(*this),t,*this);
			return *this;
		}
		template<size_t m, size_t n, class T> inline matrixMxN &postmultiply(const T &t)
		{
			scalar cp[M][n]; Somvector::multiply<m,n>(t,matrix<S,M,n>::map(cp).copy<M,n>(*this),*this);
			return *this;
		}

		//transpose matrix
		template<size_t m, size_t n> inline matrixMxN &transpose() 
		{
			scalar cp[n][m]; Somvector::transpose<m,n>(matrix<S,n,m>::map(cp).copy<n,m>(*this),*this);
			return *this;
		}

		//matrix inversion		
		template<size_t m, size_t n> inline matrixMxN &invert()
		{
			int compile[m<=M&&n<=N&&m==n]; scalar cof[m][n];
			
			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) cof[i][j] = cofactor<m,n>(i,j);
											   
			scalar rcpdet = vector<S,n>::map(cof[0]).dot<n>(vector<S,n>::map((*this)[0]));
			
			if(rcpdet) rcpdet = scalar(1)/rcpdet; //(there is a possible divide by zero here)

			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) (*this)[i][j] = cof[j][i]*rcpdet; 

			return *this;
		}
		template<size_t m, size_t n> inline scalar cofactor(size_t w, size_t z)const
		{
			int compile[m<=M&&n<=N&&m==n&&m>=2]; 

			scalar flip = w&1^z&1?-1:1; assert(w<m&&z<n);			

			if(m==2) return (*this)[!w][!z]*flip; //trivial
			
			const size_t l = m-1; //square matrix (m equals n)

			scalar sub[l][l]; size_t k = 0; //may not unroll so well

			for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) 
			{
				if(i!=w&&j!=z) sub[0][k++] = (*this)[i][j]; //sub[l*l] 
			}			

			return matrix<S,l,l>::map(sub).determinant<l,l>()*flip;
		}
		template<size_t m, size_t n> inline scalar determinant()const
		{
			int compile[m<=M&&n<=N&&m==n&&m>=1&&m<=4]; 

			switch(m) //column-major
			{
			case 1: return (*this)[0][0];
			case 2: return (*this)[0][0]*(*this)[1][1]
			              -(*this)[1][0]*(*this)[0][1];
			case 3: return (*this)[0][0]*(*this)[1][1]*(*this)[2][2]
						  -(*this)[0][0]*(*this)[2][1]*(*this)[1][2]
			              +(*this)[1][0]*(*this)[2][1]*(*this)[0][2]
						  -(*this)[1][0]*(*this)[0][1]*(*this)[2][2]
			              +(*this)[2][0]*(*this)[0][1]*(*this)[1][2]
			              -(*this)[2][0]*(*this)[1][1]*(*this)[0][2]; 
			case 4: return (*this)[3][0]*(*this)[2][1]*(*this)[1][2]*(*this)[0][3]
			              -(*this)[2][0]*(*this)[3][1]*(*this)[1][2]*(*this)[0][3]
			              -(*this)[3][0]*(*this)[1][1]*(*this)[2][2]*(*this)[0][3]
			              +(*this)[1][0]*(*this)[3][1]*(*this)[2][2]*(*this)[0][3]
			              +(*this)[2][0]*(*this)[1][1]*(*this)[3][2]*(*this)[0][3]
			              -(*this)[1][0]*(*this)[2][1]*(*this)[3][2]*(*this)[0][3]
			              -(*this)[3][0]*(*this)[2][1]*(*this)[0][2]*(*this)[1][3]
			              +(*this)[2][0]*(*this)[3][1]*(*this)[0][2]*(*this)[1][3]
			              +(*this)[3][0]*(*this)[0][1]*(*this)[2][2]*(*this)[1][3]
			              -(*this)[0][0]*(*this)[3][1]*(*this)[2][2]*(*this)[1][3]
			              -(*this)[2][0]*(*this)[0][1]*(*this)[3][2]*(*this)[1][3]
			              +(*this)[0][0]*(*this)[2][1]*(*this)[3][2]*(*this)[1][3]
			              +(*this)[3][0]*(*this)[1][1]*(*this)[0][2]*(*this)[2][3]
			              -(*this)[1][0]*(*this)[3][1]*(*this)[0][2]*(*this)[2][3]
			              -(*this)[3][0]*(*this)[0][1]*(*this)[1][2]*(*this)[2][3]
			              +(*this)[0][0]*(*this)[3][1]*(*this)[1][2]*(*this)[2][3]
			              +(*this)[1][0]*(*this)[0][1]*(*this)[3][2]*(*this)[2][3]
			              -(*this)[0][0]*(*this)[1][1]*(*this)[3][2]*(*this)[2][3]
			              -(*this)[2][0]*(*this)[1][1]*(*this)[0][2]*(*this)[3][3]
			              +(*this)[1][0]*(*this)[2][1]*(*this)[0][2]*(*this)[3][3]
			              +(*this)[2][0]*(*this)[0][1]*(*this)[1][2]*(*this)[3][3]
			              -(*this)[0][0]*(*this)[2][1]*(*this)[1][2]*(*this)[3][3]
			              -(*this)[1][0]*(*this)[0][1]*(*this)[2][2]*(*this)[3][3]
						  +(*this)[0][0]*(*this)[1][1]*(*this)[2][2]*(*this)[3][3];

			default: return 0; //compiler
			}
		}			  
		
		//glFrustum
		template<size_t m, size_t n> inline matrixMxN &frustum
		(scalar left, scalar right, scalar bottom, scalar top, scalar znear, scalar zfar) 
		{
			int compile[M==4&&N==4&&m==4&&n==4];
						
			//right handed (z-clipping is confused)
			(*this)[0][0] = (znear*2)/(right-left);
			(*this)[0][1] = 0;
			(*this)[0][2] = 0;
			(*this)[0][3] = 0;
			(*this)[1][0] = 0;
			(*this)[1][1] = (znear*2)/(top-bottom);
			(*this)[1][2] = 0;
			(*this)[1][3] = 0;
			(*this)[2][0] = (right+left)/(right-left);
			(*this)[2][1] = (top+bottom)/(top-bottom);
			//hmm: would prefer to go with OpenGL here???
			//(*this)[2][2] = -(zfar+znear)/(zfar-znear); //OpenGL
			(*this)[2][2] = zfar/(znear-zfar); //Direct3D
			(*this)[2][3] = -1;
			(*this)[3][0] = 0;
			(*this)[3][1] = 0;			
			//(*this)[3][2] = -(zfar*znear*2)/(zfar-znear); //OpenGL
			(*this)[3][2] = zfar*znear/(znear-zfar); //Direct3D
			(*this)[3][3] = 0;
			return *this;
		};
		template<size_t m, size_t n> //gluPerspective
		inline matrixMxN &perspective(scalar fovy, scalar aspect, scalar znear, scalar zfar) 
		{			
			scalar v = znear*SOMVECTOR_TAN(fovy>SOMVECTOR_PI+0.000003?fovy*SOMVECTOR_PI/360:fovy);

			scalar u = v*aspect; return frustum<m,n>(-u,u,-v,v,znear,zfar);
		};
	};
				
	//multiplication
	template<size_t m, size_t n, class T, class U, class V> 
	static inline V &multiply(const T &t, const U &u, V &v) //matrix matrix product
	{
		int conform[u.M==t.N];
		int compile[m<=t.M&&n<=u.N&&m<=v.M&&n<=v.N]; 
		for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) //row*column dot product
		{
			v[i][j] = 0; for(size_t k=0;k<u.N;k++) v[i][j]+=t[i][k]*u[k][j];
		}
		return v;
	}	
	template<size_t n, class T, class U, class V> 
	static inline V &multiply(const T &t, const U &u, V &v) //vector matrix product
	{	
		const size_t xx = n==v.N+1?v.N:n;	
		
		int conform[u.M-t.N<=1]; 
		int compile[xx<=v.N&&xx<=t.N&&u.M>=n]; 

		typedef typename T::scalar scalar;

		//homogeneous extension where t is 1 less u
		const scalar w = t.N==u.N?t[t.N-1]:scalar(1);

		for(size_t i=0;i<1;i++)
		for(size_t j=0;j<xx;j++) //row*column dot product
		{
			v[j] = xx!=n&&xx<t.N?w*u[xx][j]:0; //extended
			
			for(size_t k=0;k<t.N;k++) v[j]+=t[k]*u[k][j];
		}

		//homogeneous extension where v is 1 more t
		if(xx==n&&t.N==v.N-1) v[n-1] = t.N==u.N?0:1; 
		
		return v;
	}
	template<class T, class U, class V> 
	static inline V &multiply(const vector<T,4> &t, const U &u, V &v) //quaternions 
	{
		typedef typename U::scalar scalar;

		//u and v can be 3D vectors for implementing rotations
		int compile[u.N>=3]; const scalar u_3_ = u.N<4?0:u[3];

		if(v.N>=1) v[0] = t[3]*u[0]+t[0]*u_3_+t[1]*u[2]-t[2]*u[1];
		if(v.N>=2) v[1] = t[3]*u[1]+t[1]*u_3_+t[2]*u[0]-t[0]*u[2];
		if(v.N>=3) v[2] = t[3]*u[2]+t[2]*u_3_+t[0]*u[1]-t[1]*u[0];
		if(v.N>=4) v[3] = t[3]*u_3_-t[0]*u[0]-t[1]*u[1]-t[2]*u[2];

		return v;
	}

	//transposition 
	template<size_t m, size_t n, class T, class U> 
	static inline U &transpose(const T &t, U &u) //transpose matrix
	{
		int compile[m<=t.N&&n<=t.M&&m<=u.M&&n<=u.N]; 
		for(size_t i=0;i<m;i++) for(size_t j=0;j<n;j++) u[i][j] = t[j][i]; 
		return u;
	}

	//serialization (define as many as you need)
	template<size_t N, class T, class X, class Y, class Z> 
	static inline vector<T,3> &series(T (&t)[N], const X &m, const Y &n, const Z &z)
	{int compile[N>=3]; vector<T,3> &s = s.map(t); s[0] = m; s[1] = n; s[2] = z; return s;}
	template<size_t N, class T, class X, class Y, class Z> 
	static inline vector<T,3> &series(T t[N], const X &m, const Y &n, const Z &z)
	{int compile[N>=3]; vector<T,3> &s = s.map<N>(t); s[0] = m; s[1] = n; s[2] = z; return s;}
	template<size_t N, class T, class X, class Y, class Z, class W> 
	static inline vector<T,4> &series(T (&t)[N], const X &m, const Y &n, const Z &z, const W &w)
	{int compile[N>=4]; vector<T,4> &s = s.map(t); s[0] = m; s[1] = n; s[2] = z; s[3] = w; return s;}
	template<size_t N, class T, class X, class Y, class Z, class W> 
	static inline vector<T,4> &series(T t[N], const X &m, const Y &n, const Z &z, const W &w)
	{int compile[N>=4]; vector<T,4> &s = s.map<N>(t); s[0] = m; s[1] = n; s[2] = z; s[3] = w; return s;}
	//maybe we can recursive macro these somehow
};

#endif //SOMVECTOR_INCLUDED