#ifndef _PIUMA_LIB_MATH_H_
#define _PIUMA_LIB_MATH_H_

#include "types.h"
#include <math.h>

// Vector types
union v2
{
	struct { f32 x, y; };
	struct { f32 r, g; };
	struct { f32 u, v; };
	struct { f32 E[2]; };
};

union v2s
{
	struct { s32 x, y; };
	struct { s32 r, g; };
	struct { s32 u, v; };
	struct { s32 E[2]; };
};

union v2u
{
	struct { u32 x, y; };
	struct { u32 r, g; };
	struct { u32 u, v; };
	struct { u32 E[2]; };
};


union v3
{
	struct { f32 x, y, z;  };
	struct { f32 r, g, b;  };
	struct { f32 u, v, w_; };
	struct { f32 E[3];     };


	struct
	{
		v2 xy;
		f32 ignored0;
	};
	struct
	{
		f32 ignored1;
		v2 yz;
	};
	struct
	{
		v2 uv;
		f32 ignored2;
	};
	struct
	{
		f32 ignored3;
		v2 vw_;
	};
	struct
	{
		f32 pitch;
		f32 yaw;
		f32 roll;
	};
};

union v3s
{
	struct { s32 x, y, z;  };
	struct { s32 r, g, b;  };
	struct { s32 u, v, w_; };
	struct { s32 E[3];     };


	struct
	{
		v2s xy;
		s32 ignored0;
	};
	struct
	{
		s32 ignored1;
		v2s yz;
	};
	struct
	{
		v2s uv;
		s32 ignored2;
	};
	struct
	{
		s32 ignored3;
		v2s vw_;
	};
	struct
	{
		s32 pitch;
		s32 yaw;
		s32 roll;
	};
};

union v3u
{
	struct { u32 x, y, z;  };
	struct { u32 r, g, b;  };
	struct { u32 u, v, w_; };
	struct { u32 E[3];     };


	struct
	{
		v2u xy;
		u32 ignored0;
	};
	struct
	{
		u32 ignored1;
		v2u yz;
	};
	struct
	{
		v2u uv;
		u32 ignored2;
	};
	struct
	{
		u32 ignored3;
		v2u vw_;
	};
	struct
	{
		u32 pitch;
		u32 yaw;
		u32 roll;
	};
};


union v4
{
	struct { f32 x, y, z, w; };
	struct { f32 r, g, b, a; };
	struct { f32 E[4];       };

	struct
	{
		v2 xy;
		f32 ignored0;
		f32 ignored1;
	};
	struct
	{
		f32 ignored2;
		v2 yz;
		f32 ignored3;
	};
	struct
	{
		f32 ignored4;
		f32 ignored5;
		v2 zw;
	};
	struct
	{
		v2 uv;
		f32 ignored6;
		f32 ignored7;
	};
	struct
	{
		v3 xyz;
		f32 ignored8;
	};
	struct
	{
		v3 rgb;
		f32 ignored9;
	};
};

union v4s
{
	struct { s32 x, y, z, w; };
	struct { s32 r, g, b, a; };
	struct { s32 E[4];       };

	struct
	{
		v2s xy;
		s32 ignored0;
		s32 ignored1;
	};
	struct
	{
		s32 ignored2;
		v2s yz;
		s32 ignored3;
	};
	struct
	{
		s32 ignored4;
		s32 ignored5;
		v2s zw;
	};
	struct
	{
		v2s uv;
		s32 ignored6;
		s32 ignored7;
	};
	struct
	{
		v3s xyz;
		s32 ignored8;
	};
	struct
	{
		v3s rgb;
		s32 ignored9;
	};
};

union v4u
{
	struct { u32 x, y, z, w; };
	struct { u32 r, g, b, a; };
	struct { u32 E[4];       };

	struct
	{
		v2u xy;
		u32 ignored0;
		u32 ignored1;
	};
	struct
	{
		u32 ignored2;
		v2u yz;
		u32 ignored3;
	};
	struct
	{
		u32 ignored4;
		u32 ignored5;
		v2u zw;
	};
	struct
	{
		v2u uv;
		u32 ignored6;
		u32 ignored7;
	};
	struct
	{
		v3u xyz;
		u32 ignored8;
	};
	struct
	{
		v3u rgb;
		u32 ignored9;
	};
};

