///////////////////////////////////////////////////////////////////////////////
/// @file deMath.hpp
///
/// @brief Math lib header
///
/// @author Assassin
///
/// This file is the intellectual property of Novus Delta, LLC.. Usage of the
/// contents of this file is subject to the Destiny3D Member License which
/// can be found at http://www.destiny3d.com.  Any other usage is prohibited.
///
/// This file is distributed "AS IS" without warranty of any kind.  Novus
/// Delta, LLC. does not guarantee the fitness of the contents of this file
/// for any particular purpose.
///
/// Copyright (C) 2001-2003 Novus Delta, LLC. All Rights Reserved.
///
/// <hr>
///                                 Change History
/// <hr>
///
/// @date Jan 2002
/// @author Assassin
/// @remarks Creation
///
/// @date Jul 2002
/// @author Hootie
/// @remarks Convert to static lib
///
///////////////////////////////////////////////////////////////////////////////

#ifndef DEMATH_HPP
#define DEMATH_HPP

#if defined(USING_DESTINY3D) || defined(USING_DESTINY3D_DLL)
#ifdef _DEBUG
// no static link option
#pragma comment(lib, "deMathd")
#else
#pragma comment(lib, "deMath")
#endif //_DEBUG
#endif //USING_DESTINY3D


#include "deGlobalTypes.hpp"
#include "deStack.hpp"
#include <cmath>



struct deFractionInfo
{
protected:
	int m_Num, m_Den;
	friend class deFraction;
};
class deFraction : public deFractionInfo
{
public:
	inline deFraction()
		{ m_Num = 0; m_Den = 1; }
	inline deFraction(int Num)
		{ m_Num = Num; m_Den = 1; }
	inline deFraction(int Num, int Den)
		{ m_Num = Num; m_Den = Den; }
	inline deFraction(const deFraction & ref)
		{ m_Num = ref.m_Num; m_Den = ref.m_Den; }
	inline deFraction(const deFractionInfo & ref)
		{ m_Num = ref.m_Num; m_Den = ref.m_Den; }
	inline ~deFraction() {}

	inline operator double() const { return ((double)m_Num)/m_Den; }
	inline operator float() const { return ((float)m_Num)/m_Den; }
	inline int Num() const { return m_Num; }
	inline int Den() const { return m_Den; }

	inline deFraction operator+(const deFraction & rhs) const
		{ deFraction ret(*this); ret += rhs; return ret; }
	inline deFraction operator-(const deFraction & rhs) const
		{ deFraction ret(*this); ret -= rhs; return ret; }
	inline deFraction operator*(const deFraction & rhs) const
		{ deFraction ret(*this); ret *= rhs; return ret; }
	inline deFraction operator/(const deFraction & rhs) const
		{ deFraction ret(*this); ret /= rhs; return ret; }
	
	deFraction operator+(int rhs) const;
	deFraction operator-(int rhs) const;
	deFraction operator*(int rhs) const;
	deFraction operator/(int rhs) const;

	deFraction& operator+=(const deFraction & rhs);
	deFraction& operator-=(const deFraction & rhs);
	inline deFraction& operator*=(const deFraction & rhs)
		{ m_Num *= rhs.m_Num; m_Den *= rhs.m_Den; return *this; }
	inline deFraction& operator/=(const deFraction & rhs)
		{ m_Num *= rhs.m_Den; m_Den *= rhs.m_Num; return *this; }
	
	inline const deFraction& operator=(const deFraction & rhs)
		{ m_Num = rhs.m_Num; m_Den = rhs.m_Den; return *this; }
	inline const deFraction& operator=(int rhs)
		{ m_Num = rhs; m_Den = 1; return *this; }
	bool operator==(const deFraction & rhs);
	bool operator==(int rhs);
};

/// storage class for a Destiny3D transform
struct deTransformInfo
{
public:
	deVec3d Translation;	///< translation, in coordinate space of parent system
	deVec3d R,U,I;			///< basis vectors, hopefully orthogonal for a rigid-body transform
	deBoolean m_Identity;	///< quick bypass value used if the basis vectors form a 3x3 identity matrix. Do not alter manually, use deTransform::TestIdentity
};
/// rigid-body transformation. supports a translation and rotation, with possible scaling.
class deTransform : public deTransformInfo
{
public:
	inline ~deTransform() {}
	/// default constructor initializes to identity matrix.
	inline deTransform(void)
		{ m_Identity = deTRUE; R = U = I = Translation = deZeroVec; R.X = U.Y = I.Z = 1.0; }
	/// copy constructor.
	inline deTransform(const deTransformInfo & ref)
		{ m_Identity = ref.m_Identity; R = ref.R; U = ref.U; I = ref.I; Translation = ref.Translation; }
	/// translation constructor, sets rotation to identity.
	deTransform(const deVec3d & Trans)
		{ m_Identity = deTRUE; R = U = I = deZeroVec; R.X = U.Y = I.Z = 1.0; Translation = Trans; }
	/// rotation constructor, sets translation to identity.
	deTransform(const deVec3d & Right, const deVec3d & Up, const deVec3d & In);
	/// all data constructor.
	deTransform(const deVec3d & Right, const deVec3d & Up, const deVec3d & In, const deVec3d & Trans);
	
