//=============================================================================
//
// Copyright 1998 by Captivation Digital Laboratories, Inc.
// All Rights Reserved.
//
//=============================================================================
//
// "KRAD" -- Kamui Rendering/Animation Demonstration
//
// Note: These classes are adapted from utilities published
// in Graphics Gems IV, but they hardly resemble them anymore
//--------------------
// C++ Vector and Matrix Algebra routines
// Author: Jean-Francois DOUE
// Version 3.1 --- October 1993
//
//=============================================================================
//-----------------------------------------------------------------------------
// $Revision: 1 $
// $NoKeywords: $
//=============================================================================

#ifndef _KALGEBRA_H_
#define _KALGEBRA_H_

#include "krad.h"

//=============================================================================
//=============================================================================	FORWARD CLASS DECL
//=============================================================================

class kvector;
class kmatrix;


//=============================================================================
//=============================================================================	DEFINITIONS, et al
//=============================================================================

//--------------------------------------------------------- Useful constants

// -- D. Knuth, Art of Computer Programming, volume 1, 3rd edition, page 619
const float PI			= 3.1415926535897932384626433932795028841792f;
const float ONE_DEGREE	= 0.0174532925199432957692369076848861271344f;

//--------------------------------------------------------- Useful indices
enum {VX, VY, VZ, VW};		    // axes

//--------------------------------------------------------- Useful functions
inline float DEG2RAD (float dval)	{	return(dval*ONE_DEGREE);	}


//=============================================================================
//=============================================================================	AGGRESSIVE MATH
//=============================================================================
void InitializeAggressiveMath (void);

extern float reciprocal_table[];
extern float squareroot_table[];

inline float reciprocal (const float& d)	{	return reciprocal_table	[*(((short*)(&d))+1)];	}
inline float squareroot	(const float& d)	{	return squareroot_table	[*(((short*)(&d))+1)];	}


//=============================================================================
//=============================================================================	class kvector
//=============================================================================
// This is an implementation of a 3 dimensional vector
//-----------------------------------------------------------------------------

class kvector
{
	public:				// everyone's a 'friend'
		float	n[3];

	public:
		//----------------------------- CONSTRUCTORS
		kvector()
				{
				}

		kvector(const float d)
				{	n[0]=
					n[1]=
					n[2]=d;
				}

		kvector(const kvector& v)	// copy constructor
				{	n[0] = v.n[0];
					n[1] = v.n[1];
					n[2] = v.n[2];
				}

		kvector(const float x, const float y, const float z)
				{	n[0] = x;
					n[1] = y;
					n[2] = z;
				}

		//----------------------------- ACCESS FUNCTIONS
		float&	X (void)			{	return n[0];	}
		float&	Y (void)			{	return n[1];	}
		float&	Z (void)			{	return n[2];	}

		float	GetX(void) const	{	return n[0];	}
		float	GetY(void) const	{	return n[1];	}
		float	GetZ(void) const	{	return n[2];	}

		float&	SetX(const float d)	{	return n[0]=d;	}
		float&	SetY(const float d)	{	return n[1]=d;	}
		float&	SetZ(const float d)	{	return n[2]=d;	}

		void SetXY(const float xval, const float yval)
		{
			n[0] = xval;
			n[1] = yval;
		}

		void SetXZ(const float xval, const float zval)
		{
			n[0] = xval;
			n[2] = zval;
		}

		void SetYZ(const float yval, const float zval)
		{
			n[1] = yval;
			n[2] = zval;
		}

		kvector& SetXYZ(const float xval, const float yval, const float zval)
		{
			n[0] = xval;
			n[1] = yval;
			n[2] = zval;
			return *this;

		}

		//----------------------------- OPERATORS
		kvector& operator = (const kvector& v)
				{
					n[0] = v.n[0];
					n[1] = v.n[1];
					n[2] = v.n[2];
					return *this;
				}

		kvector& operator += (const kvector& v)
				{
					n[0] += v.n[0];
					n[1] += v.n[1];
					n[2] += v.n[2];
					return *this;
				}

		kvector& operator -= (const kvector& v)
				{
					n[0] -= v.n[0];
					n[1] -= v.n[1];
					n[2] -= v.n[2];
					return *this;
				}

		kvector& operator *= (const float d)
				{
					n[0] *= d;
					n[1] *= d;
					n[2] *= d;
					return *this;
				}

