#ifndef MATH_H
#define MATH_H

#include <cmath>
#include <limits>

#ifndef ASSERT
    #define ASSERT(condition) if ( !(condition) ) throw; else do {} while (0)
#endif

#ifndef NULL
    #define NULL 0
#endif

namespace MATH
{
    extern double g_dblTolerance;

    inline bool
    Zero(const float p_fltValue)
    {
        return ( (p_fltValue * p_fltValue) < (float)(MATH::g_dblTolerance) );
    }

    inline bool
    Zero(const double p_dblValue)
    {
        return ( (p_dblValue * p_dblValue) < MATH::g_dblTolerance );
    }

    inline bool
    Equal(const int p_intValue1, const int p_intValue2)
    {
        return ( p_intValue1 == p_intValue2 );
    }

    inline bool
    Equal(const unsigned int p_uintValue1, const unsigned int p_uintValue2)
    {
        return ( p_uintValue1 == p_uintValue2 );
    }

    inline bool
    Equal(const float p_fltValue1, const float p_fltValue2)
    {
        return Zero(p_fltValue1 - p_fltValue2);
    }

    inline bool
    Equal(const double p_dblValue1, const double p_dblValue2)
    {
        return Zero(p_dblValue1 - p_dblValue2);
    }

    inline float
    FastSqrt(const float Value)
    {
        float fltValues[3] = {Value, Value, Value};

        *((int *)(&(fltValues[1]))) = 0x1fbcf800 + (*((int *)(&(fltValues[1]))) >> 1);
        *((int *)(&(fltValues[2]))) = 0x5f3759df - (*((int *)(&(fltValues[2]))) >> 1);

        return 0.5f * (fltValues[1] + (fltValues[0] * fltValues[2]));
    }

//    inline float
//    Max(const float)
//    {
//        return 3.402823466e+38f;
//    }
//
//    inline float
//    Min(const float)
//    {
//        return -1.0f * 3.402823466e+38f;
//    }
//
//    inline int
//    Max(const int)
//    {
//        return ~((int)(1) << ((8 * sizeof(int)) - 1));
//    }
//
//    inline int
//    Min(const int)
//    {
//        return ((int)(1) << ((8 * sizeof(int)) - 1));
//    }
//
//    inline unsigned int
//    Max(const unsigned int)
//    {
//        return (unsigned int)(~0);
//    }
//
//    inline unsigned int
//    Min(const unsigned int)
//    {
//        return 0;
//    }
//
////    #pragma warning( push )
////    #pragma warning( disable:4310 )
//
//    inline char
//    Max(const char)
//    {
//        return (char)(~((char)(1) << ((8 * sizeof(char)) - 1)));
//    }
//
//    inline char
//    Min(const char)
//    {
//        return (char)(((char)(1) << ((8 * sizeof(char)) - 1)));
//    }
//
//    inline short
//    Max(const short)
//    {
//        return (short)(~((short)(1) << ((8 * sizeof(short)) - 1)));
//    }
//
//    inline short
//    Min(const short)
//    {
//        return (short)(((short)(1) << ((8 * sizeof(short)) - 1)));
//    }
//
////    #pragma warning( pop )
//
//    inline unsigned short
//    Max(const unsigned short)
//    {
//        return (unsigned short)(~0);
//    }
//
//    inline unsigned short
//    Min(const unsigned short)
//    {
//        return 0;
//    }

    inline int
    Mod(const int p_intDividend,
        const int p_intDivisor)
    {
        return p_intDividend % p_intDivisor;
    }

    inline unsigned int
    Mod(const unsigned int p_uintDividend,
        const unsigned int p_uintDivisor)
    {
        return p_uintDividend % p_uintDivisor;
    }

    inline float
    Mod(const float p_fltDividend,
        const float p_fltDivisor)
    {
        return fmodf(p_fltDividend, p_fltDivisor);
    }

    inline double
    Mod(const double p_dblDividend,
        const double p_dblDivisor)
    {
        return fmod(p_dblDividend, p_dblDivisor);
    }

    namespace DN
    {
        template <class T, unsigned int Dimensions> struct Point;
        template <class T, unsigned int Dimensions> struct Line;
        template <class T, unsigned int Dimensions> struct Triangle;
        template <class T, unsigned int Dimensions> struct Sphere;
        template <class T, unsigned int Dimensions> struct Cylinder;
        template <class T, unsigned int Dimensions> struct Box;
        template <class T, unsigned int Dimensions> struct ResultBox;
        template <class T, unsigned int Dimensions> struct Result;
        template <class T, unsigned int Dimensions> struct ResultBox;
        template <class T, unsigned int Dimensions> struct ResultLine;

        extern void * g_ptrWrapDistances;

        template <class T, unsigned int Dimensions>
        struct Point
        {
            T m_clsValues[Dimensions];

            Point() {}

            Point(const T & p_clsData)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] = p_clsData;
            }