// Matrix types
union m3x3
{
	struct { v3  row[3]; };
	struct { f32 E[3][3]; };
};
typedef m3x3 m3;

union m4x4
{
	struct { v4  row[4]; };
	struct { f32 E[4][4]; };
};
typedef m4x4 m4;


static const m3 m3_zero =
{
	0, 0, 0,
	0, 0, 0,
	0, 0, 0
};
static const m3 m3_identity =
{
	1, 0, 0,
	0, 1, 0,
	0, 0, 1
};
static const m4 m4_zero =
{
	0, 0, 0, 0,
	0, 0, 0, 0,
	0, 0, 0, 0,
	0, 0, 0, 0
};
static const m4 m4_identity =
{
	1, 0, 0, 0,
	0, 1, 0, 0,
	0, 0, 1, 0,
	0, 0, 0, 1
};


// Angles
static const f32 PI  = 3.141592653589793238462643383;
static const f32 TAU = 6.283185307179586476925286766;

#define radians(x) (x*TAU/360)
#define degrees(x) (x*360/TAU)


// Utilities
#define square(x) ((x)*(x))
#define minimum(a, b) ((a) < (b) ? (a) : (b))
#define maximum(a, b) ((a) > (b) ? (a) : (b))


inline f32 clamp(f32 min, f32 max, f32 value)
{
	if(value < min)
		return min;
	if(value > max)
		return max;
	return value;
}

inline f32 lerp(f32 start, f32 end, f32 u)
{
	//return (1 - u) * start + u * end;
	return start + (end - start) * u;
}

inline v3 lerp(v3 start, v3 end, f32 u)
{
	f32 x = lerp(start.x, end.x, u);
	f32 y = lerp(start.y, end.y, u);
	f32 z = lerp(start.z, end.z, u);
	return v3{x, y, z};
}

inline v4 lerp(v4 start, v4 end, f32 u)
{
	f32 x = lerp(start.x, end.x, u);
	f32 y = lerp(start.y, end.y, u);
	f32 z = lerp(start.z, end.z, u);
	f32 w = lerp(start.w, end.w, u);
	return v4{x, y, z, w};
}

inline f32 map(f32 new_start, f32 new_end, f32 old_start, f32 old_end, f32 u)
{
	f32 old_range = old_end - old_start;
	f32 new_range = new_end - new_start;
	return (u - old_start) / old_range * new_range + new_start;
}


// Conversions
inline v2 V2(v2s a)
{
	return {(f32)a.x, (f32)a.y};
}

inline v2 V2(v2u a)
{
	return {(f32)a.x, (f32)a.y};
}

inline v2s V2S(v2 a)
{
	return {(s32)a.x, (s32)a.y};
}

inline v2s V2S(v2u a)
{
	return {(s32)a.x, (s32)a.y};
}

inline v2u V2U(v2 a)
{
	return {(u32)a.x, (u32)a.y};
}

inline v2u V2U(v2s a)
{
	return {(u32)a.x, (u32)a.y};
}

inline v3 V3(v2 a, f32 z)
{
	return {a.x, a.y, z};
}

inline v3 V3(v3s a)
{
	return {(f32)a.x, (f32)a.y, (f32)a.z};
}

inline v3 V3(v3u a)
{
	return {(f32)a.x, (f32)a.y, (f32)a.z};
}

inline v3 V3(f32 e[3])
{
	return v3{e[0], e[1], e[2]};
}

inline v3 V3(v4 v)
{
	return v3{v.x / v.w, v.y / v.w, v.z / v.w};
}

inline v3s V3S(v2s a, s32 z)
{
	return {a.x, a.y, z};
}

inline v3s V3S(v3 a)
{
	return {(s32)a.x, (s32)a.y, (s32)a.z};
}

inline v3s V3S(v3u a)
{
	return {(s32)a.x, (s32)a.y, (s32)a.z};
}

inline v3u V3U(v2u a, u32 z)
{
	return {a.x, a.y, z};
}

inline v3u V3U(v3 a)
{
	return {(u32)a.x, (u32)a.y, (u32)a.z};
}

inline v3u V3U(v3s a)
{
	return {(u32)a.x, (u32)a.y, (u32)a.z};
}

