// takakuMatrix.h
#ifndef TAKAKU_MATRIX_H
#define TAKAKU_MATRIX_H


#include <takakuDefines.h>
#include <cmath>
#include <cstdarg>


namespace takaku
{


template <class Ty> struct __StackCast
{
   typedef Ty __Type;
};

template <> struct __StackCast<signed char>
{
   typedef signed int __Type;
};

template <> struct __StackCast<signed short>
{
   typedef signed int __Type;
};

template <> struct __StackCast<unsigned char>
{
   typedef unsigned int __Type;
};

template <> struct __StackCast<unsigned short>
{
   typedef unsigned int __Type;
};

template <> struct __StackCast<float>
{
   typedef double __Type;
};

template <class Ty,int rows,int cols> struct Matrix
{
   Ty m[rows][cols];

   typedef Ty TYPE;
   static const int ROWS = rows;
   static const int COLS = cols;

   Matrix()
   {
       for (int i=0; i<rows*cols; i++)
           ((Ty *)m)[i] = 0;
   }

   explicit Matrix(Ty * pArr)
   {
       CopyArray((Ty *)m, pArr, rows*cols);
   }

   explicit Matrix(Ty * *ppArr)
   {
       for (int i=0; i<rows; i++)
           CopyArray(m[i], ppArr[i], cols);
   }

   template <class _T> explicit Matrix(_T first, ...)
   {
       va_list args;
       va_start(args, first);
       m[0][0] = first;
       for (int i=1; i<rows*cols; i++)
           m[i/cols][i%cols] = (Ty) (va_arg(args, __StackCast<_T>::__Type));
       va_end(args);
   }

   const Ty & operator () (int i, int j)const
   {
       TAKAKU_ASSERT( i>=0 && i<rows && j>=0 && j<cols );
       return m[i][j];
   }

   Ty & operator () (int i, int j)
   {
       TAKAKU_ASSERT( i>=0 && i<rows && j>=0 && j<cols );
       return m[i][j];
   }

   Matrix<Ty,cols,rows> T()const
   {
       Matrix<Ty,cols,rows> ret;
       for(int i=0; i<rows; i++)
           for(int j=0; j<cols; j++)
               ret.m[j][i] = m[i][j];
       return ret;
   }
};

template <class Ty,int rows> struct Matrix<Ty,rows,1>
{
   Ty m[rows][1];

   typedef Ty TYPE;
   static const int ROWS = rows;
   static const int COLS = 1;

   Matrix()
   {
       for (int i=0; i<rows; i++)
           ((Ty *)m)[i] = 0;
   }

   explicit Matrix(Ty * pArr)
   {
       CopyArray((Ty *)m, pArr, rows);
   }

   explicit Matrix(Ty * *ppArr)
   {
       for (int i=0; i<rows; i++)
           CopyArray(m[i], ppArr[i], 1);
   }

   template <class _T> explicit Matrix(_T first, ...)
   {
       va_list args;
       va_start(args, first);
       m[0][0] = first;
       for (int i=1; i<rows; i++)
           m[i][0] = (Ty) (va_arg(args, __StackCast<_T>::__Type));
       va_end(args);
   }

   inline const Ty & operator () (int i, int j)const
   {
       TAKAKU_ASSERT( i>=0 && i<rows && j>=0 && j<1 );
       return m[i][0];
   }

   inline Ty & operator () (int i, int j)
   {
       TAKAKU_ASSERT( i>=0 && i<rows && j>=0 && j<1 );
       return m[i][0];
   }

   inline const Ty & operator () (int i)const
   {
       TAKAKU_ASSERT( i>=0 && i<rows );
       return m[i][0];
   }

   inline Ty & operator () (int i)
   {
       TAKAKU_ASSERT( i>=0 && i<rows );
       return m[i][0];
   }

   Matrix<Ty,1,rows> T()const
   {
       Matrix<Ty,1,rows> ret;
       for(int i=0; i<rows; i++)
           ret.m[0][i] = m[i][0];
       return ret;
   }
};

template <class Ty,int cols> struct Matrix<Ty,1,cols>
{
   Ty m[1][cols];

   typedef Ty TYPE;
   static const int ROWS = 1;
   static const int COLS = cols;

   Matrix()
   {
       for (int i=0; i<cols; i++)
           ((Ty *)m)[i] = 0;
   }

   explicit Matrix(Ty * pArr)
   {
       CopyArray((Ty *)m, pArr, cols);
   }

   explicit Matrix(Ty * *ppArr)
   {
       CopyArray(m[0], ppArr[0], cols);
   }

   template <class _T> explicit Matrix(_T first, ...)
   {
       va_list args;
       va_start(args, first);
       m[0][0] = first;
       for (int i=1; i<cols; i++)
           m[0][i] = (Ty) (va_arg(args, __StackCast<_T>::__Type));
       va_end(args);
   }

   inline const Ty & operator () (int i, int j)const
   {
       TAKAKU_ASSERT( i>=0 && i<1 && j>=0 && j<cols );
       return m[0][i];
   }

   inline Ty & operator () (int i, int j)
   {
       TAKAKU_ASSERT( i>=0 && i<1 && j>=0 && j<cols );
       return m[0][i];
   }

   inline const Ty & operator () (int i)const
   {
       TAKAKU_ASSERT( i>=0 && i<cols );
       return m[0][i];
   }

   inline Ty & operator () (int i)
   {
       TAKAKU_ASSERT( i>=0 && i<cols );
       return m[0][i];
   }

   Matrix<Ty,cols,1> T()const
   {
       Matrix<Ty,cols,1> ret;
       for(int i=0; i<cols; j++)
           ret.m[i][0] = m[0][i];
       return ret;
   }
};

template <class Ty> struct Matrix<Ty, 1, 1> // disabled
{
   Matrix() { TAKAKU_BREAK };
};

//////////////////////////////////////////////////////////////////////////

template <class Ty,int rows,int colsA, int colsB> inline Matrix<Ty, rows, colsA + colsB>
operator || (const Matrix<Ty,rows,colsA> & a, const Matrix<Ty,rows,colsB> & b)
{
   Matrix<Ty, rows, colsA + colsB> ret;
   for(int I=0; I<rows; I++)
   {
       for(int i=0; i<colsA; i++)
           ret.m[I][i] = a.m[I][i];
       for(int i=colsA, j=0; j<colsB; i++, j++)
           ret.m[I][i] = b.m[I][j];
   }
   return ret;
}

template <class Ty,int rowsA,int rowsB, int cols> inline Matrix<Ty, rowsA + rowsB, cols>
operator && (const Matrix<Ty,rowsA,cols> & a, const Matrix<Ty,rowsA,cols> & b)
{
   Matrix<Ty, rowsA + rowsB, cols> ret;
   for(int i=0; i<rowsA*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) b.m) [i];
   for(int i=rowsA*cols, j=0; j<rowsB*cols; i++, j++)
       ((Ty *) ret.m) [i] = ((Ty *) b.m) [j];
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class Tdest, class Tsrc, int rows, int cols> inline Matrix<Tdest, rows, cols>
matrix_cast(Matrix<Tsrc, rows, cols> mat)
{
   Matrix<Tdest,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Tdest *) ret.m) [i] = (Tdest)( ((Tsrc *) mat.m) [i] );
   return ret;
}