            Point(const T (& p_clsData)[Dimensions])
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] = p_clsData[i];
            }

            inline operator T * () {return &(m_clsValues[0]);}
            inline operator const T * () const {return &(m_clsValues[0]);}

            inline Point<T, Dimensions>
            operator+(const Point<T, Dimensions> & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] + Value.m_clsValues[i];
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator+(const T & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] + Value;
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator-(const Point<T, Dimensions> & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] - Value.m_clsValues[i];
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator-(const T & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] - Value;
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator*(const Point<T, Dimensions> & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] * Value.m_clsValues[i];
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator*(const T & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] * Value;
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator*(const T (& Value)[Dimensions]) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] * Value[i];
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator/(const Point<T, Dimensions> & Value) const
            {
                Point<T, Dimensions> clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] / Value.m_clsValues[i];
                return clsPoint;
            }

            inline Point<T, Dimensions>
            operator/(const T & Value) const
            {
                Point clsPoint;
                for (unsigned int i=0; i<Dimensions; i++)
                    clsPoint.m_clsValues[i] = m_clsValues[i] / Value;
                return clsPoint;
            }

            inline Point<T, Dimensions> &
            operator=(const Point<T, Dimensions> & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] = Value.m_clsValues[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator=(const T & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] = Value;
                return *this;
            }

            inline Point<T, Dimensions> &
            operator=(const T (& Value)[Dimensions])
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] = Value[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator+=(const Point<T, Dimensions> & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] += Value.m_clsValues[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator+=(const T & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] += Value;
                return *this;
            }

            inline Point<T, Dimensions> &
            operator-=(const Point<T, Dimensions> & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] -= Value.m_clsValues[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator-=(const T Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] -= Value;
                return *this;
            }

            inline Point<T, Dimensions> &
            operator*=(const Point<T, Dimensions> Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] *= Value.m_clsValues[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator*=(const T Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] *= Value;
                return *this;
            }

            inline Point<T, Dimensions> &
            operator/=(const Point<T, Dimensions> Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] /= Value.m_clsValues[i];
                return *this;
            }

            inline Point<T, Dimensions> &
            operator/=(const T Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] /= Value;
                return *this;
            }

            inline bool
            operator==(const Point<T, Dimensions> & Value) const
            {
                for (unsigned int i=0; i<Dimensions; i++)
                {
                    if ( Equal(m_clsValues[i], Value.m_clsValues[i]) == false )
                        return false;
                }
                return true;
            }

            inline bool
            operator!=(const Point<T, Dimensions> & Value) const
            {
                return !(*this == Value);
            }

            inline float
            Distance(const Point<T, Dimensions> & Value,
                     const bool p_blnFast = false) const
            {
                return ( p_blnFast )
                    ? FastSqrt(DistanceSquared(Value))
                    : sqrtf(DistanceSquared(Value));
            }

            inline T
            DistanceSquared(const Point<T, Dimensions> & Value) const
            {
                T clsDistance = (0);
                for (unsigned int i=0; i<Dimensions; i++)
                {
                    T clsDelta = m_clsValues[i] - Value.m_clsValues[i];
                    clsDistance += clsDelta * clsDelta;
                }
                return clsDistance;
            }

            inline T
            Dot(const Point<T, Dimensions> & Value) const
            {
                T clsValue = (T)(0);
                for (unsigned int i=0; i<Dimensions; i++)
                    clsValue += m_clsValues[i] * Value.m_clsValues[i];
                return clsValue;
            }

            inline float
            LengthF(const bool Fast = false) const
            {
                float fltValue = 0.0f;
                for (unsigned int i=0; i<Dimensions; i++)
                    fltValue += (float)(m_clsValues[i] * m_clsValues[i]);
                return ( Fast ) ? FastSqrt(fltValue) : sqrtf(fltValue);
            }

            inline double
            LengthD() const
            {
                double dblValue = 0.0;
                for (unsigned int i=0; i<Dimensions; i++)
                    dblValue += (double)(m_clsValues[i] * m_clsValues[i]);
                return sqrt(dblValue);
            }

            inline void
            NormalizeF(const bool Fast = false)
            {
                T clsLengthReciprocal = (T)(1.0f / LengthF(Fast));
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] *= clsValue;
            }

            inline void
            NormalizeD()
            {
                T clsLengthReciprocal = (T)(1.0 / LengthD());
                for (unsigned int i=0; i<Dimensions; i++)
                    m_clsValues[i] *= clsValue;
            }

            T    DistanceSquaredW(const Point<T, Dimensions> & Value) const;
            bool Intersect(const Point<T, Dimensions> & Value) const;
        };

        template <class T, unsigned int Dimensions>
        struct Line
        {
            Point<T, Dimensions> m_stuPoint;
            Point<T, Dimensions> m_stuDirection;

            Line() {}

            Line(const Point<T, Dimensions> & ValuePoint,
                 const Point<T, Dimensions> & ValueDirection)
                : m_stuPoint(ValuePoint),
                  m_stuDirection(ValueDirection) {}

            bool Intersect(const Point<T, Dimensions> & Point, float * Time = NULL) const;
            bool Intersect(const Line<T, Dimensions> & Line, ResultLine<T, Dimensions> * ResultRef = NULL) const;
        };

        template <class T, unsigned int Dimensions>
        struct Triangle
        {
            Point<T, Dimensions> m_stuPoints[3];
        };

        template <class T, unsigned int Dimensions>
        struct Sphere
        {
            Point<T, Dimensions> m_stuPoint;
            T                    m_clsRadius;
        };

        template <class T, unsigned int Dimensions>
        struct Cylinder
        {
            Line<T, Dimensions> m_stuOrientation;
            T                   m_clsHeight;
            T                   m_clsRadius;
        };

        template <class T, unsigned int Dimensions>
        struct Box
        {
            static void *       m_ptrWrapValues;
            static void *       m_ptrWrapDistances;
            static unsigned int m_uintWrapMask;

            T m_clsValues[2 * Dimensions];

            Box() {}
            Box(const T (& Value)[2 * Dimensions])
            {
                for (unsigned int i=0; i<(2 * Dimensions); i++)
                    m_clsValues[i] = Value[i];
            }

            inline operator T * () {return &(m_clsValues[0]);}
            inline operator const T * () const {return &(m_clsValues[0]);}

            inline Box<T, Dimensions> &
            operator+=(const Point<T, Dimensions> & Value)
            {
                for (unsigned int i=0; i<Dimensions; i++)
                {
                    if ( m_clsValues[2 * i] > Value[i] )
                        m_clsValues[2 * i] = Value[i];
                    else if ( m_clsValues[(2 * i) + 1] > Value[i] )
                        m_clsValues[(2 * i) + 1] = Value[i];
                }
                return *this;
            }

            inline void Init()
            {
                for (unsigned int i=0; i<Dimensions; i++)
                {
                    m_clsValues[2 * i] = Max((T)(0));
                    m_clsValues[(2 * i) + 1] = Min((T)(0));
                }
            }

            void          Add(const Point<T, Dimensions> & Point);
            void          Add(const Box<T, Dimensions> & Box);
            void          AdjustW();
            T             Distance(const Point<T, Dimensions> & Point) const;
            T             DistanceW(const Point<T, Dimensions> & Point) const;
            bool          Intersect(const Point<T, Dimensions> & Point) const;
            bool          Intersect(const Box<T, Dimensions> & Box, ResultBox<T, Dimensions> * Result = NULL) const;
            bool          IntersectW(const Point<T, Dimensions> & Point, const unsigned int Mask = (unsigned int)(~0)) const;
            bool          IntersectW(const Box<T, Dimensions> & Box, ResultBox<T, Dimensions> * Result = NULL) const;
            unsigned long SectorW(const Point<T, Dimensions> & Point, Box<T, Dimensions> * Result = NULL) const;
        };

        struct ResultTime
        {
            float m_fltTime;
        };

        template <class T, unsigned int Dimensions>
        struct ResultBox
        {
            int                m_intType;
            Box<T, Dimensions> m_stuResult;
        };

        template <class T, unsigned int Dimensions>
        struct ResultLine
        {
            float                m_fltTime1;
            float                m_fltTime2;
            Point<T, Dimensions> m_stuPoint;
        };
    }

    namespace D2
    {
        struct Point
        {
            float m_fltX;
            float m_fltY;

            Point() {}

            Point(const float Value)
                : m_fltX(Value),
                  m_fltY(Value) {}

            Point(const float ValueX, const float ValueY)
                : m_fltX(ValueX),
                  m_fltY(ValueY) {}

            Point(const float * Value)
            {
                m_fltX = Value[0];
                m_fltY = Value[1];
            }

            inline operator float * () {return &m_fltX;}
            inline operator const float * () const {return (const float *)(&m_fltX);}
            inline Point  operator-() const {return Point(-m_fltX, -m_fltY);}
            inline Point operator+(const Point Value) const {return Point(m_fltX + Value.m_fltX, m_fltY + Value.m_fltY);}
            inline Point operator+(const float Value) const {return Point(m_fltX + Value, m_fltY + Value);}
            inline Point operator-(const Point Value) const {return Point(m_fltX - Value.m_fltX, m_fltY - Value.m_fltY);}
            inline Point operator-(const float Value) const {return Point(m_fltX - Value, m_fltY - Value);}
            inline Point operator*(const Point Value) const {return Point(m_fltX * Value.m_fltX, m_fltY * Value.m_fltY);}
            inline Point operator*(const float Value) const {return Point(m_fltX * Value, m_fltY * Value);}
            inline Point operator/(const Point Value) const {return Point(m_fltX / Value.m_fltX, m_fltY / Value.m_fltY);}
            inline Point operator/(const float Value) const {return Point(m_fltX / Value, m_fltY / Value);}
            inline Point operator=(const float Value) {m_fltX = Value; m_fltY = Value; return (*this);}
            inline Point & operator+=(const Point Value) {m_fltX += Value.m_fltX; m_fltY += Value.m_fltY; return (*this);}
            inline Point & operator+=(const float Value) {m_fltX += Value; m_fltY += Value; return (*this);}
            inline Point & operator-=(const Point Value) {m_fltX -= Value.m_fltX; m_fltY -= Value.m_fltY; return (*this);}
            inline Point & operator-=(const float Value) {m_fltX -= Value; m_fltY -= Value; return (*this);}
            inline Point & operator*=(const Point Value) {m_fltX *= Value.m_fltX; m_fltY *= Value.m_fltY; return (*this);}
            inline Point & operator*=(const float Value) {m_fltX *= Value; m_fltY *= Value; return (*this);}
            inline Point & operator/=(const Point Value) {m_fltX /= Value.m_fltX; m_fltY /= Value.m_fltY; return (*this);}
            inline Point & operator/=(const float Value) {m_fltX /= Value; m_fltY /= Value; return (*this);}
            inline bool  operator==(const Point Value) const {return ( Equal(m_fltX, Value.m_fltX) && Equal(m_fltY, Value.m_fltY) );}
            inline float Cross(const Point & Value) const {return (m_fltX * Value.m_fltY) - (m_fltY * Value.m_fltX);}
            inline float Dot(const Point & Value) const {return (m_fltX * Value.m_fltX) + (m_fltY * Value.m_fltY);}
            inline float LengthSquared() {return (m_fltX * m_fltX) + (m_fltY * m_fltY);}
            inline float Length(const bool Fast = false) {return ( Fast ) ? FastSqrt(LengthSquared()) : sqrtf(LengthSquared());}

            template <class T>
            inline Point &
            operator=(const DN::Point<T, 2> & Value)
            {
                m_fltX = (float)(Value.m_clsValues[0]);
                m_fltY = (float)(Value.m_clsValues[1]);
                return (*this);
            }
        };

        struct Triangle
        {
            Point m_stuPoints[3];

            bool Intersect(const Point & Value) const
            {
                return ( Value.Cross(m_stuPoints[1] - m_stuPoints[0]) <= 0.0f &&
                         Value.Cross(m_stuPoints[2] - m_stuPoints[1]) <= 0.0f &&
                         Value.Cross(m_stuPoints[0] - m_stuPoints[2]) <= 0.0f );
            }
        };

        struct Box
        {
            float m_fltValues[4];

            Box() {}

            Box(const Point & Value)
            {
                m_fltValues[0] = m_fltValues[1] = Value.m_fltX;
                m_fltValues[2] = m_fltValues[3] = Value.m_fltY;
            }

            Box(const Box & Box1,
                const Box & Box2)
            {
                m_fltValues[0] = ( Box1.m_fltValues[0] < Box2.m_fltValues[0] ) ? Box1.m_fltValues[0] : Box2.m_fltValues[0];
                m_fltValues[1] = ( Box1.m_fltValues[1] > Box2.m_fltValues[1] ) ? Box1.m_fltValues[1] : Box2.m_fltValues[1];
                m_fltValues[2] = ( Box1.m_fltValues[2] < Box2.m_fltValues[2] ) ? Box1.m_fltValues[2] : Box2.m_fltValues[2];
                m_fltValues[3] = ( Box1.m_fltValues[3] > Box2.m_fltValues[3] ) ? Box1.m_fltValues[3] : Box2.m_fltValues[3];
            }

            Box(const float MinX,
                const float MaxX,
                const float MinY,
                const float MaxY)
            {
                m_fltValues[0] = MinX;
                m_fltValues[1] = MaxX;
                m_fltValues[2] = MinY;
                m_fltValues[3] = MaxY;
            }

            void Init();

            inline operator float * () {return &(m_fltValues[0]);}
            inline operator const float * () const {return &(m_fltValues[0]);}
            inline operator const MATH::DN::Box<float, 2> * () const {return (const DN::Box<float, 2> *)(this);}

            template <class T>
            inline void
            operator=(const DN::Box<T, 2> & Value)
            {
                m_fltValues[0] = (float)(Value[0]);
                m_fltValues[1] = (float)(Value[1]);
                m_fltValues[2] = (float)(Value[2]);
                m_fltValues[3] = (float)(Value[3]);
            }

            inline void operator+=(const Point & Value) {Add(Value);}
            inline bool operator==(const Box & Value) const
            {
                return ( Equal(m_fltValues[0], Value.m_fltValues[0]) &&
                         Equal(m_fltValues[1], Value.m_fltValues[1]) &&
                         Equal(m_fltValues[2], Value.m_fltValues[2]) &&
                         Equal(m_fltValues[3], Value.m_fltValues[3]) );
            }

            inline float &
            operator[](const int Index)
            {
                ASSERT(Index >= 0 && Index <= 3);
                return m_fltValues[Index];
            }

            inline const float &
            operator[](const int Index) const
            {
                ASSERT(Index >= 0 && Index <= 3);
                return m_fltValues[Index];
            }

            inline void
            Add(const Point & Value)
            {
                if ( m_fltValues[0] > Value.m_fltX )
                    m_fltValues[0] = Value.m_fltX;
                if ( m_fltValues[1] < Value.m_fltX )
                    m_fltValues[1] = Value.m_fltX;
                if ( m_fltValues[2] > Value.m_fltY )
                    m_fltValues[2] = Value.m_fltY;
                if ( m_fltValues[3] < Value.m_fltY )
                    m_fltValues[3] = Value.m_fltY;
            }

            float Height() const {return m_fltValues[3] - m_fltValues[2];}
            int   Intersect(const Point & Start, const Point & End, Point (& Result)[2]) const;
            float Width() const {return m_fltValues[1] - m_fltValues[0];}

            inline bool
            Intersect(const Point & Value) const
            {
                return (Value.m_fltX >= m_fltValues[0] && Value.m_fltX <= m_fltValues[1] &&
                        Value.m_fltY >= m_fltValues[2] && Value.m_fltY <= m_fltValues[3]);
            }
        };

        struct BoxW : public Box
        {
            static float m_fltWrapDistances[2];
            static Box   m_stuWrap;

            static inline float
            WrapValue(const float Value)
            {
                return m_stuWrap[0] + fmodf(Value - m_stuWrap[0], m_fltWrapDistances[0]);
            }

            static inline float
            ClampValue(const float Value)
            {
                if ( Value > m_stuWrap[2] )
                    return ( Value < m_stuWrap[3] ) ? Value : m_stuWrap[3];
                else
                    return m_stuWrap[2];
            }

            BoxW()
                : Box(m_stuWrap[0],
                      m_stuWrap[0],
                      m_stuWrap[1],
                      m_stuWrap[1]) {}

            BoxW(const Point & Value)
            {
                m_fltValues[0] = m_fltValues[1] = WrapValue(Value.m_fltX);
                m_fltValues[2] = m_fltValues[3] = ClampValue(Value.m_fltY);
            }

            BoxW(const Box & Value)
            {
                m_fltValues[0] = WrapValue(Value[0]);
                m_fltValues[1] = WrapValue(Value[1]);
                m_fltValues[2] = ClampValue(Value[2]);
                m_fltValues[3] = ClampValue(Value[3]);
            }

            BoxW(const BoxW & Box1,
                 const BoxW & Box2)
            {
                *this = Box1;
                Add(Box2);
            }

            BoxW(const float MinX,
                 const float MaxX,
                 const float MinY,
                 const float MaxY)
            {
                m_fltValues[0] = WrapValue(MinX);
                m_fltValues[1] = WrapValue(MaxX);
                m_fltValues[2] = ClampValue(MinY);
                float fltTemp = ClampValue(MaxY);
                if ( m_fltValues[2] < fltTemp )
                    m_fltValues[3] = fltTemp;
                else
                {
                    m_fltValues[3] = m_fltValues[2];
                    m_fltValues[2] = fltTemp;
                }
            }

            void  Add(const Point & Point);
            void  Add(const BoxW & Box);
            void  Adjust();
            bool  Intersect(const Point & Point) const;
            int   Intersect(const BoxW & Box, BoxW * Results = NULL) const;
            BoxW  Sector(const unsigned char Sector) const;
            float Width() const;
        };

        union Matrix
        {
            struct
            {
                float m_flt11;
                float m_flt12;
                float m_flt13;
                float m_flt21;
                float m_flt22;
                float m_flt23;
                float m_flt31;
                float m_flt32;
                float m_flt33;
            };
            struct
            {
                Point m_stuBasisX;
                float m_fltExtraX;
                Point m_stuBasisY;
                float m_fltExtraY;
                Point m_stuTranslation;
                float m_fltExtra;
            };
            float m_fltValues[9];

            void
            Identity()
            {
                m_flt11 = m_flt22 = m_flt33 = 1.0f;
                m_flt12 = m_flt13 = m_flt21 = m_flt23 = m_flt31 = m_flt32 = 0.0f;
            }

            Point Multiply(const Point & Value) const
            {
                return Point((m_flt11 * Value.m_fltX) + (m_flt21 * Value.m_fltY) + m_flt31,
                             (m_flt12 * Value.m_fltX) + (m_flt22 * Value.m_fltY) + m_flt32);
            }

            Point operator*(const Point & Value) const {return Multiply(Value);}
        };
    }

    namespace D3
    {
        struct Point;
        struct Line;
        struct Plane;
        struct Triangle;
        struct Sphere;
        struct Cylinder;
        union  Box;
        struct Frustum;
        union  Quaternion;
        union  Matrix;

        union Result;

        struct Point
        {
            float m_fltX;
            float m_fltY;
            float m_fltZ;

            Point() {}

            Point(const float Value)
                : m_fltX(Value),
                  m_fltY(Value),
                  m_fltZ(Value) {}

            Point(const float X,
                  const float Y,
                  const float Z)
                : m_fltX(X),
                  m_fltY(Y),
                  m_fltZ(Z) {}

            Point(const float * Values)
            {
                m_fltX = Values[0];
                m_fltY = Values[1];
                m_fltZ = Values[2];
            }

            inline operator float * () {return &m_fltX;}
            inline operator const float * () const {return (const float *)(&m_fltX);}

            inline Point
            operator-() const
            {
                return Point(-m_fltX, -m_fltY, -m_fltZ);
            }

            inline Point
            operator+(const Point & Value) const
            {
                return Point(m_fltX + Value.m_fltX,
                             m_fltY + Value.m_fltY,
                             m_fltZ + Value.m_fltZ);
            }

            inline Point
            operator+(const float Value) const
            {
                return Point(m_fltX + Value,
                             m_fltY + Value,
                             m_fltZ + Value);
            }

            inline Point
            operator-(const Point & Value) const
            {
                return Point(m_fltX - Value.m_fltX,
                             m_fltY - Value.m_fltY,
                             m_fltZ - Value.m_fltZ);
            }

            inline Point
            operator-(const float Value) const
            {
                return Point(m_fltX - Value,
                             m_fltY - Value,
                             m_fltZ - Value);
            }

            inline Point
            operator*(const Point & Value) const
            {
                return Point(m_fltX * Value.m_fltX,
                             m_fltY * Value.m_fltY,
                             m_fltZ * Value.m_fltZ);
            }

            inline Point
            operator*(const float Value) const
            {
                return Point(m_fltX * Value,
                             m_fltY * Value,
                             m_fltZ * Value);
            }

            inline Point
            operator/(const Point & Value) const
            {
                return Point(m_fltX / Value.m_fltX,
                             m_fltY / Value.m_fltY,
                             m_fltZ / Value.m_fltZ);
            }

            inline Point
            operator/(const float Value) const
            {
                return Point(m_fltX / Value,
                             m_fltY / Value,
                             m_fltZ / Value);
            }

            inline Point &
            operator=(const float Value)
            {
                m_fltX = Value;
                m_fltY = Value;
                m_fltZ = Value;
                return (*this);
            }


            template <class T>
            inline Point &
            operator=(const MATH::DN::Point<T, 3> & Value)
            {
                m_fltX = (float)(Value[0]);
                m_fltY = (float)(Value[1]);
                m_fltZ = (float)(Value[2]);
                return (*this);
            }

            inline Point &
            operator+=(const Point & Value)
            {
                m_fltX += Value.m_fltX;
                m_fltY += Value.m_fltY;
                m_fltZ += Value.m_fltZ;
                return (*this);
            }

            inline Point &
            operator+=(const float Value)
            {
                m_fltX += Value;
                m_fltY += Value;
                m_fltZ += Value;
                return (*this);
            }

            inline Point &
            operator-=(const Point & Value)
            {
                m_fltX -= Value.m_fltX;
                m_fltY -= Value.m_fltY;
                m_fltZ -= Value.m_fltZ;
                return (*this);
            }

            inline Point &
            operator-=(const float Value)
            {
                m_fltX -= Value;
                m_fltY -= Value;
                m_fltZ -= Value;
                return (*this);
            }

            inline Point &
            operator*=(const Point & Value)
            {
                m_fltX *= Value.m_fltX;
                m_fltY *= Value.m_fltY;
                m_fltZ *= Value.m_fltZ;
                return (*this);
            }

            inline Point &
            operator*=(const float Value)
            {
                m_fltX *= Value;
                m_fltY *= Value;
                m_fltZ *= Value;
                return (*this);
            }

            inline Point &
            operator/=(const Point & Value)
            {
                m_fltX /= Value.m_fltX;
                m_fltY /= Value.m_fltY;
                m_fltZ /= Value.m_fltZ;
                return (*this);
            }

            inline Point &
            operator/=(const float Value)
            {
                m_fltX /= Value;
                m_fltY /= Value;
                m_fltZ /= Value;
                return (*this);
            }

            inline bool
            operator==(const Point & Value) const
            {
                return ( Equal(m_fltX, Value.m_fltX) &&
                         Equal(m_fltY, Value.m_fltY) &&
                         Equal(m_fltZ, Value.m_fltZ) );
            }

            inline bool
            operator==(const float Value) const
            {
                return ( Equal(m_fltX, Value) &&
                         Equal(m_fltY, Value) &&
                         Equal(m_fltZ, Value) );
            }

            inline bool
            operator!=(const Point & Value) const
            {
                return !((*this) == Value);
            }

            inline bool
            operator!=(const float Value) const
            {
                return !((*this) == Value);
            }

            inline Point
            Cross(const Point & Value) const
            {
                return Point((m_fltY * Value.m_fltZ) - (m_fltZ * Value.m_fltY),
                             (m_fltZ * Value.m_fltX) - (m_fltX * Value.m_fltZ),
                             (m_fltX * Value.m_fltY) - (m_fltY * Value.m_fltX));
            }

            inline Point
            Cross(const Point & Value1,
                  const Point & Value2) const
            {
                return (Value1 - *this).Cross(Value2 - *this);
            }

            inline float Distance(const Point & Value, const bool Fast = false) const {return ((*this) - Value).Length(Fast);}
            inline float DistanceSquared(const Point & Value) const {return ((*this) - Value).LengthSquared();}
            inline float Dot(const Point & Value) const {return (m_fltX * Value.m_fltX) + (m_fltY * Value.m_fltY) + (m_fltZ * Value.m_fltZ);}
            inline float DotSelf() const {return (m_fltX * m_fltX) + (m_fltY * m_fltY) + (m_fltZ * m_fltZ);}
            inline Point Inverse() const {return -(*this);}
            inline float Length(const bool Fast = false) const {return ( Fast ) ? FastSqrt(LengthSquared()) : sqrtf(LengthSquared());}
            inline float LengthSquared() const {return (m_fltX * m_fltX) + (m_fltY * m_fltY) + (m_fltZ * m_fltZ);}

            inline Point
            Normalize(const bool Fast = false) const
            {
                float fltLengthRecip = 1.0f / Length(Fast);
                return Point(m_fltX * fltLengthRecip,
                             m_fltY * fltLengthRecip,
                             m_fltZ * fltLengthRecip);
            }

            inline Point
            Swizzle(const unsigned char p_uchrIndex1,
                    const unsigned char p_uchrIndex2,
                    const unsigned char p_uchrIndex3) const
            {
                const float * ptrValues = &m_fltX;
                return Point(ptrValues[p_uchrIndex1],
                             ptrValues[p_uchrIndex2],
                             ptrValues[p_uchrIndex3]);
            }
        };

        union Result
        {
            float m_fltTime;
            struct
            {
                float m_fltTime1;
                float m_fltTime2;
                Point m_stuPoint;
            };
            struct
            {
                int   m_intType;
                float m_fltBox[6];
            } Box;
        };

        struct Line
        {
            Point m_stuPoint;
            Point m_stuDirection;

            Line() {}

            Line(const Point & Origin,
                 const Point & Direction)
                : m_stuPoint(Origin),
                  m_stuDirection(Direction) {}

            inline float
            DistanceSquared(const Point & Value) const
            {
                Point stuTemp = Value - m_stuPoint;
                stuTemp = m_stuPoint + (m_stuDirection * (m_stuDirection.Dot(stuTemp) / m_stuDirection.Dot(m_stuDirection)));
                return Value.DistanceSquared(stuTemp);
            }

            inline Point Extend(const float p_fltAmount) const {return m_stuPoint + (m_stuDirection * p_fltAmount);}
            bool         Intersect(const Point * PointRef, Result * ResultRef) const;
            bool         Intersect(const Line & LineRef, Result * ResultRef) const;
        };

        struct Plane
        {
            Point m_stuNormal;
            float m_fltConstant;

            float
            Distance(const Point & Value) const
            {
                return m_stuNormal.Dot(Value) - m_fltConstant;
            }

            bool Intersect(const Line & Value, Result * Result) const;

            /*
            Returns:

            -1: The point lies below the plane
            0: The point lies on the plane
            1: The point lies above the plane
            */
            int
            Intersect(const Point & Value) const
            {
                float fltDistance = Distance(Value);
                if ( Zero(fltDistance) )
                    return 0;
                else if ( fltDistance < 0.0f )
                    return -1;
                else
                    return 1;
            }

            void
            Normalize(const bool Fast = false)
            {
                float fltLengthRecip = 1.0f / m_stuNormal.Length(Fast);
                m_stuNormal *= fltLengthRecip;
                m_fltConstant *= fltLengthRecip;
            }

            void
            Set(const Point & Value1,
                const Point & Value2,
                const Point & Value3)
            {
                m_stuNormal = Value1.Cross(Value2, Value3);
                m_fltConstant = m_stuNormal.Dot(Value1);
            }

            void
            Set(const Point * Values,
                const unsigned int * Indices)
            {
                Set(Values[Indices[0]],
                    Values[Indices[1]],
                    Values[Indices[2]]);
            }
        };

        struct Triangle
        {
            Point m_stuPoints[3];

            bool Intersect(const Line & Line, Result * Result = NULL) const;
        };

        struct Sphere
        {
            Point m_stuPosition;
            float m_fltRadius;

            bool Intersect(const Point & p_clsPoint) const {return ( m_stuPosition.DistanceSquared(p_clsPoint) < (m_fltRadius * m_fltRadius) );}
        };

        struct Cylinder
        {
            Line  m_stuOrientation;
            float m_fltHeight;
            float m_fltRadius;
        };

        union Quaternion
        {
            float m_fltValues[4];
            struct
            {
                float m_fltW;
                float m_fltX;
                float m_fltY;
                float m_fltZ;
            };

            Quaternion() {}
            Quaternion(const Point & EulerAngles);
            Quaternion(const Point & Axis, const float Angle);
            Quaternion(const float EulerAngleX, const float EulerAngleY, const float EulerAngleZ) {*this = Quaternion(Point(EulerAngleX, EulerAngleY, EulerAngleZ));}
            Quaternion(const float AxisX, const float AxisY, const float AxisZ, const float Angle) {*this = Quaternion(Point(AxisX, AxisY, AxisZ), Angle);}

            inline void
            Identity()
            {
                m_fltW = 1.0f;
                m_fltX = m_fltY = m_fltZ = 0.0f;
            }

            inline Quaternion
            Inverse() const
            {
                Quaternion stuReturn;
                stuReturn.m_fltW = m_fltW;
                stuReturn.m_fltX = -m_fltX;
                stuReturn.m_fltY = -m_fltY;
                stuReturn.m_fltZ = -m_fltZ;
                return stuReturn;
            }

            inline Quaternion
            Normalize(const bool Fast = false) const
            {
                Quaternion stuReturn;
                float fltLength = ( Fast )
                    ? FastSqrt((m_fltW * m_fltW) + (m_fltX * m_fltX) + (m_fltY * m_fltY) + (m_fltZ * m_fltZ))
                    : sqrtf((m_fltW * m_fltW) + (m_fltX * m_fltX) + (m_fltY * m_fltY) + (m_fltZ * m_fltZ));
                float fltLengthRecip = 1.0f / fltLength;
                stuReturn.m_fltW = m_fltW * fltLengthRecip;
                stuReturn.m_fltX = m_fltX * fltLengthRecip;
                stuReturn.m_fltY = m_fltY * fltLengthRecip;
                stuReturn.m_fltZ = m_fltZ * fltLengthRecip;

                return stuReturn;
            }

            Quaternion Multiply(const Quaternion & Value) const;
            Point      Rotate(const Point & Value) const;
            Quaternion Slerp(const Quaternion & Value, const float Amount) const;
            Point      ToAxisAngle(float & Value) const;
            Point      ToAxisAngleFast(float & Value) const;
            Matrix     ToMatrix() const;

            inline Quaternion operator*(const Quaternion & Value) const {return Multiply(Value);}
        };

        union Box
        {
            struct
            {
                float m_fltMinX;
                float m_fltMaxX;
                float m_fltMinY;
                float m_fltMaxY;
                float m_fltMinZ;
                float m_fltMaxZ;
            };
            float m_fltValues[6];

            Box() {}

            Box(const Point & Value)
            {
                m_fltValues[0] = m_fltValues[1] = Value.m_fltX;
                m_fltValues[2] = m_fltValues[3] = Value.m_fltY;
                m_fltValues[4] = m_fltValues[5] = Value.m_fltZ;
            }

            Box(const Box & Value1,
                const Box & Value2)
            {
                m_fltValues[0] = ( Value1[0] < Value2[0] ) ? Value1[0] : Value2[0];
                m_fltValues[1] = ( Value1[1] > Value2[1] ) ? Value1[1] : Value2[1];
                m_fltValues[2] = ( Value1[2] < Value2[2] ) ? Value1[2] : Value2[2];
                m_fltValues[3] = ( Value1[3] > Value2[3] ) ? Value1[3] : Value2[3];
                m_fltValues[4] = ( Value1[4] < Value2[4] ) ? Value1[4] : Value2[4];
                m_fltValues[5] = ( Value1[5] > Value2[5] ) ? Value1[5] : Value2[5];
            }

            Box(const float MinX,
                const float MaxX,
                const float MinY,
                const float MaxY,
                const float MinZ,
                const float MaxZ)
            {
                m_fltValues[0] = MinX;
                m_fltValues[1] = MaxX;
                m_fltValues[2] = MinY;
                m_fltValues[3] = MaxY;
                m_fltValues[4] = MinZ;
                m_fltValues[5] = MaxZ;
            }

            void Init();

            inline operator float * () {return &(m_fltValues[0]);}
            inline operator const float * () const {return &(m_fltValues[0]);}
            inline operator const MATH::DN::Box<float, 3> * () const {return (const DN::Box<float, 3> *)(this);}

            template <class T>
            inline void
            operator=(const MATH::DN::Box<T, 3> & Value)
            {
                m_fltValues[0] = (float)(Value[0]);
                m_fltValues[1] = (float)(Value[1]);
                m_fltValues[2] = (float)(Value[2]);
                m_fltValues[3] = (float)(Value[3]);
                m_fltValues[4] = (float)(Value[4]);
                m_fltValues[5] = (float)(Value[5]);
            }

            inline void operator+=(const Point & Value) {Add(Value);}

            inline float &
            operator[](const int Index)
            {
                #ifdef UTIL_DEBUG
                    if ( Index < 0 || Index > 5 )
                        throw;
                #endif

                return m_fltValues[Index];
            }

            inline const float &
            operator[](const int Index) const
            {
                #ifdef UTIL_DEBUG
                    if ( Index < 0 || Index > 5 )
                        throw;
                #endif

                return m_fltValues[Index];
            }

            inline void
            Add(const Point & Value)
            {
                if ( m_fltValues[0] > Value.m_fltX )
                    m_fltValues[0] = Value.m_fltX;
                if ( m_fltValues[1] < Value.m_fltX )
                    m_fltValues[1] = Value.m_fltX;
                if ( m_fltValues[2] > Value.m_fltY )
                    m_fltValues[2] = Value.m_fltY;
                if ( m_fltValues[3] < Value.m_fltY )
                    m_fltValues[3] = Value.m_fltY;
                if ( m_fltValues[4] > Value.m_fltZ )
                    m_fltValues[4] = Value.m_fltZ;
                if ( m_fltValues[5] < Value.m_fltZ )
                    m_fltValues[5] = Value.m_fltZ;
            }

            int Intersect(const Point & Value) const;
            int Intersect(const Triangle & Value) const;
            int Intersect(const Sphere & Value) const;
            int Intersect(const Box & Value) const;
        };

        struct Frustum
        {
            Plane m_stuRight;
            Plane m_stuLeft;
            Plane m_stuBottom;
            Plane m_stuTop;
            Plane m_stuFar;
            Plane m_stuNear;

            int Intersect(const Box & Value) const;
            int Intersect(const Point & Value) const;
        };

        union Matrix
        {
            struct
            {
                float m_flt11;
                float m_flt12;
                float m_flt13;
                float m_flt14;
                float m_flt21;
                float m_flt22;
                float m_flt23;
                float m_flt24;
                float m_flt31;
                float m_flt32;
                float m_flt33;
                float m_flt34;
                float m_flt41;
                float m_flt42;
                float m_flt43;
                float m_flt44;
            };
            struct
            {
                Point m_stuBasisX;
                float m_fltExtraX;
                Point m_stuBasisY;
                float m_fltExtraY;
                Point m_stuBasisZ;
                float m_fltExtraZ;
                Point m_stuTranslation;
                float m_fltExtra;
            };
            float m_fltValues[16];
            float m_fltRows[4][4];

            enum Order
            {
                ORDER_XYZ,
                ORDER_XZY,
                ORDER_YXZ,
                ORDER_YZX,
                ORDER_ZXY,
                ORDER_ZYX,

                ORDER_SRT,
                ORDER_STR,
                ORDER_RST,
                ORDER_RTS,
                ORDER_TSR,
                ORDER_TRS,
            };

        private:

            inline Matrix
            MultiplyPrivate(const Matrix & Operand) const
            {
                Matrix stuTemp;
                for (unsigned char i=0; i<4; i++)
                {
                    for (unsigned char j=0; j<4; j++)
                    {
                        stuTemp.m_fltRows[i][j] = (m_fltRows[i][0] * Operand.m_fltRows[0][j]) +
                                                  (m_fltRows[i][1] * Operand.m_fltRows[1][j]) +
                                                  (m_fltRows[i][2] * Operand.m_fltRows[2][j]) +
                                                  (m_fltRows[i][3] * Operand.m_fltRows[3][j]);
                    }
                }
                return stuTemp;
            }

            static inline void
            AxisAngleValues(const Point & Axis, const float Angle, float (& Result)[3][3])
            {
                float fltCos = cosf(Angle);
                float fltSin = sinf(Angle);
                float fltCosX = 1.0f - fltCos;

                Result[0][0] = fltCos + (Axis[0] * Axis[0] * fltCosX);
                Result[1][1] = fltCos + (Axis[1] * Axis[1] * fltCosX);
                Result[2][2] = fltCos + (Axis[2] * Axis[2] * fltCosX);

                float fltTemp1 = Axis[0] * Axis[1] * fltCosX;
                float fltTemp2 = Axis[2] * fltSin;
                Result[1][0] = fltTemp1 - fltTemp2;
                Result[0][1] = fltTemp1 + fltTemp2;

                fltTemp1 = Axis[0] * Axis[2] * fltCosX;
                fltTemp2 = Axis[1] * fltSin;
                Result[2][0] = fltTemp1 + fltTemp2;
                Result[0][2] = fltTemp1 - fltTemp2;

                fltTemp1 = Axis[1] * Axis[2] * fltCosX;
                fltTemp2 = Axis[0] * fltSin;
                Result[2][1] = fltTemp1 - fltTemp2;
                Result[1][2] = fltTemp1 + fltTemp2;
            }

            static inline void
            ValuesRotation(const float X,
                           const float Y,
                           const float Z,
                           const Matrix::Order RotationOrder,
                           float * Result)
            {
                float fltSinX = sinf(X);
                float fltCosX = cosf(X);
                float fltSinY = sinf(Y);
                float fltCosY = cosf(Y);
                float fltSinZ = sinf(Z);
                float fltCosZ = cosf(Z);
                float fltOpt1, fltOpt2;

                switch ( RotationOrder )
                {
                    case Matrix::ORDER_ZYX:

                        fltOpt1 = fltCosZ * fltSinY;
                        fltOpt2 = fltSinZ * fltSinY;
                        Result[0] = fltCosZ * fltCosY;
                        Result[1] = (fltOpt1 * fltSinX) - (fltSinZ * fltCosX);
                        Result[2] = (fltOpt1 * fltCosX) + (fltSinZ * fltSinX);
                        Result[3] = fltSinZ * fltCosY;
                        Result[4] = (fltOpt2 * fltSinX) + (fltCosZ * fltCosX);
                        Result[5] = (fltOpt2 * fltCosX) - (fltCosZ * fltSinX);
                        Result[6] = -fltSinY;
                        Result[7] = fltCosY * fltSinX;
                        Result[8] = fltCosY * fltCosX;
                        return;

                    case Matrix::ORDER_XYZ:

                        fltOpt1 = fltSinX * fltSinY;
                        fltOpt2 = fltCosX * fltSinY;
                        Result[0] = fltCosY * fltCosZ;
                        Result[1] = -(fltCosY * fltSinZ);
                        Result[2] = fltSinY;
                        Result[3] = (fltCosX * fltSinZ) + (fltOpt1 * fltCosZ);
                        Result[4] = (fltCosX * fltCosZ) - (fltOpt1 * fltSinZ);
                        Result[5] = -(fltSinX * fltCosY);
                        Result[6] = (fltSinX * fltSinZ) - (fltOpt2 * fltCosZ);
                        Result[7] = (fltSinX * fltCosZ) + (fltOpt2 * fltSinZ);
                        Result[8] = fltCosX * fltCosY;
                        return;

                    default:
                        throw;
                }
            }

        public:

            Matrix() {}
            Matrix(const Quaternion &);
            Matrix(const Point & Axis, const float Angle);
            Matrix(const Point & Scale, const Point & Rotation, const Point & Translation, const Order MatrixOrder, const Order RotationOrder);

            Matrix(const float * Scale,
                   const float * Rotation,
                   const float * Translation,
                   const Order MatrixOrder,
                   const Order RotationOrder)
            {
                *this = Matrix(*((const Point *)(Scale)),
                               *((const Point *)(Rotation)),
                               *((const Point *)(Translation)),
                               MatrixOrder,
                               RotationOrder);
            }

            Matrix(const float ScaleX,
                   const float ScaleY,
                   const float ScaleZ,
                   const float RotationX,
                   const float RotationY,
                   const float RotationZ,
                   const float TranslationX,
                   const float TranslationY,
                   const float TranslationZ,
                   const Order MatrixOrder,
                   const Order RotationOrder)
            {
                *this = Matrix(Point(ScaleX, ScaleY, ScaleZ),
                               Point(RotationX, RotationY, RotationZ),
                               Point(TranslationX, TranslationY, TranslationZ),
                               MatrixOrder,
                               RotationOrder);
            }

            void
            Identity()
            {
                m_flt11 = m_flt22 = m_flt33 = m_flt44 = 1.0f;
                m_flt12 = m_flt13 = m_flt14 = m_flt21 = m_flt23 = m_flt24 = m_flt31 = m_flt32 = m_flt34 = m_flt41 = m_flt42 = m_flt43 = 0.0f;
            }

            static Matrix IdentityMatrix() {Matrix stuReturn; stuReturn.Identity(); return stuReturn;}

            inline operator float * () {return &(m_fltValues[0]);}
            inline operator const float * () const {return (const float *)(&(m_fltValues[0]));}
            float      Determinant() const;
            Matrix     Inverse() const;
            Matrix     Multiply(const Matrix & Right) const {return (*this) * Right;}
            void       RotateX(const float p_fltAmount);
            void       RotateY(const float p_fltAmount);
            void       RotateZ(const float p_fltAmount);
            void       Rotate(const Point & Axis, const float Angle);
            void       Rotate(const Point & Point, const Order RotationOrder = ORDER_XYZ) {Rotate(Point.m_fltX, Point.m_fltY, Point.m_fltZ, RotationOrder);}
            void       Rotate(const float X, const float Y, const float Z, const Order RotationOrder = ORDER_XYZ);
            void       Scale(const float X, const float Y, const float Z);
            void       Scale(const Point & Amount) {Scale(Amount.m_fltX, Amount.m_fltY, Amount.m_fltZ);}
            Quaternion ToQuaternion() const;
            Quaternion ToQuaternionFast() const;
            void       Translate(const float X, const float Y, const float Z);
            void       Translate(const Point & Amount) {Translate(Amount.m_fltX, Amount.m_fltY, Amount.m_fltZ);}
            Matrix     Transpose() const;

            Point      operator*(const Point & Right) const {return Multiply(Right);}
            Matrix     operator*(const Matrix & Right) const {return MultiplyPrivate(Right);}
            Matrix     operator*=(const Matrix & Right) {*this = MultiplyPrivate(Right); return *this;}

            Point      Multiply(const Point & Right) const
            {
                return Point((m_flt11 * Right.m_fltX) + (m_flt21 * Right.m_fltY) + (m_flt31 * Right.m_fltZ) + m_flt41,
                             (m_flt12 * Right.m_fltX) + (m_flt22 * Right.m_fltY) + (m_flt32 * Right.m_fltZ) + m_flt42,
                             (m_flt13 * Right.m_fltX) + (m_flt23 * Right.m_fltY) + (m_flt33 * Right.m_fltZ) + m_flt43);
            }
        };
    }
}