	/// reset the transform to 0 translation and identity rotation
	void inline Reset()
		{ m_Identity = deTRUE; R = U = I = Translation = deZeroVec; R.X = U.Y = I.Z = 1.0; }
	/// convert from a 4x4 matrix, such as an opengl or direct3d one. the last column will be ignored.
	void FromFloat44(const float Mat[4][4], deBoolean d3d);
	/// convert from a 4x4 matrix, such as deMatrix44.
	void FromDouble44(const double Mat[4][4], deBoolean d3d);
	/// convert to a 4x4 matrix, such as an opengl or direct3d one. the last column will be created with (0,0,0,1).
	void ToFloat44(float Mat[4][4], deBoolean d3d) const;
	/// forces any values that are close to 0 to become exactly 0
	void SqueezeVectors();
	/// set the "right" (aka X or i) basis vector
	void SetRight(const deVec3d & Right);
	/// set the "up" (aka Y or j) basis vector
	void SetUp(const deVec3d & Up);
	/// set the "in" (aka Z or k) basis vector
	void SetIn(const deVec3d & In);
	/// set the translation of the transform. translation is in the space of the parent's coordinate system.
	void inline SetTranslation(const deVec3d & Trans) { Translation = Trans; }
	/// set all the basis vectors at once
	void SetRightUpIn(const deVec3d & Right, const deVec3d & Up, const deVec3d & In);
	/// Reset the transform and then build a rotation around the X axis. Any existing rotation or translation is lost.
	void SetXRot(deDouble r);
	/// Reset the transform and then build a rotation around the Y axis. Any existing rotation or translation is lost.
	void SetYRot(deDouble r);
	/// Reset the transform and then build a rotation around the Z axis. Any existing rotation or translation is lost.
	void SetZRot(deDouble r);
	/// Applies a rotation around the X axis to the existing data in the transform.
	void RotateX(deDouble Rads);
	/// Applies a rotation around the Y axis to the existing data in the transform.
	void RotateY(deDouble Rads);
	/// Applies a rotation around the Z axis to the existing data in the transform.
	void RotateZ(deDouble Rads);
	/// Builds a Euler angles rotation, losing any previously stored rotation or translation.
	void FreeRot(const deVec3d &Rads);
	/// Builds a Euler angles rotation, losing any previously stored rotation or translation.
	void FreeRot(deDouble x, deDouble y, deDouble z);
	/// Builds a Euler angles rotation, losing any previously stored rotation, but not translation.
	void FreeRotInPlace(const deVec3d &Rads);
	/// Builds a Euler angles rotation, losing any previously stored rotation, but not translation.
	void FreeRotInPlace(deDouble x, deDouble y, deDouble z);
	/// Scale the basis vectors by a uniform value.
	void ApplyScaling(deDouble UniformScaling);
	/// Scale the basis vectors by a distinct axis values.
	void ApplyScaling(const deVec3d & Multiples);
	/// Scale the basis vectors by a distinct axis values.
	void ApplyScaling(deDouble x, deDouble y, deDouble z);
	/// Obtain the lengths of the basis vectors (axes) of the transform
	void GetBasisScales(deDouble& x, deDouble& y, deDouble &z);
	inline void GetBasisScales(deVec3d& v) { GetBasisScales(v.x, v.y, v.z); }

	/// Tests whether the basis vectors form an identity matrix, setting the internal m_Identity member appropriately.
	deBoolean TestIdentity(deDouble Epsilon = 0.00001);
	/// Tests whether the basis vectors form an orthogonal basis (all perpendicular to one another).
	deBoolean IsOrthogonal() const;
	/// Tests the basis vectors for an orthonormality (all perpendicular to one another and magnitude=1.0).
	deBoolean IsOrthoNormal() const;
	/// Creates an orthonormal basis from the existing basis vectors.
	deBoolean OrthoNormalize();

	/// retrieve the "right" basis vector
	void GetRight(deVec3d & r) const;
	/// retrieve the "up" basis vector
	void GetUp(deVec3d & u) const;
	/// retrieve the "in" basis vector
	void GetIn(deVec3d & i) const;
	/// retrieve the "right" basis vector
	deVec3d GetRight() const;
	/// retrieve the "up" basis vector
	deVec3d GetUp() const;
	/// retrieve the "in" basis vector
	deVec3d GetIn() const;
	/// retrieve all the basis vectors at once
	void GetRightUpIn(deVec3d & Right, deVec3d & Up, deVec3d & In) const;
	/// Get a transposition of the rotation matrix, combined with a processed translation.
	/// On an orthonormal transform, this is the same operation as computing the inverse, but much quicker to compute.
	void GetTranspose(deTransformInfo & target);
	/// alternate form of GetTranspose
	deTransformInfo GetTranspose();
	/// retrieve the translation for this transform
	deVec3d inline GetTranslation() const
		{ return Translation; }