template <class Tdest, class Tsrc, int rows, int cols> inline Matrix<Tdest, rows, cols>
matrix_static_cast(Matrix<Tsrc, rows, cols> mat)
{
   Matrix<Tdest,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Tdest *) ret.m) [i] = static_cast<Tdest>( ((Tsrc *) mat.m) [i] );
   return ret;
}

template <class Tdest, class Tsrc, int rows, int cols> inline Matrix<Tdest, rows, cols>
matrix_reinterpret_cast(Matrix<Tsrc, rows, cols> mat)
{
   Matrix<Tdest,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Tdest *) ret.m) [i] = reinterpret_cast<Tdest>( ((Tsrc *) mat.m) [i] );
   return ret;
}

template <class Tdest, class Tsrc, int rows, int cols> inline Matrix<Tdest, rows, cols>
matrix_const_cast(Matrix<Tsrc, rows, cols> mat)
{
   Matrix<Tdest,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Tdest *) ret.m) [i] = const_cast<Tdest>( ((Tsrc *) mat.m) [i] );
   return ret;
}

template <class Tdest, class Tsrc, int rows, int cols> inline Matrix<Tdest, rows, cols>
matrix_dynamic_cast(Matrix<Tsrc, rows, cols> mat)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Tdest *) ret.m) [i] = dynamic_cast<Tdest>( ((Tsrc *) mat.m) [i] );
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int rows,int cols> inline bool
operator == (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   return CompareArray((Ty *) a.m, (Ty *) b.m, rows*cols);
}