template <class T, unsigned int Dimensions>
T
MATH::DN::Point<T, Dimensions>::
DistanceSquaredW(const MATH::DN::Point<T, Dimensions> & p_stuPoint) const
{
    T clsDistance = (T)(0);
    MATH::DN::Box<T, Dimensions> * ptrWrapValues = (MATH::DN::Box<T, Dimensions> *)(g_ptrWrapValues);

    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( MATH::DN::g_uintWrapMask & (1 << i) )
        {
            unsigned int uintMin = i << 1;
            unsigned int uintMax = uintMin + 1;

            T clsDistanceWrap, clsDistanceNoWrap;
            if ( m_clsValues[i] > p_stuPoint[i] )
            {
                clsDistanceWrap = (p_stuPoint[i] - ptrWrapValues->m_clsValues[uintMin]) + (ptrWrapValues->m_clsValues[uintMax] - m_clsValues[i]);
                clsDistanceNoWrap = m_clsValues[i] - p_stuPoint[i];
            }
            else
            {
                clsDistanceWrap = (m_clsValues[i] - ptrWrapValues->m_clsValues[uintMin]) + (ptrWrapValues->m_clsValues[uintMax] - p_stuPoint[i]);
                clsDistanceNoWrap = p_stuPoint[i] - m_clsValues[i];
            }

            clsDistance += ( clsDistanceNoWrap < clsDistanceWrap )
                ? (clsDistanceNoWrap * clsDistanceNoWrap)
                : (clsDistanceWrap * clsDistanceWrap);
        }
        else
            clsDistance += (m_clsValues[i] - p_stuPoint[i]) * (m_clsValues[i] - p_stuPoint[i]);
    }

    return clsDistance;
}