	/// Perform a transformation of a 3d vector using this transform.
	void Transform(deVec3d InVec, deVec3d & OutVec) const;
	/// Perform a transformation of an array of 3d vectors using this transform.
	void TransformArray(const deVec3d InVec[], deVec3d OutVec[], int num) const;
	void Transform(deVertex InVec, deVertex & OutVec) const;
	void TransformArray(const deVertex InVec[], deVertex OutVec[], int num) const;
	/// transform a 3d vector by the transposition of this transform, which is the same as the inverse for orthonormal transforms.
	void TransposeTransform(deVec3d InVec, deVec3d & OutVec) const;
	/// concatenate this transform with another transform, storing the result in "ret".
	/// ret = this * rhs
	void Multiply(const deTransformInfo & rhs, deTransformInfo & ret);
	/// return = this * rhs
	deTransform Multiply(const deTransformInfo & rhs);
};

struct deVec4d
{
	deDouble x,y,z,w;
};
class deMatrix44
{
	// uses direct3d matrix layout
public:
    union {
        struct {
            deDouble        _00, _01, _02, _03;
            deDouble        _10, _11, _12, _13;
            deDouble        _20, _21, _22, _23;
            deDouble        _30, _31, _32, _33;
        };
        deDouble m[4][4];
    };
	inline ~deMatrix44() {}
	inline deMatrix44() {}
	inline deMatrix44(const deTransformInfo & ref)
	{	
		_00 = ref.R.x; _01 = ref.R.y; _02 = ref.R.z; _03 = 0;
		_10 = ref.U.x; _11 = ref.U.y; _12 = ref.U.z; _13 = 0;
		_20 = ref.I.x; _21 = ref.I.y; _22 = ref.I.z; _23 = 0;
		_30 = ref.Translation.x; _31 = ref.Translation.y; _32 = ref.Translation.z;
		_33 = 1;
	}
	deMatrix44(const deVec3d & Right, const deVec3d & Up, const deVec3d & In, const deVec3d & Trans = deZeroVec)
	{	
		_00 = Right.x; _01 = Right.y; _02 = Right.z; _03 = 0;
		_10 = Up.x;    _11 = Up.y;    _12 = Up.z;    _13 = 0;
		_20 = In.x;    _21 = In.y;    _22 = In.z;    _23 = 0;
		_30 = Trans.x; _31 = Trans.y; _32 = Trans.z; _33 = 1;
	}
	deMatrix44(const deFloat f[4][4], deBoolean d3d);
	deMatrix44(const deDouble d[4][4], deBoolean d3d);
	void operator=(const deDouble d[4][4]);

	void Identity()
	{
		_00 = 1; _01 = 0; _02 = 0; _03 = 0;
		_10 = 0; _11 = 1; _12 = 0; _13 = 0;
		_20 = 0; _21 = 0; _22 = 1; _23 = 0;
		_30 = 0; _31 = 0; _32 = 0; _33 = 1;
	}
	void Translation(const deVec3d & T)
	{
		_00 = 1;   _01 = 0;   _02 = 0;   _03 = 0;
		_10 = 0;   _11 = 1;   _12 = 0;   _13 = 0;
		_20 = 0;   _21 = 0;   _22 = 1;   _23 = 0;
		_30 = T.x; _31 = T.y; _32 = T.z; _33 = 1;
	}
	void Scaling(deDouble x, deDouble y, deDouble z)
	{
		_00 = x; _01 = 0; _02 = 0; _03 = 0;
		_10 = 0; _11 = y; _12 = 0; _13 = 0;
		_20 = 0; _21 = 0; _22 = z; _23 = 0;
		_30 = 0; _31 = 0; _32 = 0; _33 = 1;
	}
	void Multiply(const deMatrix44 & right, deMatrix44 & out) const;
	void Transform(const deVec3d & in, deVec3d & out) const;
	void Transform(const deVec3d & in, deVec4d & out) const;
	void Transform(const deVertex & in, deVec4d & out) const;
};


/// a class for maintaining a hierarchy of transformations.
/// concatenations are performed with first-in on the left.
class deTransformStack
{
private:
	deTStack <deTransform> m_Stack;

public:
	inline deTransformStack() : m_Stack(deTStack<deTransform>::LIFO_STACK) {}
	inline ~deTransformStack() {}

	/// compute the concatenation of a transform with the previous top of the stack, and store the concatenation as the new top.
	deBoolean Push(const deTransformInfo & transform, deTransformInfo * concat = NULL);
	/// remove the last-inserted transform from the stack
	deBoolean Pop(deTransformInfo & transform);
	/// remove the last-inserted transform from the stack
	deBoolean Pop();
	/// view the concatenation of all the currently-stored transforms
	deBoolean Peek(deTransformInfo & transform);

	/// remove all the transforms from this structure
	inline deBoolean Clear() { m_Stack.Clear();}
};


/// a quaternion calculation class
class deQuaternion {
public:
	union
	{	struct
		{
			deDouble x, y, z;
		};
		deVec3d vector;
	};
	deDouble w;

	// A quaternion's length must be closer than this tolerance to 1.0 to be
	// considered a unit quaternion
	static const deDouble UNIT_TOLERANCE;

