/**
	This file is a part of the rexy/r0nk/atlas project
	Copyright (C) 2020 rexy712

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.	See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.	If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef REXY_DETAIL_MATRIX_HPP
#define REXY_DETAIL_MATRIX_HPP

#include <cstdlib> //size_t
#include <utility> //integer_sequence

namespace math::detail{

	template<size_t SW, size_t W = SW, size_t H = SW-1, size_t... Args>
	struct gen_id_tup {
		using tup = typename gen_id_tup<SW, W-1, H, Args..., 0>::tup;
	};
	template<size_t SW, size_t H, size_t... Args>
	struct gen_id_tup<SW,SW,H,Args...> {
		using tup = typename gen_id_tup<SW, SW-1, H, Args..., 1>::tup;
	};
	template<size_t SW, size_t H, size_t... Args>
	struct gen_id_tup<SW,0,H,Args...> {
		using tup = typename gen_id_tup<SW, SW, H-1, Args..., 0>::tup;
	};
	template<size_t SW, size_t... Args>
	struct gen_id_tup<SW,SW,0,Args...> {
		using tup = std::integer_sequence<size_t,Args...,1>;
	};


	template<size_t N, size_t... Args>
	struct gen_zero_tup {
		using tup = typename gen_zero_tup<N-1,Args...,0>::tup;
	};
	template<size_t... Args>
	struct gen_zero_tup<0,Args...> {
		using tup = std::integer_sequence<size_t,Args...>;
	};

	template<size_t W>
	struct id_initialization_matrix {
		using tuple = typename gen_id_tup<W>::tup;
	};


	template<size_t W, size_t H>
	struct default_initialization_matrix {
		using tuple = typename gen_zero_tup<W>::tup;
	};
	template<size_t W>
	struct default_initialization_matrix<W,W> {
		using tuple = typename id_initialization_matrix<W>::tuple;
	};

	template<typename T, size_t R>
	class mat_ref_obj
	{
	public:
		using size_type = size_t;

	protected:
		T* m_data = nullptr;
	public:
		constexpr mat_ref_obj(T* d, size_type i);
		constexpr T& operator[](size_type i);
		constexpr const T& operator[](size_type i)const;
	};
	template<typename T, size_t R>
	struct determinate_helper {
		static constexpr T perform(const matrix<T,R,R>& m){
			T sum = 0;
			T op = 1;
			for(size_t i = 0;i < R;++i){
				T item = op * m[0][i];
				matrix<T,R-1,R-1> mul(no_initialize);
				for(size_t j = 1, mj = 0;j < R;++j){
					for(size_t k = 0, mk = 0;k < R;++k){
						if(k == i)
							continue;
						mul[mj][mk] = m[j][k];
						++mk;
					}
					++mj;
				}
				sum += item * determinate_helper<T,R-1>::perform(mul);
				op = -op;
			}
			return sum;
		}
	};
	template<typename T>
	struct determinate_helper<T,3> {
		static constexpr T perform(const matrix<T,3,3>& m){
			return (m.get(0) * ((m.get(4) * m.get(8)) - (m.get(5) * m.get(7))) -
			        m.get(1) * ((m.get(3) * m.get(8)) - (m.get(5) * m.get(6))) +
			        m.get(2) * ((m.get(3) * m.get(7)) - (m.get(4) * m.get(6))));
		}
	};
	template<typename T>
	struct determinate_helper<T,2> {
		static constexpr T perform(const matrix<T,2,2>& m){
			return m.get(0) * m.get(3) - m.get(1) * m.get(2);
		}
	};

	template<typename T, size_t R>
	struct inverse_helper {
		//TODO generalized inverse
	};
	template<typename T>
	struct inverse_helper<T,2> {
		static constexpr matrix<T,2,2> perform(const matrix<T,2,2>& m){
			T det = m.determinate();
			if(!det)
				return matrix<T,2,2>(zero_initialize);
			return matrix<T,2,2>(m.get(3) / det, -(m.get(1)) / det, -(m.get(2)) / det, m.get(0) / det);
		}
	};
	template<typename T>
	struct inverse_helper<T,3> {
		static constexpr matrix<T,3,3> perform(const matrix<T,3,3>& m){
			T det = m.determinate();
			if(!det)
				return matrix<T,3,3>(zero_initialize);
			return matrix<T,3,3>(((m.get(4) * m.get(8)) - (m.get(5) * m.get(7))) / det,
			                    -((m.get(1) * m.get(8)) - (m.get(2) * m.get(7))) / det,
			                     ((m.get(1) * m.get(5)) - (m.get(2) * m.get(4))) / det,
			                    -((m.get(3) * m.get(8)) - (m.get(5) * m.get(6))) / det,
			                     ((m.get(0) * m.get(8)) - (m.get(2) * m.get(6))) / det,
			                    -((m.get(0) * m.get(5)) - (m.get(2) * m.get(3))) / det,
			                     ((m.get(3) * m.get(7)) - (m.get(4) * m.get(6))) / det,
			                    -((m.get(0) * m.get(7)) - (m.get(1) * m.get(6))) / det,
			                     ((m.get(0) * m.get(4)) - (m.get(1) * m.get(3))) / det);
		}
	};
	template<typename T>
	struct inverse_helper<T,4> {
		static constexpr matrix<T,4,4> perform(const matrix<T,4,4>& m){
			T det = m.determinate();
			if(!det)
				return matrix<T,4,4>(zero_initialize);
			//Math is power
			return matrix<T,4,4>((m.get(5) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
			                      m.get(6) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  +
			                      m.get(7) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))) / det,

			                    -(m.get(1) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
			                      m.get(2) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  +
			                      m.get(3) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))) / det,

			                     (m.get(1) * ((m.get(6)  * m.get(15)) - (m.get(7)  * m.get(14)))  -
			                      m.get(2) * ((m.get(5)  * m.get(15)) - (m.get(7)  * m.get(13)))  +
			                      m.get(3) * ((m.get(5)  * m.get(14)) - (m.get(6)  * m.get(13)))) / det,

			                    -(m.get(1) * ((m.get(6)  * m.get(11)) - (m.get(7)  * m.get(10)))  -
			                      m.get(2) * ((m.get(5)  * m.get(11)) - (m.get(7)  * m.get(9)))   +
			                      m.get(3) * ((m.get(5)  * m.get(10)) - (m.get(6)  * m.get(9))))  / det,

			                    -(m.get(4) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
			                      m.get(6) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
			                      m.get(7) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))) / det,

			                     (m.get(0) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
			                      m.get(2) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
			                      m.get(3) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))) / det,

			                    -(m.get(0) * ((m.get(6)  * m.get(15)) - (m.get(7)  * m.get(14)))  -
			                      m.get(2) * ((m.get(4)  * m.get(15)) - (m.get(7)  * m.get(12)))  +
			                      m.get(3) * ((m.get(4)  * m.get(14)) - (m.get(6)  * m.get(12)))) / det,

			                     (m.get(0) * ((m.get(6)  * m.get(11)) - (m.get(7)  * m.get(10)))  -
			                      m.get(2) * ((m.get(4)  * m.get(11)) - (m.get(7)  * m.get(8)))   +
			                      m.get(3) * ((m.get(4)  * m.get(10)) - (m.get(6)  * m.get(8))))  / det,

			                     (m.get(4) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  -
			                      m.get(5) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
			                      m.get(7) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,

			                    -(m.get(0) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  -
			                      m.get(1) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
			                      m.get(3) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,

			                     (m.get(0) * ((m.get(5)  * m.get(15)) - (m.get(7)  * m.get(13)))  -
			                      m.get(1) * ((m.get(4)  * m.get(15)) - (m.get(7)  * m.get(12)))  +
			                      m.get(3) * ((m.get(4)  * m.get(13)) - (m.get(5)  * m.get(12)))) / det,

			                    -(m.get(0) * ((m.get(5)  * m.get(11)) - (m.get(7)  * m.get(9)))   -
			                      m.get(1) * ((m.get(4)  * m.get(11)) - (m.get(7)  * m.get(8)))   +
			                      m.get(3) * ((m.get(4)  * m.get(9))  - (m.get(5)  * m.get(8))))  / det,

			                    -(m.get(4) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))  -
			                      m.get(5) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))  +
			                      m.get(6) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,

			                     (m.get(0) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))  -
			                      m.get(1) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))  +
			                      m.get(2) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,

			                    -(m.get(0) * ((m.get(5)  * m.get(14)) - (m.get(6)  * m.get(13)))  -
			                      m.get(1) * ((m.get(4)  * m.get(14)) - (m.get(6)  * m.get(12)))  +
			                      m.get(2) * ((m.get(4)  * m.get(13)) - (m.get(5)  * m.get(12)))) / det,

			                     (m.get(0) * ((m.get(5)  * m.get(10)) - (m.get(6)  * m.get(9)))   -
			                      m.get(1) * ((m.get(4)  * m.get(10)) - (m.get(6)  * m.get(8)))   +
			                      m.get(2) * ((m.get(4)  * m.get(9))  - (m.get(5)  * m.get(8))))  / det);
		}
	};

}

#include "matrix.tpp"

#endif