		kvector& operator /= (const float d)
				{
					float r = reciprocal(d);
					n[0] *= r;
					n[1] *= r;
					n[2] *= r;
					return *this;
				}

		float&	 operator [] (int i)
				{
				 	return n[i];
				}

		//----------------------------- EXPLICIT ARITHMETIC
		// *this = *this op b
		kvector& negate	(void)
				{
					n[0] = -n[0];
					n[1] = -n[1];
					n[2] = -n[2];
					return *this;
				}

		kvector& add (const kvector& b)
				{
					n[0]+=b.n[0];
					n[1]+=b.n[1];
					n[2]+=b.n[2];
					return *this;
				}

		kvector& sub (const kvector& b)
				{
					n[0]-=b.n[0];
					n[1]-=b.n[1];
					n[2]-=b.n[2];
					return *this;
				}

		kvector& mult (const float d)
				{
					n[0] *= d;
					n[1] *= d;
					n[2] *= d;
					return *this;
				}

		kvector& div (const float d)
				{
					float r = reciprocal(d);
					n[0] *= r;
					n[1] *= r;
					n[2] *= r;
					return *this;
				}

		kvector& cross (const kvector& b)
				{
					float t0	=n[1]*b.n[2] - n[2]*b.n[1];
					float t1	=n[2]*b.n[0] - n[0]*b.n[2];
					n[2]		=n[0]*b.n[1] - n[1]*b.n[0];
					n[1] = t1;
					n[0] = t0;
					return *this;
				}

		//----------------------------- EXPLICIT ARITHMETIC
		// *this = a op b
		kvector& negate	(const kvector& b)
				{
					n[0] = -b.n[0];
					n[1] = -b.n[1];
					n[2] = -b.n[2];
					return *this;
				}

		kvector& add (const kvector& a, const kvector& b)
				{
					n[0] = a.n[0]+b.n[0];
					n[1] = a.n[1]+b.n[1];
					n[2] = a.n[2]+b.n[2];
					return *this;
				}

		kvector& sub (const kvector& a, const kvector& b)
				{
					n[0] = a.n[0]-b.n[0];
					n[1] = a.n[1]-b.n[1];
					n[2] = a.n[2]-b.n[2];
					return *this;
				}

		kvector& mult (const kvector& a, const float d)
				{
					n[0] = d*a.n[0];
					n[1] = d*a.n[1];
					n[2] = d*a.n[2];
					return *this;
				}

		kvector& mult (const float d, const kvector& a)
				{
					n[0] = d*a.n[0];
					n[1] = d*a.n[1];
					n[2] = d*a.n[2];
					return *this;
				}

		kvector& div (const kvector& a, const float d)
				{
					float r=reciprocal(d);
					n[0] = r*a.n[0];
					n[1] = r*a.n[1];
					n[2] = r*a.n[2];
					return *this;
				}

		kvector& cross (const kvector& a, const kvector& b)
				{
					n[0] = a.n[1]*b.n[2] - a.n[2]*b.n[1];
					n[1] = a.n[2]*b.n[0] - a.n[0]*b.n[2];
					n[2] = a.n[0]*b.n[1] - a.n[1]*b.n[0];
					return *this;
				}

		//----------------------------- MISC FUNCTIONS
		// Special functions
		float	 length2() const
				{
					return n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
				}

		float	 length() const
				{
					return (float) squareroot(length2());
				}

		kvector& normalize()
				{
					*this /= length();
					return *this;
				}

		kvector& normalize(const kvector& a)
				{
					return div(a,a.length());
				}
};

//=============================================================================
//=============================================================================	class kmatrix
//=============================================================================
// This is an implementation of a 4x3 homogenous matrix.
//
// The fourth row is assumed to be (0,0,0,1), thus having
// no effect as long as shearing is avoided.
//
// Internally it is built up as a 3x3 matrix, with a 3x1 position vector.
//-----------------------------------------------------------------------------

class kmatrix
{
	public:				// everyone's a 'friend'
		kvector	m[3];	// matrix component
		kvector t;		// translation component

	public:
		//----------------------------- CONSTRUCTORS
		kmatrix()
				{
				}

		kmatrix(const float d)
				{
					m[0]=
					m[1]=
					m[2]=
					t	=kvector(d);
				}

		kmatrix(const kmatrix& s)	// copy constructor
				{
					m[0]=s.m[0];
					m[1]=s.m[1];
					m[2]=s.m[2];
					t	=s.t;
				}