	// Epsilon used in spherical linear interpolation
	static const deDouble SLERP_EPSILON;

	// Operators
	void operator +=(const deQuaternion&);
	void operator -=(const deQuaternion&);
	void operator *=(const deQuaternion&);
	deQuaternion operator +(const deQuaternion&) const;
	deQuaternion operator -(const deQuaternion&) const;
	deQuaternion operator *(const deQuaternion&) const;
	bool operator ==(const deQuaternion&) const;
	bool operator !=(const deQuaternion&) const;

	// Conversion operators
	operator deVec3d&(){ return vector; }
	operator const deVec3d&() const { return vector; }

	// Analytical functions
	deBoolean IsValid() const;
	deBoolean IsUnit() const;

	// Math operations
	// These are slightly faster than the +, -, * operators since they don't use the stack to
	// return the result.
	// NOTE: These functions are still slower than the +=, -= or *= operators.
	void Multiply(const deQuaternion &q, deQuaternion &result) const;
	void Add(const deQuaternion &q, deQuaternion &result) const;
	void Subtract(const deQuaternion &q, deQuaternion &result) const;

	void Set(deDouble xAxis, deDouble yAxis, deDouble zAxis, deDouble wAxis)
		{ x = xAxis;	y = yAxis;	z = zAxis;	w = wAxis; }

	// Length of quaternion (also known as 'magnitude')
	deDouble Length() const;

	// Normalize quaternion
	void Normalize();

	// Interpolate between two quaternions
	// If balance is 0, the result is 'this', if balance is 1.0 the result is 'q'
	// NOTE: The balance must >= 0 and <= 1
	void Interpolate(const deQuaternion &q, deDouble balance, deQuaternion &result);

	void RotateX(deDouble radians){ x += radians; }
	void RotateY(deDouble radians){ y += radians; }
	void RotateZ(deDouble radians){	z += radians; }

	// Conversion
	deBoolean FromMatrix(const deTransformInfo &matrix);

	// NOTE: The translation of the matrix remains unchanged.
	void ToMatrix(deTransformInfo &matrix) const;

	deQuaternion(){}
	deQuaternion(deDouble xAxis, deDouble yAxis, deDouble zAxis, deDouble wAxis){ x = xAxis;	y = yAxis;	z = zAxis;	w = wAxis; }
	deQuaternion(deVec3d &axisdeVector, deDouble www){ vector = axisdeVector; w = www;}

};



template <typename T>
static inline T SQUARE(T x) {return (x)*(x);}
/// squared length of a deVec3d
static inline deDouble Vec3d_LengthSq(const deVec3d& v)
	{return SQUARE(v.X) + SQUARE(v.Y) + SQUARE(v.Z);}
/// length of a deVec3d
static inline deDouble Vec3d_Length(const deVec3d& v)
	{return sqrt(Vec3d_LengthSq(v));}
/// multiply each component of a deVec3d by the corresponding component of another
static inline void Vec3d_ScaleByVec(deVec3d & v1, const deVec3d & v2)
	{ v1.X *= v2.X; v1.Y *= v2.Y; v1.Z *= v2.Z; }
/// compare equality of 2 floating-point values, within some Epsilon value
static inline deBoolean Compare_Epsilon(const deDouble d1, const deDouble d2, deDouble epsilon = 0.00001)
	{ return (fabs(d1 - d2) <= epsilon); }
/// compare equality of two deVec3d, within some Epsilon value
static inline deBoolean Vec3d_Compare(const deVec3d & v1, const deVec3d &v2, deDouble epsilon = 0.00001)
	{ return ((fabs(v1.x - v2.x) <= epsilon) && (fabs(v1.y - v2.y) <= epsilon) && (fabs(v1.z - v2.z) <= epsilon)); }
/// compare exact opposite equality of two deVec3d, within some Epsilon value
static inline deBoolean Vec3d_CompareOpposite(const deVec3d & v1, const deVec3d &v2, deDouble epsilon = 0.00001)
	{ return ((fabs(v1.x + v2.x) <= epsilon) && (fabs(v1.y + v2.y) <= epsilon) && (fabs(v1.z + v2.z) <= epsilon)); }
/// compute the dot-product of two deVec3d (deDouble components)
static inline deDouble Vec3d_Dot(const deVec3d& v1, const deVec3d& v2)
	{ return (v1).X*(v2).X+(v1).Y*(v2).Y+(v1).Z*(v2).Z; }
/// compute the dot-product of two deVertex (deFloat components)
static inline deFloat Vec3d_Dot(const deVertex& v1, const deVertex& v2)
	{ return (v1).X*(v2).X+(v1).Y*(v2).Y+(v1).Z*(v2).Z; }
/// compute the square of the distance between two deVec3d
static inline deDouble Vec3d_DistSq(const deVec3d & v1, const deVec3d & v2) 
	{ 
		deDouble Temp1, Temp2; 
		Temp1=(v1).x-(v2).x; Temp2 = SQUARE(Temp1); 
		Temp1=(v1).y-(v2).y; Temp2 += SQUARE(Temp1);
		Temp1=(v1).z-(v2).z; Temp2 += SQUARE(Temp1);
		return Temp2;
	}