template <class T, unsigned int Dimensions>
bool
MATH::DN::Point<T, Dimensions>::
Intersect(const MATH::DN::Point<T, Dimensions> & p_stuPoint) const
{
    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( Equal(m_clsValues[i], p_stuPoint[i]) == false )
            return false;
    }

    return true;
}

template <class T, unsigned int Dimensions>
bool
MATH::DN::Line<T, Dimensions>::
Intersect(const MATH::DN::Point<T, Dimensions> & p_stuPoint,
          float * p_ptrTime) const
{
    float fltValues[2];
    if ( m_stuDirection.m_fltX != 0.0f )
    {
        if ( p_ptrTime )
            *p_ptrTime = (p_ptrPoint->m_fltX - m_stuPoint.m_fltX) / m_stuDirection.m_fltX;
        fltValues[0] = (m_stuPoint.m_fltY + (p_ptrResult->m_fltTime * m_stuDirection.m_fltY)) - p_ptrPoint->m_fltY;
        fltValues[1] = (m_stuPoint.m_fltZ + (p_ptrResult->m_fltTime * m_stuDirection.m_fltZ)) - p_ptrPoint->m_fltZ;
    }
    else if ( m_stuDirection.m_fltY != 0.0f )
    {
        if ( p_ptrTime )
            *p_ptrTime = (p_ptrPoint->m_fltY - m_stuPoint.m_fltY) / m_stuDirection.m_fltY;
        fltValues[0] = (m_stuPoint.m_fltX + (p_ptrResult->m_fltTime * m_stuDirection.m_fltX)) - p_ptrPoint->m_fltX;
        fltValues[1] = (m_stuPoint.m_fltZ + (p_ptrResult->m_fltTime * m_stuDirection.m_fltZ)) - p_ptrPoint->m_fltZ;
    }
    else if ( m_stuDirection.m_fltZ != 0.0f )
    {
        if ( p_ptrTime )
            *p_ptrTime = (p_ptrPoint->m_fltZ - m_stuPoint.m_fltZ) / m_stuDirection.m_fltZ;
        fltValues[0] = (m_stuPoint.m_fltX + (p_ptrResult->m_fltTime * m_stuDirection.m_fltX)) - p_ptrPoint->m_fltX;
        fltValues[1] = (m_stuPoint.m_fltY + (p_ptrResult->m_fltTime * m_stuDirection.m_fltY)) - p_ptrPoint->m_fltY;
    }
    else
        return false;

    return ( Zero(fltValues[0]) && Zero(fltValues[1]) );
}