		kmatrix(const kvector& v0, const kvector& v1, const kvector& v2, const kvector& pos)
				{
					m[0]=v0;
					m[1]=v1;
					m[2]=v2;
					t	=pos;
				}

		//----------------------------- ACCESS FUNCTIONS
		kvector& Position(void)
				{
					return t;
				}

		kvector& SetPosition(const kvector& a)
				{
					return t = a;
				}
		kvector& SetPosition (const float xval, const float yval, const float zval)
				{
					return (t.SetXYZ (xval, yval, zval));
				}

		kvector& Row(int i)
				{
					return m[i];
				}

		kvector& SetRow(int i, const kvector& a)
				{
					return m[i] = a;
				}
		kvector& SetRow (int i, const float xval, const float yval, const float zval)
				{
					return (m[i].SetXYZ (xval, yval, zval));
				}

		//----------------------------- OPERATORS
		kmatrix& operator = (const kmatrix& a)
				{
					m[0] = a.m[0];
					m[1] = a.m[1];
					m[2] = a.m[2];
					t	 = a.t;
					return *this;
				}

		kmatrix& operator += (const kmatrix& a)
				{
					m[0] += a.m[0];
					m[1] += a.m[1];
					m[2] += a.m[2];
					t    += a.t;
					return *this;
				}

		kmatrix& operator -= (const kmatrix& a)
				{
					m[0] -= a.m[0];
					m[1] -= a.m[1];
					m[2] -= a.m[2];
					t    -= a.t;
					return *this;
				}

		kmatrix& operator *= (const float d)
				{
					m[0] *= d;
					m[1] *= d;
					m[2] *= d;
					t    *= d;
					return *this;
				}

		kmatrix& operator /= (const float d)
				{
					float r = reciprocal(d);
					m[0] *= r;
					m[1] *= r;
					m[2] *= r;
					t    *= r;
					return *this;
				}

		//----------------------------- EXPLICIT ARITHMETIC
		// *this = *this op b
		kmatrix& negate (void)
				{
					m[0].negate();
					m[1].negate();
					m[2].negate();
					t   .negate();
					return *this;
				}

		kmatrix& add (const kmatrix& b)
				{
					m[0].add(b.m[0]);
					m[1].add(b.m[1]);
					m[2].add(b.m[2]);
					t   .add(b.t   );
					return *this;
				}

		kmatrix& sub (const kmatrix& b)
				{
					m[0].sub(b.m[0]);
					m[1].sub(b.m[1]);
					m[2].sub(b.m[2]);
					t   .sub(b.t   );
					return *this;
				}

		kmatrix& mult (const float d)
				{
					m[0].mult(d);
					m[1].mult(d);
					m[2].mult(d);
					t   .mult(d);
					return *this;
				}

		kmatrix& div (const float d)
				{
					float r = reciprocal(d);
					m[0].mult(r);
					m[1].mult(r);
					m[2].mult(r);
					t   .mult(r);
					return *this;
				}

		//----------------------------- EXPLICIT ARITHMETIC
		// *this = a op b
		kmatrix& negate (const kmatrix& b)
				{
					m[0].negate(b.m[0]);
					m[1].negate(b.m[1]);
					m[2].negate(b.m[2]);
					t   .negate(b.t   );
					return *this;
				}

		kmatrix& add (const kmatrix& a, const kmatrix& b)
				{
					m[0].add(a.m[0],b.m[0]);
					m[1].add(a.m[1],b.m[1]);
					m[2].add(a.m[2],b.m[2]);
					t   .add(a.t   ,b.t   );
					return *this;
				}

		kmatrix& sub (const kmatrix& a, const kmatrix& b)
				{
					m[0].sub(a.m[0],b.m[0]);
					m[1].sub(a.m[1],b.m[1]);
					m[2].sub(a.m[2],b.m[2]);
					t   .sub(a.t   ,b.t   );
					return *this;
				}

		kmatrix& mult (const kmatrix& a, const float d)
				{
					m[0].mult(a.m[0],d);
					m[1].mult(a.m[1],d);
					m[2].mult(a.m[2],d);
					t   .mult(a.t   ,d);
					return *this;
				}

		kmatrix& mult (const float d, const kmatrix& a)
				{
					m[0].mult(d,a.m[0]);
					m[1].mult(d,a.m[1]);
					m[2].mult(d,a.m[2]);
					t   .mult(d,a.t   );
					return *this;
				}

