1

public/mathlib/IceKey.H

@@ -0,0 +1,62 @@
+// Purpose: Header file for the C++ ICE encryption class.
+//			Taken from public domain code, as written by Matthew Kwan - July 1996
+//			http://www.darkside.com.au/ice/
+
+#ifndef _IceKey_H
+#define _IceKey_H
+
+/*
+The IceKey class is used for encrypting and decrypting 64-bit blocks of data 
+with the ICE (Information Concealment Engine) encryption algorithm. 
+
+The constructor creates a new IceKey object that can be used to encrypt and decrypt data. 
+The level of encryption determines the size of the key, and hence its speed. 
+Level 0 uses the Thin-ICE variant, which is an 8-round cipher taking an 8-byte key. 
+This is the fastest option, and is generally considered to be at least as secure as DES, 
+although it is not yet certain whether it is as secure as its key size. 
+
+For levels n greater than zero, a 16n-round cipher is used, taking 8n-byte keys. 
+Although not as fast as level 0, these are very very secure. 
+
+Before an IceKey can be used to encrypt data, its key schedule must be set with the set() member function. 
+The length of the key required is determined by the level, as described above. 
+
+The member functions encrypt() and decrypt() encrypt and decrypt respectively data 
+in blocks of eight chracters, using the specified key. 
+
+Two functions keySize() and blockSize() are provided 
+which return the key and block size respectively, measured in bytes. 
+The key size is determined by the level, while the block size is always 8. 
+
+The destructor zeroes out and frees up all memory associated with the key. 
+*/
+
+class IceSubkey;
+
+class IceKey {
+    public:
+	IceKey (int n);
+	~IceKey ();
+
+	void		set (const unsigned char *key);
+
+	void		encrypt (const unsigned char *plaintext,
+					unsigned char *ciphertext) const;
+
+	void		decrypt (const unsigned char *ciphertext,
+					unsigned char *plaintext) const;
+
+	int		keySize () const;
+
+	int		blockSize () const;
+
+    private:
+	void		scheduleBuild (unsigned short *k, int n,
+							const int *keyrot);
+
+	int		_size;
+	int		_rounds;
+	IceSubkey	*_keysched;
+};
+
+#endif

public/mathlib/anorms.h

@@ -0,0 +1,25 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//=============================================================================//
+
+#ifndef ANORMS_H
+#define ANORMS_H
+#ifdef _WIN32
+#pragma once
+#endif
+
+
+#include "mathlib/vector.h"
+
+
+#define NUMVERTEXNORMALS	162
+
+// the angle between consecutive g_anorms[] vectors is ~14.55 degrees
+#define VERTEXNORMAL_CONE_INNER_ANGLE	DEG2RAD(7.275)
+
+extern Vector g_anorms[NUMVERTEXNORMALS];
+
+
+#endif // ANORMS_H

public/mathlib/bumpvects.h

@@ -0,0 +1,37 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $Workfile:     $
+// $Date:         $
+// $NoKeywords: $
+//=============================================================================//
+
+#ifndef BUMPVECTS_H
+#define BUMPVECTS_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include "mathlib/mathlib.h"
+
+#define OO_SQRT_2 0.70710676908493042f
+#define OO_SQRT_3 0.57735025882720947f
+#define OO_SQRT_6 0.40824821591377258f
+// sqrt( 2 / 3 )
+#define OO_SQRT_2_OVER_3 0.81649661064147949f
+
+#define NUM_BUMP_VECTS 3
+
+const TableVector g_localBumpBasis[NUM_BUMP_VECTS] = 
+{
+	{	OO_SQRT_2_OVER_3, 0.0f, OO_SQRT_3 },
+	{  -OO_SQRT_6, OO_SQRT_2, OO_SQRT_3 },
+	{  -OO_SQRT_6, -OO_SQRT_2, OO_SQRT_3 }
+};
+
+void GetBumpNormals( const Vector& sVect, const Vector& tVect, const Vector& flatNormal, 
+					 const Vector& phongNormal, Vector bumpNormals[NUM_BUMP_VECTS] );
+
+#endif // BUMPVECTS_H

public/mathlib/compressed_3d_unitvec.h

@@ -0,0 +1,284 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+#ifndef _3D_UNITVEC_H
+#define _3D_UNITVEC_H
+
+
+#define UNITVEC_DECLARE_STATICS \
+   float cUnitVector::mUVAdjustment[0x2000]; \
+   Vector cUnitVector::mTmpVec;
+
+// upper 3 bits
+#define SIGN_MASK  0xe000
+#define XSIGN_MASK 0x8000
+#define YSIGN_MASK 0x4000
+#define ZSIGN_MASK 0x2000
+
+// middle 6 bits - xbits
+#define TOP_MASK  0x1f80
+
+// lower 7 bits - ybits
+#define BOTTOM_MASK  0x007f
+
+// unitcomp.cpp : A Unit Vector to 16-bit word conversion
+// algorithm based on work of Rafael Baptista (rafael@oroboro.com)
+// Accuracy improved by O.D. (punkfloyd@rocketmail.com)
+// Used with Permission.
+
+// a compressed unit vector. reasonable fidelty for unit
+// vectors in a 16 bit package. Good enough for surface normals
+// we hope.
+class cUnitVector // : public c3dMathObject
+{
+public:
+   cUnitVector() { mVec = 0; }
+   cUnitVector( const Vector& vec )
+   {
+      packVector( vec );
+   }
+   cUnitVector( unsigned short val ) { mVec = val; }
+
+   cUnitVector& operator=( const Vector& vec )
+   { packVector( vec ); return *this; }
+
+   operator Vector()
+   {
+      unpackVector( mTmpVec );
+      return mTmpVec;
+   }
+
+   void packVector( const Vector& vec )
+   {
+      // convert from Vector to cUnitVector
+
+      Assert( vec.IsValid());
+      Vector tmp = vec;
+
+      // input vector does not have to be unit length
+      // Assert( tmp.length() <= 1.001f );
+
+      mVec = 0;
+      if ( tmp.x < 0 ) { mVec |= XSIGN_MASK; tmp.x = -tmp.x; }
+      if ( tmp.y < 0 ) { mVec |= YSIGN_MASK; tmp.y = -tmp.y; }
+      if ( tmp.z < 0 ) { mVec |= ZSIGN_MASK; tmp.z = -tmp.z; }
+
+      // project the normal onto the plane that goes through
+      // X0=(1,0,0),Y0=(0,1,0),Z0=(0,0,1).
+      // on that plane we choose an (projective!) coordinate system
+      // such that X0->(0,0), Y0->(126,0), Z0->(0,126),(0,0,0)->Infinity
+
+      // a little slower... old pack was 4 multiplies and 2 adds.
+      // This is 2 multiplies, 2 adds, and a divide....
+      float w = 126.0f / ( tmp.x + tmp.y + tmp.z );
+      long xbits = (long)( tmp.x * w );
+      long ybits = (long)( tmp.y * w );
+
+      Assert( xbits <  127 );
+      Assert( xbits >= 0   );
+      Assert( ybits <  127 );
+      Assert( ybits >= 0   );
+
+      // Now we can be sure that 0<=xp<=126, 0<=yp<=126, 0<=xp+yp<=126
+      // however for the sampling we want to transform this triangle
+      // into a rectangle.
+      if ( xbits >= 64 )
+      {
+         xbits = 127 - xbits;
+         ybits = 127 - ybits;
+      }
+
+      // now we that have xp in the range (0,127) and yp in
+      // the range (0,63), we can pack all the bits together
+      mVec |= ( xbits << 7 );
+      mVec |= ybits;
+   }
+
+   void unpackVector( Vector& vec )
+   {
+      // if we do a straightforward backward transform
+      // we will get points on the plane X0,Y0,Z0
+      // however we need points on a sphere that goes through
+      // these points. Therefore we need to adjust x,y,z so
+      // that x^2+y^2+z^2=1 by normalizing the vector. We have
+      // already precalculated the amount by which we need to
+      // scale, so all we do is a table lookup and a
+      // multiplication
+
+      // get the x and y bits
+      long xbits = (( mVec & TOP_MASK ) >> 7 );
+      long ybits = ( mVec & BOTTOM_MASK );
+
+      // map the numbers back to the triangle (0,0)-(0,126)-(126,0)
+      if (( xbits + ybits ) >= 127 )
+      {
+         xbits = 127 - xbits;
+         ybits = 127 - ybits;
+      }
+
+      // do the inverse transform and normalization
+      // costs 3 extra multiplies and 2 subtracts. No big deal.
+      float uvadj = mUVAdjustment[mVec & ~SIGN_MASK];
+      vec.x = uvadj * (float) xbits;
+      vec.y = uvadj * (float) ybits;
+      vec.z = uvadj * (float)( 126 - xbits - ybits );
+
+      // set all the sign bits
+      if ( mVec & XSIGN_MASK ) vec.x = -vec.x;
+      if ( mVec & YSIGN_MASK ) vec.y = -vec.y;
+      if ( mVec & ZSIGN_MASK ) vec.z = -vec.z;
+
+      Assert( vec.IsValid());
+   }
+
+   static void initializeStatics()
+   {
+      for ( int idx = 0; idx < 0x2000; idx++ )
+      {
+         long xbits = idx >> 7;
+         long ybits = idx & BOTTOM_MASK;
+
+         // map the numbers back to the triangle (0,0)-(0,127)-(127,0)
+         if (( xbits + ybits ) >= 127 )
+         {
+            xbits = 127 - xbits;
+            ybits = 127 - ybits;
+         }
+
+         // convert to 3D vectors
+         float x = (float)xbits;
+         float y = (float)ybits;
+         float z = (float)( 126 - xbits - ybits );
+		
+         // calculate the amount of normalization required
+         mUVAdjustment[idx] = 1.0f / sqrtf( y*y + z*z + x*x );
+         Assert( _finite( mUVAdjustment[idx]));
+
+         //cerr << mUVAdjustment[idx] << "\t";
+         //if ( xbits == 0 ) cerr << "\n";
+      }
+   }
+
+#if 0
+   void test()
+   {
+      #define TEST_RANGE 4
+      #define TEST_RANDOM 100
+      #define TEST_ANGERROR 1.0
+
+      float maxError = 0;
+      float avgError = 0;
+      int numVecs = 0;
+
+      {for ( int x = -TEST_RANGE; x < TEST_RANGE; x++ )
+      {
+         for ( int y = -TEST_RANGE; y < TEST_RANGE; y++ )
+         {
+            for ( int z = -TEST_RANGE; z < TEST_RANGE; z++ )
+            {
+               if (( x + y + z ) == 0 ) continue;
+
+               Vector vec( (float)x, (float)y, (float)z );
+               Vector vec2;
+
+               vec.normalize();
+               packVector( vec );
+               unpackVector( vec2 );
+
+               float ang = vec.dot( vec2 );
+               ang = (( fabs( ang ) > 0.99999f ) ? 0 : (float)acos(ang));
+
+               if (( ang > TEST_ANGERROR ) | ( !_finite( ang )))
+               {
+                  cerr << "error: " << ang << endl;
+                  cerr << "orig vec:       " << vec.x << ",\t"
+                       << vec.y << ",\t" << vec.z << "\tmVec: "
+                       << mVec << endl;
+                  cerr << "quantized vec2: " << vec2.x
+                       << ",\t" << vec2.y << ",\t"
+                       << vec2.z << endl << endl;
+               }
+               avgError += ang;
+               numVecs++;
+               if ( maxError < ang ) maxError = ang;
+            }
+         }
+      }}
+
+      for ( int w = 0; w < TEST_RANDOM; w++ )
+      {
+         Vector vec( genRandom(), genRandom(), genRandom());
+         Vector vec2;
+         vec.normalize();
+
+         packVector( vec );
+         unpackVector( vec2 );
+
+         float ang =vec.dot( vec2 );
+         ang = (( ang > 0.999f ) ? 0 : (float)acos(ang));
+
+         if (( ang > TEST_ANGERROR ) | ( !_finite( ang )))
+         {
+            cerr << "error: " << ang << endl;
+            cerr << "orig vec:       " << vec.x << ",\t"
+                 << vec.y << ",\t" << vec.z << "\tmVec: "
+                 << mVec << endl;
+            cerr << "quantized vec2: " << vec2.x << ",\t"
+                 << vec2.y << ",\t"
+                 << vec2.z << endl << endl;
+         }
+         avgError += ang;
+         numVecs++;
+         if ( maxError < ang ) maxError = ang;
+      }
+
+      { for ( int x = 0; x < 50; x++ )
+      {
+         Vector vec( (float)x, 25.0f, 0.0f );
+         Vector vec2;
+
+         vec.normalize();
+         packVector( vec );
+         unpackVector( vec2 );
+
+         float ang = vec.dot( vec2 );
+         ang = (( fabs( ang ) > 0.999f ) ? 0 : (float)acos(ang));
+
+         if (( ang > TEST_ANGERROR ) | ( !_finite( ang )))
+         {
+            cerr << "error: " << ang << endl;
+            cerr << "orig vec:       " << vec.x << ",\t"
+                 << vec.y << ",\t" << vec.z << "\tmVec: "
+                 << mVec << endl;
+            cerr << "   quantized vec2: " << vec2.x << ",\t"
+                 << vec2.y << ",\t" << vec2.z << endl << endl;
+         }
+
+         avgError += ang;
+         numVecs++;
+         if ( maxError < ang ) maxError = ang;
+      }}
+
+      cerr << "max angle error: " << maxError
+           << ", average error: " << avgError / numVecs
+           << ", num tested vecs: " << numVecs << endl;
+   }
+
+   friend ostream& operator<< ( ostream& os, const cUnitVector& vec )
+   { os << vec.mVec; return os; }
+#endif
+
+//protected: // !!!!
+
+   unsigned short mVec;
+   static float mUVAdjustment[0x2000];
+   static Vector mTmpVec;
+};
+
+#endif // _3D_VECTOR_H
+
+

public/mathlib/compressed_light_cube.h

@@ -0,0 +1,24 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//=============================================================================//
+
+#ifndef COMPRESSED_LIGHT_CUBE_H
+#define COMPRESSED_LIGHT_CUBE_H
+#ifdef _WIN32
+#pragma once
+#endif
+
+
+#include "mathlib/mathlib.h"
+
+
+struct CompressedLightCube
+{
+	DECLARE_BYTESWAP_DATADESC();
+	ColorRGBExp32 m_Color[6];
+};
+
+
+#endif // COMPRESSED_LIGHT_CUBE_H

public/mathlib/compressed_vector.h