template <class Ty,int rows,int cols> inline bool
operator != (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   return ! CompareArray((Ty *) a.m, (Ty *) b.m, rows*cols);
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator += (Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   for(int i = 0; i<rows*cols; i++)
       ((Ty *) a.m) [i] += ((Ty *) b.m) [i];
   return a;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator + (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) a.m) [i] + ((Ty *) b.m) [i];
   return ret;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator + (const Matrix<Ty,rows,cols> & mat)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = + ((Ty *) mat.m) [i];
   return ret;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator -= (Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   for(int i=0; i <rows*cols; i++)
       ((Ty *) a.m) [i] -= ((Ty *) b.m) [i];
   return a;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator - (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) a.m) [i] - ((Ty *) b.m) [i];
   return ret;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator - (const Matrix<Ty,rows,cols> & mat)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = - ((Ty *) mat.m) [i];
   return ret;
}

//////////////////////////////////////////////////////////////////////////
/*
template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator &= (Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   for(int i=0; i<rows*cols; i++)
       ((Ty *) a.m) [i] *= ((Ty *) b.m) [i];
   return a;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator & (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) a.m) [i] * ((Ty *) b.m) [i];
   return ret;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator |= (Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   for(int i=0; i<rows*cols; i++)
       ((Ty *) a.m) [i] /= ((Ty *) b.m) [i];
   return a;
}

template <class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator | (const Matrix<Ty,rows,cols> & a, const Matrix<Ty,rows,cols> & b)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) a.m) [i] / ((Ty *) b.m) [i];
   return ret;
}
*/
//////////////////////////////////////////////////////////////////////////

template <class Ty,int M,int S,int N> inline Matrix<Ty,M,N>
operator * (const Matrix<Ty,M,S> & a, const Matrix<Ty,S,N> & b)
{
   Matrix<Ty,M,N> ret; // fill all with 0 at first
   for(int i=0; i<M; i++)
       for(int j=0; j<N; j++)
           for(int k=0; k<S; k++)
               ret.m[i][j] += a.m[i][k] * b.m[k][j];
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class K,class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator *= (Matrix<Ty,rows,cols> & mat, K k)
{
   for(int i=0; i<rows*cols; i++)
       ((Ty *) mat.m) [i] *= k;
   return mat;
}

template <class K,class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator * (const Matrix<Ty,rows,cols> & mat, K k)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) mat.m) [i] * k;
   return ret;
}

template <class K,class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator * (K k, const Matrix<Ty,rows,cols> & mat)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = k * ((Ty *) mat.m) [i];
   return ret;
}

template <class K,class Ty,int rows,int cols> inline Matrix<Ty,rows,cols> &
operator /= (Matrix<Ty,rows,cols> & mat, K k)
{
   for(int i=0; i<rows*cols; i++)
       ((Ty *) mat.m) [i] /= k;
   return mat;
}

template <class K,class Ty,int rows,int cols> inline Matrix<Ty,rows,cols>
operator / (const Matrix<Ty,rows,cols> & mat, K k)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = ((Ty *) mat.m) [i] / k;
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int n> inline Matrix<Ty,n,n>
Identity()
{
   Matrix<Ty,n,n> ret; // fill all with 0
   for(int i=0; i<n; i++)
       ret.m[i][i] = 1;
   return ret;
}