inline v4 V4(v2 a, f32 z, f32 w)
{
	v4 result;

	result.xy = a;
	result.z  = z;
	result.w  = w;

	return result;
}

inline v4 V4(v3 a, f32 w)
{
	v4 result;

	result.xyz = a;
	result.w   = w;

	return result;
}

inline v4 V4(f32 e[4])
{
	return v4{e[0], e[1], e[2], e[3]};
}

inline m3 M3(m4 m)
{
	return m3{ m.row[0].xyz, m.row[1].xyz, m.row[2].xyz };
}

inline m4 M4(m3 m)
{
	return m4{ V4(m.row[0], 0), V4(m.row[1], 0), V4(m.row[2], 0), v4{0,0,0,1} };
}


// Operators
inline v2 operator+(v2 a)
{
	return a;
}
inline v2 operator-(v2 a)
{
	v2 result = { -a.x, -a.y };
	return result;
}

inline bool operator==(v2 a, v2 b)
{
	return a.x == b.x && a.y == b.y;
}

inline bool operator!=(v2 a, v2 b)
{
	return !(a == b);
}

inline v2 operator+(v2 a, v2 b)
{
	v2 result;

	result.x = a.x + b.x;
	result.y = a.y + b.y;

	return result;
}
inline v2 operator-(v2 a, v2 b)
{
	v2 result;

	result.x = a.x - b.x;
	result.y = a.y - b.y;

	return result;
}
inline v2 operator*(v2 a, v2 b)
{
	v2 result;

	result.x = a.x * b.x;
	result.y = a.y * b.y;

	return result;
}
inline v2 operator/(v2 a, v2 b)
{
	v2 result;

	result.x = a.x / b.x;
	result.y = a.y / b.y;

	return result;
}
inline v2 operator*(v2 a, f32 b)
{
	v2 result;

	result.x = a.x * b;
	result.y = a.y * b;

	return result;
}
inline v2 operator/(v2 a, f32 b)
{
	v2 result;

	result.x = a.x / b;
	result.y = a.y / b;

	return result;
}
inline v2 operator*(f32 a, v2 b)
{
	v2 result;

	result.x = a * b.x;
	result.y = a * b.y;

	return result;
}

inline v2 & operator+=(v2 &a, v2 b)
{
	a = a + b;
	return a;
}
inline v2 & operator-=(v2 &a, v2 b)
{
	a = a - b;
	return a;
}
inline v2 & operator*=(v2 &a, v2 b)
{
	a = a * b;
	return a;
}
inline v2 & operator/=(v2 &a, v2 b)
{
	a = a / b;
	return a;
}
inline v2 & operator*=(v2 &a, f32 b)
{
	a = a * b;
	return a;
}
inline v2 & operator/=(v2 &a, f32 b)
{
	a = a / b;
	return a;
}

inline v2s operator+(v2s a)
{
	return a;
}
inline v2s operator-(v2s a)
{
	v2s result = { -a.x, -a.y };
	return result;
}

inline bool operator==(v2s a, v2s b)
{
	return a.x == b.x && a.y == b.y;
}

inline bool operator!=(v2s a, v2s b)
{
	return !(a == b);
}

inline v2s operator+(v2s a, v2s b)
{
	v2s result;

	result.x = a.x + b.x;
	result.y = a.y + b.y;

	return result;
}
inline v2s operator-(v2s a, v2s b)
{
	v2s result;

	result.x = a.x - b.x;
	result.y = a.y - b.y;

	return result;
}
inline v2s operator*(v2s a, v2s b)
{
	v2s result;

	result.x = a.x * b.x;
	result.y = a.y * b.y;

	return result;
}
inline v2s operator/(v2s a, v2s b)
{
	v2s result;

	result.x = a.x / b.x;
	result.y = a.y / b.y;

	return result;
}
inline v2s operator*(v2s a, f32 b)
{
	v2s result;

	result.x = a.x * b;
	result.y = a.y * b;

	return result;
}
inline v2s operator/(v2s a, f32 b)
{
	v2s result;

	result.x = a.x / b;
	result.y = a.y / b;

	return result;
}
inline v2s operator*(f32 a, v2s b)
{
	v2s result;

	result.x = a * b.x;
	result.y = a * b.y;

	return result;
}