template <class T, unsigned int Dimensions>
bool
MATH::DN::Line<T, Dimensions>::
Intersect(const MATH::DN::Line<T, Dimensions> & p_stuLine,
          MATH::DN::ResultLine<T, Dimensions> * p_ptrResult) const
{
    for (unsigned int i=0; i<Dimensions; i++)
    {
        for (unsigned int j=0; j<Dimensions; j++)
        {
            if ( i == j )
                continue;

            T clsDenominator = (m_stuDirection[i] * p_ptrLine->m_stuDirection[j]) - (m_stuDirection[j] * p_ptrLine->m_stuDirection[i]);
            if ( clsDenominator != (T)(0) )
            {
                float fltNumerator = (float)((p_ptrLine->m_stuDirection[i] * (m_stuPoint[j] - p_ptrLine->m_stuPoint[j])) +
                                             (p_ptrLine->m_stuDirection[j] * (p_ptrLine->m_stuPoint[i] - m_stuPoint[i])));
                if ( p_ptrResult )
                    p_ptrResult->m_fltTime1 = fltNumerator / ((float)(clsDenominator));

                for (i=0; i<Dimensions; i++)
                {
                    if ( p_ptrLine->m_stuDirection[i] != (T)(0) )
                        break;
                }

                if ( p_ptrResult )
                    p_ptrResult->m_fltTime2 = ((float)(m_stuPoint[i] + (p_ptrResult->m_fltTime1 * m_stuDirection[i]) - p_ptrLine->m_stuPoint[i])) / ((float)(p_ptrLine->m_stuDirection[i]));
                Point<T, Dimensions> stuIntersection;
                for (i=0; i<Dimensions; i++)
                {
                    float fltValue1 = (float)(m_stuPoint[i]) + (p_ptrResult->m_fltTime1 * (float)(m_stuDirection[i]));
                    float fltValue2 = (float)(p_ptrLine->m_stuPoint[i]) + (p_ptrResult->m_fltTime2 * (float)(p_ptrLine->m_stuDirection[i]));
                    if ( Zero(fltValue1 - fltValue2) == false )
                        return false;
                    stuIntersection[i] = (T)(fltValue1);
                }
                if ( p_ptrResult )
                    p_ptrResult->m_stuPoint = stuIntersection;

                return true;
            }
        }
    }

    return false;
}