template <class Ty,int rows,int cols> inline bool
IsIdentity(const Matrix<Ty,rows,cols> &mat)
{
   if(rows != cols) return false;
   for(int i=0; i<rows; i++)
       for(int j=0; j<cols; j++)
           if(i==j)
           {
               if(mat.m[i][j] != 1) return false;
           }
           Else
           {
               if(mat.m[i][j] != 0) return false;
           }
   return true;
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int dim> inline Ty
Norm(const Matrix<Ty,dim,1> & vec)
{
   Ty tmp = 0;
   for(int i=0; i<dim; i++)
       tmp += vec.m[i][0] * vec.m[i][0];
   return std::sqrt(tmp);
}

template <class F,class Ty,int dim> inline F    // floatでない場合、このoverloadを呼び出さないと警告される
Norm(const Matrix<Ty,dim,1> & vec)
{
   Ty tmp = 0;
   for(int i=0; i<dim; i++)
       tmp += vec.m[i][0] * vec.m[i][0];
   return (F)std::sqrt((F)tmp);
}

template <class Ty,int dim> inline Ty
Norm(const Matrix<Ty,1,dim> & vec)
{
   Ty tmp = 0;
   for(int i=0; i<dim; i++)
       tmp += vec.m[0][i] * vec.m[0][i];
   return std::sqrt(tmp);
}

template <class F,class Ty,int dim> inline F    // floatでない場合、このoverloadを呼び出さないと警告される
Norm(const Matrix<Ty,1,dim> & vec)
{
   Ty tmp = 0;
   for(int i=0; i<dim; i++)
       tmp += vec.m[0][i] * vec.m[0][i];
   return (F)std::sqrt((F)tmp);
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int dim> inline Matrix<Ty,dim,1>
Normalize(const Matrix<Ty,dim,1> & vec)
{
   return vec / Norm(vec);
}

template <class F,class Ty,int dim> inline Matrix<F,dim,1>  // floatでない場合、このoverloadを呼び出さないと警告される
Normalize(const Matrix<Ty,dim,1> & vec)
{
   return vec / Norm<F>(vec);
}

template <class Ty,int dim> inline Matrix<Ty,dim,1>
Normalize(const Matrix<Ty,1,dim> & vec)
{
   return vec / Norm(vec);
}

template <class F,class Ty,int dim> inline Matrix<F,dim,1>  // floatでない場合、このoverloadを呼び出さないと警告される
Normalize(const Matrix<Ty,1,dim> & vec)
{
   return vec / Norm<F>(vec);
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int dim> inline Ty
Dot(const Matrix<Ty,dim,1> & a, const Matrix<Ty,dim,1> & b)
{
   Ty ret = 0;
   for(int i=0; i<dim; i++)
       ret += a.m[i][0] * b.m[i][0];
   return ret;
}

template <class Ty,int dim> inline Ty
Dot(const Matrix<Ty,1,dim> & a, const Matrix<Ty,1,dim> & b)
{
   Ty ret = 0;
   for(int i=0; i<dim; i++)
       ret += a.m[0][i] * b.m[0][i];
   return ret;
}

template <class Ty,int dim> inline Ty
Dot(const Matrix<Ty,dim,1> & a, const Matrix<Ty,1,dim> & b)
{
   Ty ret = 0;
   for(int i=0; i<dim; i++)
       ret += a.m[i][0] * b.m[0][i];
   return ret;
}

template <class Ty,int dim> inline Ty
Dot(const Matrix<Ty,1,dim> & a, const Matrix<Ty,dim,1> & b)
{
   Ty ret = 0;
   for(int i=0; i<dim; i++)
       ret += a.m[0][i] * b.m[0][i];
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,int dim> inline bool
operator < (const Matrix<Ty,dim,1> & a, const Matrix<Ty,dim,1> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[i][0] < b.m[i][0];
   return ret;
}

template <class Ty,int dim> inline bool
operator < (const Matrix<Ty,1,dim> & a, const Matrix<Ty,1,dim> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[0][i] < b.m[0][i];
   return ret;
}

template <class Ty,int dim> inline bool
operator <= (const Matrix<Ty,dim,1> & a, const Matrix<Ty,dim,1> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[i][0] <= b.m[i][0];
   return ret;
}

template <class Ty,int dim> inline bool
operator <= (const Matrix<Ty,1,dim> & a, const Matrix<Ty,1,dim> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[0][i] <= b.m[0][i];
   return ret;
}

template <class Ty,int dim> inline bool
operator > (const Matrix<Ty,dim,1> & a, const Matrix<Ty,dim,1> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[i][0] > b.m[i][0];
   return ret;
}

template <class Ty,int dim> inline bool
operator > (const Matrix<Ty,1,dim> & a, const Matrix<Ty,1,dim> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[0][i] > b.m[0][i];
   return ret;
}

template <class Ty,int dim> inline bool
operator >= (const Matrix<Ty,dim,1> & a, const Matrix<Ty,dim,1> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[i][0] >= b.m[i][0];
   return ret;
}

template <class Ty,int dim> inline bool
operator >= (const Matrix<Ty,1,dim> & a, const Matrix<Ty,1,dim> & b)
{
   bool ret = false;
   for(int i=0; i<dim; i++)
       ret &&= a.m[0][i] >= b.m[0][i];
   return ret;
}

//////////////////////////////////////////////////////////////////////////

template <class Ty,class R,int rows,int cols,class F> inline Matrix<R,rows,cols>
Each(const Matrix<Ty,rows,cols> & mat, F func)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((R *) ret.m) [i] = func( ((Ty *) mat.m) [i] );
   return ret;
}

template <class Ty,int rows,int cols,class F> inline Matrix<Ty,rows,cols>
Each(const Matrix<Ty,rows,cols> & mat, F func)
{
   Matrix<Ty,rows,cols> ret;
   for(int i=0; i<rows*cols; i++)
       ((Ty *) ret.m) [i] = func( ((Ty *) mat.m) [i] );
   return ret;
}

//////////////////////////////////////////////////////////////////////////

typedef Matrix<double,2,1>              double2x1, double2;
typedef Matrix<double,3,1>              double3x1, double3;
typedef Matrix<double,4,1>              double4x1, double4;

typedef Matrix<double,1,2>              double1x2;
typedef Matrix<double,1,3>              double1x3;
typedef Matrix<double,1,4>              double1x4;

typedef Matrix<double,2,2>              double2x2;
typedef Matrix<double,3,3>              double3x3;
typedef Matrix<double,4,4>              double4x4;

typedef Matrix<float,2,1>               float2x1, float2;
typedef Matrix<float,3,1>               float3x1, float3;
typedef Matrix<float,4,1>               float4x1, float4;

typedef Matrix<float,1,2>               float1x2;
typedef Matrix<float,1,3>               float1x3;
typedef Matrix<float,1,4>               float1x4;

typedef Matrix<float,2,2>               float2x2;
typedef Matrix<float,3,3>               float3x3;
typedef Matrix<float,4,4>               float4x4;

//////////////////////////////////////////////////////////////////////////

typedef Matrix<long,2,1>                long2x1, long2;
typedef Matrix<long,3,1>                long3x1, long3;
typedef Matrix<long,4,1>                long4x1, long4;

typedef Matrix<long,1,2>                long1x2;
typedef Matrix<long,1,3>                long1x3;
typedef Matrix<long,1,4>                long1x4;

typedef Matrix<long,2,2>                long2x2;
typedef Matrix<long,3,3>                long3x3;
typedef Matrix<long,4,4>                long4x4;

typedef Matrix<int,2,1>                 int2x1, int2;
typedef Matrix<int,3,1>                 int3x1, int3;
typedef Matrix<int,4,1>                 int4x1, int4;

typedef Matrix<int,1,2>                 int1x2;
typedef Matrix<int,1,3>                 int1x3;
typedef Matrix<int,1,4>                 int1x4;

typedef Matrix<int,2,2>                 int2x2;
typedef Matrix<int,3,3>                 int3x3;
typedef Matrix<int,4,4>                 int4x4;

typedef Matrix<short,2,1>               short2x1, short2;
typedef Matrix<short,3,1>               short3x1, short3;
typedef Matrix<short,4,1>               short4x1, short4;

typedef Matrix<short,1,2>               short1x2;
typedef Matrix<short,1,3>               short1x3;
typedef Matrix<short,1,4>               short1x4;

typedef Matrix<short,2,2>               short2x2;
typedef Matrix<short,3,3>               short3x3;
typedef Matrix<short,4,4>               short4x4;

//////////////////////////////////////////////////////////////////////////

typedef Matrix<byte,2,1>                byte2x1, byte2;
typedef Matrix<byte,3,1>                byte3x1, byte3;
typedef Matrix<byte,4,1>                byte4x1, byte4;

typedef Matrix<byte,1,2>                byte1x2;
typedef Matrix<byte,1,3>                byte1x3;
typedef Matrix<byte,1,4>                byte1x4;

typedef Matrix<byte,2,2>                byte2x2;
typedef Matrix<byte,3,3>                byte3x3;
typedef Matrix<byte,4,4>                byte4x4;

typedef Matrix<word,2,1>                word2x1, word2;
typedef Matrix<word,3,1>                word3x1, word3;
typedef Matrix<word,4,1>                word4x1, word4;

typedef Matrix<word,1,2>                word1x2;
typedef Matrix<word,1,3>                word1x3;
typedef Matrix<word,1,4>                word1x4;

typedef Matrix<word,2,2>                word2x2;
typedef Matrix<word,3,3>                word3x3;
typedef Matrix<word,4,4>                word4x4;

typedef Matrix<dword,2,1>               dword2x1, dword2;
typedef Matrix<dword,3,1>               dword3x1, dword3;
typedef Matrix<dword,4,1>               dword4x1, dword4;

typedef Matrix<dword,1,2>               dword1x2;
typedef Matrix<dword,1,3>               dword1x3;
typedef Matrix<dword,1,4>               dword1x4;

typedef Matrix<dword,2,2>               dword2x2;
typedef Matrix<dword,3,3>               dword3x3;
typedef Matrix<dword,4,4>               dword4x4;

//////////////////////////////////////////////////////////////////////////

typedef Matrix<char,2,1>                char2x1, char2;
typedef Matrix<char,3,1>                char3x1, char3;
typedef Matrix<char,4,1>                char4x1, char4;

typedef Matrix<char,1,2>                char1x2;
typedef Matrix<char,1,3>                char1x3;
typedef Matrix<char,1,4>                char1x4;

typedef Matrix<char,2,2>                char2x2;
typedef Matrix<char,3,3>                char3x3;
typedef Matrix<char,4,4>                char4x4;

typedef Matrix<wchar,2,1>               wchar2x1, wchar2;
typedef Matrix<wchar,3,1>               wchar3x1, wchar3;
typedef Matrix<wchar,4,1>               wchar4x1, wchar4;

typedef Matrix<wchar,1,2>               wchar1x2;
typedef Matrix<wchar,1,3>               wchar1x3;
typedef Matrix<wchar,1,4>               wchar1x4;

typedef Matrix<wchar,2,2>               wchar2x2;
typedef Matrix<wchar,3,3>               wchar3x3;
typedef Matrix<wchar,4,4>               wchar4x4;

//////////////////////////////////////////////////////////////////////////

typedef Matrix<bool,2,1>                bool2x1, bool2;
typedef Matrix<bool,3,1>                bool3x1, bool3;
typedef Matrix<bool,4,1>                bool4x1, bool4;

typedef Matrix<bool,1,2>                bool1x2;
typedef Matrix<bool,1,3>                bool1x3;
typedef Matrix<bool,1,4>                bool1x4;

typedef Matrix<bool,2,2>                bool2x2;
typedef Matrix<bool,3,3>                bool3x3;
typedef Matrix<bool,4,4>                bool4x4;


}//end of ns takaku


#End If