@@ -0,0 +1,608 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef COMPRESSED_VECTOR_H
+#define COMPRESSED_VECTOR_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include <math.h>
+#include <float.h>
+
+// For vec_t, put this somewhere else?
+#include "basetypes.h"
+
+// For rand(). We really need a library!
+#include <stdlib.h>
+
+#include "tier0/dbg.h"
+#include "mathlib/vector.h"
+
+#include "mathlib/mathlib.h"
+
+#if defined( _X360 )
+#pragma bitfield_order( push, lsb_to_msb )
+#endif
+//=========================================================
+// fit a 3D vector into 32 bits
+//=========================================================
+
+class Vector32
+{
+public:
+	// Construction/destruction:
+	Vector32(void); 
+	Vector32(vec_t X, vec_t Y, vec_t Z);
+
+	// assignment
+	Vector32& operator=(const Vector &vOther);
+	operator Vector ();
+
+private:
+	unsigned short x:10;
+	unsigned short y:10;
+	unsigned short z:10;
+	unsigned short exp:2;
+};
+
+inline Vector32& Vector32::operator=(const Vector &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	static float expScale[4] = { 4.0f, 16.0f, 32.f, 64.f };
+
+	float fmax = Max( fabs( vOther.x ), fabs( vOther.y ) );
+	fmax = Max( fmax, (float)fabs( vOther.z ) );
+
+	for (exp = 0; exp < 3; exp++)
+	{
+		if (fmax < expScale[exp])
+			break;
+	}
+	Assert( fmax < expScale[exp] );
+
+	float fexp = 512.0f / expScale[exp];
+
+	x = Clamp( (int)(vOther.x * fexp) + 512, 0, 1023 );
+	y = Clamp( (int)(vOther.y * fexp) + 512, 0, 1023 );
+	z = Clamp( (int)(vOther.z * fexp) + 512, 0, 1023 );
+	return *this; 
+}
+
+
+inline Vector32::operator Vector ()
+{
+	Vector tmp;
+
+	static float expScale[4] = { 4.0f, 16.0f, 32.f, 64.f };
+
+	float fexp = expScale[exp] / 512.0f;
+
+	tmp.x = (((int)x) - 512) * fexp;
+	tmp.y = (((int)y) - 512) * fexp;
+	tmp.z = (((int)z) - 512) * fexp; 
+	return tmp; 
+}
+
+
+//=========================================================
+// Fit a unit vector into 32 bits
+//=========================================================
+
+class Normal32
+{
+public:
+	// Construction/destruction:
+	Normal32(void); 
+	Normal32(vec_t X, vec_t Y, vec_t Z);
+
+	// assignment
+	Normal32& operator=(const Vector &vOther);
+	operator Vector ();
+
+private:
+	unsigned short x:15;
+	unsigned short y:15;
+	unsigned short zneg:1;
+};
+
+
+inline Normal32& Normal32::operator=(const Vector &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	x = Clamp( (int)(vOther.x * 16384) + 16384, 0, 32767 );
+	y = Clamp( (int)(vOther.y * 16384) + 16384, 0, 32767 );
+	zneg = (vOther.z < 0);
+	//x = vOther.x; 
+	//y = vOther.y; 
+	//z = vOther.z; 
+	return *this; 
+}
+
+
+inline Normal32::operator Vector ()
+{
+	Vector tmp;
+
+	tmp.x = ((int)x - 16384) * (1 / 16384.0);
+	tmp.y = ((int)y - 16384) * (1 / 16384.0);
+	tmp.z = sqrt( 1 - tmp.x * tmp.x - tmp.y * tmp.y );
+	if (zneg)
+		tmp.z = -tmp.z;
+	return tmp; 
+}
+
+
+//=========================================================
+// 64 bit Quaternion
+//=========================================================
+
+class Quaternion64
+{
+public:
+	// Construction/destruction:
+	Quaternion64(void); 
+	Quaternion64(vec_t X, vec_t Y, vec_t Z);
+
+	// assignment
+	// Quaternion& operator=(const Quaternion64 &vOther);
+	Quaternion64& operator=(const Quaternion &vOther);
+	operator Quaternion ();
+private:
+	uint64 x:21;
+	uint64 y:21;
+	uint64 z:21;
+	uint64 wneg:1;
+};
+
+
+inline Quaternion64::operator Quaternion ()	
+{
+	Quaternion tmp;
+
+	// shift to -1048576, + 1048575, then round down slightly to -1.0 < x < 1.0
+	tmp.x = ((int)x - 1048576) * (1 / 1048576.5f);
+	tmp.y = ((int)y - 1048576) * (1 / 1048576.5f);
+	tmp.z = ((int)z - 1048576) * (1 / 1048576.5f);
+	tmp.w = sqrt( 1 - tmp.x * tmp.x - tmp.y * tmp.y - tmp.z * tmp.z );
+	if (wneg)
+		tmp.w = -tmp.w;
+	return tmp; 
+}
+
+inline Quaternion64& Quaternion64::operator=(const Quaternion &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	x = Clamp( (int)(vOther.x * 1048576) + 1048576, 0, 2097151 );
+	y = Clamp( (int)(vOther.y * 1048576) + 1048576, 0, 2097151 );
+	z = Clamp( (int)(vOther.z * 1048576) + 1048576, 0, 2097151 );
+	wneg = (vOther.w < 0);
+	return *this; 
+}
+
+//=========================================================
+// 48 bit Quaternion
+//=========================================================
+
+class Quaternion48
+{
+public:
+	// Construction/destruction:
+	Quaternion48(void); 
+	Quaternion48(vec_t X, vec_t Y, vec_t Z);
+
+	// assignment
+	// Quaternion& operator=(const Quaternion48 &vOther);
+	Quaternion48& operator=(const Quaternion &vOther);
+	operator Quaternion ();
+private:
+	unsigned short x:16;
+	unsigned short y:16;
+	unsigned short z:15;
+	unsigned short wneg:1;
+};
+
+
+inline Quaternion48::operator Quaternion ()	
+{
+	Quaternion tmp;
+
+	tmp.x = ((int)x - 32768) * (1 / 32768.0);
+	tmp.y = ((int)y - 32768) * (1 / 32768.0);
+	tmp.z = ((int)z - 16384) * (1 / 16384.0);
+	tmp.w = sqrt( 1 - tmp.x * tmp.x - tmp.y * tmp.y - tmp.z * tmp.z );
+	if (wneg)
+		tmp.w = -tmp.w;
+	return tmp; 
+}
+
+inline Quaternion48& Quaternion48::operator=(const Quaternion &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	x = Clamp( (int)(vOther.x * 32768) + 32768, 0, 65535 );
+	y = Clamp( (int)(vOther.y * 32768) + 32768, 0, 65535 );
+	z = Clamp( (int)(vOther.z * 16384) + 16384, 0, 32767 );
+	wneg = (vOther.w < 0);
+	return *this; 
+}
+
+//=========================================================
+// 32 bit Quaternion
+//=========================================================
+
+class Quaternion32
+{
+public:
+	// Construction/destruction:
+	Quaternion32(void); 
+	Quaternion32(vec_t X, vec_t Y, vec_t Z);
+
+	// assignment
+	// Quaternion& operator=(const Quaternion48 &vOther);
+	Quaternion32& operator=(const Quaternion &vOther);
+	operator Quaternion ();
+private:
+	unsigned int x:11;
+	unsigned int y:10;
+	unsigned int z:10;
+	unsigned int wneg:1;
+};
+
+
+inline Quaternion32::operator Quaternion ()	
+{
+	Quaternion tmp;
+
+	tmp.x = ((int)x - 1024) * (1 / 1024.0);
+	tmp.y = ((int)y - 512) * (1 / 512.0);
+	tmp.z = ((int)z - 512) * (1 / 512.0);
+	tmp.w = sqrt( 1 - tmp.x * tmp.x - tmp.y * tmp.y - tmp.z * tmp.z );
+	if (wneg)
+		tmp.w = -tmp.w;
+	return tmp; 
+}
+
+inline Quaternion32& Quaternion32::operator=(const Quaternion &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	x = Clamp( (int)(vOther.x * 1024) + 1024, 0, 2047 );
+	y = Clamp( (int)(vOther.y * 512) + 512, 0, 1023 );
+	z = Clamp( (int)(vOther.z * 512) + 512, 0, 1023 );
+	wneg = (vOther.w < 0);
+	return *this; 
+}
+
+//=========================================================
+// 16 bit float
+//=========================================================
+
+
+const int float32bias = 127;
+const int float16bias = 15;
+
+const float maxfloat16bits = 65504.0f;
+
+class float16
+{
+public:
+	//float16() {}
+	//float16( float f ) { m_storage.rawWord = ConvertFloatTo16bits(f); }
+
+	void Init() { m_storage.rawWord = 0; }
+//	float16& operator=(const float16 &other) { m_storage.rawWord = other.m_storage.rawWord; return *this; }
+//	float16& operator=(const float &other) { m_storage.rawWord = ConvertFloatTo16bits(other); return *this; }
+//	operator unsigned short () { return m_storage.rawWord; }
+//	operator float () { return Convert16bitFloatTo32bits( m_storage.rawWord ); }
+	unsigned short GetBits() const 
+	{ 
+		return m_storage.rawWord; 
+	}
+	float GetFloat() const 
+	{ 
+		return Convert16bitFloatTo32bits( m_storage.rawWord ); 
+	}
+	void SetFloat( float in ) 
+	{ 
+		m_storage.rawWord = ConvertFloatTo16bits( in ); 
+	}
+
+	bool IsInfinity() const
+	{
+		return m_storage.bits.biased_exponent == 31 && m_storage.bits.mantissa == 0;
+	}
+	bool IsNaN() const
+	{
+		return m_storage.bits.biased_exponent == 31 && m_storage.bits.mantissa != 0;
+	}
+
+	bool operator==(const float16 other) const { return m_storage.rawWord == other.m_storage.rawWord; }
+	bool operator!=(const float16 other) const { return m_storage.rawWord != other.m_storage.rawWord; }
+	
+//	bool operator< (const float other) const	   { return GetFloat() < other; }
+//	bool operator> (const float other) const	   { return GetFloat() > other; }
+
+protected:
+	union float32bits
+	{
+		float rawFloat;
+		struct 
+		{
+			unsigned int mantissa : 23;
+			unsigned int biased_exponent : 8;
+			unsigned int sign : 1;
+		} bits;
+	};
+
+	union float16bits
+	{
+		unsigned short rawWord;
+		struct
+		{
+			unsigned short mantissa : 10;
+			unsigned short biased_exponent : 5;
+			unsigned short sign : 1;
+		} bits;
+	};
+
+	static bool IsNaN( float16bits in )
+	{
+		return in.bits.biased_exponent == 31 && in.bits.mantissa != 0;
+	}
+	static bool IsInfinity( float16bits in )
+	{
+		return in.bits.biased_exponent == 31 && in.bits.mantissa == 0;
+	}
+
+	// 0x0001 - 0x03ff
+	static unsigned short ConvertFloatTo16bits( float input )
+	{
+		if ( input > maxfloat16bits )
+			input = maxfloat16bits;
+		else if ( input < -maxfloat16bits )
+			input = -maxfloat16bits;
+
+		float16bits output;
+		float32bits inFloat;
+
+		inFloat.rawFloat = input;
+
+		output.bits.sign = inFloat.bits.sign;
+
+		if ( (inFloat.bits.biased_exponent==0) && (inFloat.bits.mantissa==0) ) 
+		{ 
+			// zero
+			output.bits.mantissa = 0;
+			output.bits.biased_exponent = 0;
+		}
+		else if ( (inFloat.bits.biased_exponent==0) && (inFloat.bits.mantissa!=0) ) 
+		{  
+			// denorm -- denorm float maps to 0 half
+			output.bits.mantissa = 0;
+			output.bits.biased_exponent = 0;
+		}
+		else if ( (inFloat.bits.biased_exponent==0xff) && (inFloat.bits.mantissa==0) ) 
+		{ 
+#if 0
+			// infinity
+			output.bits.mantissa = 0;
+			output.bits.biased_exponent = 31;
+#else
+			// infinity maps to maxfloat
+			output.bits.mantissa = 0x3ff;
+			output.bits.biased_exponent = 0x1e;
+#endif
+		}
+		else if ( (inFloat.bits.biased_exponent==0xff) && (inFloat.bits.mantissa!=0) ) 
+		{ 
+#if 0
+			// NaN
+			output.bits.mantissa = 1;
+			output.bits.biased_exponent = 31;
+#else
+			// NaN maps to zero
+			output.bits.mantissa = 0;
+			output.bits.biased_exponent = 0;
+#endif
+		}
+		else 
+		{ 
+			// regular number
+			int new_exp = inFloat.bits.biased_exponent-127;
+
+			if (new_exp<-24) 
+			{ 
+				// this maps to 0
+				output.bits.mantissa = 0;
+				output.bits.biased_exponent = 0;
+			}
+
+			if (new_exp<-14) 
+			{
+				// this maps to a denorm
+				output.bits.biased_exponent = 0;
+				unsigned int exp_val = ( unsigned int )( -14 - ( inFloat.bits.biased_exponent - float32bias ) );
+				if( exp_val > 0 && exp_val < 11 )
+				{
+					output.bits.mantissa = ( 1 << ( 10 - exp_val ) ) + ( inFloat.bits.mantissa >> ( 13 + exp_val ) );
+				}
+			}
+			else if (new_exp>15) 
+			{ 
+#if 0
+				// map this value to infinity
+				output.bits.mantissa = 0;
+				output.bits.biased_exponent = 31;
+#else
+				// to big. . . maps to maxfloat
+				output.bits.mantissa = 0x3ff;
+				output.bits.biased_exponent = 0x1e;
+#endif
+			}
+			else 
+			{
+				output.bits.biased_exponent = new_exp+15;
+				output.bits.mantissa = (inFloat.bits.mantissa >> 13);
+			}
+		}
+		return output.rawWord;
+	}
+
+	static float Convert16bitFloatTo32bits( unsigned short input )
+	{
+		float32bits output;
+		const float16bits &inFloat = *((float16bits *)&input);
+
+		if( IsInfinity( inFloat ) )
+		{
+			return maxfloat16bits * ( ( inFloat.bits.sign == 1 ) ? -1.0f : 1.0f );
+		}
+		if( IsNaN( inFloat ) )
+		{
+			return 0.0;
+		}
+		if( inFloat.bits.biased_exponent == 0 && inFloat.bits.mantissa != 0 )
+		{
+			// denorm
+			const float half_denorm = (1.0f/16384.0f); // 2^-14
+			float mantissa = ((float)(inFloat.bits.mantissa)) / 1024.0f;
+			float sgn = (inFloat.bits.sign)? -1.0f :1.0f;
+			output.rawFloat = sgn*mantissa*half_denorm;
+		}
+		else
+		{
+			// regular number
+			unsigned mantissa = inFloat.bits.mantissa;
+			unsigned biased_exponent = inFloat.bits.biased_exponent;
+			unsigned sign = ((unsigned)inFloat.bits.sign) << 31;
+			biased_exponent = ( (biased_exponent - float16bias + float32bias) * (biased_exponent != 0) ) << 23;
+			mantissa <<= (23-10);
+
+			*((unsigned *)&output) = ( mantissa | biased_exponent | sign );
+		}
+		
+		return output.rawFloat;
+	}
+
+
+	float16bits m_storage;
+};
+
+class float16_with_assign : public float16
+{
+public:
+	float16_with_assign() {}
+	float16_with_assign( float f ) { m_storage.rawWord = ConvertFloatTo16bits(f); }
+
+	float16& operator=(const float16 &other) { m_storage.rawWord = ((float16_with_assign &)other).m_storage.rawWord; return *this; }
+	float16& operator=(const float &other) { m_storage.rawWord = ConvertFloatTo16bits(other); return *this; }
+//	operator unsigned short () const { return m_storage.rawWord; }
+	operator float () const { return Convert16bitFloatTo32bits( m_storage.rawWord ); }
+};
+
+//=========================================================
+// Fit a 3D vector in 48 bits
+//=========================================================
+
+class Vector48
+{
+public:
+	// Construction/destruction:
+	Vector48(void) {}
+	Vector48(vec_t X, vec_t Y, vec_t Z) { x.SetFloat( X ); y.SetFloat( Y ); z.SetFloat( Z ); }
+
+	// assignment
+	Vector48& operator=(const Vector &vOther);
+	operator Vector ();
+
+	const float operator[]( int i ) const { return (((float16 *)this)[i]).GetFloat(); }
+
+	float16 x;
+	float16 y;
+	float16 z;
+};
+
+inline Vector48& Vector48::operator=(const Vector &vOther)	
+{
+	CHECK_VALID(vOther);
+
+	x.SetFloat( vOther.x );
+	y.SetFloat( vOther.y );
+	z.SetFloat( vOther.z );
+	return *this; 
+}
+
+
+inline Vector48::operator Vector ()
+{
+	Vector tmp;
+
+	tmp.x = x.GetFloat();
+	tmp.y = y.GetFloat();
+	tmp.z = z.GetFloat(); 
+
+	return tmp;
+}
+
+//=========================================================
+// Fit a 2D vector in 32 bits
+//=========================================================
+
+class Vector2d32
+{
+public:
+	// Construction/destruction:
+	Vector2d32(void) {}
+	Vector2d32(vec_t X, vec_t Y) { x.SetFloat( X ); y.SetFloat( Y ); }
+
+	// assignment
+	Vector2d32& operator=(const Vector &vOther);
+	Vector2d32& operator=(const Vector2D &vOther);
+
+	operator Vector2D ();
+
+	void Init( vec_t ix = 0.f, vec_t iy = 0.f);
+
+	float16_with_assign x;
+	float16_with_assign y;
+};
+
+inline Vector2d32& Vector2d32::operator=(const Vector2D &vOther)	
+{
+	x.SetFloat( vOther.x );
+	y.SetFloat( vOther.y );
+	return *this; 
+}
+
+inline Vector2d32::operator Vector2D ()
+{
+	Vector2D tmp;
+
+	tmp.x = x.GetFloat();
+	tmp.y = y.GetFloat();
+
+	return tmp;
+}
+
+inline void Vector2d32::Init( vec_t ix, vec_t iy )
+{
+	x.SetFloat(ix);
+	y.SetFloat(iy);
+}
+
+#if defined( _X360 )
+#pragma bitfield_order( pop )
+#endif
+
+#endif
+

public/mathlib/halton.h

@@ -0,0 +1,71 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+// $Id$
+
+// halton.h - classes, etc for generating numbers using the Halton pseudo-random sequence.  See
+// http://halton-sequences.wikiverse.org/.
+//
+// what this function is useful for is any sort of sampling/integration problem where
+// you want to solve it by random sampling. Each call the NextValue() generates
+// a random number between 0 and 1, in an unclumped manner, so that the space can be more
+// or less evenly sampled with a minimum number of samples.
+//
+// It is NOT useful for generating random numbers dynamically, since the outputs aren't
+// particularly random.
+//
+// To generate multidimensional sample values (points in a plane, etc), use two
+// HaltonSequenceGenerator_t's, with different (primes) bases.
+
+#ifndef HALTON_H
+#define HALTON_H
+
+#include <tier0/platform.h>
+#include <mathlib/vector.h>
+
+class HaltonSequenceGenerator_t
+{
+	int seed;
+	int base;
+	float fbase;											//< base as a float
+
+public:
+	HaltonSequenceGenerator_t(int base);					//< base MUST be prime, >=2
+
+	float GetElement(int element);
+
+	inline float NextValue(void)
+	{
+		return GetElement(seed++);
+	}
+
+};
+
+
+class DirectionalSampler_t									//< pseudo-random sphere sampling
+{
+	HaltonSequenceGenerator_t zdot;
+	HaltonSequenceGenerator_t vrot;
+public:
+	DirectionalSampler_t(void)
+		: zdot(2),vrot(3)
+	{
+	}
+
+	Vector NextValue(void)
+	{
+		float zvalue=zdot.NextValue();
+		zvalue=2*zvalue-1.0;								// map from 0..1 to -1..1
+		float phi=acos(zvalue);
+		// now, generate a random rotation angle for x/y
+		float theta=2.0*M_PI*vrot.NextValue();
+		float sin_p=sin(phi);
+		return Vector(cos(theta)*sin_p,
+					  sin(theta)*sin_p,
+					  zvalue);
+
+	}
+};
+
+
+
+
+#endif // halton_h

public/mathlib/lightdesc.h

@@ -0,0 +1,173 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//===========================================================================//
+
+// light structure definitions.
+#ifndef LIGHTDESC_H
+#define LIGHTDESC_H
+
+#include <mathlib/ssemath.h>
+#include <mathlib/vector.h>
+
+//-----------------------------------------------------------------------------
+// Light structure
+//-----------------------------------------------------------------------------
+
+enum LightType_t
+{
+	MATERIAL_LIGHT_DISABLE = 0,
+	MATERIAL_LIGHT_POINT,
+	MATERIAL_LIGHT_DIRECTIONAL,
+	MATERIAL_LIGHT_SPOT,
+};
+
+enum LightType_OptimizationFlags_t
+{
+	LIGHTTYPE_OPTIMIZATIONFLAGS_HAS_ATTENUATION0 = 1,
+	LIGHTTYPE_OPTIMIZATIONFLAGS_HAS_ATTENUATION1 = 2,
+	LIGHTTYPE_OPTIMIZATIONFLAGS_HAS_ATTENUATION2 = 4,
+	LIGHTTYPE_OPTIMIZATIONFLAGS_DERIVED_VALUES_CALCED = 8,
+};
+
+struct LightDesc_t 
+{
+    LightType_t m_Type;										//< MATERIAL_LIGHT_xxx
+	Vector m_Color;											//< color+intensity 
+    Vector m_Position;										//< light source center position
+    Vector m_Direction;										//< for SPOT, direction it is pointing
+    float  m_Range;											//< distance range for light.0=infinite
+    float m_Falloff;										//< angular falloff exponent for spot lights
+    float m_Attenuation0;									//< constant distance falloff term
+    float m_Attenuation1;									//< linear term of falloff
+    float m_Attenuation2;									//< quadatic term of falloff
+    float m_Theta;											//< inner cone angle. no angular falloff 
+															//< within this cone
+    float m_Phi;											//< outer cone angle
+
+	// the values below are derived from the above settings for optimizations
+	// These aren't used by DX8. . used for software lighting.
+	float m_ThetaDot;
+	float m_PhiDot;
+	unsigned int m_Flags;
+protected:
+	float OneOver_ThetaDot_Minus_PhiDot;
+	float m_RangeSquared;
+public:
+
+	void RecalculateDerivedValues(void);			 // calculate m_xxDot, m_Type for changed parms
+
+	LightDesc_t(void)
+	{
+	}
+
+	// constructors for various useful subtypes
+
+	// a point light with infinite range
+	LightDesc_t( const Vector &pos, const Vector &color )
+	{
+		InitPoint( pos, color );
+	}
+	
+	/// a simple light. cone boundaries in radians. you pass a look_at point and the
+	/// direciton is derived from that.
+	LightDesc_t( const Vector &pos, const Vector &color, const Vector &point_at,
+				float inner_cone_boundary, float outer_cone_boundary )
+	{
+		InitSpot( pos, color, point_at, inner_cone_boundary, outer_cone_boundary );
+	}
+
+	void InitPoint( const Vector &pos, const Vector &color );
+	void InitDirectional( const Vector &dir, const Vector &color );
+	void InitSpot(const Vector &pos, const Vector &color, const Vector &point_at,
+		float inner_cone_boundary, float outer_cone_boundary );
+
+	/// Given 4 points and 4 normals, ADD lighting from this light into "color".
+	void ComputeLightAtPoints( const FourVectors &pos, const FourVectors &normal,
+							   FourVectors &color, bool DoHalfLambert=false ) const;
+	void ComputeNonincidenceLightAtPoints( const FourVectors &pos, FourVectors &color ) const;
+	void ComputeLightAtPointsForDirectional( const FourVectors &pos,
+											 const FourVectors &normal,
+											 FourVectors &color, bool DoHalfLambert=false ) const;
+
+	// warning - modifies color!!! set color first!!
+	void SetupOldStyleAttenuation( float fQuadatricAttn, float fLinearAttn, float fConstantAttn );
+
+	void SetupNewStyleAttenuation( float fFiftyPercentDistance, float fZeroPercentDistance );
+
+
+/// given a direction relative to the light source position, is this ray within the
+	/// light cone (for spotlights..non spots consider all rays to be within their cone)
+	bool IsDirectionWithinLightCone(const Vector &rdir) const
+	{
+		return ((m_Type!=MATERIAL_LIGHT_SPOT) || (rdir.Dot(m_Direction)>=m_PhiDot));
+	}
+
+	float OneOverThetaDotMinusPhiDot() const
+	{
+		return OneOver_ThetaDot_Minus_PhiDot;
+	}
+};
+
+
+//-----------------------------------------------------------------------------
+// a point light with infinite range
+//-----------------------------------------------------------------------------
+inline void LightDesc_t::InitPoint( const Vector &pos, const Vector &color )
+{
+	m_Type=MATERIAL_LIGHT_POINT;
+	m_Color=color;
+	m_Position=pos;
+	m_Range=0.0;									// infinite
+	m_Attenuation0=1.0;
+	m_Attenuation1=0;
+	m_Attenuation2=0;
+	RecalculateDerivedValues();
+}
+
+
+//-----------------------------------------------------------------------------
+// a directional light with infinite range
+//-----------------------------------------------------------------------------
+inline void LightDesc_t::InitDirectional( const Vector &dir, const Vector &color )
+{
+	m_Type=MATERIAL_LIGHT_DIRECTIONAL;
+	m_Color=color;
+	m_Direction=dir;
+	m_Range=0.0;									// infinite
+	m_Attenuation0=1.0;
+	m_Attenuation1=0;
+	m_Attenuation2=0;
+	RecalculateDerivedValues();
+}
+
+
+//-----------------------------------------------------------------------------
+// a simple light. cone boundaries in radians. you pass a look_at point and the
+// direciton is derived from that.
+//-----------------------------------------------------------------------------
+inline void LightDesc_t::InitSpot(const Vector &pos, const Vector &color, const Vector &point_at,
+	float inner_cone_boundary, float outer_cone_boundary)
+{
+	m_Type=MATERIAL_LIGHT_SPOT;
+	m_Color=color;
+	m_Position=pos;
+	m_Direction=point_at;
+	m_Direction-=pos;
+	VectorNormalizeFast(m_Direction);
+	m_Falloff=5.0;										// linear angle falloff
+	m_Theta=inner_cone_boundary;
+	m_Phi=outer_cone_boundary;
+
+	m_Range=0.0;										// infinite
+
+	m_Attenuation0=1.0;
+	m_Attenuation1=0;
+	m_Attenuation2=0;
+	RecalculateDerivedValues();
+}
+
+
+#endif
+

public/mathlib/math_pfns.h

@@ -0,0 +1,80 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//=====================================================================================//
+
+#ifndef _MATH_PFNS_H_
+#define _MATH_PFNS_H_
+
+#if defined( _X360 )
+#include <xboxmath.h>
+#endif
+
+#if !defined( _X360 )
+
+// These globals are initialized by mathlib and redirected based on available fpu features
+extern float (*pfSqrt)(float x);
+extern float (*pfRSqrt)(float x);
+extern float (*pfRSqrtFast)(float x);
+extern void  (*pfFastSinCos)(float x, float *s, float *c);
+extern float (*pfFastCos)(float x);
+
+// The following are not declared as macros because they are often used in limiting situations,
+// and sometimes the compiler simply refuses to inline them for some reason
+#define FastSqrt(x)			(*pfSqrt)(x)
+#define	FastRSqrt(x)		(*pfRSqrt)(x)
+#define FastRSqrtFast(x)    (*pfRSqrtFast)(x)
+#define FastSinCos(x,s,c)   (*pfFastSinCos)(x,s,c)
+#define FastCos(x)			(*pfFastCos)(x)
+
+#if defined(__i386__) || defined(_M_IX86)
+// On x86, the inline FPU or SSE sqrt instruction is faster than
+// the overhead of setting up a function call and saving/restoring
+// the FPU or SSE register state and can be scheduled better, too.
+#undef FastSqrt
+#define FastSqrt(x)			::sqrtf(x)
+#endif
+
+#endif // !_X360
+
+#if defined( _X360 )
+
+FORCEINLINE float _VMX_Sqrt( float x )
+{
+	return __fsqrts( x );
+}
+
+FORCEINLINE float _VMX_RSqrt( float x )
+{
+	float rroot = __frsqrte( x );
+
+	// Single iteration NewtonRaphson on reciprocal square root estimate
+	return (0.5f * rroot) * (3.0f - (x * rroot) * rroot);
+}
+
+FORCEINLINE float _VMX_RSqrtFast( float x )
+{
+	return __frsqrte( x );
+}
+
+FORCEINLINE void _VMX_SinCos( float a, float *pS, float *pC )
+{
+	XMScalarSinCos( pS, pC, a );
+}
+
+FORCEINLINE float _VMX_Cos( float a )
+{
+	return XMScalarCos( a );
+}
+
+// the 360 has fixed hw and calls directly
+#define FastSqrt(x)			_VMX_Sqrt(x)
+#define	FastRSqrt(x)		_VMX_RSqrt(x)
+#define FastRSqrtFast(x)	_VMX_RSqrtFast(x)
+#define FastSinCos(x,s,c)	_VMX_SinCos(x,s,c)
+#define FastCos(x)			_VMX_Cos(x)
+
+#endif // _X360
+
+#endif // _MATH_PFNS_H_