template <class T, unsigned int Dimensions>
void
MATH::DN::Box<T, Dimensions>::
Add(const MATH::DN::Point<T, Dimensions> & p_stuPoint)
{
    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( m_clsValues[2 * i] > p_stuPoint[i] )
            m_clsValues[2 * i] = p_stuPoint[i];
        if ( m_clsValues[(2 * i) + 1] < p_stuPoint[i] )
            m_clsValues[(2 * i) + 1] = p_stuPoint[i];
    }
}

template <class T, unsigned int Dimensions>
void
MATH::DN::Box<T, Dimensions>::
AdjustW()
{
    ASSERT(g_ptrWrapDistances != NULL && g_ptrWrapValues != NULL);

    T * ptrWrapDistances = (T *)(g_ptrWrapDistances);
    MATH::DN::Box<T, Dimensions> & clsWrap = *((MATH::DN::Box<T, Dimensions> *)(g_ptrWrapValues));
    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( MATH::DN::g_uintWrapMask & (1 << i) )
        {
            unsigned int uintMin = i << 1;
            unsigned int uintMax = uintMin + 1;
            m_clsValues[uintMin] = clsWrap[uintMin] + Mod(m_clsValues[uintMin], ptrWrapDistances[i]);
            m_clsValues[uintMax] = clsWrap[uintMax] + Mod(m_clsValues[uintMax], ptrWrapDistances[i]);
        }
    }
}