		kmatrix& div (const kmatrix& a, const float d)
				{
					float r = reciprocal(d);
					m[0].mult(a.m[0],r);
					m[1].mult(a.m[1],r);
					m[2].mult(a.m[2],r);
					t   .mult(a.t   ,r);
					return *this;
				}

		//----------------------------- MISC FUNCTIONS
		kmatrix& inverse();
		kmatrix& inverse(const kmatrix& a);

		//----------------------------- MATRIX MANIPULATION FUNCTIONS
		void SetPos (const kvector& v)	{	t = v;	};
#if 0
	void ApplyRotX();
	void ApplyRotY();
	void ApplyRotZ();
	void ApplyRotXY();
	void ApplyRotXZ();
	void ApplyRotYX();
	void ApplyRotYZ();
	void ApplyRotZX();
	void ApplyRotZY();
	void ApplyRotXYZ();
	void ApplyRotXZY();
	void ApplyRotYXZ();
	void ApplyRotYZX();
	void ApplyRotZXY();
	void ApplyRotZYX();

	void ApplyScaleX();
	void ApplyScaleY();
	void ApplyScaleZ();
	void ApplyScaleXY();
	void ApplyScaleXZ();
	void ApplyScaleYZ();
	void ApplyScaleXYZ();

	void ApplyTransX();
	void ApplyTransY();
	void ApplyTransZ();
	void ApplyTransXY();
	void ApplyTransXZ();
	void ApplyTransYZ();
	void ApplyTransXYZ();

	void ApplyMatrix();

	void SetPosX()
	void SetPosY()
	void SetPosZ()
	void SetPosXY()
	void SetPosXZ()
	void SetPosYZ()
	void SetPosXYZ()

	// GetMatrix();
	// GetPos();
#endif
};

//=============================================================================
//=============================================================================	Global space operators on math types
//=============================================================================

// -- kvector Operator Arithmetic operation --
inline kvector  operator - (const kvector& a)
		{
			return kvector(-a.n[0], -a.n[1], -a.n[2]);
		}

inline kvector operator + (const kvector& a, const kvector& b)
		{
			return kvector(a.n[0]+b.n[0], a.n[1]+b.n[1], a.n[2]+b.n[2]);
		}

inline kvector operator - (const kvector& a, const kvector& b)
		{
			return kvector(a.n[0]-b.n[0], a.n[1]-b.n[1], a.n[2]-b.n[2]);
		}

inline kvector operator * (const kvector& a, const float d)
		{
			return kvector(d*a.n[0], d*a.n[1], d*a.n[2]);
		}

inline kvector operator * (const float d, const kvector& a)
		{
			return a*d;
		}

inline float operator * (const kvector& a, const kvector& b)	// DOT PRODUCT
		{
			return (a.n[0]*b.n[0] + a.n[1]*b.n[1] + a.n[2]*b.n[2]);
		}

inline kvector operator / (const kvector& a, const float d)
		{
		 	float r=reciprocal(d);
			return kvector(a.n[0]*r, a.n[1]*r, a.n[2]*r);
		}

inline kvector operator ^ (const kvector& a, const kvector& b)	// CROSS PRODUCT
		{
			return kvector(	a.n[1]*b.n[2] - a.n[2]*b.n[1],
							a.n[2]*b.n[0] - a.n[0]*b.n[2],
							a.n[0]*b.n[1] - a.n[1]*b.n[0]);
		}

inline int operator == (const kvector& a, const kvector& b)
		{
			return (a.n[0]==b.n[0]) && (a.n[1]==b.n[1]) && (a.n[2]==b.n[2]);
		}

inline int operator != (const kvector& a, const kvector& b)
		{
			return !(a==b);
		}

// -- kmatrix Operator Arithmetic operation --
//		friend mat4 operator - (const mat4& a);			    // -m1
//		friend mat4 operator + (const mat4& a, const mat4& b);	    // m1 + m2
//		friend mat4 operator - (const mat4& a, const mat4& b);	    // m1 - m2
//		friend mat4 operator * (const mat4& a, const float d)	    // m1 * 4.0
//		friend mat4 operator * (const float d, const mat4& a)	    // 4.0 * m1
//		friend mat4 operator / (const mat4& a, const float d);	    // m1 / 3.0
//		friend int operator == (const mat4& a, const mat4& b);	    // m1 == m2 ?
//		friend int operator != (const mat4& a, const mat4& b);	    // m1 != m2 ?