public/mathlib/matrixmath.h

@@ -0,0 +1,385 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//  A set of generic, template-based matrix functions.
+//===========================================================================//
+
+#ifndef MATRIXMATH_H
+#define MATRIXMATH_H
+
+#include <stdarg.h>
+
+// The operations in this file can perform basic matrix operations on matrices represented
+// using any class that supports the necessary operations:
+//
+//  .Element( row, col )  - return the element at a given matrox position
+//  .SetElement( row, col, val ) - modify an element
+//  .Width(), .Height() - get dimensions
+//  .SetDimensions( nrows, ncols) - set a matrix to be un-initted and the appropriate size
+//
+// Generally, vectors can be used with these functions by using N x 1 matrices to represent them.
+//  Matrices are addressed as row, column, and indices are 0-based
+//
+//
+// Note that the template versions of these routines are defined for generality - it is expected
+// that template specialization is used for common high performance cases.
+
+namespace MatrixMath
+{
+	/// M *= flScaleValue
+	template<class MATRIXCLASS>
+	void ScaleMatrix( MATRIXCLASS &matrix, float flScaleValue )
+	{
+		for( int i = 0; i < matrix.Height(); i++ )
+		{
+			for( int j = 0; j < matrix.Width(); j++ )
+			{
+				matrix.SetElement( i, j, flScaleValue * matrix.Element( i, j ) );
+			}
+		}
+	}
+
+	/// AppendElementToMatrix - same as setting the element, except only works when all calls
+	/// happen in top to bottom left to right order, end you have to call FinishedAppending when
+	/// done. For normal matrix classes this is not different then SetElement, but for
+	/// CSparseMatrix, it is an accelerated way to fill a matrix from scratch.
+	template<class MATRIXCLASS>
+	FORCEINLINE void AppendElement( MATRIXCLASS &matrix, int nRow, int nCol, float flValue )
+	{
+		matrix.SetElement( nRow, nCol, flValue );			// default implementation
+	}
+
+	template<class MATRIXCLASS>
+	FORCEINLINE void FinishedAppending( MATRIXCLASS &matrix ) {} // default implementation
+
+	/// M += fl
+	template<class MATRIXCLASS>
+	void AddToMatrix( MATRIXCLASS &matrix, float flAddend )
+	{
+		for( int i = 0; i < matrix.Height(); i++ )
+		{
+			for( int j = 0; j < matrix.Width(); j++ )
+			{
+				matrix.SetElement( i, j, flAddend + matrix.Element( i, j ) );
+			}
+		}
+	}
+
+	/// transpose
+	template<class MATRIXCLASSIN, class MATRIXCLASSOUT>
+	void TransposeMatrix( MATRIXCLASSIN const &matrixIn, MATRIXCLASSOUT *pMatrixOut )
+	{
+		pMatrixOut->SetDimensions( matrixIn.Width(), matrixIn.Height() );
+		for( int i = 0; i < pMatrixOut->Height(); i++ )
+		{
+			for( int j = 0; j < pMatrixOut->Width(); j++ )
+			{
+				AppendElement( *pMatrixOut, i, j, matrixIn.Element( j, i ) );
+			}
+		}
+		FinishedAppending( *pMatrixOut );
+	}
+
+	/// copy
+	template<class MATRIXCLASSIN, class MATRIXCLASSOUT>
+	void CopyMatrix( MATRIXCLASSIN const &matrixIn, MATRIXCLASSOUT *pMatrixOut )
+	{
+		pMatrixOut->SetDimensions( matrixIn.Height(), matrixIn.Width() );
+		for( int i = 0; i < matrixIn.Height(); i++ )
+		{
+			for( int j = 0; j < matrixIn.Width(); j++ )
+			{
+				AppendElement( *pMatrixOut, i, j, matrixIn.Element( i, j ) );
+			}
+		}
+		FinishedAppending( *pMatrixOut );
+	}
+
+
+
+	/// M+=M
+	template<class MATRIXCLASSIN, class MATRIXCLASSOUT>
+	void AddMatrixToMatrix( MATRIXCLASSIN const &matrixIn, MATRIXCLASSOUT *pMatrixOut )
+	{
+		for( int i = 0; i < matrixIn.Height(); i++ )
+		{
+			for( int j = 0; j < matrixIn.Width(); j++ )
+			{
+				pMatrixOut->SetElement( i, j, pMatrixOut->Element( i, j ) + matrixIn.Element( i, j ) );
+			}
+		}
+	}
+
+	// M += scale * M
+	template<class MATRIXCLASSIN, class MATRIXCLASSOUT>
+	void AddScaledMatrixToMatrix( float flScale, MATRIXCLASSIN const &matrixIn, MATRIXCLASSOUT *pMatrixOut )
+	{
+		for( int i = 0; i < matrixIn.Height(); i++ )
+		{
+			for( int j = 0; j < matrixIn.Width(); j++ )
+			{
+				pMatrixOut->SetElement( i, j, pMatrixOut->Element( i, j ) + flScale * matrixIn.Element( i, j ) );
+			}
+		}
+	}
+
+
+	// simple way to initialize a matrix with constants from code.
+	template<class MATRIXCLASSOUT> 
+	void SetMatrixToIdentity( MATRIXCLASSOUT *pMatrixOut, float flDiagonalValue = 1.0 )
+	{
+		for( int i = 0; i < pMatrixOut->Height(); i++ )
+		{
+			for( int j = 0; j < pMatrixOut->Width(); j++ )
+			{
+				AppendElement( *pMatrixOut, i, j, ( i == j ) ? flDiagonalValue : 0 );
+			}
+		}
+		FinishedAppending( *pMatrixOut );
+	}
+
+	//// simple way to initialize a matrix with constants from code
+	template<class MATRIXCLASSOUT> 
+	void SetMatrixValues( MATRIXCLASSOUT *pMatrix, int nRows, int nCols, ... )
+	{
+		va_list argPtr;
+		va_start( argPtr, nCols );
+
+		pMatrix->SetDimensions( nRows, nCols );
+		for( int nRow = 0; nRow < nRows; nRow++ )
+		{
+			for( int nCol = 0; nCol < nCols; nCol++ )
+			{
+				double flNewValue = va_arg( argPtr, double );
+				pMatrix->SetElement( nRow, nCol, flNewValue );
+			}
+		}
+		va_end( argPtr );
+	}
+
+
+	/// row and colum accessors. treat a row or a column as a column vector
+	template<class MATRIXTYPE> class MatrixRowAccessor
+	{
+	public:
+		FORCEINLINE MatrixRowAccessor( MATRIXTYPE const &matrix, int nRow )
+		{
+			m_pMatrix = &matrix;
+			m_nRow = nRow;
+		}
+
+		FORCEINLINE float Element( int nRow, int nCol ) const
+		{
+			Assert( nCol == 0 );
+			return m_pMatrix->Element( m_nRow, nRow );
+		}
+
+		FORCEINLINE int Width( void ) const { return 1; };
+		FORCEINLINE int Height( void ) const { return m_pMatrix->Width(); }
+
+	private:
+		MATRIXTYPE const *m_pMatrix;
+		int m_nRow;
+	};
+
+	template<class MATRIXTYPE> class MatrixColumnAccessor
+	{
+	public:
+		FORCEINLINE MatrixColumnAccessor( MATRIXTYPE const &matrix, int nColumn )
+		{
+			m_pMatrix = &matrix;
+			m_nColumn = nColumn;
+		}
+
+		FORCEINLINE float Element( int nRow, int nColumn ) const
+		{
+			Assert( nColumn == 0 );
+			return m_pMatrix->Element( nRow, m_nColumn );
+		}
+
+		FORCEINLINE int Width( void ) const { return 1; }
+		FORCEINLINE int Height( void ) const { return m_pMatrix->Height(); }
+	private:
+		MATRIXTYPE const *m_pMatrix;
+		int m_nColumn;
+	};
+
+	/// this translator acts as a proxy for the transposed matrix
+	template<class MATRIXTYPE> class MatrixTransposeAccessor
+	{
+	public:
+		FORCEINLINE MatrixTransposeAccessor( MATRIXTYPE const & matrix )
+		{
+			m_pMatrix = &matrix;
+		}
+
+		FORCEINLINE float Element( int nRow, int nColumn ) const
+		{
+			return m_pMatrix->Element( nColumn, nRow );
+		}
+
+		FORCEINLINE int Width( void ) const { return m_pMatrix->Height(); }
+		FORCEINLINE int Height( void ) const { return m_pMatrix->Width(); }
+	private:
+		MATRIXTYPE const *m_pMatrix;
+	};
+
+	/// this tranpose returns a wrapper around it's argument, allowing things like AddMatrixToMatrix( Transpose( matA ), &matB ) without an extra copy
+	template<class MATRIXCLASSIN>
+	MatrixTransposeAccessor<MATRIXCLASSIN> TransposeMatrix( MATRIXCLASSIN const &matrixIn )
+	{
+		return MatrixTransposeAccessor<MATRIXCLASSIN>( matrixIn );
+	}
+
+
+	/// retrieve rows and columns
+	template<class MATRIXTYPE>
+	FORCEINLINE MatrixColumnAccessor<MATRIXTYPE> MatrixColumn( MATRIXTYPE const &matrix, int nColumn )
+	{
+		return MatrixColumnAccessor<MATRIXTYPE>( matrix, nColumn );
+	}
+
+	template<class MATRIXTYPE>
+	FORCEINLINE MatrixRowAccessor<MATRIXTYPE> MatrixRow( MATRIXTYPE const &matrix, int nRow )
+	{
+		return MatrixRowAccessor<MATRIXTYPE>( matrix, nRow );
+	}
+
+	//// dot product between vectors (or rows and/or columns via accessors)
+	template<class MATRIXACCESSORATYPE, class MATRIXACCESSORBTYPE >
+	float InnerProduct( MATRIXACCESSORATYPE const &vecA, MATRIXACCESSORBTYPE const &vecB )
+	{
+		Assert( vecA.Width() == 1 );
+		Assert( vecB.Width() == 1 );
+		Assert( vecA.Height() == vecB.Height() );
+		double flResult = 0;
+		for( int i = 0; i < vecA.Height(); i++ )
+		{
+			flResult += vecA.Element( i, 0 ) * vecB.Element( i, 0 );
+		}
+		return flResult;
+	}
+
+
+
+	/// matrix x matrix multiplication
+	template<class MATRIXATYPE, class MATRIXBTYPE, class MATRIXOUTTYPE>
+	void MatrixMultiply( MATRIXATYPE const &matA, MATRIXBTYPE const &matB, MATRIXOUTTYPE *pMatrixOut )
+	{
+		Assert( matA.Width() == matB.Height() );
+		pMatrixOut->SetDimensions( matA.Height(), matB.Width() );
+		for( int i = 0; i < matA.Height(); i++ )
+		{
+			for( int j = 0; j < matB.Width(); j++ )
+			{
+				pMatrixOut->SetElement( i, j, InnerProduct( MatrixRow( matA, i ), MatrixColumn( matB, j ) ) );
+			}
+		}
+	}
+
+	/// solve Ax=B via the conjugate graident method. Code and naming conventions based on the
+	/// wikipedia article.
+	template<class ATYPE, class XTYPE, class BTYPE>
+	void ConjugateGradient( ATYPE const &matA, BTYPE const &vecB, XTYPE &vecX, float flTolerance = 1.0e-20 )
+	{
+		XTYPE vecR;
+		vecR.SetDimensions( vecX.Height(), 1 );
+		MatrixMultiply( matA, vecX, &vecR );
+		ScaleMatrix( vecR, -1 );
+		AddMatrixToMatrix( vecB, &vecR );
+		XTYPE vecP;
+		CopyMatrix( vecR, &vecP );
+		float flRsOld = InnerProduct( vecR, vecR );
+		for( int nIter = 0; nIter < 100; nIter++ )
+		{
+			XTYPE vecAp;
+			MatrixMultiply( matA, vecP, &vecAp );
+			float flDivisor = InnerProduct( vecAp, vecP );
+			float flAlpha = flRsOld / flDivisor;
+			AddScaledMatrixToMatrix( flAlpha, vecP, &vecX );
+			AddScaledMatrixToMatrix( -flAlpha, vecAp, &vecR );
+			float flRsNew = InnerProduct( vecR, vecR );
+			if ( flRsNew < flTolerance )
+			{
+				break;
+			}
+			ScaleMatrix( vecP, flRsNew / flRsOld );
+			AddMatrixToMatrix( vecR, &vecP );
+			flRsOld = flRsNew;
+		}
+	}
+
+	/// solve (A'*A) x=B via the conjugate gradient method. Code and naming conventions based on
+	/// the wikipedia article. Same as Conjugate gradient but allows passing in two matrices whose
+	/// product is used as the A matrix (in order to preserve sparsity)
+	template<class ATYPE, class APRIMETYPE, class XTYPE, class BTYPE>
+	void ConjugateGradient( ATYPE const &matA, APRIMETYPE const &matAPrime, BTYPE const &vecB, XTYPE &vecX, float flTolerance = 1.0e-20 )
+	{
+		XTYPE vecR1;
+		vecR1.SetDimensions( vecX.Height(), 1 );
+		MatrixMultiply( matA, vecX, &vecR1 );
+		XTYPE vecR;
+		vecR.SetDimensions( vecR1.Height(), 1 );
+		MatrixMultiply( matAPrime, vecR1, &vecR );
+		ScaleMatrix( vecR, -1 );
+		AddMatrixToMatrix( vecB, &vecR );
+		XTYPE vecP;
+		CopyMatrix( vecR, &vecP );
+		float flRsOld = InnerProduct( vecR, vecR );
+		for( int nIter = 0; nIter < 100; nIter++ )
+		{
+			XTYPE vecAp1;
+			MatrixMultiply( matA, vecP, &vecAp1 );
+			XTYPE vecAp;
+			MatrixMultiply( matAPrime, vecAp1, &vecAp );
+			float flDivisor = InnerProduct( vecAp, vecP );
+			float flAlpha = flRsOld / flDivisor;
+			AddScaledMatrixToMatrix( flAlpha, vecP, &vecX );
+			AddScaledMatrixToMatrix( -flAlpha, vecAp, &vecR );
+			float flRsNew = InnerProduct( vecR, vecR );
+			if ( flRsNew < flTolerance )
+			{
+				break;
+			}
+			ScaleMatrix( vecP, flRsNew / flRsOld );
+			AddMatrixToMatrix( vecR, &vecP );
+			flRsOld = flRsNew;
+		}
+	}
+
+	
+	template<class ATYPE,  class XTYPE, class BTYPE>
+	void LeastSquaresFit( ATYPE const &matA, BTYPE const &vecB, XTYPE &vecX )
+	{
+		// now, generate the normal equations
+		BTYPE vecBeta;
+		MatrixMath::MatrixMultiply( MatrixMath::TransposeMatrix( matA ), vecB, &vecBeta );
+
+		vecX.SetDimensions( matA.Width(), 1 );
+		MatrixMath::SetMatrixToIdentity( &vecX );
+
+		ATYPE matATransposed;
+		TransposeMatrix( matA, &matATransposed );
+		ConjugateGradient( matA, matATransposed, vecBeta, vecX, 1.0e-20 );
+	}
+
+};
+
+/// a simple fixed-size matrix class
+template<int NUMROWS, int NUMCOLS> class CFixedMatrix
+{
+public:
+	FORCEINLINE int Width( void ) const { return NUMCOLS; }
+	FORCEINLINE int Height( void ) const { return NUMROWS; }
+	FORCEINLINE float Element( int nRow, int nCol ) const { return m_flValues[nRow][nCol]; }
+	FORCEINLINE void SetElement( int nRow, int nCol, float flValue ) { m_flValues[nRow][nCol] = flValue; }
+	FORCEINLINE void SetDimensions( int nNumRows, int nNumCols ) { Assert( ( nNumRows == NUMROWS ) && ( nNumCols == NUMCOLS ) ); }
+
+private:
+	float m_flValues[NUMROWS][NUMCOLS];
+};
+
+
+
+#endif //matrixmath_h

public/mathlib/noise.h

@@ -0,0 +1,35 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+//=====================================================================================//
+
+#ifndef NOISE_H
+#define NOISE_H
+
+#include <math.h>
+#include "basetypes.h"
+#include "mathlib/vector.h"
+#include "tier0/dbg.h"
+
+
+// The following code is the c-ification of Ken Perlin's new noise algorithm
+// "JAVA REFERENCE IMPLEMENTATION OF IMPROVED NOISE - COPYRIGHT 2002 KEN PERLIN"
+// as available here: http://mrl.nyu.edu/~perlin/noise/
+// it generates a single octave of noise in the -1..1 range
+// this should at some point probably replace SparseConvolutionNoise - jd
+float ImprovedPerlinNoise( Vector const &pnt );
+
+// get the noise value at a point. Output range is 0..1.
+float SparseConvolutionNoise( Vector const &pnt );
+
+// get the noise value at a point, passing a custom noise shaping function. The noise shaping
+// function should map the domain 0..1 to 0..1.
+float SparseConvolutionNoise(Vector const &pnt, float (*pNoiseShapeFunction)(float) );
+
+// returns a 1/f noise. more octaves take longer
+float FractalNoise( Vector const &pnt, int n_octaves );
+
+// returns a abs(f)*1/f noise i.e. turbulence
+float Turbulence( Vector const &pnt, int n_octaves );
+#endif // NOISE_H

public/mathlib/polyhedron.h

@@ -0,0 +1,73 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef POLYHEDRON_H_
+#define	POLYHEDRON_H_
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include "mathlib/mathlib.h"
+
+
+
+struct Polyhedron_IndexedLine_t
+{
+	unsigned short iPointIndices[2];
+};
+
+struct Polyhedron_IndexedLineReference_t
+{
+	unsigned short iLineIndex;
+	unsigned char iEndPointIndex; //since two polygons reference any one line, one needs to traverse the line backwards, this flags that behavior
+};
+
+struct Polyhedron_IndexedPolygon_t
+{
+	unsigned short iFirstIndex;
+	unsigned short iIndexCount;
+	Vector polyNormal;
+};
+
+class CPolyhedron //made into a class because it's going virtual to support distinctions between temp and permanent versions
+{
+public:
+	Vector *pVertices;
+	Polyhedron_IndexedLine_t *pLines;
+	Polyhedron_IndexedLineReference_t *pIndices;
+	Polyhedron_IndexedPolygon_t *pPolygons;
+	
+	unsigned short iVertexCount;
+	unsigned short iLineCount;
+	unsigned short iIndexCount;
+	unsigned short iPolygonCount;
+
+	virtual ~CPolyhedron( void ) {};
+	virtual void Release( void ) = 0;
+	Vector Center( void );
+};
+
+class CPolyhedron_AllocByNew : public CPolyhedron
+{
+public:
+	virtual void Release( void );
+	static CPolyhedron_AllocByNew *Allocate( unsigned short iVertices, unsigned short iLines, unsigned short iIndices, unsigned short iPolygons ); //creates the polyhedron along with enough memory to hold all it's data in a single allocation
+
+private:
+	CPolyhedron_AllocByNew( void ) { }; //CPolyhedron_AllocByNew::Allocate() is the only way to create one of these.
+};
+
+CPolyhedron *GeneratePolyhedronFromPlanes( const float *pOutwardFacingPlanes, int iPlaneCount, float fOnPlaneEpsilon, bool bUseTemporaryMemory = false ); //be sure to polyhedron->Release()
+CPolyhedron *ClipPolyhedron( const CPolyhedron *pExistingPolyhedron, const float *pOutwardFacingPlanes, int iPlaneCount, float fOnPlaneEpsilon, bool bUseTemporaryMemory = false ); //this does NOT modify/delete the existing polyhedron
+
+CPolyhedron *GetTempPolyhedron( unsigned short iVertices, unsigned short iLines, unsigned short iIndices, unsigned short iPolygons ); //grab the temporary polyhedron. Avoids new/delete for quick work. Can only be in use by one chunk of code at a time
+
+
+#endif //#ifndef POLYHEDRON_H_
+

public/mathlib/quantize.h