/// compute the distance between two deVec3d
static inline deDouble Vec3d_Dist(const deVec3d & v1, const deVec3d & v2) 
	{ return (deDouble)sqrt(Vec3d_DistSq(v1,v2));}
/// compute the square of the distance between two deVertex
static inline deFloat Vec3d_DistSq(const deVertex & v1, const deVertex & v2) 
	{ 
		deFloat Temp1, Temp2; 
		Temp1=(v1).x-(v2).x; Temp2 = SQUARE(Temp1); 
		Temp1=(v1).y-(v2).y; Temp2 += SQUARE(Temp1);
		Temp1=(v1).z-(v2).z; Temp2 += SQUARE(Temp1);
		return Temp2;
	}
/// compute the distance between two deVertex
static inline deFloat Vec3d_Dist(const deVertex & v1, const deVertex & v2) 
	{ return (deFloat)sqrtf(Vec3d_DistSq(v1, v2));}
/// compute the inverse of a deVec3d (overloaded operator-)
static inline deVec3d operator-(const deVec3d &v)
	{deVec3d ret={-v.x,-v.y,-v.z};return ret;}
/// compute the component difference between two deVec3d (overloaded operator-)
static inline deVec3d operator-(const deVec3d &v1, const deVec3d &v2)
	{deVec3d ret={v1.x-v2.x,v1.y-v2.y,v1.z-v2.z};return ret;}
/// compute the component sum of two deVec3d (overloaded operator+)
static inline deVec3d operator+(const deVec3d &v1, const deVec3d &v2)
	{deVec3d ret={v1.x+v2.x,v1.y+v2.y,v1.z+v2.z};return ret;}
/// compute the scalar product of a deVec3d and a floating-point value (overloaded operator*)
static inline deVec3d operator*(const deVec3d &v1, deDouble d)
	{deVec3d ret={v1.x*d,v1.y*d,v1.z*d};return ret;}
/// subtract a deVec3d from another one (overloaded operator-=)
static inline void operator-=(deVec3d &v1, const deVec3d &v2)
	{v1.x -= v2.x; v1.y -= v2.y; v1.z -= v2.z; }
/// add a deVec3d to another one (overloaded operator+=)
static inline void operator+=(deVec3d &v1, const deVec3d &v2)
	{v1.x += v2.x; v1.y += v2.y; v1.z += v2.z; }
/// multiply a deVec3d by a scalar value (overloaded operator*=)
static inline void operator*=(deVec3d &v1, deDouble d)
	{v1.x *= d; v1.y *= d; v1.z *= d; }
/// divide a deVec3d by a scalar value (overloaded operator/=)
static inline void operator/=(deVec3d &v1, deDouble d)
	{deDouble d2 = 1/d; v1.x *= d2; v1.y *= d2; v1.z *= d2; }
#ifdef _DEBUG
// this should work even on packing of 8 due to a double's size
static inline deDouble & Vec3d_GetElement(const deVec3d& v, long i)
	{if (i < 0 || i >= 3) DEBUGRETURN(*((deDouble*)&v),"Vec3d_GetElement: Invalid index"); return *( ((deDouble*)&v) + i ); }
#else
/// obtain a reference to an indexed value inside a 3d vector (valid index values: {0,1,2}, map to {x,y,z})
static inline deDouble & Vec3d_GetElement(const deVec3d& v, long i)
	{return *(((deDouble*)&v) + i); }
#endif

/// compute the inverse of a deVertex (overloaded operator-)
static inline deVertex operator-(const deVertex &v)
	{deVertex ret={-v.x,-v.y,-v.z};return ret;}
/// compute the component difference between two deVertex (overloaded operator-)
static inline deVertex operator-(const deVertex &v1, const deVertex &v2)
	{deVertex ret={v1.x-v2.x,v1.y-v2.y,v1.z-v2.z};return ret;}
/// compute the component sum of two deVertex (overloaded operator+)
static inline deVertex operator+(const deVertex &v1, const deVertex &v2)
	{deVertex ret={v1.x+v2.x,v1.y+v2.y,v1.z+v2.z};return ret;}
/// compute the scalar product of a deVertex and a floating-point value (overloaded operator*)
static inline deVertex operator*(const deVertex &v1, deFloat f)
	{deVertex ret={v1.x*f,v1.y*f,v1.z*f};return ret;}
/// subtract a deVertex from another one (overloaded operator-=)
static inline void operator-=(deVertex &v1, const deVertex &v2)
	{v1.x -= v2.x; v1.y -= v2.y; v1.z -= v2.z; }
/// add a deVertex to another one (overloaded operator+=)
static inline void operator+=(deVertex &v1, const deVertex &v2)
	{v1.x += v2.x; v1.y += v2.y; v1.z += v2.z; }
/// multiply a deVertex by a scalar value (overloaded operator*=)
static inline void operator*=(deVertex &v1, deFloat d)
	{v1.x *= d; v1.y *= d; v1.z *= d; }