inline kmatrix operator - (const kmatrix& a)
		{
			return kmatrix(-a.m[0], -a.m[1], -a.m[2], -a.t);
		}

inline kmatrix operator + (const kmatrix& a, const kmatrix& b)
		{
			return kmatrix(	a.m[0] + b.m[0],
							a.m[1] + b.m[1],
							a.m[2] + b.m[2],
							a.t    + b.t   );
		}

inline kmatrix operator - (const kmatrix& a, const kmatrix& b)
		{
			return kmatrix(	a.m[0] - b.m[0],
							a.m[1] - b.m[1],
							a.m[2] - b.m[2],
							a.t    - b.t   );
		}

inline kmatrix operator * (kmatrix& a, kmatrix& b)    		// m1 * m2
		{
		    #define ROWCOL(i,j) a.m[i].n[0]*b.m[0].n[j] + a.m[i].n[1]*b.m[1].n[j] + a.m[i].n[2]*b.m[2].n[j]
		    #define ROWT(i)     a.m[i].n[0]*b.t   .n[0] + a.m[i].n[1]*b.t   .n[1] + a.m[i].n[2]*b.t   .n[2] + a.t.n[i]
		    return kmatrix	(	kvector(ROWCOL(0,0), ROWCOL(0,1), ROWCOL(0,2)),
							    kvector(ROWCOL(1,0), ROWCOL(1,1), ROWCOL(1,2)),
							    kvector(ROWCOL(2,0), ROWCOL(2,1), ROWCOL(2,2)),
							    kvector(ROWT	(0),	// position
										ROWT	(1),
										ROWT	(2)	)
							);
			#undef ROWCOL
			#undef ROWT
		}

inline kmatrix operator * (const kmatrix& a, const float d)
		{
			return kmatrix(	a.m[0] * d,
							a.m[1] * d,
							a.m[2] * d,
							a.t    * d);
		}

inline kmatrix operator * (const float d, const kmatrix& a)
		{
			return kmatrix(	a.m[0] * d,
							a.m[1] * d,
							a.m[2] * d,
							a.t    * d);
		}

inline kmatrix operator / (const kmatrix& a, const float d)
		{
			float r = reciprocal(d);
			return kmatrix(	a.m[0] * r,
							a.m[1] * r,
							a.m[2] * r,
							a.t    * r);
		}

inline int operator == (const kmatrix& a, const kmatrix& b)
		{
			return (	(a.m[0] == b.m[0])
					&&	(a.m[1] == b.m[1])
					&&	(a.m[2] == b.m[2])
					&&  (a.t    == b.t   ));
		}

inline int operator != (const kmatrix& a, const kmatrix& b)
		{
			return !(a == b);
		}


// -- kvector-kmatrix Operator Arithmetic operation --
inline kvector operator * (const kmatrix& a, const kvector& v)    // linear transform
		{
		    #define ROWCOL(i) a.m[i].n[0]*v.n[0] + a.m[i].n[1]*v.n[1] + a.m[i].n[2]*v.n[2] + a.t.n[i]
		    return kvector (ROWCOL(0), ROWCOL(1), ROWCOL(2));
		    #undef ROWCOL
		}

inline kvector& ApplyMatrix (kvector& dst, const kmatrix& a, const kvector& src)    // linear transform
		{
		    dst.n[0] = a.m[0].n[0]*src.n[0] + a.m[0].n[1]*src.n[1] + a.m[0].n[2]*src.n[2] + a.t.n[0];
		    dst.n[1] = a.m[1].n[0]*src.n[0] + a.m[1].n[1]*src.n[1] + a.m[1].n[2]*src.n[2] + a.t.n[1];
		    dst.n[2] = a.m[2].n[0]*src.n[0] + a.m[2].n[1]*src.n[1] + a.m[2].n[2]*src.n[2] + a.t.n[2];
		    return dst;
		}

// -- Printing --
ostream& operator << (ostream& s, kvector& v);
ostream& operator << (ostream& s, kmatrix& m);

//=============================================================================
const kmatrix Identity	(	kvector	(1.0f, 0.0f, 0.0f),
							kvector (0.0f, 1.0f, 0.0f),
							kvector (0.0f, 0.0f, 1.0f),
							kvector (0.0f, 0.0f, 0.0f)	);

//=============================================================================
//=============================================================================
//=============================================================================

#endif  // Headerfile Guard

//=============================================================================	END-OF-INCLUDE