@@ -0,0 +1,141 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+#ifndef QUANTIZE_H
+#define QUANTIZE_H
+
+#ifndef STRING_H
+#include <string.h>
+#endif
+
+#define MAXDIMS 768
+#define MAXQUANT 16000
+
+
+#include <tier0/platform.h>
+
+struct Sample;
+
+struct QuantizedValue {
+	double MinError;											// minimum possible error. used
+	// for neighbor searches.
+	struct QuantizedValue *Children[2];						// splits
+	int32 value;											// only exists for leaf nodes
+	struct Sample *Samples;									// every sample quantized into this
+	// entry
+	int32 NSamples;											// how many were quantized to this.
+	int32 TotSamples;
+	double *ErrorMeasure;									// variance measure for each dimension
+	double TotalError;										// sum of errors
+	uint8 *Mean;											// average value of each dimension
+	uint8 *Mins;											// min box for children and this
+	uint8 *Maxs;											// max box for children and this
+	int NQuant;												// the number of samples which were
+															// quantzied to this node since the
+															// last time OptimizeQuantizer()
+															// was called.
+	int *Sums;												// sum used by OptimizeQuantizer
+	int sortdim;											// dimension currently sorted along.
+};
+
+struct Sample {
+	int32 ID;												// identifier of this sample. can
+															// be used for any purpose.
+	int32 Count;											// number of samples this sample
+															// represents
+	int32 QNum;										   // what value this sample ended up quantized
+															// to.
+	struct QuantizedValue *qptr;							// ptr to what this was quantized to.
+	uint8 Value[1];										   // array of values for multi-dimensional
+	// variables.
+};
+
+void FreeQuantization(struct QuantizedValue *t);
+
+struct QuantizedValue *Quantize(struct Sample *s, int nsamples, int ndims,
+								int nvalues, uint8 *weights, int value0=0);
+
+int CompressSamples(struct Sample *s, int nsamples, int ndims);
+
+struct QuantizedValue *FindMatch(uint8 const *sample,
+								 int ndims,uint8 *weights,
+								 struct QuantizedValue *QTable);
+void PrintSamples(struct Sample const *s, int nsamples, int ndims);
+
+struct QuantizedValue *FindQNode(struct QuantizedValue const *q, int32 code);
+
+inline struct Sample *NthSample(struct Sample *s, int i, int nd)
+{
+	uint8 *r=(uint8 *) s;
+	r+=i*(sizeof(*s)+(nd-1));
+	return (struct Sample *) r;
+}
+
+inline struct Sample *AllocSamples(int ns, int nd)
+{
+	size_t size5=(sizeof(struct Sample)+(nd-1))*ns;
+	void *ret=new uint8[size5];
+	memset(ret,0,size5);
+	for(int i=0;i<ns;i++)
+		NthSample((struct Sample *)ret,i,nd)->Count=1;
+	return (struct Sample *) ret;
+}
+
+
+// MinimumError: what is the min error which will occur if quantizing
+// a sample to the given qnode? This is just the error if the qnode
+// is a leaf.
+double MinimumError(struct QuantizedValue const *q, uint8 const *sample,
+					int ndims, uint8 const *weights);
+double MaximumError(struct QuantizedValue const *q, uint8 const *sample,
+					int ndims, uint8 const *weights);
+
+void PrintQTree(struct QuantizedValue const *p,int idlevel=0);
+void OptimizeQuantizer(struct QuantizedValue *q, int ndims);
+
+// RecalculateVelues: update the means in a sample tree, based upon
+// the samples. can be used to reoptimize when samples are deleted,
+// for instance.
+
+void RecalculateValues(struct QuantizedValue *q, int ndims);
+
+extern double SquaredError;	// may be reset and examined. updated by
+															// FindMatch()
+
+
+
+
+// the routines below can be used for uniform quantization via dart-throwing.
+typedef void (*GENERATOR)(void *);    // generate a random sample
+typedef double (*COMPARER)(void const *a, void const *b);
+
+void *DartThrow(int NResults, int NTries, size_t itemsize, GENERATOR gen,
+				COMPARER cmp);
+void *FindClosestDart(void *items,int NResults, size_t itemsize,
+					  COMPARER cmp, void *lookfor, int *idx);
+
+
+
+
+// color quantization of 24 bit images
+#define QUANTFLAGS_NODITHER 1	// don't do Floyd-steinberg dither
+
+extern void ColorQuantize(
+uint8 const	*pImage,			// 4 byte pixels ARGB
+int			nWidth,
+int			nHeight,
+int			nFlags, 			// QUANTFLAGS_xxx
+int			nColors,			// # of colors to fill in in palette
+uint8		*pOutPixels,		// where to store resulting 8 bit pixels
+uint8		*pOutPalette,		// where to store resulting 768-byte palette
+int			nFirstColor);		// first color to use in mapping
+
+
+
+
+
+#endif

public/mathlib/simdvectormatrix.h

@@ -0,0 +1,142 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: Provide a class (SSE/SIMD only) holding a 2d matrix of class FourVectors,
+// for high speed processing in tools.
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef SIMDVECTORMATRIX_H
+#define SIMDVECTORMATRIX_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+
+#include <string.h>
+#include "tier0/platform.h"
+#include "tier0/dbg.h"
+#include "tier1/utlsoacontainer.h"
+#include "mathlib/ssemath.h"
+
+class CSIMDVectorMatrix
+{
+public:
+	int m_nWidth;											// in actual vectors
+	int m_nHeight;
+
+	int m_nPaddedWidth;										// # of 4x wide elements
+
+	FourVectors *m_pData;
+
+protected:
+	void Init( void )
+	{
+		m_pData = NULL;
+		m_nWidth = 0;
+		m_nHeight = 0;
+		m_nPaddedWidth = 0;
+	}
+
+	int NVectors( void ) const
+	{
+		return m_nHeight * m_nPaddedWidth;
+	}
+
+public:
+	// constructors and destructors
+	CSIMDVectorMatrix( void )
+	{
+		Init();
+	}
+
+	~CSIMDVectorMatrix( void )
+	{
+		if ( m_pData )
+			delete[] m_pData;
+	}
+
+	// set up storage and fields for m x n matrix. destroys old data
+	void SetSize( int width, int height )
+	{
+		if ( ( ! m_pData ) || ( width != m_nWidth ) || ( height != m_nHeight ) )
+		{
+			if ( m_pData )
+				delete[] m_pData;
+			
+			m_nWidth = width;
+			m_nHeight = height;
+			
+			m_nPaddedWidth = ( m_nWidth + 3) >> 2;
+			m_pData = NULL;
+			if ( width && height )
+				m_pData = new FourVectors[ m_nPaddedWidth * m_nHeight ];
+		}
+	}
+
+	CSIMDVectorMatrix( int width, int height )
+	{
+		Init();
+		SetSize( width, height );
+	}
+
+	CSIMDVectorMatrix &operator=( CSIMDVectorMatrix const &src )
+	{
+		SetSize( src.m_nWidth, src.m_nHeight );
+		if ( m_pData )
+			memcpy( m_pData, src.m_pData, m_nHeight*m_nPaddedWidth*sizeof(m_pData[0]) ); 
+		return *this;
+	}
+
+	CSIMDVectorMatrix &operator+=( CSIMDVectorMatrix const &src );
+
+	CSIMDVectorMatrix &operator*=( Vector const &src );
+
+	// create from an RGBA float bitmap. alpha ignored.
+	void CreateFromRGBA_FloatImageData(int srcwidth, int srcheight, float const *srcdata );
+
+	// create from 3 fields in a csoa
+	void CreateFromCSOAAttributes( CSOAContainer const *pSrc,
+								   int nAttrIdx0, int nAttrIdx1, int nAttrIdx2 );
+
+	// Element access. If you are calling this a lot, you don't want to use this class, because
+	// you're not getting the sse advantage
+	Vector Element(int x, int y) const
+	{
+		Assert( m_pData );
+		Assert( x < m_nWidth );
+		Assert( y < m_nHeight );
+		Vector ret;
+		FourVectors const *pData=m_pData+y*m_nPaddedWidth+(x >> 2);
+
+		int xo=(x & 3);
+		ret.x=pData->X( xo );
+		ret.y=pData->Y( xo );
+		ret.z=pData->Z( xo );
+		return ret;
+	}
+
+	//addressing the individual fourvectors elements
+	FourVectors &CompoundElement(int x, int y)
+	{
+		Assert( m_pData );
+		Assert( y < m_nHeight );
+		Assert( x < m_nPaddedWidth );
+		return m_pData[x + m_nPaddedWidth*y ];
+	}
+
+	// math operations on the whole image
+	void Clear( void )
+	{
+		Assert( m_pData );
+		memset( m_pData, 0, m_nHeight*m_nPaddedWidth*sizeof(m_pData[0]) );
+	}
+
+	void RaiseToPower( float power );
+};
+
+
+
+#endif

public/mathlib/spherical_geometry.h

@@ -0,0 +1,73 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: Functions for spherical geometry.
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef SPHERICAL_GEOMETRY_H
+#define SPHERICAL_GEOMETRY_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include <math.h>
+#include <float.h>
+
+// see http://mathworld.wolfram.com/SphericalTrigonometry.html
+
+// return the spherical distance, in radians, between 2 points on the unit sphere.
+FORCEINLINE float UnitSphereLineSegmentLength( Vector const &a, Vector const &b )
+{
+	// check unit length
+	Assert( fabs( VectorLength( a ) - 1.0 ) < 1.0e-3 );
+	Assert( fabs( VectorLength( b ) - 1.0 ) < 1.0e-3 );
+	return acos( DotProduct( a, b ) );
+}
+
+
+// given 3 points on the unit sphere, return the spherical area (in radians) of the triangle they form.
+// valid for "small" triangles.
+FORCEINLINE float UnitSphereTriangleArea( Vector const &a, Vector const &b , Vector const &c )
+{
+	float flLengthA = UnitSphereLineSegmentLength( b, c );
+	float flLengthB = UnitSphereLineSegmentLength( c, a );
+	float flLengthC = UnitSphereLineSegmentLength( a, b );
+	
+	if ( ( flLengthA == 0. ) || ( flLengthB == 0. ) || ( flLengthC == 0. ) )
+		return 0.;											// zero area triangle
+			
+	// now, find the 3 incribed angles for the triangle
+	float flHalfSumLens = 0.5 * ( flLengthA + flLengthB + flLengthC );
+	float flSinSums = sin( flHalfSumLens );
+	float flSinSMinusA= sin( flHalfSumLens - flLengthA );
+	float flSinSMinusB= sin( flHalfSumLens - flLengthB );
+	float flSinSMinusC= sin( flHalfSumLens - flLengthC );
+	
+	float flTanAOver2 = sqrt ( ( flSinSMinusB * flSinSMinusC ) / ( flSinSums * flSinSMinusA ) );
+	float flTanBOver2 = sqrt ( ( flSinSMinusA * flSinSMinusC ) / ( flSinSums * flSinSMinusB ) );
+	float flTanCOver2 = sqrt ( ( flSinSMinusA * flSinSMinusB ) / ( flSinSums * flSinSMinusC ) );
+
+	// Girards formula : area = sum of angles - pi.
+	return 2.0 * ( atan( flTanAOver2 ) + atan( flTanBOver2 ) + atan( flTanCOver2 ) ) - M_PI;
+}
+
+// spherical harmonics-related functions. Best explanation at http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.pdf
+
+// Evaluate associated legendre polynomial P( l, m ) at flX, using recurrence relation
+float AssociatedLegendrePolynomial( int nL, int nM, float flX );
+
+// Evaluate order N spherical harmonic with spherical coordinates
+// nL = band, 0..N
+// nM = -nL .. nL
+// theta = 0..M_PI
+// phi = 0.. 2 * M_PHI
+float SphericalHarmonic( int nL, int nM, float flTheta, float flPhi );
+
+// evaluate spherical harmonic with normalized vector direction
+float SphericalHarmonic( int nL, int nM, Vector const &vecDirection );
+
+
+#endif // SPHERICAL_GEOMETRY_H

public/mathlib/ssequaternion.h

@@ -0,0 +1,367 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: - defines SIMD "structure of arrays" classes and functions.
+//
+//===========================================================================//
+#ifndef SSEQUATMATH_H
+#define SSEQUATMATH_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+
+#include "mathlib/ssemath.h"
+
+// Use this #define to allow SSE versions of Quaternion math
+// to exist on PC.
+// On PC, certain horizontal vector operations are not supported.
+// This causes the SSE implementation of quaternion math to mix the
+// vector and scalar floating point units, which is extremely 
+// performance negative if you don't compile to native SSE2 (which 
+// we don't as of Sept 1, 2007). So, it's best not to allow these
+// functions to exist at all. It's not good enough to simply replace
+// the contents of the functions with scalar math, because each call
+// to LoadAligned and StoreAligned will result in an unnecssary copy
+// of the quaternion, and several moves to and from the XMM registers.
+//
+// Basically, the problem you run into is that for efficient SIMD code,
+// you need to load the quaternions and vectors into SIMD registers and
+// keep them there as long as possible while doing only SIMD math,
+// whereas for efficient scalar code, each time you copy onto or ever
+// use a fltx4, it hoses your pipeline. So the difference has to be
+// in the management of temporary variables in the calling function,
+// not inside the math functions.
+//
+// If you compile assuming the presence of SSE2, the MSVC will abandon
+// the traditional x87 FPU operations altogether and make everything use
+// the SSE2 registers, which lessens this problem a little.
+
+// permitted only on 360, as we've done careful tuning on its Altivec math:
+#ifdef _X360
+#define ALLOW_SIMD_QUATERNION_MATH 1  // not on PC!
+#endif
+
+
+
+//---------------------------------------------------------------------
+// Load/store quaternions
+//---------------------------------------------------------------------
+#ifndef _X360
+#if ALLOW_SIMD_QUATERNION_MATH
+// Using STDC or SSE
+FORCEINLINE fltx4 LoadAlignedSIMD( const QuaternionAligned & pSIMD )
+{
+	fltx4 retval = LoadAlignedSIMD( pSIMD.Base() );
+	return retval;
+}
+
+FORCEINLINE fltx4 LoadAlignedSIMD( const QuaternionAligned * RESTRICT pSIMD )
+{
+	fltx4 retval = LoadAlignedSIMD( pSIMD );
+	return retval;
+}
+
+FORCEINLINE void StoreAlignedSIMD( QuaternionAligned * RESTRICT pSIMD, const fltx4 & a )
+{
+	StoreAlignedSIMD( pSIMD->Base(), a );
+}
+#endif
+#else
+
+// for the transitional class -- load a QuaternionAligned
+FORCEINLINE fltx4 LoadAlignedSIMD( const QuaternionAligned & pSIMD )
+{
+	fltx4 retval = XMLoadVector4A( pSIMD.Base() );
+	return retval;
+}
+
+FORCEINLINE fltx4 LoadAlignedSIMD( const QuaternionAligned * RESTRICT pSIMD )
+{
+	fltx4 retval = XMLoadVector4A( pSIMD );
+	return retval;
+}
+
+FORCEINLINE void StoreAlignedSIMD( QuaternionAligned * RESTRICT pSIMD, const fltx4 & a )
+{
+	XMStoreVector4A( pSIMD->Base(), a );
+}
+
+#endif
+
+
+#if ALLOW_SIMD_QUATERNION_MATH
+//---------------------------------------------------------------------
+// Make sure quaternions are within 180 degrees of one another, if not, reverse q
+//---------------------------------------------------------------------
+FORCEINLINE fltx4 QuaternionAlignSIMD( const fltx4 &p, const fltx4 &q )
+{
+	// decide if one of the quaternions is backwards
+	fltx4 a = SubSIMD( p, q );
+	fltx4 b = AddSIMD( p, q );
+	a = Dot4SIMD( a, a );
+	b = Dot4SIMD( b, b );
+	fltx4 cmp = CmpGtSIMD( a, b );
+	fltx4 result = MaskedAssign( cmp, NegSIMD(q), q );
+	return result;
+}
+
+//---------------------------------------------------------------------
+// Normalize Quaternion
+//---------------------------------------------------------------------
+#if USE_STDC_FOR_SIMD
+
+FORCEINLINE fltx4 QuaternionNormalizeSIMD( const fltx4 &q )
+{
+	fltx4 radius, result;
+	radius = Dot4SIMD( q, q );
+
+	if ( SubFloat( radius, 0 ) ) // > FLT_EPSILON && ((radius < 1.0f - 4*FLT_EPSILON) || (radius > 1.0f + 4*FLT_EPSILON))
+	{
+		float iradius = 1.0f / sqrt( SubFloat( radius, 0 ) );
+		result = ReplicateX4( iradius );
+		result = MulSIMD( result, q );
+		return result;
+	}
+	return q;
+}
+
+#else
+
+// SSE + X360 implementation
+FORCEINLINE fltx4 QuaternionNormalizeSIMD( const fltx4 &q )
+{
+	fltx4 radius, result, mask;
+	radius = Dot4SIMD( q, q );
+	mask = CmpEqSIMD( radius, Four_Zeros ); // all ones iff radius = 0
+	result = ReciprocalSqrtSIMD( radius );
+	result = MulSIMD( result, q );
+	return MaskedAssign( mask, q, result );	// if radius was 0, just return q
+}
+
+#endif
+
+
+//---------------------------------------------------------------------
+// 0.0 returns p, 1.0 return q.
+//---------------------------------------------------------------------
+FORCEINLINE fltx4 QuaternionBlendNoAlignSIMD( const fltx4 &p, const fltx4 &q, float t )
+{
+	fltx4 sclp, sclq, result;
+	sclq = ReplicateX4( t );
+	sclp = SubSIMD( Four_Ones, sclq );
+	result = MulSIMD( sclp, p );
+	result = MaddSIMD( sclq, q, result );
+	return QuaternionNormalizeSIMD( result );
+}
+
+
+//---------------------------------------------------------------------
+// Blend Quaternions
+//---------------------------------------------------------------------
+FORCEINLINE fltx4 QuaternionBlendSIMD( const fltx4 &p, const fltx4 &q, float t )
+{
+	// decide if one of the quaternions is backwards
+	fltx4 q2, result;
+	q2 = QuaternionAlignSIMD( p, q );
+	result = QuaternionBlendNoAlignSIMD( p, q2, t );
+	return result;
+}
+
+
+//---------------------------------------------------------------------
+// Multiply Quaternions
+//---------------------------------------------------------------------
+#ifndef _X360
+
+// SSE and STDC
+FORCEINLINE fltx4 QuaternionMultSIMD( const fltx4 &p, const fltx4 &q )
+{
+	// decide if one of the quaternions is backwards
+	fltx4 q2, result;
+	q2 = QuaternionAlignSIMD( p, q );
+	SubFloat( result, 0 ) =  SubFloat( p, 0 ) * SubFloat( q2, 3 ) + SubFloat( p, 1 ) * SubFloat( q2, 2 ) - SubFloat( p, 2 ) * SubFloat( q2, 1 ) + SubFloat( p, 3 ) * SubFloat( q2, 0 );
+	SubFloat( result, 1 ) = -SubFloat( p, 0 ) * SubFloat( q2, 2 ) + SubFloat( p, 1 ) * SubFloat( q2, 3 ) + SubFloat( p, 2 ) * SubFloat( q2, 0 ) + SubFloat( p, 3 ) * SubFloat( q2, 1 );
+	SubFloat( result, 2 ) =  SubFloat( p, 0 ) * SubFloat( q2, 1 ) - SubFloat( p, 1 ) * SubFloat( q2, 0 ) + SubFloat( p, 2 ) * SubFloat( q2, 3 ) + SubFloat( p, 3 ) * SubFloat( q2, 2 );
+	SubFloat( result, 3 ) = -SubFloat( p, 0 ) * SubFloat( q2, 0 ) - SubFloat( p, 1 ) * SubFloat( q2, 1 ) - SubFloat( p, 2 ) * SubFloat( q2, 2 ) + SubFloat( p, 3 ) * SubFloat( q2, 3 );
+	return result;
+}
+
+#else 
+
+// X360
+extern const fltx4 g_QuatMultRowSign[4];
+FORCEINLINE fltx4 QuaternionMultSIMD( const fltx4 &p, const fltx4 &q )
+{
+	fltx4 q2, row, result;
+	q2 = QuaternionAlignSIMD( p, q );
+
+	row = XMVectorSwizzle( q2, 3, 2, 1, 0 );
+	row = MulSIMD( row, g_QuatMultRowSign[0] );
+	result = Dot4SIMD( row, p );
+
+	row = XMVectorSwizzle( q2, 2, 3, 0, 1 );
+	row = MulSIMD( row, g_QuatMultRowSign[1] );
+	row = Dot4SIMD( row, p );
+	result = __vrlimi( result, row, 4, 0 );
+	
+	row = XMVectorSwizzle( q2, 1, 0, 3, 2 );
+	row = MulSIMD( row, g_QuatMultRowSign[2] );
+	row = Dot4SIMD( row, p );
+	result = __vrlimi( result, row, 2, 0 );
+	
+	row = MulSIMD( q2, g_QuatMultRowSign[3] );
+	row = Dot4SIMD( row, p );
+	result = __vrlimi( result, row, 1, 0 );
+	return result;
+}
+
+#endif
+
+
+//---------------------------------------------------------------------
+// Quaternion scale
+//---------------------------------------------------------------------
+#ifndef _X360
+
+// SSE and STDC
+FORCEINLINE fltx4 QuaternionScaleSIMD( const fltx4 &p, float t )
+{
+	float r;
+	fltx4 q;
+
+	// FIXME: nick, this isn't overly sensitive to accuracy, and it may be faster to 
+	// use the cos part (w) of the quaternion (sin(omega)*N,cos(omega)) to figure the new scale.
+	float sinom = sqrt( SubFloat( p, 0 ) * SubFloat( p, 0 ) + SubFloat( p, 1 ) * SubFloat( p, 1 ) + SubFloat( p, 2 ) * SubFloat( p, 2 ) );
+	sinom = min( sinom, 1.f );
+
+	float sinsom = sin( asin( sinom ) * t );
+
+	t = sinsom / (sinom + FLT_EPSILON);
+	SubFloat( q, 0 ) = t * SubFloat( p, 0 );
+	SubFloat( q, 1 ) = t * SubFloat( p, 1 );
+	SubFloat( q, 2 ) = t * SubFloat( p, 2 );
+
+	// rescale rotation
+	r = 1.0f - sinsom * sinsom;
+
+	// Assert( r >= 0 );
+	if (r < 0.0f) 
+		r = 0.0f;
+	r = sqrt( r );
+
+	// keep sign of rotation
+	SubFloat( q, 3 ) = fsel( SubFloat( p, 3 ), r, -r );
+	return q;
+}
+
+#else
+
+// X360
+FORCEINLINE fltx4 QuaternionScaleSIMD( const fltx4 &p, float t )
+{
+	fltx4 sinom = Dot3SIMD( p, p );
+	sinom = SqrtSIMD( sinom );
+	sinom = MinSIMD( sinom, Four_Ones );
+	fltx4 sinsom = ArcSinSIMD( sinom );
+	fltx4 t4 = ReplicateX4( t );
+	sinsom = MulSIMD( sinsom, t4 );
+	sinsom = SinSIMD( sinsom );
+	sinom = AddSIMD( sinom, Four_Epsilons );
+	sinom = ReciprocalSIMD( sinom );
+	t4 = MulSIMD( sinsom, sinom );
+	fltx4 result = MulSIMD( p, t4 );
+
+	// rescale rotation
+	sinsom = MulSIMD( sinsom, sinsom );
+	fltx4 r = SubSIMD( Four_Ones, sinsom );
+	r = MaxSIMD( r, Four_Zeros );
+	r = SqrtSIMD( r );
+
+	// keep sign of rotation
+	fltx4 cmp = CmpGeSIMD( p, Four_Zeros );
+	r = MaskedAssign( cmp, r, NegSIMD( r ) );
+
+	result = __vrlimi(result, r, 1, 0);
+	return result;
+}
+
+#endif
+
+
+//-----------------------------------------------------------------------------
+// Quaternion sphereical linear interpolation
+//-----------------------------------------------------------------------------
+#ifndef _X360
+
+// SSE and STDC
+FORCEINLINE fltx4 QuaternionSlerpNoAlignSIMD( const fltx4 &p, const fltx4 &q, float t )
+{
+	float omega, cosom, sinom, sclp, sclq;
+
+	fltx4 result;
+
+	// 0.0 returns p, 1.0 return q.
+	cosom = SubFloat( p, 0 ) * SubFloat( q, 0 ) + SubFloat( p, 1 ) * SubFloat( q, 1 ) + 
+		SubFloat( p, 2 ) * SubFloat( q, 2 ) + SubFloat( p, 3 ) * SubFloat( q, 3 );
+
+	if ( (1.0f + cosom ) > 0.000001f ) 
+	{
+		if ( (1.0f - cosom ) > 0.000001f ) 
+		{
+			omega = acos( cosom );
+			sinom = sin( omega );
+			sclp = sin( (1.0f - t)*omega) / sinom;
+			sclq = sin( t*omega ) / sinom;
+		}
+		else 
+		{
+			// TODO: add short circuit for cosom == 1.0f?
+			sclp = 1.0f - t;
+			sclq = t;
+		}
+		SubFloat( result, 0 ) = sclp * SubFloat( p, 0 ) + sclq * SubFloat( q, 0 );
+		SubFloat( result, 1 ) = sclp * SubFloat( p, 1 ) + sclq * SubFloat( q, 1 );
+		SubFloat( result, 2 ) = sclp * SubFloat( p, 2 ) + sclq * SubFloat( q, 2 );
+		SubFloat( result, 3 ) = sclp * SubFloat( p, 3 ) + sclq * SubFloat( q, 3 );
+	}
+	else 
+	{
+		SubFloat( result, 0 ) = -SubFloat( q, 1 );
+		SubFloat( result, 1 ) =  SubFloat( q, 0 );
+		SubFloat( result, 2 ) = -SubFloat( q, 3 );
+		SubFloat( result, 3 ) =  SubFloat( q, 2 );
+		sclp = sin( (1.0f - t) * (0.5f * M_PI));
+		sclq = sin( t * (0.5f * M_PI));
+		SubFloat( result, 0 ) = sclp * SubFloat( p, 0 ) + sclq * SubFloat( result, 0 );
+		SubFloat( result, 1 ) = sclp * SubFloat( p, 1 ) + sclq * SubFloat( result, 1 );
+		SubFloat( result, 2 ) = sclp * SubFloat( p, 2 ) + sclq * SubFloat( result, 2 );
+	}
+
+	return result;
+}
+
+#else
+
+// X360
+FORCEINLINE fltx4 QuaternionSlerpNoAlignSIMD( const fltx4 &p, const fltx4 &q, float t )
+{
+	return XMQuaternionSlerp( p, q, t );
+}
+
+#endif
+
+
+FORCEINLINE fltx4 QuaternionSlerpSIMD( const fltx4 &p, const fltx4 &q, float t )
+{
+	fltx4 q2, result;
+	q2 = QuaternionAlignSIMD( p, q );
+	result = QuaternionSlerpNoAlignSIMD( p, q2, t );
+	return result;
+}
+
+
+#endif // ALLOW_SIMD_QUATERNION_MATH
+
+#endif // SSEQUATMATH_H
+