inline v2s & operator+=(v2s &a, v2s b)
{
	a = a + b;
	return a;
}
inline v2s & operator-=(v2s &a, v2s b)
{
	a = a - b;
	return a;
}
inline v2s & operator*=(v2s &a, v2s b)
{
	a = a * b;
	return a;
}
inline v2s & operator/=(v2s &a, v2s b)
{
	a = a / b;
	return a;
}
inline v2s & operator*=(v2s &a, f32 b)
{
	a = a * b;
	return a;
}
inline v2s & operator/=(v2s &a, f32 b)
{
	a = a / b;
	return a;
}

inline v2u operator+(v2u a)
{
	return a;
}
inline v2u operator-(v2u a)
{
	v2u result = { -a.x, -a.y };
	return result;
}

inline bool operator==(v2u a, v2u b)
{
	return a.x == b.x && a.y == b.y;
}

inline bool operator!=(v2u a, v2u b)
{
	return !(a == b);
}

inline v2u operator+(v2u a, v2u b)
{
	v2u result;

	result.x = a.x + b.x;
	result.y = a.y + b.y;

	return result;
}
inline v2u operator-(v2u a, v2u b)
{
	v2u result;

	result.x = a.x - b.x;
	result.y = a.y - b.y;

	return result;
}
inline v2u operator*(v2u a, v2u b)
{
	v2u result;

	result.x = a.x * b.x;
	result.y = a.y * b.y;

	return result;
}
inline v2u operator/(v2u a, v2u b)
{
	v2u result;

	result.x = a.x / b.x;
	result.y = a.y / b.y;

	return result;
}
inline v2u operator*(v2u a, f32 b)
{
	v2u result;

	result.x = a.x * b;
	result.y = a.y * b;

	return result;
}
inline v2u operator/(v2u a, f32 b)
{
	v2u result;

	result.x = a.x / b;
	result.y = a.y / b;

	return result;
}
inline v2u operator*(f32 a, v2u b)
{
	v2u result;

	result.x = a * b.x;
	result.y = a * b.y;

	return result;
}

inline v2u & operator+=(v2u &a, v2u b)
{
	a = a + b;
	return a;
}
inline v2u & operator-=(v2u &a, v2u b)
{
	a = a - b;
	return a;
}
inline v2u & operator*=(v2u &a, v2u b)
{
	a = a * b;
	return a;
}
inline v2u & operator/=(v2u &a, v2u b)
{
	a = a / b;
	return a;
}
inline v2u & operator*=(v2u &a, f32 b)
{
	a = a * b;
	return a;
}
inline v2u & operator/=(v2u &a, f32 b)
{
	a = a / b;
	return a;
}



inline v3 operator+(v3 a)
{
	return a;
}
inline v3 operator-(v3 a)
{
	v3 result = { -a.x, -a.y, -a.z };
	return result;
}

inline bool operator==(v3 a, v3 b)
{
	return a.x == b.x && a.y == b.y && a.z == b.z;
}

inline bool operator!=(v3 a, v3 b)
{
	return !(a == b);
}

inline v3 operator+(v3 a, v3 b)
{
	v3 result;

	result.x = a.x + b.x;
	result.y = a.y + b.y;
	result.z = a.z + b.z;

	return result;
}
inline v3 operator-(v3 a, v3 b)
{
	v3 result;

	result.x = a.x - b.x;
	result.y = a.y - b.y;
	result.z = a.z - b.z;

	return result;
}
inline v3 operator*(v3 a, v3 b)
{
	v3 result;

	result.x = a.x * b.x;
	result.y = a.y * b.y;
	result.z = a.z * b.z;

	return result;
}
inline v3 operator/(v3 a, v3 b)
{
	v3 result;

	result.x = a.x / b.x;
	result.y = a.y / b.y;
	result.z = a.z / b.z;

	return result;
}
inline v3 operator*(v3 a, f32 b)
{
	v3 result;

	result.x = a.x * b;
	result.y = a.y * b;
	result.z = a.z * b;

	return result;
}
inline v3 operator/(v3 a, f32 b)
{
	v3 result;

	result.x = a.x / b;
	result.y = a.y / b;
	result.z = a.z / b;

	return result;
}
inline v3 operator*(f32 a, v3 b)
{
	v3 result;

	result.x = a * b.x;
	result.y = a * b.y;
	result.z = a * b.z;

	return result;
}