/// divide a deVertex by a scalar value (overloaded operator/=)
static inline void operator/=(deVertex &v1, deFloat d)
	{deFloat d2 = 1/d; v1.x *= d2; v1.y *= d2; v1.z *= d2; }

static inline deTexCoord operator-(const deTexCoord &v1, const deTexCoord &v2)
	{deTexCoord ret={v1.u-v2.u,v1.v-v2.v};return ret;}
static inline deTexCoord operator+(const deTexCoord &v1, const deTexCoord &v2)
	{deTexCoord ret={v1.u+v2.u,v1.v+v2.v};return ret;}
static inline deTexCoord operator*(const deTexCoord &v1, deFloat f)
	{deTexCoord ret={v1.u*f,v1.v*f};return ret;}
static inline void operator-=(deTexCoord &v1, const deTexCoord &v2)
	{v1.u -= v2.u; v1.v -= v2.v; }
static inline void operator+=(deTexCoord &v1, const deTexCoord &v2)
	{v1.u += v2.u; v1.v += v2.v; }
static inline void operator*=(deTexCoord &v1, deFloat d)
	{v1.u *= d; v1.v *= d; }
static inline void operator/=(deTexCoord &v1, deFloat d)
	{deFloat d2 = 1/d; v1.u *= d2; v1.v *= d2; }

/// multiply each component of a deColor by the corresponding component of another
static inline void operator*=(deColor& lhs, const deColor& rhs)
	{lhs.a *= rhs.a; lhs.r *= rhs.r; lhs.g *= rhs.g; lhs.b *= rhs.b;}

/// "squeeze" the components of a deVec3d: if they are within Epsilon units of zero, set them to zero.
/// Useful for getting "perfect" orthographic orientations after moving around a bit
void Vec3d_Squeeze(deVec3d &v, deDouble Epsilon = 0.00001);
/// compute the cross product of two deVec3d
void Vec3d_Cross(const deVec3d &V1, const deVec3d &V2, deVec3d &Result);
/// compute the cross product of two deVertex
void Vec3d_Cross(const deVertex &V1, const deVertex &V2, deVertex &Result);
static inline deVec3d Vec3d_Cross(const deVec3d &V1, const deVec3d &V2)
	{deVec3d ret; Vec3d_Cross(V1,V2,ret); return ret;}
/// normalize a deVec3d
/// @return the length of the original vector
deDouble Vec3d_Normalize(deVec3d &V);
/// normalize a deVertex
/// @return the length of the original vector
deFloat Vertex_Normalize(deVertex &V);
/// compute the 32-bit deARGB value for a normalized deVec3d (for dot3 bumpmapping)
deARGB Vec3d_ConvertToDot3(const deVec3d &V1, deFloat Height);
deARGB Vec3d_ConvertToDot3Clamp(const deVec3d &V1, deFloat Height);