template <class T, unsigned int Dimensions>
T
MATH::DN::Box<T, Dimensions>::
DistanceW(const MATH::DN::Point<T, Dimensions> & p_stuPoint) const
{
    ASSERT(g_ptrWrapDistances != NULL && g_ptrWrapValues != NULL);

    T                          clsDistance = (T)(0);
    const T *                  ptrWrapDistances = (const T *)(g_ptrWrapDistances);
    const Box<T, Dimensions> & clsWrap = *((Box<T, Dimensions> *)(g_ptrWrapValues));
    for (unsigned int i=0; i<Dimensions; i++)
    {
        unsigned int uintMin = i << 1;
        unsigned int uintMax = uintMin + 1;

        if ( MATH::DN::g_uintWrapMask & (1 << i) )
        {
            T clsDistanceWrap, clsDistanceNoWrap;
            if ( m_clsValues[uintMin] <= m_clsValues[uintMax] )
            {
                if ( p_stuPoint[i] < m_clsValues[uintMin] )
                {
                    clsDistanceWrap = (p_stuPoint[i] - clsWrap[uintMin]) + (clsWrap[uintMax] - m_clsValues[uintMax]);
                    clsDistanceNoWrap = m_clsValues[uintMin] - p_stuPoint[i];
                    clsDistance += ( clsDistanceNoWrap < clsDistanceWrap )
                        ? (clsDistanceNoWrap * clsDistanceNoWrap)
                        : (clsDistanceWrap * clsDistanceWrap);
                }
                else if ( p_stuPoint[i] > m_clsValues[uintMax] )
                {
                    clsDistanceWrap = (clsWrap[uintMax] - p_stuPoint[i]) + (m_clsValues[uintMin] - clsWrap[uintMin]);
                    clsDistanceNoWrap = p_stuPoint[i] - m_clsValues[uintMax];
                    clsDistance += ( clsDistanceNoWrap < clsDistanceWrap )
                        ? (clsDistanceNoWrap * clsDistanceNoWrap)
                        : (clsDistanceWrap * clsDistanceWrap);
                }
            }
            else
            {
                if ( p_stuPoint[i] > m_clsValues[uintMax] &&
                     p_stuPoint[i] < m_clsValues[uintMin] )
                {
                    clsDistanceWrap = p_stuPoint[i] - m_clsValues[uintMax];
                    clsDistanceNoWrap = m_clsValues[uintMin] - p_stuPoint[i];
                    clsDistance += ( clsDistanceNoWrap < clsDistanceWrap )
                        ? (clsDistanceNoWrap * clsDistanceNoWrap)
                        : (clsDistanceWrap * clsDistanceWrap);
                }
            }
        }
        else
        {
            if ( p_stuPoint[i] < m_clsValues[uintMin] )
                clsDistance += (m_clsValues[uintMin] - p_stuPoint[i]) * (m_clsValues[uintMin] - p_stuPoint[i]);
            else if ( p_stuPoint[i] > m_clsValues[uintMax] )
                clsDistance += (p_stuPoint[i] - m_clsValues[uintMax]) * (p_stuPoint[i] - m_clsValues[uintMax]);
        }
    }

    return clsDistance;
}


template <class T, unsigned int Dimensions>
bool
MATH::DN::Box<T, Dimensions>::
Intersect(const MATH::DN::Point<T, Dimensions> & p_stuPoint) const
{
    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( p_stuPoint[i] < m_clsValues[2 * i] || p_stuPoint[i] > m_clsValues[(2 * i) + 1] )
            return false;
    }

    return true;
}

//
// Return Values:
//
// 1: No Box is completely within the other
// 2: Parameter Box is contained within this Box
// 3: This Box is contained within parameter Box
// 4: The Boxes are the same (contain each other)
//
template <class T, unsigned int Dimensions>
bool
MATH::DN::Box<T, Dimensions>::
Intersect(const MATH::DN::Box<T, Dimensions> & p_stuBox,
          MATH::DN::ResultBox<T, Dimensions> * p_ptrResult) const
{
    unsigned char uchrThisInsideParam = 1;
    unsigned char uchrParamInsideThis = 1;
    for (unsigned int i=0; i<Dimensions; i++)
    {
        unsigned int uintMin = i << 1;
        unsigned int uintMax = uintMin + 1;
        if ( m_clsValues[uintMin] > p_stuBox[uintMax] || m_clsValues[uintMax] < p_stuBox[uintMin] )
            return false;

        if ( m_clsValues[uintMin] < p_stuBox[uintMin] || m_clsValues[uintMax] > p_stuBox[uintMax] )
            uchrThisInsideParam = 0;
        else if ( p_stuBox[uintMin] < m_clsValues[uintMin] || p_stuBox[uintMax] > m_clsValues[uintMax] )
            uchrParamInsideThis = 0;

        if ( p_ptrResult )
        {
            p_ptrResult->m_stuResult[uintMin] = ( p_stuBox[uintMin] < m_clsValues[uintMin] )
                ? p_stuBox[uintMin]
                : m_clsValues[uintMin];

            p_ptrResult->m_stuResult[uintMax] = ( p_stuBox[uintMax] < m_clsValues[uintMax] )
                ? p_stuBox[uintMax]
                : m_clsValues[uintMax];
        }
    }

    if ( p_ptrResult )
    {
        char chrResults[] = {1, 2, 3, 4};
        p_ptrResult->m_intType = chrResults[uchrParamInsideThis | (uchrThisInsideParam << 1)];
    }

    return true;
}

