/**
	This file is a part of our_dick
	Copyright (C) 2020 rexy712

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU Affero 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 Affero General Public License for more details.

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

#ifndef REXY_MAT_TPP
#define REXY_MAT_TPP

#include <cstdlib> //size_t
#include <cmath> //sin, cos
#include <type_traits> //decay_t, declval
#include "detail/matrix.hpp"
#include "quat.hpp"

namespace math{

	template<typename T, size_t R, size_t C>
  template<size_t... Ss>
  constexpr matrix_base<T,R,C>::matrix_base(std::integer_sequence<size_type,Ss...>):
		m_data{Ss...}{}

	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::matrix_base():
		matrix_base(typename detail::default_initialization_matrix<Columns,Rows>::tuple{}){}

	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::matrix_base(detail::zero_initialize_t):
		m_data{}{}
	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::matrix_base(detail::no_initialize_t){}

	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::matrix_base(value_type v){
		for(size_type i = 0; i < Columns*Rows; ++i)
			m_data[i] = v;
	}
	template<typename T, size_t R, size_t C>
	template<typename... Args>
	constexpr matrix_base<T,R,C>::matrix_base(Args&&... args):
		m_data{std::forward<Args>(args)...}{}

	template<typename T, size_t R, size_t C>
	template<typename U>
	constexpr matrix_base<T,R,C>::matrix_base(const matrix_base<U,Columns,Rows>& m){
		using mat = decltype(m);
		for(typename mat::size_type i = 0; i < mat::Columns*mat::Rows; ++i)
			m_data[i] = m.get(i);
	}

	template<typename T, size_t R, size_t C>
	template<typename U>
	constexpr matrix_base<T,R,C>& matrix_base<T,R,C>::operator=(const matrix_base<U,Columns,Rows>& m){
		using mat = decltype(m);
		for(typename mat::size_type i = 0; i < mat::Columns*mat::Rows; ++i)
			m_data[i] = m.get(i);
		return *this;
	}


	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::operator[](size_type x){
		return detail::mat_ref_obj<value_type,Rows>{m_data, x};
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::operator[](size_type x)const{
		return detail::mat_ref_obj<const value_type,Rows>{m_data, x};
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::get(size_type x, size_type y) -> reference{
		return m_data[x+(y*Rows)];
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::get(size_type x, size_type y)const -> const_reference{
		return m_data[x+(y*Rows)];
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::get(size_type i) -> reference{
		return m_data[i];
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::get(size_type i)const -> const_reference{
		return m_data[i];
	}

	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::columns()const -> size_type{
		return Columns;
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::rows()const -> size_type{
		return Rows;
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::size()const -> size_type{
		return Columns*Rows;
	}

	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::raw() -> pointer{
		return m_data;
	}
	template<typename T, size_t R, size_t C>
	constexpr auto matrix_base<T,R,C>::raw()const -> const_pointer{
		return m_data;
	}
	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::operator pointer(){
		return m_data;
	}
	template<typename T, size_t R, size_t C>
	constexpr matrix_base<T,R,C>::operator const_pointer()const{
		return m_data;
	}

	template<typename T, size_t R, size_t C>
	template<typename U>
	constexpr matrix<T,R,C>& matrix<T,R,C>::operator=(const matrix<U,R,C>& m){
		base::operator=(m);
		return *this;
	}

	template<typename T, size_t R>
	constexpr matrix<T,R,R>::matrix(detail::id_initialize_t):
		base(){}
	template<typename T, size_t R>
	template<typename U>
	constexpr matrix<T,R,R>& matrix<T,R,R>::operator=(const matrix<U,R,R>& m){
		base::operator=(m);
		return *this;
	}

	template<typename T, size_t R>
	constexpr auto matrix<T,R,R>::determinate()const -> value_type{
		return math::determinate(*this);
	}
	template<typename T, size_t R>
	constexpr auto matrix<T,R,R>::trace()const -> value_type{
		value_type sum = 0;
		for(size_type i = 0; i < R; ++i){
			sum += get(i, i);
		}
	}
	template<typename T, size_t R>
	constexpr matrix<T,R,R> matrix<T,R,R>::transpose()const{
		matrix m(no_initialize);
		for(size_type i = 0; i < R; ++i){
			for(size_type j = 0; j < R; ++j){
				m.get(j, i) = get(i, j);
			}
		}
		return m;
	}
	template<typename T, size_t R>
	constexpr matrix<T,R,R> matrix<T,R,R>::inverse()const{
		return math::inverse(*this);
	}


	template<typename T, size_t R>
	constexpr T determinate(const matrix<T,R,R>& m){
		return detail::determinate_helper<T,R>::perform(m);
	}
	template<typename T, size_t R>
	constexpr matrix<T,R,R> inverse(const matrix<T,R,R>& m){
		return detail::inverse_helper<T,R>::perform(m);
	}

	template<typename T>
	matrix<T,2,2> rotation2d_pure(T angle){
		return rotation2d_pure(std::sin(angle), std::cos(angle));
	}
	template<typename T>
	constexpr matrix<T,2,2> rotation2d_pure(T sin, T cos){
		return matrix<T,2,2>(cos, sin, -sin, cos);
	}
	template<typename T>
	constexpr matrix<T,2,2> scale2d(T x, T y){
		return matrix<T,2,2>(x, 0, 0, y);
	}

	template<typename T>
	matrix<T,3,3> rotation2d(T angle){
		return rotation2d(std::sin(angle), std::cos(angle));
	}
	template<typename T>
	constexpr matrix<T,3,3> rotation2d(T sin, T cos){
		return matrix<T,3,3>(cos, -sin, 0,
		                     sin,  cos, 0,
		                     0,    0,   1);
	}
	template<typename T>
	matrix<T,3,3> rotation2d(T x, T y, T z){
		quaternion<T> q(x, y, z);
		return q.to_mat3();
	}

	template<typename T>
	matrix<T,4,4> fov_projection(T fov, T asp, T near, T far){
		T r = near * std::tan(fov / T{2.0});
		return matrix<T,4,4>((near / r) / asp, 0,          0,                               0,
		                     0,                (near / r), 0,                               0,
		                     0,                0,          (far + near) / (near - far),    -1,
		                     0,                0,          (2 * near * far) / (near - far), 0);
	}
	template<typename T>
	matrix<T,4,4> fov_asymetric_projection(T fovl, T fovr, T fovb, T fovt, T asp, T n, T f){
		T l = n * std::tan(fovl);
		T r = n * std::tan(fovr);
		T b = n * std::tan(fovb);
		T t = n * std::tan(fovt);

		return matrix<T,4,4>(((2 * n) / (r - l)) * asp, 0,                 0,                     0,
		                     0,                         (2 * n) / (t - b), 0,                     0,
		                     (r + l) / (r - l),         (t + b) / (t - b), (f + n) / (n - f),    -1,
		                     0,                         0,                 (2 * n * f) / (n - f), 0);
	}
	template<typename T>
	matrix<T,4,4> ortho_projection(T w, T h, T n, T f){
		return matrix<T,4,4>(2 / w, 0,     0,                 0,
		                     0,     2 / h, 0,                 0,
		                     0,     0,     2 / (n - f),       0,
		                     0,     0,     (n + f) / (n - f), 1);
	}
	template<typename T>
	matrix<T,4,4> ortho_asymetric_projection(T l, T r, T b, T t, T n, T f){
		return matrix<T,4,4>(2 / (r - l),       0,                 0,                 0,
		                     0,                 2 / (t - b),       0,                 0,
		                     0,                 0,                 2 / (n - f),       0,
		                     (r + l) / (l - r), (t + b) / (b - t), (n + f) / (n - f), 1);
	}
	template<typename T>
	constexpr matrix<T,4,4> rotation3d(T angle_x, T angle_y, T angle_z){
		quaternion<T> q(angle_x, angle_y, angle_z);
		return q.to_mat4();
	}
	template<typename T>
	constexpr matrix<T,4,4> translation3d(T x, T y, T z){
		return matrix<T,4,4>(1, 0, 0, 0,
		                     0, 1, 0, 0,
		                     0, 0, 1, 0,
		                     x, y, z, 1);
	}
	template<typename T>
	constexpr matrix<T,4,4> scale3d(T x, T y, T z){
		return matrix<T,4,4>(x, 0, 0, 0,
		                     0, y, 0, 0,
		                     0, 0, z, 0,
		                     0, 0, 0, 1);
	}


	template<typename T, typename U, size_t R, size_t C>
	constexpr bool operator==(const matrix_base<T,R,C>& left, const matrix_base<U,R,C> right){
		for(size_t i = 0; i < left.size(); ++i){
			if(left.get(i) != right.get(i))
				return false;
		}
		return true;
	}
	template<typename T, typename U, size_t R, size_t C>
	constexpr bool operator!=(const matrix_base<T,R,C>& left, const matrix_base<U,R,C> right){
		return !(left == right);
	}

	template<typename T, typename U, size_t R1, size_t C1, size_t C2>
	constexpr auto operator*(const matrix<T,R1,C1>& left, const matrix<U,C1,C2>& right){
		using res_t = decltype(std::declval<T>() * std::declval<U>());
		matrix<res_t,R1,C2> res(zero_initialize);
		size_t index = 0;
		for(size_t i = 0; i < left.rows(); ++i){
			for(size_t j = 0; j < right.columns(); ++j){
				for(size_t k = 0; k < right.rows(); ++k){
					res.get(index) += left[i][k] * right[k][j];
				}
				++index;
			}
		}
	}
	template<typename T, typename U, size_t C, size_t R, std::enable_if_t<!is_matrix<U>::value,int>>
	constexpr auto operator*(const matrix<T,R,C>& left, U&& right){
		using res_t = decltype(std::declval<T>() * std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < left.size(); ++i){
			res.get(i) = left.get(i) * std::forward<U>(right);
		}
		return res;
	}
	template<typename T, typename U, size_t C, size_t R, std::enable_if_t<!is_matrix<U>::value,int>>
	constexpr auto operator*(U&& left, const matrix<T,R,C>& right){
		using res_t = decltype(std::declval<T>() * std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < right.size(); ++i){
			res.get(i) = std::forward<U>(left) * right.get(i);
		}
		return res;
	}
	template<typename T, typename U, size_t C, size_t R, std::enable_if_t<!is_matrix<U>::value,int>>
	constexpr auto operator/(const matrix<T,R,C>& left, U&& right){
		using res_t = decltype(std::declval<T>() / std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < left.size(); ++i){
			res.get(i) = left.get(i) / std::forward<U>(right);
		}
		return res;
	}
	template<typename T, typename U, size_t C, size_t R>
	constexpr auto operator+(const matrix<T,R,C>& left, const matrix<U,R,C>& right){
		using res_t = decltype(std::declval<T>() + std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < left.size(); ++i){
			res.get(i) = left.get(i) + right.get(i);
		}
		return res;
	}
	template<typename T, typename U, size_t C, size_t R>
	constexpr auto operator-(const matrix<T,R,C>& left, const matrix<U,R,C>& right){
		using res_t = decltype(std::declval<T>() - std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < left.size(); ++i){
			res.get(i) = left.get(i) - right.get(i);
		}
		return res;
	}
	template<typename T, typename U, size_t C, size_t R>
	constexpr auto operator-(const matrix<T,R,C>& left){
		using res_t = decltype(-std::declval<U>());
		matrix<res_t,R,C> res(no_initialize);
		for(size_t i = 0; i < left.size(); ++i){
			res.get(i) = -left.get(i);
		}
		return res;
	}


	template<typename T, typename U, size_t R>
	constexpr decltype(auto) operator*=(matrix<T,R,R>& left, const matrix<U,R,R>& right){
		//have to evaluate entire expression first since matrix multiplication depends on reusing many elements
		//cannot be expression templatized, TODO
		return (left = (left * right));
	}
	template<typename T, typename U, size_t C, size_t R, std::enable_if_t<!is_matrix<U>::value,int>>
	constexpr decltype(auto) operator*=(matrix<T,R,C>& left, U&& right){
		for(size_t i = 0; i < left.size(); ++i){
			left.get(i) = left.get(i) * std::forward<U>(right);
		}
		return left;
	}
	template<typename T, typename U, size_t C, size_t R, std::enable_if_t<!is_matrix<U>::value,int>>
	constexpr decltype(auto) operator/=(matrix<T,R,C>& left, U&& right){
		for(size_t i = 0; i < left.size(); ++i){
			left.get(i) = left.get(i) / std::forward<U>(right);
		}
		return left;
	}
	template<typename T, typename U, size_t C, size_t R>
	constexpr decltype(auto) operator+=(matrix<T,R,C>& left, const matrix<U,R,C>& right){
		for(size_t i = 0; i < left.size(); ++i){
			left.get(i) = left.get(i) + right.get(i);
		}
		return left;
	}
	template<typename T, typename U, size_t C, size_t R>
	constexpr decltype(auto) operator-=(matrix<T,R,C>& left, const matrix<U,R,C>& right){
		for(size_t i = 0; i < left.size(); ++i){
			left.get(i) = left.get(i) - right.get(i);
		}
		return left;
	}

}

#endif