inline v3 & operator+=(v3 &a, v3 b)
{
	a = a + b;
	return a;
}
inline v3 & operator-=(v3 &a, v3 b)
{
	a = a - b;
	return a;
}
inline v3 & operator*=(v3 &a, v3 b)
{
	a = a * b;
	return a;
}
inline v3 & operator/=(v3 &a, v3 b)
{
	a = a / b;
	return a;
}
inline v3 & operator*=(v3 &a, f32 b)
{
	a = a * b;
	return a;
}
inline v3 & operator/=(v3 &a, f32 b)
{
	a = a / b;
	return a;
}



inline v4 operator+(v4 a)
{
	return a;
}
inline v4 operator-(v4 a)
{
	v4 result = { -a.x, -a.y, -a.z, -a.w };
	return result;
}

inline bool operator==(v4 a, v4 b)
{
	return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w;
}

inline bool operator!=(v4 a, v4 b)
{
	return !(a == b);
}

inline v4 operator+(v4 a, v4 b)
{
	v4 result;

	result.x = a.x + b.x;
	result.y = a.y + b.y;
	result.z = a.z + b.z;
	result.w = a.w + b.w;

	return result;
}
inline v4 operator-(v4 a, v4 b)
{
	v4 result;

	result.x = a.x - b.x;
	result.y = a.y - b.y;
	result.z = a.z - b.z;
	result.w = a.w - b.w;

	return result;
}
inline v4 operator*(v4 a, v4 b)
{
	v4 result;

	result.x = a.x * b.x;
	result.y = a.y * b.y;
	result.z = a.z * b.z;
	result.w = a.w * b.w;

	return result;
}
inline v4 operator/(v4 a, v4 b)
{
	v4 result;

	result.x = a.x / b.x;
	result.y = a.y / b.y;
	result.z = a.z / b.z;
	result.w = a.w / b.w;

	return result;
}
inline v4 operator*(v4 a, f32 b)
{
	v4 result;

	result.x = a.x * b;
	result.y = a.y * b;
	result.z = a.z * b;
	result.w = a.w * b;

	return result;
}
inline v4 operator/(v4 a, f32 b)
{
	v4 result;

	result.x = a.x / b;
	result.y = a.y / b;
	result.z = a.z / b;
	result.w = a.w / b;

	return result;
}
inline v4 operator*(f32 a, v4 b)
{
	v4 result;

	result.x = a * b.x;
	result.y = a * b.y;
	result.z = a * b.z;
	result.w = a * b.w;

	return result;
}

inline v4 & operator+=(v4 &a, v4 b)
{
	a = a + b;
	return a;
}
inline v4 & operator-=(v4 &a, v4 b)
{
	a = a - b;
	return a;
}
inline v4 & operator*=(v4 &a, v4 b)
{
	a = a * b;
	return a;
}
inline v4 & operator/=(v4 &a, v4 b)
{
	a = a / b;
	return a;
}
inline v4 & operator*=(v4 &a, f32 b)
{
	a = a * b;
	return a;
}
inline v4 & operator/=(v4 &a, f32 b)
{
	a = a / b;
	return a;
}


// Vector functions
inline f32 dot(v2 a, v2 b)
{
	return a.x*b.x + a.y*b.y;
}

inline f32 dot(v3 a, v3 b)
{
	return a.x*b.x + a.y*b.y + a.z*b.z;
}

inline f32 dot(v4 a, v4 b)
{
	return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
}

inline v3 cross(v3 a, v3 b)
{
	v3 result;

	result.x = a.y*b.z - a.z*b.y;
	result.y = a.z*b.x - a.x*b.z;
	result.z = a.x*b.y - a.y*b.x;

	return result;
}


inline f32 length(v2 a)
{
	return sqrt(dot(a,a));
}

inline f32 length(v3 a)
{
	return sqrt(dot(a,a));
}

inline f32 length(v4 a)
{
	return sqrt(dot(a,a));
}

inline f32 length2(v2 a)
{
	return dot(a,a);
}

inline f32 length2(v3 a)
{
	return dot(a,a);
}