/// Some useful geometry methods
namespace deGeometry
{
	// standard true/false test functions
	/// tests whether the inner AABB is actually inside the outer AABB
	deBoolean BoxInBox(const deVec3d & OuterMin, const deVec3d & OuterMax, const deVec3d & InnerMin, const deVec3d & InnerMax);
	/// tests whether the inner box is actually inside the outer box
	deBoolean BoxInBox(const deAABB & OuterBox, const deAABB & InnerBox);
	/// tests whether two AABBs intersect
	deBoolean BoxBoxIntersect(const deVec3d & Min1, const deVec3d & Max1, const deVec3d & Min2, const deVec3d & Max2);
	/// tests whether two AABBs intersect
	deBoolean BoxBoxIntersect(const deAABB & Box1, const deAABB & Box2);
	/// tests whether a point lies within an AABB
	deBoolean PointInBox(const deVec3d & pt, const deVec3d &Min, const deVec3d &Max);
	/// tests whether a point lies within an AABB
	deBoolean PointInBox(const deVec3d & pt, const deAABB & Box);
	/// returns the distance from a plane to a point, along the plane normal
	static inline deDouble PlaneDist(const deVec3d & pt, const deVec3d & norm, const deDouble & dist)
	{
		return Vec3d_Dot(pt, norm) - dist;
	}
	/// tests whether a point lies on the front of a plane (positive half-space)
	static inline deBoolean PointOnPlaneFront(const deVec3d & pt, const deVec3d & norm, const deDouble & dist)
	{ return (Vec3d_Dot(pt,norm) >= dist) ? deTRUE : deFALSE; }	
	/// tests whether a point lies on the back of a plane (negative half-space)
	static inline deBoolean PointOnPlaneBack(const deVec3d & pt, const deVec3d & norm, const deDouble & dist)
	{ return (Vec3d_Dot(pt,norm) <= dist) ? deTRUE : deFALSE; }	
	/// tests whether a point lies within an AABB
	static inline deBoolean PointInAABB(const deVec3d &Point, const deVec3d &Min, const deVec3d &Max)
	{
		if (Point.x < Min.x || Point.y < Min.y || Point.z < Min.z) return deFALSE;
		if (Point.x > Max.x || Point.y > Max.y || Point.z > Max.z) return deFALSE;
		return deTRUE;
	}
	/// calculates the geometric center of an AABB
	static inline void AABBCenter(const deAABB & bbox, deVec3d & vec)
	{
		vec = bbox.Min; vec += bbox.Max; vec *= 0.5;
	}
	/*
	needs fixing, radius is conservative
	/// calculates circumscribed radius of an AABB, from the origin
	static inline deDouble AABBRadius(const deAABB & BBox)
	{
		deDouble len1, len2;
		len1 = Vec3d_LengthSq(BBox.Min);
		len2 = Vec3d_LengthSq(BBox.Max);
		if (len1 < 1 && len2 < 1)
			len1 = sqrt(min(len1, len2));
		else
			len1 = sqrt(max(len1, len2));
		return len1;
	}
	*/
	/// calculates diagonal size of an AABB, from min to max points
	static inline deDouble AABBDiagonal(const deAABB & BBox)
	{
		deVec3d len;
		len = BBox.Max;
		len -= BBox.Min;
		return Vec3d_Length(len);
	}
	/// tests whether a point lies within an AABB
	static inline deBoolean PointInAABB(const deVec3d &Point, const deAABB & Box)
	{ return PointInAABB(Point, Box.Min, Box.Max); }
	/// expand a given AABB by adding a point. Result: BBox is enlarged to hold input BBox space and pt
	static inline void AABBAddPoint(deAABB & BBox, const deVec3d & pt)
	{
		BBox.Min.x = min(BBox.Min.x, pt.x);
		BBox.Min.y = min(BBox.Min.y, pt.y);
		BBox.Min.z = min(BBox.Min.z, pt.z);
		BBox.Max.x = max(BBox.Max.x, pt.x);
		BBox.Max.y = max(BBox.Max.y, pt.y);
		BBox.Max.z = max(BBox.Max.z, pt.z);
	}
	/// compute the union of two AABBs, and store the resulting AABB in the first argument (In_Out)
	static inline void AABBUnion(deAABB & In_Out, const deAABB & Input)
	{
		In_Out.Min.x = min(In_Out.Min.x, Input.Min.x);
		In_Out.Min.y = min(In_Out.Min.y, Input.Min.y);
		In_Out.Min.z = min(In_Out.Min.z, Input.Min.z);
		In_Out.Max.x = max(In_Out.Max.x, Input.Max.x);
		In_Out.Max.y = max(In_Out.Max.y, Input.Max.y);
		In_Out.Max.z = max(In_Out.Max.z, Input.Max.z);
	}
	/// test whether two AABBs intersect, using a half-open interval
	static inline deBoolean AABBIntersect(const deVec3d &Min1, const deVec3d &Max1, const deVec3d &Min2, const deVec3d &Max2)
	{
		// if one box's max is less than the other's min, they don't intersect
		if (Max1.x < Min2.x || Max1.y < Min2.y || Max1.z < Min2.z) return deFALSE;
		// and vice-versa
		if (Max2.x < Min1.x || Max2.y < Min1.y || Max2.z < Min1.z) return deFALSE;
		// otherwise, they do (I think)
		return deTRUE;
	}
	/// test whether two AABBs intersect, using a half-open interval
	static inline deBoolean AABBIntersect(const deAABB & Box1, const deAABB & Box2)
	{ return AABBIntersect(Box1.Min, Box1.Max, Box2.Min, Box2.Max); }
	/// tests whether the inner AABB is actually inside the outer AABB
	static inline deBoolean AABBInAABB(const deVec3d &InMin, const deVec3d &InMax, const deVec3d &OutMin, const deVec3d &OutMax)
	{
		// if the inner box's min is less than the outer box's min, it's not inside
		if (InMin.x < OutMin.x || InMin.y < OutMin.y || InMin.z < OutMin.z) return deFALSE;
		// and opposite relation with max
		if (InMax.x > OutMax.x || InMax.y > OutMax.y || InMax.z > OutMax.z) return deFALSE;
		// otherwise, they do (I think)
		return deTRUE;
	}
	/// tests whether the inner AABB is actually inside the outer AABB
	static inline deBoolean AABBInAABB(const deAABB & InBox, const deAABB & OutBox)
	{ return AABBInAABB(InBox.Min, InBox.Max, OutBox.Min, OutBox.Max); }
	/// classify where an AABB lies relative to a plane.
	/// -1 => all negative, 0 => intersects, 1 => all positive
	int AABBPlaneClassify(const deAABB & BBox, const deVec3d & Norm, const deDouble & Dist);
	