template <class T, unsigned int Dimensions>
bool
MATH::DN::Box<T, Dimensions>::
IntersectW(const MATH::DN::Point<T, Dimensions> & p_stuPoint,
           const unsigned int p_uintMask) const
{
    for (unsigned int i=0; i<Dimensions; i++)
    {
        if ( (p_uintMask & (1 << i)) )
        {
            unsigned int uintMin = i << 1;
            unsigned int uintMax = uintMin + 1;
            if ( (MATH::DN::g_uintWrapMask & (1 << i)) == 0 || m_clsValues[uintMin] < m_clsValues[uintMax] )
            {
                if ( p_stuPoint[i] < m_clsValues[uintMin] ||
                     p_stuPoint[i] > m_clsValues[uintMax] )
                {
                    return false;
                }
            }
            else if ( p_stuPoint[i] < m_clsValues[uintMin] &&
                      p_stuPoint[i] > m_clsValues[uintMax] )
            {
                return false;
            }
        }
    }

    return true;
}

//
// Return Values:
//
// 1: Intersection, no Box is completely within the other
// 2: Intersection, no Box is completely within the other, but it includes a wrapped axis
// 3: Intersection, parameter Box is contained within this Box
// 4: Intersection, this Box is contained within parameter Box
// 5: Intersection, the Boxes are the same (contain each other)
//
template <class T, unsigned int Dimensions>
bool
MATH::DN::Box<T, Dimensions>::
IntersectW(const MATH::DN::Box<T, Dimensions> & p_stuBox,
           MATH::DN::ResultBox<T, Dimensions> * p_ptrResult) const
{
    bool blnWrapped = false;
    bool blnOneInsideTwo = true;
    bool blnTwoInsideOne = true;

    for (unsigned int i=0; i<Dimensions; i++)
    {
        unsigned int uintMin = i << 1;
        unsigned int uintMax = uintMin + 1;
        if ( MATH::DN::g_uintWrapMask & (1 << i) )
        {
            if ( m_clsValues[uintMin] < m_clsValues[uintMax] )
            {
                if ( p_stuBox[uintMin] < p_stuBox[uintMax] )
                {
                    if ( !((m_clsValues[uintMin] >= p_stuBox[uintMin] && m_clsValues[uintMin] <= p_stuBox[uintMax]) ||
                           (m_clsValues[uintMax] >= p_stuBox[uintMin] && m_clsValues[uintMax] <= p_stuBox[uintMax]) ||
                           (m_clsValues[uintMin] < p_stuBox[uintMin]  && m_clsValues[uintMax] >  p_stuBox[uintMax])) )
                    {
                        return false;
                    }
                    else if ( m_clsValues[uintMin] < p_stuBox[uintMin] || m_clsValues[uintMax] > p_stuBox[uintMax] )
                        blnOneInsideTwo = false;
                    else if ( p_stuBox[uintMin] < m_clsValues[uintMin] || p_stuBox[uintMax] > m_clsValues[uintMax] )
                        blnTwoInsideOne = false;
                }
                else
                {
                    blnTwoInsideOne = false;
                    if ( !((m_clsValues[uintMin] <= p_stuBox[uintMax]) || (m_clsValues[uintMax] >= p_stuBox[uintMin])) )
                        return 0;
                    else if ( !((m_clsValues[uintMin] >= p_stuBox[uintMin]) || (m_clsValues[uintMax] < p_stuBox[uintMax])) )
                        blnOneInsideTwo = false;

                    if ( p_stuBox[uintMax] >= m_clsValues[uintMin] && p_stuBox[uintMin] <= m_clsValues[uintMax] )
                        blnWrapped = true;
                }
            }
            else
            {
                if ( p_stuBox[uintMin] < p_stuBox[uintMax] )
                {
                    blnOneInsideTwo = false;
                    if ( !((p_stuBox[uintMin] <= m_clsValues[uintMax]) || (p_stuBox[uintMax] >= m_clsValues[uintMin])) )
                        return false;
                    else if ( !((p_stuBox[uintMin] >= m_clsValues[uintMin]) || (p_stuBox[uintMax] < m_clsValues[uintMax])) )
                        blnTwoInsideOne = false;

                    if ( m_clsValues[uintMax] >= p_stuBox[uintMin] && m_clsValues[uintMin] <= p_stuBox[uintMax] )
                        blnWrapped = true;
                }
                else
                {
                    if ( m_clsValues[uintMin] < p_stuBox[uintMin] || m_clsValues[uintMax] > p_stuBox[uintMax] )
                        blnOneInsideTwo = false;
                    if ( p_stuBox[uintMin] < m_clsValues[uintMin] || p_stuBox[uintMax] > m_clsValues[uintMax] )
                        blnTwoInsideOne = false;

                    if ( m_clsValues[uintMin] <= p_stuBox[uintMax] || p_stuBox[uintMin] <= m_clsValues[uintMax] )
                        blnWrapped = true;
                }
            }
        }
        else
        {
            if ( m_clsValues[uintMax] < p_stuBox[uintMin] ||
                 m_clsValues[uintMin] > p_stuBox[uintMax] )
            {
                return false;
            }
            if ( m_clsValues[uintMin] < p_stuBox[uintMin] ||
                 m_clsValues[uintMax] > p_stuBox[uintMax] )
            {
                blnOneInsideTwo = false;
            }
            if ( p_stuBox[uintMin] < m_clsValues[uintMin] ||
                 p_stuBox[uintMax] > m_clsValues[uintMax] )
            {
                blnTwoInsideOne = false;
            }
        }
    }

    if ( p_ptrResult )
    {
        if ( blnOneInsideTwo && blnTwoInsideOne )
            p_ptrResult->m_intType = 5;
        else if ( blnOneInsideTwo )
            p_ptrResult->m_intType = 4;
        else if ( blnTwoInsideOne )
            p_ptrResult->m_intType = 3;
        else
            p_ptrResult->m_intType = ( blnWrapped == false ) ? 1 : 2;
    }

    return true;
}

template <class T, unsigned int Dimensions>
unsigned long
MATH::DN::Box<T, Dimensions>::
SectorW(const MATH::DN::Point<T, Dimensions> & p_stuPoint,
        MATH::DN::Box<T, Dimensions> * p_ptrBox) const
{
    unsigned long ulngSector = 0;
    MATH::DN::Box<T, Dimensions> * ptrWrapValues = (MATH::DN::Box<T, Dimensions> *)(g_ptrWrapValues);
    MATH::DN::Box<T, Dimensions>   stuTemp;
    if ( p_ptrBox == NULL )
        p_ptrBox = &stuTemp;

    for (unsigned int i=0; i<Dimensions; i++)
    {
        unsigned int uintMin = i << 1;
        unsigned int uintMax = uintMin + 1;
        if ( (MATH::DN::g_uintWrapMask & (1 << i)) == 0 || m_clsValues[uintMin] < m_clsValues[uintMax] )
        {
            //
            // Calculate the middle of the area in this axis
            //
            T clsValue = (m_clsValues[uintMin] + m_clsValues[uintMax]) / ((T)(2));

            //
            // Pick the sector based on where the point is
            //
            if ( p_stuPoint[i] <= clsValue )
            {
                p_ptrBox->m_clsValues[uintMin] = m_clsValues[uintMin];
                p_ptrBox->m_clsValues[uintMax] = clsValue;
            }
            else
            {
                p_ptrBox->m_clsValues[uintMin] = clsValue;
                p_ptrBox->m_clsValues[uintMax] = m_clsValues[uintMax];
                ulngSector |= (1 << i);
            }
        }
        else
        {
            //
            // Calculate the middle of the node in this axis
            //
            T clsValue = UTIL::FixValueWrap(m_clsValues[uintMin] + (((ptrWrapValues->m_clsValues[uintMax] - m_clsValues[uintMin]) + (m_clsValues[uintMax] - ptrWrapValues->m_clsValues[uintMin])) / (T)(2)),
                                            ptrWrapValues->m_clsValues[uintMin],
                                            ptrWrapValues->m_clsValues[uintMax]);

            //
            // Pick the sector based on where the point is
            //
            bool blnLeft = ( ((p_stuPoint[i] >= m_clsValues[uintMin]) && (clsValue >= m_clsValues[uintMin])) ||
                             ((p_stuPoint[i] <= m_clsValues[uintMax]) && (clsValue <= m_clsValues[uintMax])) )
                ? ( p_stuPoint[i] <= clsValue )
                : ( p_stuPoint[i] > clsValue );

            if ( blnLeft )
            {
                p_ptrBox->m_clsValues[uintMin] = m_clsValues[uintMin];
                p_ptrBox->m_clsValues[uintMax] = clsValue;
            }
            else
            {
                p_ptrBox->m_clsValues[uintMin] = clsValue;
                p_ptrBox->m_clsValues[uintMax] = m_clsValues[uintMax];
                ulngSector |= (1 << i);
            }
        }
    }

    return ulngSector;
}

#endif