inline f32 length2(v4 a)
{
	return dot(a,a);
}


inline f32 distance(v2 a, v2 b)
{
	return length(a - b);
}

inline f32 distance(v3 a, v3 b)
{
	return length(a - b);
}

inline f32 distance(v4 a, v4 b)
{
	return length(a - b);
}


inline f32 distance2(v2 a, v2 b)
{
	return length2(a - b);
}

inline f32 distance2(v3 a, v3 b)
{
	return length2(a - b);
}

inline f32 distance2(v4 a, v4 b)
{
	return length2(a - b);
}


inline v2 normalize(v2 a)
{
	return a / length(a);
}

inline v3 normalize(v3 a)
{
	return a / length(a);
}

inline v4 normalize(v4 a)
{
	return a / length(a);
}


// Matrix operators
inline v3 extract_column(m3 m, u8 column_index)
{
	v3 result;

	result.x = m.E[0][column_index];
	result.y = m.E[1][column_index];
	result.z = m.E[2][column_index];

	return result;
}

inline v4 extract_column(m4 m, u8 column_index)
{
	v4 result;

	result.x = m.E[0][column_index];
	result.y = m.E[1][column_index];
	result.z = m.E[2][column_index];
	result.w = m.E[3][column_index];

	return result;
}

inline v3 operator*(m3 m, v3 v)
{
	v3 result;

	result.x = dot(m.row[0], v);
	result.y = dot(m.row[1], v);
	result.z = dot(m.row[2], v);

	return result;
}

inline m3 operator*(m3 a, m3 b)
{
	m3 result;

	v3 c0 = extract_column(b, 0);
	v3 c1 = extract_column(b, 1);
	v3 c2 = extract_column(b, 2);

	result.E[0][0] = dot(a.row[0], c0);
	result.E[0][1] = dot(a.row[0], c1);
	result.E[0][2] = dot(a.row[0], c2);

	result.E[1][0] = dot(a.row[1], c0);
	result.E[1][1] = dot(a.row[1], c1);
	result.E[1][2] = dot(a.row[1], c2);

	result.E[2][0] = dot(a.row[2], c0);
	result.E[2][1] = dot(a.row[2], c1);
	result.E[2][2] = dot(a.row[2], c2);

	return result;
}

inline v4 operator*(m4 m, v4 v)
{
	v4 result;

	result.x = dot(m.row[0], v);
	result.y = dot(m.row[1], v);
	result.z = dot(m.row[2], v);
	result.w = dot(m.row[3], v);

	return result;
}

inline m4 operator*(m4 a, m4 b)
{
	m4 result;

	v4 c0 = extract_column(b, 0);
	v4 c1 = extract_column(b, 1);
	v4 c2 = extract_column(b, 2);
	v4 c3 = extract_column(b, 3);

	result.E[0][0] = dot(a.row[0], c0);
	result.E[0][1] = dot(a.row[0], c1);
	result.E[0][2] = dot(a.row[0], c2);
	result.E[0][3] = dot(a.row[0], c3);

	result.E[1][0] = dot(a.row[1], c0);
	result.E[1][1] = dot(a.row[1], c1);
	result.E[1][2] = dot(a.row[1], c2);
	result.E[1][3] = dot(a.row[1], c3);

	result.E[2][0] = dot(a.row[2], c0);
	result.E[2][1] = dot(a.row[2], c1);
	result.E[2][2] = dot(a.row[2], c2);
	result.E[2][3] = dot(a.row[2], c3);

	result.E[3][0] = dot(a.row[3], c0);
	result.E[3][1] = dot(a.row[3], c1);
	result.E[3][2] = dot(a.row[3], c2);
	result.E[3][3] = dot(a.row[3], c3);

	return result;
}


// Matrix functions
inline m3 m3_make_diagonal_matrix(f32 value)
{
	m3 result = m3_zero;
	result.E[0][0] = value;
	result.E[1][1] = value;
	result.E[2][2] = value;
	return result;
}

inline m4 m4_make_diagonal_matrix(f32 value)
{
	m4 result = m4_zero;
	result.E[0][0] = value;
	result.E[1][1] = value;
	result.E[2][2] = value;
	result.E[3][3] = value;
	return result;
}

