/************************************************************************
 * Matrix<int D1, int D2, class T_FLOAT>
 *
 * 要素の数が静的に決まる２次元行列クラス。
 * 同じ要素数の２次元配列との演算もサポート。
 * ただし、掛け算とInverse()は縦横同じ要素数でなければ爆弾。
 * 基本的に例外なし、実行時チェックなしの速度重視クラス。
 ************************************************************************/

#ifndef ATS_MATRIXXXX_H
#define ATS_MATRIXXXX_H

// ats_MatrixXxX.h

#include <iostream>
#include <cmath>
#include <stdexcept>

namespace ats
{
	template <class T_FLOAT>
	class Matrix3D; // 循環参照解決

	// D = Dimension
	template <int D1, int D2, class T_FLOAT>
	class Matrix
	{
	friend Matrix3D<T_FLOAT>;
	private:
		T_FLOAT m_matrix[D1][D2];

	public:
		// Constructor / Destructor -----------------------------------------------------
		Matrix() 
		{
			Initialize();
		}

		// boolをとるコンストラクタについて ----------------------------
		// 普通に使えば、初期化は保証されているが、初期化しないという
		// 選択肢を提供する。これは派生を重ねる時に生じるコンストラク
		// タのオーバーヘッドを抑えるためである。
		// -------------------------------------------------------------
		Matrix(const bool init) 
		{
			// 初期化するかしないかを選択する。
			if(init){
				Initialize();
			}
		}

		Matrix(const T_FLOAT a[D1][D2]) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = a[i][j];
			}}
		}

		// Copy Constructor
		Matrix(const Matrix& a) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = a(i,j);
			}}
		}

		virtual ~Matrix() 
		{}

		virtual void Initialize() 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = 0;
			}}
		}

		virtual Matrix<D1,D2,T_FLOAT>& Inverse() 
		{
			Matrix<D1,D2,T_FLOAT> tmp(false);
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					tmp(i,j) = m_matrix[j][i];
			}}
			return (*this = tmp);
		}

		// Getter / Setter ---------------------------------------------------------
		inline T_FLOAT& operator()(int a, int b)
		{
			return m_matrix[a][b];
		}

		inline T_FLOAT operator()(int a, int b) const
		{
			return m_matrix[a][b];
		}

		// 代入 ---------------------------------------------------------------------
		Matrix<D1,D2,T_FLOAT>& operator=(const Matrix<D1,D2,T_FLOAT>& a) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = a(i,j);
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator+=(const Matrix<D1,D2,T_FLOAT>& a) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] += a(i,j);
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator-=(const Matrix<D1,D2,T_FLOAT>& a) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] -= a(i,j);
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator*=(const Matrix<D1,D2,T_FLOAT>& a) 
		{
			Matrix<D1,D2,T_FLOAT> tmp;
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					for(int m=0; m < D1; m++){
							tmp(i,j) += m_matrix[i][m] * a(m,j);
			}}}
			return (*this = tmp);
		}

		// 演算 ----------------------------------------------------------------------
		Matrix<D1,D2,T_FLOAT> operator+(const Matrix<D1,D2,T_FLOAT>& a) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp(false);
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					tmp[i][j] = m_matrix[i][j] + a(i,j);
			}}
			return tmp;
		}

		Matrix<D1,D2,T_FLOAT> operator-(const Matrix<D1,D2,T_FLOAT>& a) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp(false);
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					tmp(i,j) = m_matrix[i][j] - a(i,j);
			}}
			return tmp;
		}

		Matrix<D1,D2,T_FLOAT> operator*(const Matrix<D1,D2,T_FLOAT>& a) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp;
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					for(int m=0; m < D1; m++){
							tmp(i,j) += m_matrix[i][m] * a(m,j);
			}}}
			return tmp;
		}

		// 代入(値) ----------------------------------------------------------------
		Matrix<D1,D2,T_FLOAT>& operator=(const T_FLOAT a) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = a
			}}
			return *this;
		}

		// 代入(配列) -------------------------------------------------------------
		Matrix<D1,D2,T_FLOAT>& operator=(const T_FLOAT a[D1][D2]) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] = a[i][j];
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator+=(const T_FLOAT a[D1][D2]) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] += a[i][j];
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator-=(const T_FLOAT a[D1][D2]) 
		{
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					m_matrix[i][j] -= a[i][j];
			}}
			return *this;
		}

		Matrix<D1,D2,T_FLOAT>& operator*=(const T_FLOAT a[D1][D2]) 
		{
			Matrix<D1,D2,T_FLOAT> tmp;
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					for(int m=0; m < D1; m++){
							tmp(i,j) += m_matrix[i][m] * a[m][j];
			}}}
			return (*this = tmp);
		}

		// 演算(配列) --------------------------------------------------------------
		Matrix<D1,D2,T_FLOAT> operator+(const T_FLOAT a[D1][D2]) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp(false);
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					tmp(i,j) = m_matrix[i][j] + a[i][j];
			}}
			return tmp;
		}

		Matrix<D1,D2,T_FLOAT> operator-(const T_FLOAT a[D1][D2]) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp(false);
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					tmp(i,j) = m_matrix[i][j] - a[i][j];
			}}
			return tmp;
		}

		Matrix<D1,D2,T_FLOAT> operator*(const T_FLOAT a[D1][D2]) const 
		{
			Matrix<D1,D2,T_FLOAT> tmp;
			for(int i=0; i < D1; i++){
				for(int j=0; j < D2; j++){
					for(int m=0; m < D1; m++){
							tmp(i,j) += m_matrix[i][m] * a[m][j];
			}}}
			return tmp;
		}
	};

	template <int D1, int D2, class T_FLOAT>
	std::ostream& operator<<(std::ostream& out, const Matrix<D1, D2, T_FLOAT>& a)
	{
		for(int i=0; i < D1; i++){
			for(int j=0; j < D2; j++){
				out << a(i, j) << ' ';
			}
			out << '\n';
		}
		return out;
	}
}

#endif // #ifndef ATS_MATRIXXXX_H