	// these take a non-const reference to fill in values
	/// intersect a plane and a line segment
	deBoolean PlaneLineIntersect(const deVec3d & Norm, const deDouble & Dist, const deVec3d & Pt1, const deVec3d & Pt2, deVec3d & Res, deDouble & Percent);
	/// intersect a plane and a ray
	deBoolean PlaneRayIntersect(const deVec3d & Norm, const deDouble & Dist, const deVec3d & Start, const deVec3d & Dir, deDouble & Multiple);
	/// intersect a sphere and a ray
	deBoolean SphereRayIntersect(const deVec3d & Center, const deDouble & Radius, const deVec3d & Start, const deVec3d & Dir, deDouble & Multiple);
	/// intersect a sphere and a line segment
	deBoolean SphereLineIntersect(const deVec3d & Center, const deDouble & Radius, const deVec3d & Start, const deVec3d & End, deVec3d & Res, deDouble & Percent);
	/// perform a linear-sweep AABB collision test
	deBoolean AABBCollision(const deAABB & MovingBox, const deVec3d & MoveStart, const deVec3d & MoveEnd, 
		const deAABB & StaticBox, const deVec3d & StaticPos, deVec3d & MovingColPt, deVec3d & ColNormal, deDouble & Percent);
	/// create an AABB that contains a set of points
	deBoolean ConstructAABB(const deVec3d * VecArray, long NumVecs, deVec3d & Min, deVec3d & Max);
	/// create an OBB using an AABB and a transformation
	deBoolean MakeOBBFromAABB(const deAABB & LocalBox, const deTransform & WorldMatrix, deOBB & WorldBox);
	/// get the 8 corners of an AABB
	deBoolean GetPointsOfAABB(const deAABB & BBox, deVec3d * ArrayOf8);
	/// get the 8 corners of an OBB
	deBoolean GetPointsOfOBB(const deOBB & WorldBox, deVec3d * ArrayOf8);
	/// create an AABB that contains a set of points
	static inline deBoolean ConstructAABB(const deVec3d * VecArray, long NumVecs, deAABB & Box)
	{ return ConstructAABB(VecArray, NumVecs, Box.Min, Box.Max); }
	/// create an AABB that contains an OBB
	static inline deBoolean MakeAABBFromOBB(const deOBB & InBox, deAABB & OutBox)
	{
		deVec3d vecarray[8];
		if (!GetPointsOfOBB(InBox, vecarray)) return deFALSE;
		if (!ConstructAABB(vecarray, 8, OutBox)) return deFALSE;
		return deTRUE;
	}
	/// create a BSphere from an AABB
	static inline void MakeBoundSphereFromAABB(const deAABB & BBox, deBoundSphere & OutSphere)
	{
		OutSphere.Radius = AABBDiagonal(BBox) * 0.5;
		AABBCenter(BBox, OutSphere.Center);
	}
	/// intersect an OBB and a ray
	deBoolean RayOBBIntersect(const deVec3d & RayStart, const deVec3d & RayDir, const deOBB & Box, deVec3d & Result, deVec3d & Normal, deDouble & Percent);
	
	/// intersect a plane and a line segment
	static inline deBoolean PlaneLineIntersect(const dePlane & Plane, const deVec3d & Pt1, const deVec3d & Pt2, deVec3d & Res, deDouble & Percent)
	{ return PlaneLineIntersect(Plane.Norm, Plane.Dist, Pt1, Pt2, Res, Percent); }
	/// intersect a plane and a ray
	static inline deBoolean PlaneRayIntersect(const dePlane & Plane, const deVec3d & Start, const deVec3d & Dir, deDouble & Multiple)
	{ return PlaneRayIntersect(Plane.Norm, Plane.Dist, Start, Dir, Multiple); }

	/// interpolate a Vec3d over a range [0,1]
	static inline void Vec3d_Interpolate(deVec3d & Output, const deVec3d & In1, const deVec3d& In2, deDouble Percent)
	{
		Output.x = In1.x*(1-Percent) + In2.x*Percent;
		Output.y = In1.y*(1-Percent) + In2.y*Percent;
		Output.z = In1.z*(1-Percent) + In2.z*Percent;
	}
	/// interpolate a Vertex over a range [0,1]
	static inline void Vertex_Interpolate(deVertex & Output, const deVertex & In1, const deVertex& In2, deDouble Percent)
	{
		Output.x = deFloat(In1.x*(1-Percent) + In2.x*Percent);
		Output.y = deFloat(In1.y*(1-Percent) + In2.y*Percent);
		Output.z = deFloat(In1.z*(1-Percent) + In2.z*Percent);
	}
	/// interpolate a TexCoord2 over a range [0,1]
	static inline void TexCoord_Interpolate(deTexCoord2 & Output, const deTexCoord2 & In1, const deTexCoord2& In2, deDouble Percent)
	{
		Output.u = deFloat(In1.u*(1-Percent) + In2.u*Percent);
		Output.v = deFloat(In1.v*(1-Percent) + In2.v*Percent);
	}
	/// interpolate an ARGB over a range [0,1]
	static inline void ARGB_Interpolate(deARGB& Output, deARGB Color1, deARGB Color2, deFloat Percent)
	{
		Output = 
			(deARGB((1-Percent)*((Color1>>24)&0xff) + (Percent)*((Color2>>24)&0xff)) << 24) |
			(deARGB((1-Percent)*((Color1>>16)&0xff) + (Percent)*((Color2>>16)&0xff)) << 16) |
			(deARGB((1-Percent)*((Color1>> 8)&0xff) + (Percent)*((Color2>> 8)&0xff)) << 8)  |
			(deARGB((1-Percent)*((Color1    )&0xff) + (Percent)*((Color2    )&0xff)));
	}
};


#endif
