#ifndef LIN_ALG_H
#define LIN_ALG_H

#include <algorithm>
#include <iostream>
#include "orient_point.h"


template <class T>
class matrix_t
{
public:
	matrix_t() {}
	matrix_t(const point_t<T>& x, const point_t<T>& y, const point_t<T>& z)
		{a=x; b=y; c=z;}
	matrix_t(const matrix_t<T>& M) {a=M.a; b=M.b; c=M.c;}

	matrix_t<T>& operator= (const matrix_t<T>& M)
		{a=M.a; b=M.b; c=M.c; return *this;}

	void MakeTranspose()
		{
			std::swap(a.y, b.x);
			std::swap(c.y, b.z);
			std::swap(a.z, c.x);
		}

	matrix_t<T>& operator *= (const matrix_t<T>& M)
		{return *this = *this * M;}

	matrix_t<T>& operator += (const matrix_t<T>& M)
		{a+=M.a; b+=M.b; c+=M.c; return *this;}
	
	matrix_t<T>& operator -= (const matrix_t<T>& M)
		{a-=M.a; b-=M.b; c-=M.c; return *this;}
	
	matrix_t<T> operator + (const matrix_t<T>& M) const
		{return matrix_t<T>(*this) += M;}
	
	matrix_t<T> operator - (const matrix_t<T>& M) const
		{return matrix_t<T>(*this) -= M;}
	
	matrix_t<T> operator * (const matrix_t<T>& M) const
		{
			matrix_t<T> tmp(M), R;
			tmp.MakeTranspose();
			point_t<T> A(a * tmp.a, a * tmp.b, a * tmp.c);
			point_t<T> B(b * tmp.a, b * tmp.b, b * tmp.c);
			point_t<T> C(c * tmp.a, c * tmp.b, c * tmp.c);
			return matrix_t<T>(A, B, C);
		}

	point_t<T> operator * (const point_t<T>& V) const
		{return point_t<T>(a * V, b * V, c * V);}

	double Det() const {return a.cross(b) * c;}

	static double Det(point_t<T>& x, point_t<T>& y, point_t<T>& z)
		{return x.cross(y) * z;}

	matrix_t<T> Inverse() const
		{
			point_t<T> A(b.cross(c));
			point_t<T> B(c.cross(a));
			point_t<T> C(a.cross(b));
			double D = 1.0 / (C * c);
			matrix_t<T> M(A * D, B * D, C * D);
			M.MakeTranspose();
			return M;
		}

	point_t<T> a, b, c;
};

template <class T>
std::ostream& operator << (std::ostream& o, const point_t<T>& V)
{
	o << "(" << V.x << ", " << V.y << ", " << V.z << ")";
	return o;
}

template <class T>
std::ostream& operator << (std::ostream& o, const matrix_t<T>& V)
{
	o << "[" << V.a << ", " << V.b << ", " << V.c << "]";
	return o;
}

template <class T>
void update_min(point_t<T>&mn, const point_t<T>& v)
{
	mn.x = std::min(mn.x, v.x);
	mn.y = std::min(mn.y, v.y);
	mn.z = std::min(mn.z, v.z);
}

template <class T>
void update_max(point_t<T>&mx, const point_t<T>& v)
{
	mx.x = std::max(mx.x, v.x);
	mx.y = std::max(mx.y, v.y);
	mx.z = std::max(mx.z, v.z);
}

template <class T>
void update_min_max(point_t<T>&mn, point_t<T>&mx,
						   const point_t<T>& v)
{
	update_min(mn, v);
	update_max(mx, v);
}


static point_t<double> unit_x(1,0,0), unit_y(0,1,0), unit_z(0,0,1);
static point_t<double> vec_zero(0,0,0), vec_one(1,1,1);
static matrix_t<double> mat_zero(vec_zero, vec_zero, vec_zero);
static matrix_t<double> mat_ident(unit_x, unit_y, unit_z);

typedef point_t<float> vector3f;
typedef point_t<double> vector3d;
typedef matrix_t<float> matrix3f;
typedef matrix_t<double> matrix3d;

#endif