public/mathlib/vector2d.h

@@ -0,0 +1,670 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef VECTOR2D_H
+#define VECTOR2D_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include <math.h>
+#include <float.h>
+
+// For vec_t, put this somewhere else?
+#include "tier0/basetypes.h"
+
+// For rand(). We really need a library!
+#include <stdlib.h>
+
+#include "tier0/dbg.h"
+#include "mathlib/math_pfns.h"
+
+//=========================================================
+// 2D Vector2D
+//=========================================================
+
+class Vector2D					
+{
+public:
+	// Members
+	vec_t x, y;
+
+	// Construction/destruction
+	Vector2D(void);
+	Vector2D(vec_t X, vec_t Y);
+	Vector2D(const float *pFloat);
+
+	// Initialization
+	void Init(vec_t ix=0.0f, vec_t iy=0.0f);
+
+	// Got any nasty NAN's?
+	bool IsValid() const;
+
+	// array access...
+	vec_t operator[](int i) const;
+	vec_t& operator[](int i);
+
+	// Base address...
+	vec_t* Base();
+	vec_t const* Base() const;
+
+	// Initialization methods
+	void Random( float minVal, float maxVal );
+
+	// equality
+	bool operator==(const Vector2D& v) const;
+	bool operator!=(const Vector2D& v) const;	
+
+	// arithmetic operations
+	Vector2D&	operator+=(const Vector2D &v);			
+	Vector2D&	operator-=(const Vector2D &v);		
+	Vector2D&	operator*=(const Vector2D &v);			
+	Vector2D&	operator*=(float s);
+	Vector2D&	operator/=(const Vector2D &v);		
+	Vector2D&	operator/=(float s);					
+
+	// negate the Vector2D components
+	void	Negate(); 
+
+	// Get the Vector2D's magnitude.
+	vec_t	Length() const;
+
+	// Get the Vector2D's magnitude squared.
+	vec_t	LengthSqr(void) const;
+
+	// return true if this vector is (0,0) within tolerance
+	bool IsZero( float tolerance = 0.01f ) const
+	{
+		return (x > -tolerance && x < tolerance &&
+				y > -tolerance && y < tolerance);
+	}
+
+	// Normalize in place and return the old length.
+	vec_t	NormalizeInPlace();
+
+	// Compare length.
+	bool	IsLengthGreaterThan( float val ) const;
+	bool	IsLengthLessThan( float val ) const;
+
+	// Get the distance from this Vector2D to the other one.
+	vec_t	DistTo(const Vector2D &vOther) const;
+
+	// Get the distance from this Vector2D to the other one squared.
+	vec_t	DistToSqr(const Vector2D &vOther) const;		
+
+	// Copy
+	void	CopyToArray(float* rgfl) const;	
+
+	// Multiply, add, and assign to this (ie: *this = a + b * scalar). This
+	// is about 12% faster than the actual Vector2D equation (because it's done per-component
+	// rather than per-Vector2D).
+	void	MulAdd(const Vector2D& a, const Vector2D& b, float scalar);	
+
+	// Dot product.
+	vec_t	Dot(const Vector2D& vOther) const;			
+
+	// assignment
+	Vector2D& operator=(const Vector2D &vOther);
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// copy constructors
+	Vector2D(const Vector2D &vOther);
+
+	// arithmetic operations
+	Vector2D	operator-(void) const;
+				
+	Vector2D	operator+(const Vector2D& v) const;	
+	Vector2D	operator-(const Vector2D& v) const;	
+	Vector2D	operator*(const Vector2D& v) const;	
+	Vector2D	operator/(const Vector2D& v) const;	
+	Vector2D	operator*(float fl) const;
+	Vector2D	operator/(float fl) const;			
+	
+	// Cross product between two vectors.
+	Vector2D	Cross(const Vector2D &vOther) const;		
+
+	// Returns a Vector2D with the min or max in X, Y, and Z.
+	Vector2D	Min(const Vector2D &vOther) const;
+	Vector2D	Max(const Vector2D &vOther) const;
+
+#else
+
+private:
+	// No copy constructors allowed if we're in optimal mode
+	Vector2D(const Vector2D& vOther);
+#endif
+};
+
+//-----------------------------------------------------------------------------
+
+const Vector2D vec2_origin(0,0);
+const Vector2D vec2_invalid( FLT_MAX, FLT_MAX );
+
+//-----------------------------------------------------------------------------
+// Vector2D related operations
+//-----------------------------------------------------------------------------
+
+// Vector2D clear
+void Vector2DClear( Vector2D& a );
+
+// Copy
+void Vector2DCopy( const Vector2D& src, Vector2D& dst );
+
+// Vector2D arithmetic
+void Vector2DAdd( const Vector2D& a, const Vector2D& b, Vector2D& result );
+void Vector2DSubtract( const Vector2D& a, const Vector2D& b, Vector2D& result );
+void Vector2DMultiply( const Vector2D& a, vec_t b, Vector2D& result );
+void Vector2DMultiply( const Vector2D& a, const Vector2D& b, Vector2D& result );
+void Vector2DDivide( const Vector2D& a, vec_t b, Vector2D& result );
+void Vector2DDivide( const Vector2D& a, const Vector2D& b, Vector2D& result );
+void Vector2DMA( const Vector2D& start, float s, const Vector2D& dir, Vector2D& result );
+
+// Store the min or max of each of x, y, and z into the result.
+void Vector2DMin( const Vector2D &a, const Vector2D &b, Vector2D &result );
+void Vector2DMax( const Vector2D &a, const Vector2D &b, Vector2D &result );
+
+#define Vector2DExpand( v ) (v).x, (v).y
+
+// Normalization
+vec_t Vector2DNormalize( Vector2D& v );
+
+// Length
+vec_t Vector2DLength( const Vector2D& v );
+
+// Dot Product
+vec_t DotProduct2D(const Vector2D& a, const Vector2D& b);
+
+// Linearly interpolate between two vectors
+void Vector2DLerp(const Vector2D& src1, const Vector2D& src2, vec_t t, Vector2D& dest );
+
+
+//-----------------------------------------------------------------------------
+//
+// Inlined Vector2D methods
+//
+//-----------------------------------------------------------------------------
+
+
+//-----------------------------------------------------------------------------
+// constructors
+//-----------------------------------------------------------------------------
+
+inline Vector2D::Vector2D(void)									
+{ 
+#ifdef _DEBUG
+	// Initialize to NAN to catch errors
+	x = y = VEC_T_NAN;
+#endif
+}
+
+inline Vector2D::Vector2D(vec_t X, vec_t Y)						
+{ 
+	x = X; y = Y;
+	Assert( IsValid() );
+}
+
+inline Vector2D::Vector2D(const float *pFloat)					
+{
+	Assert( pFloat );
+	x = pFloat[0]; y = pFloat[1];	
+	Assert( IsValid() );
+}
+
+
+//-----------------------------------------------------------------------------
+// copy constructor
+//-----------------------------------------------------------------------------
+
+inline Vector2D::Vector2D(const Vector2D &vOther)					
+{ 
+	Assert( vOther.IsValid() );
+	x = vOther.x; y = vOther.y;
+}
+
+//-----------------------------------------------------------------------------
+// initialization
+//-----------------------------------------------------------------------------
+
+inline void Vector2D::Init( vec_t ix, vec_t iy )    
+{ 
+	x = ix; y = iy;
+	Assert( IsValid() );
+}
+
+inline void Vector2D::Random( float minVal, float maxVal )
+{
+	x = minVal + ((float)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+	y = minVal + ((float)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+}
+
+inline void Vector2DClear( Vector2D& a )
+{
+	a.x = a.y = 0.0f;
+}
+
+//-----------------------------------------------------------------------------
+// assignment
+//-----------------------------------------------------------------------------
+
+inline Vector2D& Vector2D::operator=(const Vector2D &vOther)	
+{
+	Assert( vOther.IsValid() );
+	x=vOther.x; y=vOther.y;
+	return *this; 
+}
+
+//-----------------------------------------------------------------------------
+// Array access
+//-----------------------------------------------------------------------------
+
+inline vec_t& Vector2D::operator[](int i)
+{
+	Assert( (i >= 0) && (i < 2) );
+	return ((vec_t*)this)[i];
+}
+
+inline vec_t Vector2D::operator[](int i) const
+{
+	Assert( (i >= 0) && (i < 2) );
+	return ((vec_t*)this)[i];
+}
+
+//-----------------------------------------------------------------------------
+// Base address...
+//-----------------------------------------------------------------------------
+
+inline vec_t* Vector2D::Base()
+{
+	return (vec_t*)this;
+}
+
+inline vec_t const* Vector2D::Base() const
+{
+	return (vec_t const*)this;
+}
+
+//-----------------------------------------------------------------------------
+// IsValid?
+//-----------------------------------------------------------------------------
+
+inline bool Vector2D::IsValid() const
+{
+	return IsFinite(x) && IsFinite(y);
+}
+
+//-----------------------------------------------------------------------------
+// comparison
+//-----------------------------------------------------------------------------
+
+inline bool Vector2D::operator==( const Vector2D& src ) const
+{
+	Assert( src.IsValid() && IsValid() );
+	return (src.x == x) && (src.y == y);
+}
+
+inline bool Vector2D::operator!=( const Vector2D& src ) const
+{
+	Assert( src.IsValid() && IsValid() );
+	return (src.x != x) || (src.y != y);
+}
+
+
+//-----------------------------------------------------------------------------
+// Copy
+//-----------------------------------------------------------------------------
+
+inline void Vector2DCopy( const Vector2D& src, Vector2D& dst )
+{
+	Assert( src.IsValid() );
+	dst.x = src.x;
+	dst.y = src.y;
+}
+
+inline void	Vector2D::CopyToArray(float* rgfl) const		
+{ 
+	Assert( IsValid() );
+	Assert( rgfl );
+	rgfl[0] = x; rgfl[1] = y; 
+}
+
+//-----------------------------------------------------------------------------
+// standard math operations
+//-----------------------------------------------------------------------------
+
+inline void Vector2D::Negate()
+{ 
+	Assert( IsValid() );
+	x = -x; y = -y;
+} 
+
+inline Vector2D& Vector2D::operator+=(const Vector2D& v)	
+{ 
+	Assert( IsValid() && v.IsValid() );
+	x+=v.x; y+=v.y;	
+	return *this;
+}
+
+inline Vector2D& Vector2D::operator-=(const Vector2D& v)	
+{ 
+	Assert( IsValid() && v.IsValid() );
+	x-=v.x; y-=v.y;	
+	return *this;
+}
+
+inline Vector2D& Vector2D::operator*=(float fl)	
+{
+	x *= fl;
+	y *= fl;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector2D& Vector2D::operator*=(const Vector2D& v)	
+{ 
+	x *= v.x;
+	y *= v.y;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector2D& Vector2D::operator/=(float fl)	
+{
+	Assert( fl != 0.0f );
+	float oofl = 1.0f / fl;
+	x *= oofl;
+	y *= oofl;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector2D& Vector2D::operator/=(const Vector2D& v)	
+{ 
+	Assert( v.x != 0.0f && v.y != 0.0f );
+	x /= v.x;
+	y /= v.y;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline void Vector2DAdd( const Vector2D& a, const Vector2D& b, Vector2D& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x + b.x;
+	c.y = a.y + b.y;
+}
+
+inline void Vector2DSubtract( const Vector2D& a, const Vector2D& b, Vector2D& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x - b.x;
+	c.y = a.y - b.y;
+}
+
+inline void Vector2DMultiply( const Vector2D& a, vec_t b, Vector2D& c )
+{
+	Assert( a.IsValid() && IsFinite(b) );
+	c.x = a.x * b;
+	c.y = a.y * b;
+}
+
+inline void Vector2DMultiply( const Vector2D& a, const Vector2D& b, Vector2D& c )
+{				  
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x * b.x;
+	c.y = a.y * b.y;
+}
+
+
+inline void Vector2DDivide( const Vector2D& a, vec_t b, Vector2D& c )
+{
+	Assert( a.IsValid() );
+	Assert( b != 0.0f );
+	vec_t oob = 1.0f / b;
+	c.x = a.x * oob;
+	c.y = a.y * oob;
+}
+
+inline void Vector2DDivide( const Vector2D& a, const Vector2D& b, Vector2D& c )
+{
+	Assert( a.IsValid() );
+	Assert( (b.x != 0.0f) && (b.y != 0.0f) );
+	c.x = a.x / b.x;
+	c.y = a.y / b.y;
+}
+
+inline void Vector2DMA( const Vector2D& start, float s, const Vector2D& dir, Vector2D& result )
+{
+	Assert( start.IsValid() && IsFinite(s) && dir.IsValid() );
+	result.x = start.x + s*dir.x;
+	result.y = start.y + s*dir.y;
+}
+
+// FIXME: Remove
+// For backwards compatability
+inline void	Vector2D::MulAdd(const Vector2D& a, const Vector2D& b, float scalar)
+{
+	x = a.x + b.x * scalar;
+	y = a.y + b.y * scalar;
+}
+
+inline void Vector2DLerp(const Vector2D& src1, const Vector2D& src2, vec_t t, Vector2D& dest )
+{
+	dest[0] = src1[0] + (src2[0] - src1[0]) * t;
+	dest[1] = src1[1] + (src2[1] - src1[1]) * t;
+}
+
+//-----------------------------------------------------------------------------
+// dot, cross
+//-----------------------------------------------------------------------------
+inline vec_t DotProduct2D(const Vector2D& a, const Vector2D& b) 
+{ 
+	Assert( a.IsValid() && b.IsValid() );
+	return( a.x*b.x + a.y*b.y ); 
+}
+
+// for backwards compatability
+inline vec_t Vector2D::Dot( const Vector2D& vOther ) const
+{
+	return DotProduct2D( *this, vOther );
+}
+
+
+//-----------------------------------------------------------------------------
+// length
+//-----------------------------------------------------------------------------
+inline vec_t Vector2DLength( const Vector2D& v )
+{
+	Assert( v.IsValid() );
+	return (vec_t)FastSqrt(v.x*v.x + v.y*v.y);		
+}
+
+inline vec_t Vector2D::LengthSqr(void) const	
+{ 
+	Assert( IsValid() );
+	return (x*x + y*y);		
+}
+
+inline vec_t Vector2D::NormalizeInPlace()
+{
+	return Vector2DNormalize( *this );
+}
+
+inline bool Vector2D::IsLengthGreaterThan( float val ) const
+{
+	return LengthSqr() > val*val;
+}
+
+inline bool Vector2D::IsLengthLessThan( float val ) const
+{
+	return LengthSqr() < val*val;
+}
+
+inline vec_t Vector2D::Length(void) const	
+{
+	return Vector2DLength( *this );
+}
+
+
+inline void Vector2DMin( const Vector2D &a, const Vector2D &b, Vector2D &result )
+{
+	result.x = (a.x < b.x) ? a.x : b.x;
+	result.y = (a.y < b.y) ? a.y : b.y;
+}
+
+
+inline void Vector2DMax( const Vector2D &a, const Vector2D &b, Vector2D &result )
+{
+	result.x = (a.x > b.x) ? a.x : b.x;
+	result.y = (a.y > b.y) ? a.y : b.y;
+}
+
+
+//-----------------------------------------------------------------------------
+// Normalization
+//-----------------------------------------------------------------------------
+inline vec_t Vector2DNormalize( Vector2D& v )
+{
+	Assert( v.IsValid() );
+	vec_t l = v.Length();
+	if (l != 0.0f)
+	{
+		v /= l;
+	}
+	else
+	{
+		v.x = v.y = 0.0f; 
+	}
+	return l;
+}
+
+
+//-----------------------------------------------------------------------------
+// Get the distance from this Vector2D to the other one 
+//-----------------------------------------------------------------------------
+inline vec_t Vector2D::DistTo(const Vector2D &vOther) const
+{
+	Vector2D delta;
+	Vector2DSubtract( *this, vOther, delta );
+	return delta.Length();
+}
+
+inline vec_t Vector2D::DistToSqr(const Vector2D &vOther) const
+{
+	Vector2D delta;
+	Vector2DSubtract( *this, vOther, delta );
+	return delta.LengthSqr();
+}
+
+
+//-----------------------------------------------------------------------------
+// Computes the closest point to vecTarget no farther than flMaxDist from vecStart
+//-----------------------------------------------------------------------------
+inline void ComputeClosestPoint2D( const Vector2D& vecStart, float flMaxDist, const Vector2D& vecTarget, Vector2D *pResult )
+{
+	Vector2D vecDelta;
+	Vector2DSubtract( vecTarget, vecStart, vecDelta );
+	float flDistSqr = vecDelta.LengthSqr();
+	if ( flDistSqr <= flMaxDist * flMaxDist )
+	{
+		*pResult = vecTarget;
+	}
+	else
+	{
+		vecDelta /= FastSqrt( flDistSqr );
+		Vector2DMA( vecStart, flMaxDist, vecDelta, *pResult );
+	}
+}
+
+
+
+//-----------------------------------------------------------------------------
+//
+// Slow methods
+//
+//-----------------------------------------------------------------------------
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+//-----------------------------------------------------------------------------
+// Returns a Vector2D with the min or max in X, Y, and Z.
+//-----------------------------------------------------------------------------
+
+inline Vector2D Vector2D::Min(const Vector2D &vOther) const
+{
+	return Vector2D(x < vOther.x ? x : vOther.x, 
+		y < vOther.y ? y : vOther.y);
+}
+
+inline Vector2D Vector2D::Max(const Vector2D &vOther) const
+{
+	return Vector2D(x > vOther.x ? x : vOther.x, 
+		y > vOther.y ? y : vOther.y);
+}
+
+
+//-----------------------------------------------------------------------------
+// arithmetic operations
+//-----------------------------------------------------------------------------
+
+inline Vector2D Vector2D::operator-(void) const
+{ 
+	return Vector2D(-x,-y);				
+}
+
+inline Vector2D Vector2D::operator+(const Vector2D& v) const	
+{ 
+	Vector2D res;
+	Vector2DAdd( *this, v, res );
+	return res;	
+}
+
+inline Vector2D Vector2D::operator-(const Vector2D& v) const	
+{ 
+	Vector2D res;
+	Vector2DSubtract( *this, v, res );
+	return res;	
+}
+
+inline Vector2D Vector2D::operator*(float fl) const	
+{ 
+	Vector2D res;
+	Vector2DMultiply( *this, fl, res );
+	return res;	
+}
+
+inline Vector2D Vector2D::operator*(const Vector2D& v) const	
+{ 
+	Vector2D res;
+	Vector2DMultiply( *this, v, res );
+	return res;	
+}
+
+inline Vector2D Vector2D::operator/(float fl) const	
+{ 
+	Vector2D res;
+	Vector2DDivide( *this, fl, res );
+	return res;	
+}
+
+inline Vector2D Vector2D::operator/(const Vector2D& v) const	
+{ 
+	Vector2D res;
+	Vector2DDivide( *this, v, res );
+	return res;	
+}
+
+inline Vector2D operator*(float fl, const Vector2D& v)	
+{ 
+	return v * fl; 
+}
+
+#endif //slow
+
+#endif // VECTOR2D_H
+

public/mathlib/vector4d.h

@@ -0,0 +1,686 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+
+#ifndef VECTOR4D_H
+#define VECTOR4D_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include <math.h>
+#include <stdlib.h>		// For rand(). We really need a library!
+#include <float.h>
+#if !defined( _X360 )
+#include <xmmintrin.h>	// For SSE
+#endif
+#include "basetypes.h"	// For vec_t, put this somewhere else?
+#include "tier0/dbg.h"
+#include "mathlib/math_pfns.h"
+
+// forward declarations
+class Vector;
+class Vector2D;
+
+//=========================================================
+// 4D Vector4D
+//=========================================================
+
+class Vector4D					
+{
+public:
+	// Members
+	vec_t x, y, z, w;
+
+	// Construction/destruction
+	Vector4D(void);
+	Vector4D(vec_t X, vec_t Y, vec_t Z, vec_t W);
+	Vector4D(const float *pFloat);
+
+	// Initialization
+	void Init(vec_t ix=0.0f, vec_t iy=0.0f, vec_t iz=0.0f, vec_t iw=0.0f);
+
+	// Got any nasty NAN's?
+	bool IsValid() const;
+
+	// array access...
+	vec_t operator[](int i) const;
+	vec_t& operator[](int i);
+
+	// Base address...
+	inline vec_t* Base();
+	inline vec_t const* Base() const;
+
+	// Cast to Vector and Vector2D...
+	Vector& AsVector3D();
+	Vector const& AsVector3D() const;
+
+	Vector2D& AsVector2D();
+	Vector2D const& AsVector2D() const;
+
+	// Initialization methods
+	void Random( vec_t minVal, vec_t maxVal );
+
+	// equality
+	bool operator==(const Vector4D& v) const;
+	bool operator!=(const Vector4D& v) const;	
+
+	// arithmetic operations
+	Vector4D&	operator+=(const Vector4D &v);			
+	Vector4D&	operator-=(const Vector4D &v);		
+	Vector4D&	operator*=(const Vector4D &v);			
+	Vector4D&	operator*=(float s);
+	Vector4D&	operator/=(const Vector4D &v);		
+	Vector4D&	operator/=(float s);					
+
+	// negate the Vector4D components
+	void	Negate(); 
+
+	// Get the Vector4D's magnitude.
+	vec_t	Length() const;
+
+	// Get the Vector4D's magnitude squared.
+	vec_t	LengthSqr(void) const;
+
+	// return true if this vector is (0,0,0,0) within tolerance
+	bool IsZero( float tolerance = 0.01f ) const
+	{
+		return (x > -tolerance && x < tolerance &&
+				y > -tolerance && y < tolerance &&
+				z > -tolerance && z < tolerance &&
+				w > -tolerance && w < tolerance);
+	}
+
+	// Get the distance from this Vector4D to the other one.
+	vec_t	DistTo(const Vector4D &vOther) const;
+
+	// Get the distance from this Vector4D to the other one squared.
+	vec_t	DistToSqr(const Vector4D &vOther) const;		
+
+	// Copy
+	void	CopyToArray(float* rgfl) const;	
+
+	// Multiply, add, and assign to this (ie: *this = a + b * scalar). This
+	// is about 12% faster than the actual Vector4D equation (because it's done per-component
+	// rather than per-Vector4D).
+	void	MulAdd(Vector4D const& a, Vector4D const& b, float scalar);	
+
+	// Dot product.
+	vec_t	Dot(Vector4D const& vOther) const;			
+
+	// No copy constructors allowed if we're in optimal mode
+#ifdef VECTOR_NO_SLOW_OPERATIONS
+private:
+#else
+public:
+#endif
+	Vector4D(Vector4D const& vOther);
+
+	// No assignment operators either...
+	Vector4D& operator=( Vector4D const& src );
+};
+
+const Vector4D vec4_origin( 0.0f, 0.0f, 0.0f, 0.0f );
+const Vector4D vec4_invalid( FLT_MAX, FLT_MAX, FLT_MAX, FLT_MAX );
+
+//-----------------------------------------------------------------------------
+// SSE optimized routines
+//-----------------------------------------------------------------------------
+
+class ALIGN16 Vector4DAligned : public Vector4D
+{
+public:
+	Vector4DAligned(void) {}
+	Vector4DAligned( vec_t X, vec_t Y, vec_t Z, vec_t W );
+
+	inline void Set( vec_t X, vec_t Y, vec_t Z, vec_t W );
+	inline void InitZero( void );
+
+	inline __m128 &AsM128() { return *(__m128*)&x; }
+	inline const __m128 &AsM128() const { return *(const __m128*)&x; } 
+
+private:
+	// No copy constructors allowed if we're in optimal mode
+	Vector4DAligned( Vector4DAligned const& vOther );
+
+	// No assignment operators either...
+	Vector4DAligned& operator=( Vector4DAligned const& src );
+} ALIGN16_POST;
+
+//-----------------------------------------------------------------------------
+// Vector4D related operations
+//-----------------------------------------------------------------------------
+
+// Vector4D clear
+void Vector4DClear( Vector4D& a );
+
+// Copy
+void Vector4DCopy( Vector4D const& src, Vector4D& dst );
+
+// Vector4D arithmetic
+void Vector4DAdd( Vector4D const& a, Vector4D const& b, Vector4D& result );
+void Vector4DSubtract( Vector4D const& a, Vector4D const& b, Vector4D& result );
+void Vector4DMultiply( Vector4D const& a, vec_t b, Vector4D& result );
+void Vector4DMultiply( Vector4D const& a, Vector4D const& b, Vector4D& result );
+void Vector4DDivide( Vector4D const& a, vec_t b, Vector4D& result );
+void Vector4DDivide( Vector4D const& a, Vector4D const& b, Vector4D& result );
+void Vector4DMA( Vector4D const& start, float s, Vector4D const& dir, Vector4D& result );
+
+// Vector4DAligned arithmetic
+void Vector4DMultiplyAligned( Vector4DAligned const& a, vec_t b, Vector4DAligned& result );
+
+
+#define Vector4DExpand( v ) (v).x, (v).y, (v).z, (v).w
+
+// Normalization
+vec_t Vector4DNormalize( Vector4D& v );
+
+// Length
+vec_t Vector4DLength( Vector4D const& v );
+
+// Dot Product
+vec_t DotProduct4D(Vector4D const& a, Vector4D const& b);
+
+// Linearly interpolate between two vectors
+void Vector4DLerp(Vector4D const& src1, Vector4D const& src2, vec_t t, Vector4D& dest );
+
+
+//-----------------------------------------------------------------------------
+//
+// Inlined Vector4D methods
+//
+//-----------------------------------------------------------------------------
+
+
+//-----------------------------------------------------------------------------
+// constructors
+//-----------------------------------------------------------------------------
+
+inline Vector4D::Vector4D(void)									
+{ 
+#ifdef _DEBUG
+	// Initialize to NAN to catch errors
+	x = y = z = w = VEC_T_NAN;
+#endif
+}
+
+inline Vector4D::Vector4D(vec_t X, vec_t Y, vec_t Z, vec_t W )
+{ 
+	x = X; y = Y; z = Z; w = W;
+	Assert( IsValid() );
+}
+
+inline Vector4D::Vector4D(const float *pFloat)					
+{
+	Assert( pFloat );
+	x = pFloat[0]; y = pFloat[1]; z = pFloat[2]; w = pFloat[3];	
+	Assert( IsValid() );
+}
+
+
+//-----------------------------------------------------------------------------
+// copy constructor
+//-----------------------------------------------------------------------------
+
+inline Vector4D::Vector4D(const Vector4D &vOther)					
+{ 
+	Assert( vOther.IsValid() );
+	x = vOther.x; y = vOther.y; z = vOther.z; w = vOther.w;
+}
+
+//-----------------------------------------------------------------------------
+// initialization
+//-----------------------------------------------------------------------------
+
+inline void Vector4D::Init( vec_t ix, vec_t iy, vec_t iz, vec_t iw )
+{ 
+	x = ix; y = iy; z = iz;	w = iw;
+	Assert( IsValid() );
+}
+
+inline void Vector4D::Random( vec_t minVal, vec_t maxVal )
+{
+	x = minVal + ((vec_t)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+	y = minVal + ((vec_t)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+	z = minVal + ((vec_t)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+	w = minVal + ((vec_t)rand() / VALVE_RAND_MAX) * (maxVal - minVal);
+}
+
+inline void Vector4DClear( Vector4D& a )
+{
+	a.x = a.y = a.z = a.w = 0.0f;
+}
+
+//-----------------------------------------------------------------------------
+// assignment
+//-----------------------------------------------------------------------------
+
+inline Vector4D& Vector4D::operator=(const Vector4D &vOther)	
+{
+	Assert( vOther.IsValid() );
+	x=vOther.x; y=vOther.y; z=vOther.z; w=vOther.w;
+	return *this; 
+}
+
+//-----------------------------------------------------------------------------
+// Array access
+//-----------------------------------------------------------------------------
+
+inline vec_t& Vector4D::operator[](int i)
+{
+	Assert( (i >= 0) && (i < 4) );
+	return ((vec_t*)this)[i];
+}
+
+inline vec_t Vector4D::operator[](int i) const
+{
+	Assert( (i >= 0) && (i < 4) );
+	return ((vec_t*)this)[i];
+}
+
+//-----------------------------------------------------------------------------
+// Cast to Vector and Vector2D...
+//-----------------------------------------------------------------------------
+
+inline Vector& Vector4D::AsVector3D()
+{
+	return *(Vector*)this;
+}
+
+inline Vector const& Vector4D::AsVector3D() const
+{
+	return *(Vector const*)this;
+}
+
+inline Vector2D& Vector4D::AsVector2D()
+{
+	return *(Vector2D*)this;
+}
+
+inline Vector2D const& Vector4D::AsVector2D() const
+{
+	return *(Vector2D const*)this;
+}
+
+//-----------------------------------------------------------------------------
+// Base address...
+//-----------------------------------------------------------------------------
+
+inline vec_t* Vector4D::Base()
+{
+	return (vec_t*)this;
+}
+
+inline vec_t const* Vector4D::Base() const
+{
+	return (vec_t const*)this;
+}
+
+//-----------------------------------------------------------------------------
+// IsValid?
+//-----------------------------------------------------------------------------
+
+inline bool Vector4D::IsValid() const
+{
+	return IsFinite(x) && IsFinite(y) && IsFinite(z) && IsFinite(w);
+}
+
+//-----------------------------------------------------------------------------
+// comparison
+//-----------------------------------------------------------------------------
+
+inline bool Vector4D::operator==( Vector4D const& src ) const
+{
+	Assert( src.IsValid() && IsValid() );
+	return (src.x == x) && (src.y == y) && (src.z == z) && (src.w == w);
+}
+
+inline bool Vector4D::operator!=( Vector4D const& src ) const
+{
+	Assert( src.IsValid() && IsValid() );
+	return (src.x != x) || (src.y != y) || (src.z != z) || (src.w != w);
+}
+
+
+//-----------------------------------------------------------------------------
+// Copy
+//-----------------------------------------------------------------------------
+
+inline void Vector4DCopy( Vector4D const& src, Vector4D& dst )
+{
+	Assert( src.IsValid() );
+	dst.x = src.x;
+	dst.y = src.y;
+	dst.z = src.z;
+	dst.w = src.w;
+}
+
+inline void	Vector4D::CopyToArray(float* rgfl) const		
+{ 
+	Assert( IsValid() );
+	Assert( rgfl );
+	rgfl[0] = x; rgfl[1] = y; rgfl[2] = z; rgfl[3] = w;
+}
+
+//-----------------------------------------------------------------------------
+// standard math operations
+//-----------------------------------------------------------------------------
+
+inline void Vector4D::Negate()
+{ 
+	Assert( IsValid() );
+	x = -x; y = -y; z = -z; w = -w;
+} 
+
+inline Vector4D& Vector4D::operator+=(const Vector4D& v)	
+{ 
+	Assert( IsValid() && v.IsValid() );
+	x+=v.x; y+=v.y; z += v.z; w += v.w;	
+	return *this;
+}
+
+inline Vector4D& Vector4D::operator-=(const Vector4D& v)	
+{ 
+	Assert( IsValid() && v.IsValid() );
+	x-=v.x; y-=v.y; z -= v.z; w -= v.w;
+	return *this;
+}
+
+inline Vector4D& Vector4D::operator*=(float fl)	
+{
+	x *= fl;
+	y *= fl;
+	z *= fl;
+	w *= fl;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector4D& Vector4D::operator*=(Vector4D const& v)	
+{ 
+	x *= v.x;
+	y *= v.y;
+	z *= v.z;
+	w *= v.w;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector4D& Vector4D::operator/=(float fl)	
+{
+	Assert( fl != 0.0f );
+	float oofl = 1.0f / fl;
+	x *= oofl;
+	y *= oofl;
+	z *= oofl;
+	w *= oofl;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline Vector4D& Vector4D::operator/=(Vector4D const& v)	
+{ 
+	Assert( v.x != 0.0f && v.y != 0.0f && v.z != 0.0f && v.w != 0.0f );
+	x /= v.x;
+	y /= v.y;
+	z /= v.z;
+	w /= v.w;
+	Assert( IsValid() );
+	return *this;
+}
+
+inline void Vector4DAdd( Vector4D const& a, Vector4D const& b, Vector4D& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x + b.x;
+	c.y = a.y + b.y;
+	c.z = a.z + b.z;
+	c.w = a.w + b.w;
+}
+
+inline void Vector4DSubtract( Vector4D const& a, Vector4D const& b, Vector4D& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x - b.x;
+	c.y = a.y - b.y;
+	c.z = a.z - b.z;
+	c.w = a.w - b.w;
+}
+
+inline void Vector4DMultiply( Vector4D const& a, vec_t b, Vector4D& c )
+{
+	Assert( a.IsValid() && IsFinite(b) );
+	c.x = a.x * b;
+	c.y = a.y * b;
+	c.z = a.z * b;
+	c.w = a.w * b;
+}
+
+inline void Vector4DMultiply( Vector4D const& a, Vector4D const& b, Vector4D& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+	c.x = a.x * b.x;
+	c.y = a.y * b.y;
+	c.z = a.z * b.z;
+	c.w = a.w * b.w;
+}
+
+inline void Vector4DDivide( Vector4D const& a, vec_t b, Vector4D& c )
+{
+	Assert( a.IsValid() );
+	Assert( b != 0.0f );
+	vec_t oob = 1.0f / b;
+	c.x = a.x * oob;
+	c.y = a.y * oob;
+	c.z = a.z * oob;
+	c.w = a.w * oob;
+}
+
+inline void Vector4DDivide( Vector4D const& a, Vector4D const& b, Vector4D& c )
+{
+	Assert( a.IsValid() );
+	Assert( (b.x != 0.0f) && (b.y != 0.0f) && (b.z != 0.0f) && (b.w != 0.0f) );
+	c.x = a.x / b.x;
+	c.y = a.y / b.y;
+	c.z = a.z / b.z;
+	c.w = a.w / b.w;
+}
+
+inline void Vector4DMA( Vector4D const& start, float s, Vector4D const& dir, Vector4D& result )
+{
+	Assert( start.IsValid() && IsFinite(s) && dir.IsValid() );
+	result.x = start.x + s*dir.x;
+	result.y = start.y + s*dir.y;
+	result.z = start.z + s*dir.z;
+	result.w = start.w + s*dir.w;
+}
+
+// FIXME: Remove
+// For backwards compatability
+inline void	Vector4D::MulAdd(Vector4D const& a, Vector4D const& b, float scalar)
+{
+	x = a.x + b.x * scalar;
+	y = a.y + b.y * scalar;
+	z = a.z + b.z * scalar;
+	w = a.w + b.w * scalar;
+}
+
+inline void Vector4DLerp(const Vector4D& src1, const Vector4D& src2, vec_t t, Vector4D& dest )
+{
+	dest[0] = src1[0] + (src2[0] - src1[0]) * t;
+	dest[1] = src1[1] + (src2[1] - src1[1]) * t;
+	dest[2] = src1[2] + (src2[2] - src1[2]) * t;
+	dest[3] = src1[3] + (src2[3] - src1[3]) * t;
+}
+
+//-----------------------------------------------------------------------------
+// dot, cross
+//-----------------------------------------------------------------------------
+
+inline vec_t DotProduct4D(const Vector4D& a, const Vector4D& b) 
+{ 
+	Assert( a.IsValid() && b.IsValid() );
+	return( a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w ); 
+}
+
+// for backwards compatability
+inline vec_t Vector4D::Dot( Vector4D const& vOther ) const
+{
+	return DotProduct4D( *this, vOther );
+}
+
+
+//-----------------------------------------------------------------------------
+// length
+//-----------------------------------------------------------------------------
+
+inline vec_t Vector4DLength( Vector4D const& v )
+{				   
+	Assert( v.IsValid() );
+	return (vec_t)FastSqrt(v.x*v.x + v.y*v.y + v.z*v.z + v.w*v.w);		
+}
+
+inline vec_t Vector4D::LengthSqr(void) const	
+{ 
+	Assert( IsValid() );
+	return (x*x + y*y + z*z + w*w);		
+}
+
+inline vec_t Vector4D::Length(void) const	
+{
+	return Vector4DLength( *this );
+}
+
+
+//-----------------------------------------------------------------------------
+// Normalization
+//-----------------------------------------------------------------------------
+
+// FIXME: Can't use until we're un-macroed in mathlib.h
+inline vec_t Vector4DNormalize( Vector4D& v )
+{
+	Assert( v.IsValid() );
+	vec_t l = v.Length();
+	if (l != 0.0f)
+	{
+		v /= l;
+	}
+	else
+	{
+		v.x = v.y = v.z = v.w = 0.0f;
+	}
+	return l;
+}
+
+//-----------------------------------------------------------------------------
+// Get the distance from this Vector4D to the other one 
+//-----------------------------------------------------------------------------
+
+inline vec_t Vector4D::DistTo(const Vector4D &vOther) const
+{
+	Vector4D delta;
+	Vector4DSubtract( *this, vOther, delta );
+	return delta.Length();
+}
+
+inline vec_t Vector4D::DistToSqr(const Vector4D &vOther) const
+{
+	Vector4D delta;
+	Vector4DSubtract( *this, vOther, delta );
+	return delta.LengthSqr();
+}
+
+
+//-----------------------------------------------------------------------------
+// Vector4DAligned routines
+//-----------------------------------------------------------------------------
+
+inline Vector4DAligned::Vector4DAligned( vec_t X, vec_t Y, vec_t Z, vec_t W )
+{ 
+	x = X; y = Y; z = Z; w = W;
+	Assert( IsValid() );
+}
+
+inline void Vector4DAligned::Set( vec_t X, vec_t Y, vec_t Z, vec_t W )
+{ 
+	x = X; y = Y; z = Z; w = W;
+	Assert( IsValid() );
+}
+
+inline void Vector4DAligned::InitZero( void )
+{ 
+#if !defined( _X360 )
+	this->AsM128() = _mm_set1_ps( 0.0f );
+#else
+	this->AsM128() = __vspltisw( 0 );
+#endif
+	Assert( IsValid() );
+}
+
+inline void Vector4DMultiplyAligned( Vector4DAligned const& a, Vector4DAligned const& b, Vector4DAligned& c )
+{
+	Assert( a.IsValid() && b.IsValid() );
+#if !defined( _X360 )
+	c.x = a.x * b.x;
+	c.y = a.y * b.y;
+	c.z = a.z * b.z;
+	c.w = a.w * b.w;
+#else
+	c.AsM128() = __vmulfp( a.AsM128(), b.AsM128() );
+#endif
+}
+
+inline void Vector4DWeightMAD( vec_t w, Vector4DAligned const& vInA, Vector4DAligned& vOutA, Vector4DAligned const& vInB, Vector4DAligned& vOutB )
+{
+	Assert( vInA.IsValid() && vInB.IsValid() && IsFinite(w) );
+
+#if !defined( _X360 )
+	vOutA.x += vInA.x * w;
+	vOutA.y += vInA.y * w;
+	vOutA.z += vInA.z * w;
+	vOutA.w += vInA.w * w;
+
+	vOutB.x += vInB.x * w;
+	vOutB.y += vInB.y * w;
+	vOutB.z += vInB.z * w;
+	vOutB.w += vInB.w * w;
+#else
+    __vector4 temp;
+
+    temp = __lvlx( &w, 0 );
+    temp = __vspltw( temp, 0 );
+
+	vOutA.AsM128() = __vmaddfp( vInA.AsM128(), temp, vOutA.AsM128() );
+	vOutB.AsM128() = __vmaddfp( vInB.AsM128(), temp, vOutB.AsM128() );
+#endif
+}
+
+inline void Vector4DWeightMADSSE( vec_t w, Vector4DAligned const& vInA, Vector4DAligned& vOutA, Vector4DAligned const& vInB, Vector4DAligned& vOutB )
+{
+	Assert( vInA.IsValid() && vInB.IsValid() && IsFinite(w) );
+
+#if !defined( _X360 )
+	// Replicate scalar float out to 4 components
+    __m128 packed = _mm_set1_ps( w );
+
+	// 4D SSE Vector MAD
+	vOutA.AsM128() = _mm_add_ps( vOutA.AsM128(), _mm_mul_ps( vInA.AsM128(), packed ) );
+	vOutB.AsM128() = _mm_add_ps( vOutB.AsM128(), _mm_mul_ps( vInB.AsM128(), packed ) );
+#else
+    __vector4 temp;
+
+    temp = __lvlx( &w, 0 );
+    temp = __vspltw( temp, 0 );
+
+	vOutA.AsM128() = __vmaddfp( vInA.AsM128(), temp, vOutA.AsM128() );
+	vOutB.AsM128() = __vmaddfp( vInB.AsM128(), temp, vOutB.AsM128() );
+#endif
+}
+
+#endif // VECTOR4D_H
+

public/mathlib/vmatrix.h

@@ -0,0 +1,947 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $NoKeywords: $
+//
+//=============================================================================//
+//
+// VMatrix always postmultiply vectors as in Ax = b.
+// Given a set of basis vectors ((F)orward, (L)eft, (U)p), and a (T)ranslation, 
+// a matrix to transform a vector into that space looks like this:
+// Fx Lx Ux Tx
+// Fy Ly Uy Ty
+// Fz Lz Uz Tz
+// 0   0  0  1
+
+// Note that concatenating matrices needs to multiply them in reverse order.
+// ie: if I want to apply matrix A, B, then C, the equation needs to look like this:
+// C * B * A * v
+// ie:
+// v = A * v;
+// v = B * v;
+// v = C * v;
+//=============================================================================
+
+#ifndef VMATRIX_H
+#define VMATRIX_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include <string.h>
+#include "mathlib/vector.h"
+#include "mathlib/vplane.h"
+#include "mathlib/vector4d.h"
+#include "mathlib/mathlib.h"
+
+struct cplane_t;
+
+
+class VMatrix
+{
+public:
+
+	VMatrix();
+	VMatrix(
+		vec_t m00, vec_t m01, vec_t m02, vec_t m03,
+		vec_t m10, vec_t m11, vec_t m12, vec_t m13,
+		vec_t m20, vec_t m21, vec_t m22, vec_t m23,
+		vec_t m30, vec_t m31, vec_t m32, vec_t m33
+		);
+
+	// Creates a matrix where the X axis = forward
+	// the Y axis = left, and the Z axis = up
+	VMatrix( const Vector& forward, const Vector& left, const Vector& up );
+	VMatrix( const Vector& forward, const Vector& left, const Vector& up, const Vector& translation );
+	
+	// Construct from a 3x4 matrix
+	VMatrix( const matrix3x4_t& matrix3x4 );
+
+	// Set the values in the matrix.
+	void		Init( 
+		vec_t m00, vec_t m01, vec_t m02, vec_t m03,
+		vec_t m10, vec_t m11, vec_t m12, vec_t m13,
+		vec_t m20, vec_t m21, vec_t m22, vec_t m23,
+		vec_t m30, vec_t m31, vec_t m32, vec_t m33 
+		);
+
+
+	// Initialize from a 3x4
+	void		Init( const matrix3x4_t& matrix3x4 );
+
+	// array access
+	inline float* operator[](int i)
+	{ 
+		return m[i]; 
+	}
+
+	inline const float* operator[](int i) const
+	{ 
+		return m[i]; 
+	}
+
+	// Get a pointer to m[0][0]
+	inline float *Base()
+	{
+		return &m[0][0];
+	}
+
+	inline const float *Base() const
+	{
+		return &m[0][0];
+	}
+
+	void		SetLeft(const Vector &vLeft);
+	void		SetUp(const Vector &vUp);
+	void		SetForward(const Vector &vForward);
+
+	void		GetBasisVectors(Vector &vForward, Vector &vLeft, Vector &vUp) const;
+	void		SetBasisVectors(const Vector &vForward, const Vector &vLeft, const Vector &vUp);
+
+	// Get/set the translation.
+	Vector &	GetTranslation( Vector &vTrans ) const;
+	void		SetTranslation(const Vector &vTrans);
+
+	void		PreTranslate(const Vector &vTrans);
+	void		PostTranslate(const Vector &vTrans);
+
+	const matrix3x4_t& As3x4() const;
+	void		CopyFrom3x4( const matrix3x4_t &m3x4 );
+	void		Set3x4( matrix3x4_t& matrix3x4 ) const;
+
+	bool		operator==( const VMatrix& src ) const;
+	bool		operator!=( const VMatrix& src ) const { return !( *this == src ); }
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// Access the basis vectors.
+	Vector		GetLeft() const;
+	Vector		GetUp() const;
+	Vector		GetForward() const;
+	Vector		GetTranslation() const;
+#endif
+
+
+// Matrix->vector operations.
+public:
+	// Multiply by a 3D vector (same as operator*).
+	void		V3Mul(const Vector &vIn, Vector &vOut) const;
+
+	// Multiply by a 4D vector.
+	void		V4Mul(const Vector4D &vIn, Vector4D &vOut) const;
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// Applies the rotation (ignores translation in the matrix). (This just calls VMul3x3).
+	Vector		ApplyRotation(const Vector &vVec) const;
+
+	// Multiply by a vector (divides by w, assumes input w is 1).
+	Vector		operator*(const Vector &vVec) const;
+
+	// Multiply by the upper 3x3 part of the matrix (ie: only apply rotation).
+	Vector		VMul3x3(const Vector &vVec) const;
+
+	// Apply the inverse (transposed) rotation (only works on pure rotation matrix)
+	Vector		VMul3x3Transpose(const Vector &vVec) const;
+
+	// Multiply by the upper 3 rows.
+	Vector		VMul4x3(const Vector &vVec) const;
+
+	// Apply the inverse (transposed) transformation (only works on pure rotation/translation)
+	Vector		VMul4x3Transpose(const Vector &vVec) const;
+#endif
+
+
+// Matrix->plane operations.
+public:
+	// Transform the plane. The matrix can only contain translation and rotation.
+	void		TransformPlane( const VPlane &inPlane, VPlane &outPlane ) const;
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// Just calls TransformPlane and returns the result.
+	VPlane		operator*(const VPlane &thePlane) const;
+#endif
+
+// Matrix->matrix operations.
+public:
+
+	VMatrix&	operator=(const VMatrix &mOther);
+	
+	// Multiply two matrices (out = this * vm).
+	void		MatrixMul( const VMatrix &vm, VMatrix &out ) const;
+
+	// Add two matrices.
+	const VMatrix& operator+=(const VMatrix &other);
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// Just calls MatrixMul and returns the result.	
+	VMatrix		operator*(const VMatrix &mOther) const;
+
+	// Add/Subtract two matrices.
+	VMatrix		operator+(const VMatrix &other) const;
+	VMatrix		operator-(const VMatrix &other) const;
+
+	// Negation.
+	VMatrix		operator-() const;
+
+	// Return inverse matrix. Be careful because the results are undefined 
+	// if the matrix doesn't have an inverse (ie: InverseGeneral returns false).
+	VMatrix		operator~() const;
+#endif
+
+// Matrix operations.
+public:
+	// Set to identity.
+	void		Identity();
+
+	bool		IsIdentity() const;
+
+	// Setup a matrix for origin and angles.
+	void		SetupMatrixOrgAngles( const Vector &origin, const QAngle &vAngles );
+	
+	// Setup a matrix for angles and no translation.
+	void		SetupMatrixAngles( const QAngle &vAngles );
+
+	// General inverse. This may fail so check the return!
+	bool		InverseGeneral(VMatrix &vInverse) const;
+	
+	// Does a fast inverse, assuming the matrix only contains translation and rotation.
+	void		InverseTR( VMatrix &mRet ) const;
+
+	// Usually used for debug checks. Returns true if the upper 3x3 contains
+	// unit vectors and they are all orthogonal.
+	bool		IsRotationMatrix() const;
+	
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// This calls the other InverseTR and returns the result.
+	VMatrix		InverseTR() const;
+
+	// Get the scale of the matrix's basis vectors.
+	Vector		GetScale() const;
+
+	// (Fast) multiply by a scaling matrix setup from vScale.
+	VMatrix		Scale(const Vector &vScale);	
+
+	// Normalize the basis vectors.
+	VMatrix		NormalizeBasisVectors() const;
+
+	// Transpose.
+	VMatrix		Transpose() const;
+
+	// Transpose upper-left 3x3.
+	VMatrix		Transpose3x3() const;
+#endif
+
+public:
+	// The matrix.
+	vec_t		m[4][4];
+};
+
+
+
+//-----------------------------------------------------------------------------
+// Helper functions.
+//-----------------------------------------------------------------------------
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+// Setup an identity matrix.
+VMatrix		SetupMatrixIdentity();
+
+// Setup as a scaling matrix.
+VMatrix		SetupMatrixScale(const Vector &vScale);
+
+// Setup a translation matrix.
+VMatrix		SetupMatrixTranslation(const Vector &vTranslation);
+
+// Setup a matrix to reflect around the plane.
+VMatrix		SetupMatrixReflection(const VPlane &thePlane);
+
+// Setup a matrix to project from vOrigin onto thePlane.
+VMatrix		SetupMatrixProjection(const Vector &vOrigin, const VPlane &thePlane);
+
+// Setup a matrix to rotate the specified amount around the specified axis.
+VMatrix		SetupMatrixAxisRot(const Vector &vAxis, vec_t fDegrees);
+
+// Setup a matrix from euler angles. Just sets identity and calls MatrixAngles.
+VMatrix		SetupMatrixAngles(const QAngle &vAngles);
+
+// Setup a matrix for origin and angles.
+VMatrix		SetupMatrixOrgAngles(const Vector &origin, const QAngle &vAngles);
+
+#endif
+
+#define VMatToString(mat)	(static_cast<const char *>(CFmtStr("[ (%f, %f, %f), (%f, %f, %f), (%f, %f, %f), (%f, %f, %f) ]", mat.m[0][0], mat.m[0][1], mat.m[0][2], mat.m[0][3], mat.m[1][0], mat.m[1][1], mat.m[1][2], mat.m[1][3], mat.m[2][0], mat.m[2][1], mat.m[2][2], mat.m[2][3], mat.m[3][0], mat.m[3][1], mat.m[3][2], mat.m[3][3] ))) // ** Note: this generates a temporary, don't hold reference!
+
+//-----------------------------------------------------------------------------
+// Returns the point at the intersection on the 3 planes.
+// Returns false if it can't be solved (2 or more planes are parallel).
+//-----------------------------------------------------------------------------
+bool PlaneIntersection( const VPlane &vp1, const VPlane &vp2, const VPlane &vp3, Vector &vOut );
+
+
+//-----------------------------------------------------------------------------
+// These methods are faster. Use them if you want faster code
+//-----------------------------------------------------------------------------
+void MatrixSetIdentity( VMatrix &dst );
+void MatrixTranspose( const VMatrix& src, VMatrix& dst );
+void MatrixCopy( const VMatrix& src, VMatrix& dst );
+void MatrixMultiply( const VMatrix& src1, const VMatrix& src2, VMatrix& dst );
+
+// Accessors
+void MatrixGetColumn( const VMatrix &src, int nCol, Vector *pColumn );
+void MatrixSetColumn( VMatrix &src, int nCol, const Vector &column );
+void MatrixGetRow( const VMatrix &src, int nCol, Vector *pColumn );
+void MatrixSetRow( VMatrix &src, int nCol, const Vector &column );
+
+// Vector3DMultiply treats src2 as if it's a direction vector
+void Vector3DMultiply( const VMatrix& src1, const Vector& src2, Vector& dst );
+
+// Vector3DMultiplyPosition treats src2 as if it's a point (adds the translation)
+inline void Vector3DMultiplyPosition( const VMatrix& src1, const VectorByValue src2, Vector& dst );
+
+// Vector3DMultiplyPositionProjective treats src2 as if it's a point 
+// and does the perspective divide at the end
+void Vector3DMultiplyPositionProjective( const VMatrix& src1, const Vector &src2, Vector& dst );
+
+// Vector3DMultiplyPosition treats src2 as if it's a direction 
+// and does the perspective divide at the end
+// NOTE: src1 had better be an inverse transpose to use this correctly
+void Vector3DMultiplyProjective( const VMatrix& src1, const Vector &src2, Vector& dst );
+
+void Vector4DMultiply( const VMatrix& src1, const Vector4D& src2, Vector4D& dst );
+
+// Same as Vector4DMultiply except that src2 has an implicit W of 1
+void Vector4DMultiplyPosition( const VMatrix& src1, const Vector &src2, Vector4D& dst );
+
+// Multiplies the vector by the transpose of the matrix
+void Vector3DMultiplyTranspose( const VMatrix& src1, const Vector& src2, Vector& dst );
+void Vector4DMultiplyTranspose( const VMatrix& src1, const Vector4D& src2, Vector4D& dst );
+
+// Transform a plane
+void MatrixTransformPlane( const VMatrix &src, const cplane_t &inPlane, cplane_t &outPlane );
+
+// Transform a plane that has an axis-aligned normal
+void MatrixTransformAxisAlignedPlane( const VMatrix &src, int nDim, float flSign, float flDist, cplane_t &outPlane );
+
+void MatrixBuildTranslation( VMatrix& dst, float x, float y, float z );
+void MatrixBuildTranslation( VMatrix& dst, const Vector &translation );
+
+inline void MatrixTranslate( VMatrix& dst, const Vector &translation )
+{
+	VMatrix matTranslation, temp;
+	MatrixBuildTranslation( matTranslation, translation );
+	MatrixMultiply( dst, matTranslation, temp );
+	dst = temp;
+}
+
+
+void MatrixBuildRotationAboutAxis( VMatrix& dst, const Vector& vAxisOfRot, float angleDegrees );
+void MatrixBuildRotateZ( VMatrix& dst, float angleDegrees );
+
+inline void MatrixRotate( VMatrix& dst, const Vector& vAxisOfRot, float angleDegrees )
+{
+	VMatrix rotation, temp;
+	MatrixBuildRotationAboutAxis( rotation, vAxisOfRot, angleDegrees );
+	MatrixMultiply( dst, rotation, temp );
+	dst = temp;
+}
+
+// Builds a rotation matrix that rotates one direction vector into another
+void MatrixBuildRotation( VMatrix &dst, const Vector& initialDirection, const Vector& finalDirection );
+
+// Builds a scale matrix
+void MatrixBuildScale( VMatrix &dst, float x, float y, float z );
+void MatrixBuildScale( VMatrix &dst, const Vector& scale );
+
+// Build a perspective matrix.
+// zNear and zFar are assumed to be positive.
+// You end up looking down positive Z, X is to the right, Y is up.
+// X range: [0..1]
+// Y range: [0..1]
+// Z range: [0..1]
+void MatrixBuildPerspective( VMatrix &dst, float fovX, float fovY, float zNear, float zFar );
+
+//-----------------------------------------------------------------------------
+// Given a projection matrix, take the extremes of the space in transformed into world space and
+// get a bounding box.
+//-----------------------------------------------------------------------------
+void CalculateAABBFromProjectionMatrix( const VMatrix &worldToVolume, Vector *pMins, Vector *pMaxs );
+
+//-----------------------------------------------------------------------------
+// Given a projection matrix, take the extremes of the space in transformed into world space and
+// get a bounding sphere.
+//-----------------------------------------------------------------------------
+void CalculateSphereFromProjectionMatrix( const VMatrix &worldToVolume, Vector *pCenter, float *pflRadius );
+
+//-----------------------------------------------------------------------------
+// Given an inverse projection matrix, take the extremes of the space in transformed into world space and
+// get a bounding box.
+//-----------------------------------------------------------------------------
+void CalculateAABBFromProjectionMatrixInverse( const VMatrix &volumeToWorld, Vector *pMins, Vector *pMaxs );
+
+//-----------------------------------------------------------------------------
+// Given an inverse projection matrix, take the extremes of the space in transformed into world space and
+// get a bounding sphere.
+//-----------------------------------------------------------------------------
+void CalculateSphereFromProjectionMatrixInverse( const VMatrix &volumeToWorld, Vector *pCenter, float *pflRadius );
+
+//-----------------------------------------------------------------------------
+// Calculate frustum planes given a clip->world space transform.
+//-----------------------------------------------------------------------------
+void FrustumPlanesFromMatrix( const VMatrix &clipToWorld, Frustum_t &frustum );
+
+//-----------------------------------------------------------------------------
+// Setup a matrix from euler angles. 
+//-----------------------------------------------------------------------------
+void MatrixFromAngles( const QAngle& vAngles, VMatrix& dst );
+
+//-----------------------------------------------------------------------------
+// Creates euler angles from a matrix 
+//-----------------------------------------------------------------------------
+void MatrixToAngles( const VMatrix& src, QAngle& vAngles );
+
+//-----------------------------------------------------------------------------
+// Does a fast inverse, assuming the matrix only contains translation and rotation.
+//-----------------------------------------------------------------------------
+void MatrixInverseTR( const VMatrix& src, VMatrix &dst );
+
+//-----------------------------------------------------------------------------
+// Inverts any matrix at all
+//-----------------------------------------------------------------------------
+bool MatrixInverseGeneral(const VMatrix& src, VMatrix& dst);
+
+//-----------------------------------------------------------------------------
+// Computes the inverse transpose
+//-----------------------------------------------------------------------------
+void MatrixInverseTranspose( const VMatrix& src, VMatrix& dst );
+
+
+
+//-----------------------------------------------------------------------------
+// VMatrix inlines.
+//-----------------------------------------------------------------------------
+inline VMatrix::VMatrix()
+{
+}
+
+inline VMatrix::VMatrix(
+	vec_t m00, vec_t m01, vec_t m02, vec_t m03,
+	vec_t m10, vec_t m11, vec_t m12, vec_t m13,
+	vec_t m20, vec_t m21, vec_t m22, vec_t m23,
+	vec_t m30, vec_t m31, vec_t m32, vec_t m33)
+{
+	Init(
+		m00, m01, m02, m03,
+		m10, m11, m12, m13,
+		m20, m21, m22, m23,
+		m30, m31, m32, m33
+		);
+}
+
+
+inline VMatrix::VMatrix( const matrix3x4_t& matrix3x4 )
+{
+	Init( matrix3x4 );
+}
+
+
+//-----------------------------------------------------------------------------
+// Creates a matrix where the X axis = forward
+// the Y axis = left, and the Z axis = up
+//-----------------------------------------------------------------------------
+inline VMatrix::VMatrix( const Vector& xAxis, const Vector& yAxis, const Vector& zAxis )
+{
+	Init(
+		xAxis.x, yAxis.x, zAxis.x, 0.0f,
+		xAxis.y, yAxis.y, zAxis.y, 0.0f,
+		xAxis.z, yAxis.z, zAxis.z, 0.0f,
+		0.0f, 0.0f, 0.0f, 1.0f
+		);
+}
+
+inline VMatrix::VMatrix( const Vector& xAxis, const Vector& yAxis, const Vector& zAxis, const Vector& translation )
+{
+	Init(
+		xAxis.x, yAxis.x, zAxis.x, translation.x,
+		xAxis.y, yAxis.y, zAxis.y, translation.y,
+		xAxis.z, yAxis.z, zAxis.z, translation.z,
+		0.0f, 0.0f, 0.0f, 1.0f
+		);
+}
+
+
+inline void VMatrix::Init(
+	vec_t m00, vec_t m01, vec_t m02, vec_t m03,
+	vec_t m10, vec_t m11, vec_t m12, vec_t m13,
+	vec_t m20, vec_t m21, vec_t m22, vec_t m23,
+	vec_t m30, vec_t m31, vec_t m32, vec_t m33
+	)
+{
+	m[0][0] = m00;
+	m[0][1] = m01;
+	m[0][2] = m02;
+	m[0][3] = m03;
+
+	m[1][0] = m10;
+	m[1][1] = m11;
+	m[1][2] = m12;
+	m[1][3] = m13;
+
+	m[2][0] = m20;
+	m[2][1] = m21;
+	m[2][2] = m22;
+	m[2][3] = m23;
+
+	m[3][0] = m30;
+	m[3][1] = m31;
+	m[3][2] = m32;
+	m[3][3] = m33;
+}
+
+
+//-----------------------------------------------------------------------------
+// Initialize from a 3x4
+//-----------------------------------------------------------------------------
+inline void VMatrix::Init( const matrix3x4_t& matrix3x4 )
+{
+	memcpy(m, matrix3x4.Base(), sizeof( matrix3x4_t ) );
+
+	m[3][0] = 0.0f;
+	m[3][1] = 0.0f;
+	m[3][2] = 0.0f;
+	m[3][3] = 1.0f;	
+}
+
+
+//-----------------------------------------------------------------------------
+// Methods related to the basis vectors of the matrix
+//-----------------------------------------------------------------------------
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline Vector VMatrix::GetForward() const
+{
+	return Vector(m[0][0], m[1][0], m[2][0]);
+}
+
+inline Vector VMatrix::GetLeft() const
+{
+	return Vector(m[0][1], m[1][1], m[2][1]);
+}
+
+inline Vector VMatrix::GetUp() const
+{
+	return Vector(m[0][2], m[1][2], m[2][2]);
+}
+
+#endif
+
+inline void VMatrix::SetForward(const Vector &vForward)
+{
+	m[0][0] = vForward.x;
+	m[1][0] = vForward.y;
+	m[2][0] = vForward.z;
+}
+
+inline void VMatrix::SetLeft(const Vector &vLeft)
+{
+	m[0][1] = vLeft.x;
+	m[1][1] = vLeft.y;
+	m[2][1] = vLeft.z;
+}
+
+inline void VMatrix::SetUp(const Vector &vUp)
+{
+	m[0][2] = vUp.x;
+	m[1][2] = vUp.y;
+	m[2][2] = vUp.z;
+}
+
+inline void VMatrix::GetBasisVectors(Vector &vForward, Vector &vLeft, Vector &vUp) const
+{
+	vForward.Init( m[0][0], m[1][0], m[2][0] );
+	vLeft.Init( m[0][1], m[1][1], m[2][1] );
+	vUp.Init( m[0][2], m[1][2], m[2][2] );
+}
+
+inline void VMatrix::SetBasisVectors(const Vector &vForward, const Vector &vLeft, const Vector &vUp)
+{
+	SetForward(vForward);
+	SetLeft(vLeft);
+	SetUp(vUp);
+}
+
+
+//-----------------------------------------------------------------------------
+// Methods related to the translation component of the matrix
+//-----------------------------------------------------------------------------
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline Vector VMatrix::GetTranslation() const
+{
+	return Vector(m[0][3], m[1][3], m[2][3]);
+}
+
+#endif
+
+inline Vector& VMatrix::GetTranslation( Vector &vTrans ) const
+{
+	vTrans.x = m[0][3];
+	vTrans.y = m[1][3];
+	vTrans.z = m[2][3];
+	return vTrans;
+}
+
+inline void VMatrix::SetTranslation(const Vector &vTrans)
+{
+	m[0][3] = vTrans.x;
+	m[1][3] = vTrans.y;
+	m[2][3] = vTrans.z;
+}
+
+		  
+//-----------------------------------------------------------------------------
+// appply translation to this matrix in the input space
+//-----------------------------------------------------------------------------
+inline void VMatrix::PreTranslate(const Vector &vTrans)
+{
+	Vector tmp;
+	Vector3DMultiplyPosition( *this, vTrans, tmp );
+	m[0][3] = tmp.x;
+	m[1][3] = tmp.y;
+	m[2][3] = tmp.z;
+}
+
+
+//-----------------------------------------------------------------------------
+// appply translation to this matrix in the output space
+//-----------------------------------------------------------------------------
+inline void VMatrix::PostTranslate(const Vector &vTrans)
+{
+	m[0][3] += vTrans.x;
+	m[1][3] += vTrans.y;
+	m[2][3] += vTrans.z;
+}
+
+inline const matrix3x4_t& VMatrix::As3x4() const
+{
+	return *((const matrix3x4_t*)this);
+}
+
+inline void VMatrix::CopyFrom3x4( const matrix3x4_t &m3x4 )
+{
+	memcpy( m, m3x4.Base(), sizeof( matrix3x4_t ) );
+	m[3][0] = m[3][1] = m[3][2] = 0;
+	m[3][3] = 1;
+}
+
+inline void	VMatrix::Set3x4( matrix3x4_t& matrix3x4 ) const
+{
+	memcpy(matrix3x4.Base(), m, sizeof( matrix3x4_t ) );
+}
+
+
+//-----------------------------------------------------------------------------
+// Matrix math operations
+//-----------------------------------------------------------------------------
+inline const VMatrix& VMatrix::operator+=(const VMatrix &other)
+{
+	for(int i=0; i < 4; i++)
+	{
+		for(int j=0; j < 4; j++)
+		{
+			m[i][j] += other.m[i][j];
+		}
+	}
+
+	return *this;
+}
+
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline VMatrix VMatrix::operator+(const VMatrix &other) const
+{
+	VMatrix ret;
+	for(int i=0; i < 16; i++)
+	{
+		((float*)ret.m)[i] = ((float*)m)[i] + ((float*)other.m)[i];
+	}
+	return ret;
+}
+
+inline VMatrix VMatrix::operator-(const VMatrix &other) const
+{
+	VMatrix ret;
+
+	for(int i=0; i < 4; i++)
+	{
+		for(int j=0; j < 4; j++)
+		{
+			ret.m[i][j] = m[i][j] - other.m[i][j];
+		}
+	}
+
+	return ret;
+}
+
+inline VMatrix VMatrix::operator-() const
+{
+	VMatrix ret;
+	for( int i=0; i < 16; i++ )
+	{
+		((float*)ret.m)[i] = ((float*)m)[i];
+	}
+	return ret;
+}
+
+#endif // VECTOR_NO_SLOW_OPERATIONS
+
+
+//-----------------------------------------------------------------------------
+// Vector transformation
+//-----------------------------------------------------------------------------
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline Vector VMatrix::operator*(const Vector &vVec) const
+{
+	Vector vRet;
+	vRet.x = m[0][0]*vVec.x + m[0][1]*vVec.y + m[0][2]*vVec.z + m[0][3];
+	vRet.y = m[1][0]*vVec.x + m[1][1]*vVec.y + m[1][2]*vVec.z + m[1][3];
+	vRet.z = m[2][0]*vVec.x + m[2][1]*vVec.y + m[2][2]*vVec.z + m[2][3];
+
+	return vRet;
+}
+
+inline Vector VMatrix::VMul4x3(const Vector &vVec) const
+{
+	Vector vResult;
+	Vector3DMultiplyPosition( *this, vVec, vResult );
+	return vResult;
+}
+
+
+inline Vector VMatrix::VMul4x3Transpose(const Vector &vVec) const
+{
+	Vector tmp = vVec;
+	tmp.x -= m[0][3];
+	tmp.y -= m[1][3];
+	tmp.z -= m[2][3];
+
+	return Vector(
+		m[0][0]*tmp.x + m[1][0]*tmp.y + m[2][0]*tmp.z,
+		m[0][1]*tmp.x + m[1][1]*tmp.y + m[2][1]*tmp.z,
+		m[0][2]*tmp.x + m[1][2]*tmp.y + m[2][2]*tmp.z
+		);
+}
+
+inline Vector VMatrix::VMul3x3(const Vector &vVec) const
+{
+	return Vector(
+		m[0][0]*vVec.x + m[0][1]*vVec.y + m[0][2]*vVec.z,
+		m[1][0]*vVec.x + m[1][1]*vVec.y + m[1][2]*vVec.z,
+		m[2][0]*vVec.x + m[2][1]*vVec.y + m[2][2]*vVec.z
+		);
+}
+
+inline Vector VMatrix::VMul3x3Transpose(const Vector &vVec) const
+{
+	return Vector(
+		m[0][0]*vVec.x + m[1][0]*vVec.y + m[2][0]*vVec.z,
+		m[0][1]*vVec.x + m[1][1]*vVec.y + m[2][1]*vVec.z,
+		m[0][2]*vVec.x + m[1][2]*vVec.y + m[2][2]*vVec.z
+		);
+}
+
+#endif // VECTOR_NO_SLOW_OPERATIONS
+
+
+inline void VMatrix::V3Mul(const Vector &vIn, Vector &vOut) const
+{
+	vec_t rw;
+
+	rw = 1.0f / (m[3][0]*vIn.x + m[3][1]*vIn.y + m[3][2]*vIn.z + m[3][3]);
+	vOut.x = (m[0][0]*vIn.x + m[0][1]*vIn.y + m[0][2]*vIn.z + m[0][3]) * rw;
+	vOut.y = (m[1][0]*vIn.x + m[1][1]*vIn.y + m[1][2]*vIn.z + m[1][3]) * rw;
+	vOut.z = (m[2][0]*vIn.x + m[2][1]*vIn.y + m[2][2]*vIn.z + m[2][3]) * rw;
+}
+
+inline void VMatrix::V4Mul(const Vector4D &vIn, Vector4D &vOut) const
+{
+	vOut[0] = m[0][0]*vIn[0] + m[0][1]*vIn[1] + m[0][2]*vIn[2] + m[0][3]*vIn[3];
+	vOut[1] = m[1][0]*vIn[0] + m[1][1]*vIn[1] + m[1][2]*vIn[2] + m[1][3]*vIn[3];
+	vOut[2] = m[2][0]*vIn[0] + m[2][1]*vIn[1] + m[2][2]*vIn[2] + m[2][3]*vIn[3];
+	vOut[3] = m[3][0]*vIn[0] + m[3][1]*vIn[1] + m[3][2]*vIn[2] + m[3][3]*vIn[3];
+}
+
+
+//-----------------------------------------------------------------------------
+// Plane transformation
+//-----------------------------------------------------------------------------
+inline void	VMatrix::TransformPlane( const VPlane &inPlane, VPlane &outPlane ) const
+{
+	Vector vTrans;
+	Vector3DMultiply( *this, inPlane.m_Normal, outPlane.m_Normal );
+	outPlane.m_Dist = inPlane.m_Dist * DotProduct( outPlane.m_Normal, outPlane.m_Normal );
+	outPlane.m_Dist += DotProduct( outPlane.m_Normal, GetTranslation( vTrans ) );
+}
+
+
+//-----------------------------------------------------------------------------
+// Other random stuff
+//-----------------------------------------------------------------------------
+inline void VMatrix::Identity()
+{
+	MatrixSetIdentity( *this );
+}
+
+
+inline bool VMatrix::IsIdentity() const
+{
+	return 
+		m[0][0] == 1.0f && m[0][1] == 0.0f && m[0][2] == 0.0f && m[0][3] == 0.0f &&
+		m[1][0] == 0.0f && m[1][1] == 1.0f && m[1][2] == 0.0f && m[1][3] == 0.0f &&
+		m[2][0] == 0.0f && m[2][1] == 0.0f && m[2][2] == 1.0f && m[2][3] == 0.0f &&
+		m[3][0] == 0.0f && m[3][1] == 0.0f && m[3][2] == 0.0f && m[3][3] == 1.0f;
+}
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline Vector VMatrix::ApplyRotation(const Vector &vVec) const
+{
+	return VMul3x3(vVec);
+}
+
+inline VMatrix VMatrix::operator~() const
+{
+	VMatrix mRet;
+	InverseGeneral(mRet);
+	return mRet;
+}
+
+#endif
+
+
+//-----------------------------------------------------------------------------
+// Accessors
+//-----------------------------------------------------------------------------
+inline void MatrixGetColumn( const VMatrix &src, int nCol, Vector *pColumn )
+{
+	Assert( (nCol >= 0) && (nCol <= 3) );
+
+	pColumn->x = src[0][nCol];
+	pColumn->y = src[1][nCol];
+	pColumn->z = src[2][nCol];
+}
+
+inline void MatrixSetColumn( VMatrix &src, int nCol, const Vector &column )
+{
+	Assert( (nCol >= 0) && (nCol <= 3) );
+
+	src.m[0][nCol] = column.x;
+	src.m[1][nCol] = column.y;
+	src.m[2][nCol] = column.z;
+}
+
+inline void MatrixGetRow( const VMatrix &src, int nRow, Vector *pRow )
+{
+	Assert( (nRow >= 0) && (nRow <= 3) );
+	*pRow = *(Vector*)src[nRow];
+}
+
+inline void MatrixSetRow( VMatrix &dst, int nRow, const Vector &row )
+{
+	Assert( (nRow >= 0) && (nRow <= 3) );
+	*(Vector*)dst[nRow] = row;
+}
+
+
+//-----------------------------------------------------------------------------
+// Vector3DMultiplyPosition treats src2 as if it's a point (adds the translation)
+//-----------------------------------------------------------------------------
+// NJS: src2 is passed in as a full vector rather than a reference to prevent the need
+// for 2 branches and a potential copy in the body.  (ie, handling the case when the src2
+// reference is the same as the dst reference ).
+inline void Vector3DMultiplyPosition( const VMatrix& src1, const VectorByValue src2, Vector& dst )
+{
+	dst[0] = src1[0][0] * src2.x + src1[0][1] * src2.y + src1[0][2] * src2.z + src1[0][3];
+	dst[1] = src1[1][0] * src2.x + src1[1][1] * src2.y + src1[1][2] * src2.z + src1[1][3];
+	dst[2] = src1[2][0] * src2.x + src1[2][1] * src2.y + src1[2][2] * src2.z + src1[2][3];
+}
+
+
+//-----------------------------------------------------------------------------
+// Transform a plane that has an axis-aligned normal
+//-----------------------------------------------------------------------------
+inline void MatrixTransformAxisAlignedPlane( const VMatrix &src, int nDim, float flSign, float flDist, cplane_t &outPlane )
+{
+	// See MatrixTransformPlane in the .cpp file for an explanation of the algorithm.
+	MatrixGetColumn( src, nDim, &outPlane.normal );
+	outPlane.normal *= flSign;
+	outPlane.dist = flDist * DotProduct( outPlane.normal, outPlane.normal );
+
+	// NOTE: Writing this out by hand because it doesn't inline (inline depth isn't large enough)
+	// This should read outPlane.dist += DotProduct( outPlane.normal, src.GetTranslation );
+	outPlane.dist += outPlane.normal.x * src.m[0][3] + outPlane.normal.y * src.m[1][3] + outPlane.normal.z * src.m[2][3];
+}
+
+
+//-----------------------------------------------------------------------------
+// Matrix equality test
+//-----------------------------------------------------------------------------
+inline bool MatricesAreEqual( const VMatrix &src1, const VMatrix &src2, float flTolerance )
+{
+	for ( int i = 0; i < 3; ++i )
+	{
+		for ( int j = 0; j < 3; ++j )
+		{
+			if ( fabs( src1[i][j] - src2[i][j] ) > flTolerance )
+				return false;
+		}
+	}
+	return true;
+}
+
+//-----------------------------------------------------------------------------
+//
+//-----------------------------------------------------------------------------
+void MatrixBuildOrtho( VMatrix& dst, double left, double top, double right, double bottom, double zNear, double zFar );
+void MatrixBuildPerspectiveX( VMatrix& dst, double flFovX, double flAspect, double flZNear, double flZFar );
+void MatrixBuildPerspectiveOffCenterX( VMatrix& dst, double flFovX, double flAspect, double flZNear, double flZFar, double bottom, double top, double left, double right );
+void MatrixBuildPerspectiveZRange( VMatrix& dst, double flZNear, double flZFar );
+
+inline void MatrixOrtho( VMatrix& dst, double left, double top, double right, double bottom, double zNear, double zFar )
+{
+	VMatrix mat;
+	MatrixBuildOrtho( mat, left, top, right, bottom, zNear, zFar );
+
+	VMatrix temp;
+	MatrixMultiply( dst, mat, temp );
+	dst = temp;
+}
+
+inline void MatrixPerspectiveX( VMatrix& dst, double flFovX, double flAspect, double flZNear, double flZFar )
+{
+	VMatrix mat;
+	MatrixBuildPerspectiveX( mat, flFovX, flAspect, flZNear, flZFar );
+
+	VMatrix temp;
+	MatrixMultiply( dst, mat, temp );
+	dst = temp;
+}
+
+inline void MatrixPerspectiveOffCenterX( VMatrix& dst, double flFovX, double flAspect, double flZNear, double flZFar, double bottom, double top, double left, double right )
+{
+	VMatrix mat;
+	MatrixBuildPerspectiveOffCenterX( mat, flFovX, flAspect, flZNear, flZFar, bottom, top, left, right );
+
+	VMatrix temp;
+	MatrixMultiply( dst, mat, temp );
+	dst = temp;
+}
+
+#endif
+
+

public/mathlib/vplane.h

@@ -0,0 +1,182 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+//
+// Purpose: 
+//
+// $Workfile:     $
+// $Date:         $
+// $NoKeywords: $
+//=============================================================================//
+
+#ifndef VPLANE_H
+#define VPLANE_H
+
+#ifdef _WIN32
+#pragma once
+#endif
+
+#include "mathlib/vector.h"
+
+typedef int SideType;
+
+// Used to represent sides of things like planes.
+#define	SIDE_FRONT	0
+#define	SIDE_BACK	1
+#define	SIDE_ON		2
+
+#define VP_EPSILON	0.01f
+
+
+class VPlane
+{
+public:
+				VPlane();
+				VPlane(const Vector &vNormal, vec_t dist);
+
+	void		Init(const Vector &vNormal, vec_t dist);
+
+	// Return the distance from the point to the plane.
+	vec_t		DistTo(const Vector &vVec) const;
+
+	// Copy.
+	VPlane&		operator=(const VPlane &thePlane);
+
+	// Returns SIDE_ON, SIDE_FRONT, or SIDE_BACK.
+	// The epsilon for SIDE_ON can be passed in.
+	SideType	GetPointSide(const Vector &vPoint, vec_t sideEpsilon=VP_EPSILON) const;
+
+	// Returns SIDE_FRONT or SIDE_BACK.
+	SideType	GetPointSideExact(const Vector &vPoint) const;
+
+	// Classify the box with respect to the plane.
+	// Returns SIDE_ON, SIDE_FRONT, or SIDE_BACK
+	SideType	BoxOnPlaneSide(const Vector &vMin, const Vector &vMax) const;
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+	// Flip the plane.
+	VPlane		Flip();
+
+	// Get a point on the plane (normal*dist).
+	Vector		GetPointOnPlane() const;
+
+	// Snap the specified point to the plane (along the plane's normal).
+	Vector		SnapPointToPlane(const Vector &vPoint) const;
+#endif
+
+public:
+	Vector		m_Normal;
+	vec_t		m_Dist;
+
+#ifdef VECTOR_NO_SLOW_OPERATIONS
+private:
+	// No copy constructors allowed if we're in optimal mode
+	VPlane(const VPlane& vOther);
+#endif
+};
+
+
+//-----------------------------------------------------------------------------
+// Inlines.
+//-----------------------------------------------------------------------------
+inline VPlane::VPlane()
+{
+}
+
+inline VPlane::VPlane(const Vector &vNormal, vec_t dist)
+{
+	m_Normal = vNormal;
+	m_Dist = dist;
+}
+
+inline void	VPlane::Init(const Vector &vNormal, vec_t dist)
+{
+	m_Normal = vNormal;
+	m_Dist = dist;
+}
+
+inline vec_t VPlane::DistTo(const Vector &vVec) const
+{
+	return vVec.Dot(m_Normal) - m_Dist;
+}
+
+inline VPlane& VPlane::operator=(const VPlane &thePlane)
+{
+	m_Normal = thePlane.m_Normal;
+	m_Dist = thePlane.m_Dist;
+	return *this;
+}
+
+#ifndef VECTOR_NO_SLOW_OPERATIONS
+
+inline VPlane VPlane::Flip()
+{
+	return VPlane(-m_Normal, -m_Dist);
+}
+
+inline Vector VPlane::GetPointOnPlane() const
+{
+	return m_Normal * m_Dist;
+}
+
+inline Vector VPlane::SnapPointToPlane(const Vector &vPoint) const
+{
+	return vPoint - m_Normal * DistTo(vPoint);
+}
+
+#endif
+
+inline SideType VPlane::GetPointSide(const Vector &vPoint, vec_t sideEpsilon) const
+{
+	vec_t fDist;
+
+	fDist = DistTo(vPoint);
+	if(fDist >= sideEpsilon)
+		return SIDE_FRONT;
+	else if(fDist <= -sideEpsilon)
+		return SIDE_BACK;
+	else
+		return SIDE_ON;
+}
+
+inline SideType VPlane::GetPointSideExact(const Vector &vPoint) const
+{
+	return DistTo(vPoint) > 0.0f ? SIDE_FRONT : SIDE_BACK;
+}
+
+
+// BUGBUG: This should either simply use the implementation in mathlib or cease to exist.
+// mathlib implementation is much more efficient.  Check to see that VPlane isn't used in
+// performance critical code.
+inline SideType VPlane::BoxOnPlaneSide(const Vector &vMin, const Vector &vMax) const
+{
+	int i, firstSide, side;
+	TableVector vPoints[8] = 
+	{
+		{ vMin.x, vMin.y, vMin.z },
+		{ vMin.x, vMin.y, vMax.z },
+		{ vMin.x, vMax.y, vMax.z },
+		{ vMin.x, vMax.y, vMin.z },
+
+		{ vMax.x, vMin.y, vMin.z },
+		{ vMax.x, vMin.y, vMax.z },
+		{ vMax.x, vMax.y, vMax.z },
+		{ vMax.x, vMax.y, vMin.z },
+	};
+
+	firstSide = GetPointSideExact(vPoints[0]);
+	for(i=1; i < 8; i++)
+	{
+		side = GetPointSideExact(vPoints[i]);
+
+		// Does the box cross the plane?
+		if(side != firstSide)
+			return SIDE_ON;
+	}
+
+	// Ok, they're all on the same side, return that.
+	return firstSide;
+}
+
+
+
+
+#endif // VPLANE_H