inline f32 determinant_of_minor(m3 m, int skip_r, int skip_c)
{
	f32 minor[2][2];
	for(int r = 0, source_r = 0; r < 2; r++, source_r++)
	{
		if(source_r == skip_r)
			source_r++;
		for(int c = 0, source_c = 0; c < 2; c++, source_c++)
		{
			if(source_c == skip_c)
				source_c++;
			minor[r][c] = m.E[source_r][source_c];
		}
	}

	return minor[0][0] * minor[1][1] - minor[0][1] * minor[1][0];
}

inline f32 determinant(m3 m)
{
	return
		+ m.E[0][0] * m.E[1][1] * m.E[2][2]
		+ m.E[0][1] * m.E[1][2] * m.E[2][0]
		+ m.E[0][2] * m.E[1][0] * m.E[2][1]
		- m.E[0][2] * m.E[1][1] * m.E[2][0]
		- m.E[0][1] * m.E[1][0] * m.E[2][2]
		- m.E[0][0] * m.E[1][2] * m.E[2][1];
}

inline m3 transpose(m3 m)
{
	m3 result;

	result.E[0][0] = m.E[0][0];
	result.E[0][1] = m.E[1][0];
	result.E[0][2] = m.E[2][0];

	result.E[1][0] = m.E[0][1];
	result.E[1][1] = m.E[1][1];
	result.E[1][2] = m.E[2][1];

	result.E[2][0] = m.E[0][2];
	result.E[2][1] = m.E[1][2];
	result.E[2][2] = m.E[2][2];

	result.E[3][0] = m.E[0][3];
	result.E[3][1] = m.E[1][3];
	result.E[3][2] = m.E[2][3];

	return result;
}

inline m3 inverse(m3 m)
{
	m3 result;

	f32 one_over_det = 1.0 / determinant(m);
	f32 sign = 1;
	for(int r = 0; r < 3; r++)
	{
		for(int c = 0; c < 3; c++)
		{
			result.E[c][r] = sign * determinant_of_minor(m, r, c) * one_over_det;
			sign = -sign;
		}
	}

	return result;
}

inline f32 determinant_of_minor(m4 m, int skip_r, int skip_c)
{
	m3 minor;
	for(int r = 0, source_r = 0; r < 3; r++, source_r++)
	{
		if(source_r == skip_r)
			source_r++;
		for(int c = 0, source_c = 0; c < 3; c++, source_c++)
		{
			if(source_c == skip_c)
				source_c++;
			minor.E[r][c] = m.E[source_r][source_c];
		}
	}
	return determinant(minor);
}

inline f32 determinant(m4 m)
{
	f32 result = 0;
	f32 sign = 1;
	int r = 0;
	// Laplace theorem
	for(int c = 0; c < 4; c++)
	{
		result += m.E[r][c] * sign * determinant_of_minor(m, r, c);
		sign = -sign;
	}
	return result;
}

inline m4 transpose(m4 m)
{
	m4 result;

	result.E[0][0] = m.E[0][0];
	result.E[0][1] = m.E[1][0];
	result.E[0][2] = m.E[2][0];
	result.E[0][3] = m.E[3][0];

	result.E[1][0] = m.E[0][1];
	result.E[1][1] = m.E[1][1];
	result.E[1][2] = m.E[2][1];
	result.E[1][3] = m.E[3][1];

	result.E[2][0] = m.E[0][2];
	result.E[2][1] = m.E[1][2];
	result.E[2][2] = m.E[2][2];
	result.E[2][3] = m.E[3][2];

	result.E[3][0] = m.E[0][3];
	result.E[3][1] = m.E[1][3];
	result.E[3][2] = m.E[2][3];
	result.E[3][3] = m.E[3][3];

	return result;
}

inline m4 inverse(m4 m)
{
	m4 result;

	f32 one_over_det = 1.0 / determinant(m);
	f32 sign_a = 1;
	for(int r = 0; r < 4; r++)
	{
		f32 sign_b = sign_a;
		for(int c = 0; c < 4; c++)
		{
			result.E[c][r] = sign_b * determinant_of_minor(m, r, c) * one_over_det;
			sign_b = -sign_b;
		}
		sign_a = -sign_a;
	}

	return result;
}

#endif