#ifndef MATRIX_H

#define MATRIX_H

#include <cmath>

#include "vector.h"

class Matrix {

public:

	union {

		Real m[3][3];

		struct {

			Real m00, m01, m02;

			Real m10, m11, m12;

			Real m20, m21, m22;

		};

	};

	Matrix(

		Real m00, Real m01, Real m02,

		Real m10, Real m11, Real m12,

		Real m20, Real m21, Real m22

	):

		m00(m00), m01(m01), m02(m02),

		m10(m10), m11(m11), m12(m12),

		m20(m20), m21(m21), m22(m22)

	{ }

	Matrix(): Matrix(

		1, 0, 0,

		0, 1, 0,

		0, 0, 1 ) { }

	Matrix(

		const Vector &axis_x,

		const Vector &axis_y,

		const Vector &offset = Vector()

	): Matrix(

		axis_x.x, axis_x.y, 0,

		axis_y.x, axis_y.y, 0,

		offset.x, offset.y, 1 ) { }

	const Vector& axis(int index) const { return *(const Vector*)(m[index]); }

	const Vector& axis_x() const { return axis(0); }

	const Vector& axis_y() const { return axis(1); }

	const Vector& offset() const { return axis(2); }

	Vector& axis(int index) { return *(Vector*)(m[index]); }

	Vector& axis_x() { return axis(0); }

	Vector& axis_y() { return axis(1); }

	Vector& offset() { return axis(2); }

	Matrix& set_identity()

		{ return *this = Matrix(); }

	bool is_identity() const

		{ return *this == Matrix(); }

	Matrix& set_scale(const Real &x, const Real &y) {

		return *this = Matrix( x, 0, 0,

							   0, y, 0,

							   0, 0, 1 );

	}

	Matrix& set_scale(const Real &s)

		{ return set_scale(s, s); }

	Matrix& set_scale(const Vector &s)

		{ return set_scale(s.x, s.y); }

	Matrix& set_rotate(const Real &angle) {

		Real s = sin(angle), c = cos(angle);

		return *this = Matrix( c, s, 0,

							  -s, c, 0,

							   0, 0, 1 );

	}

	Matrix& set_translate(const Real &x, const Real &y) {

		return *this = Matrix( 1, 0, 0,

							   0, 1, 0,

							   x, y, 1 );

	}

	Matrix& set_translate(const Vector &t)

		{ return set_translate(t.x, t.y); }

	void transform(

		Real &out_x, Real &out_y, Real &out_z,

		const Real &x, const Real &y, const Real &z ) const

	{

		out_x = x*m00 + y*m10 + z*m20;

		out_y = x*m01 + y*m11 + z*m21;

		out_z = x*m02 + y*m12 + z*m22;

	}

	void transform(Real &out_x, Real &out_y, const Real &x, const Real &y, bool translate = true) const

		{ Real z; transform(out_x, out_y, z, x, y, translate ? 1.0 : 0.0); }

	Vector transform(const Vector &v, bool translate = true)const

		{ Vector vv; transform(vv.x, vv.y, v.x, v.y, translate); return vv; }

	bool operator==(const Matrix &other) const {

		return equal(m00, other.m00) && equal(m01, other.m01) && equal(m02, other.m02)

			&& equal(m10, other.m10) && equal(m11, other.m11) && equal(m12, other.m12)

			&& equal(m20, other.m20) && equal(m21, other.m21) && equal(m22, other.m22);

	}

	bool operator!=(const Matrix &other) const

		{ return !(*this == other); }

	Matrix& operator*=(const Matrix &other)

		{ return *this = *this * other; }

	Matrix operator*(const Matrix &other) const {

		Matrix m;

		transform(m.m00, m.m01, m.m02, other.m00, other.m01, other.m02);

		transform(m.m10, m.m11, m.m12, other.m10, other.m11, other.m12);

		transform(m.m20, m.m21, m.m22, other.m20, other.m21, other.m22);

		return m;

	}

	Real det() const {

		return m00*(m11*m22 - m12*m21)

			 - m01*(m10*m22 - m12*m20)

			 + m02*(m10*m21 - m11*m20);

	}

	bool is_invertible() const

		{ return !equal(det(), 0); }

	Matrix get_inverted() const;

	Matrix& invert()

		{ return *this = get_inverted(); }

};

#endif