# HG changeset patch
# User Eris Caffee <discordia@eldalin.com>
# Date 1315326399 18000
# Node ID 11e216148d1cf293eff7f7ddb4f1e1526255d9e2
# Parent  f6e6f3c8f7eb463287b355aa78d6e0222221b44a
Making progress in Matrix.  Reorganized code a bit.  Reworked Vector routines to use += -= *=

diff -r f6e6f3c8f7eb -r 11e216148d1c .hgignore
--- a/.hgignore	Sun Sep 04 12:46:08 2011 -0500
+++ b/.hgignore	Tue Sep 06 11:26:39 2011 -0500
@@ -2,4 +2,4 @@
 
 build*/*
 *~
-
+tmp/*
diff -r f6e6f3c8f7eb -r 11e216148d1c CMakeLists.txt
--- a/CMakeLists.txt	Sun Sep 04 12:46:08 2011 -0500
+++ b/CMakeLists.txt	Tue Sep 06 11:26:39 2011 -0500
@@ -138,8 +138,11 @@
 # test driver.  Also, it just seems good practice to explicitly list things since
 # it avoids accidental inclusion of things that ought not to be included.
 set (SRCS_${App_Name}
-  src/Math.cpp
+  src/Vector.cpp
+  src/Matrix.cpp
   include/Math.h
+  include/Vector.h
+  include/Matrix.h
 )
 
 include_directories (
diff -r f6e6f3c8f7eb -r 11e216148d1c include/Math.h
--- a/include/Math.h	Sun Sep 04 12:46:08 2011 -0500
+++ b/include/Math.h	Tue Sep 06 11:26:39 2011 -0500
@@ -2,6 +2,7 @@
 #define MATH_H_
 
 #include "Vector.h"
+#include "Matrix.h"
 
 namespace arda 
    {
diff -r f6e6f3c8f7eb -r 11e216148d1c include/Matrix.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/include/Matrix.h	Tue Sep 06 11:26:39 2011 -0500
@@ -0,0 +1,412 @@
+#ifndef MATRIX_H_
+#define MATRIX_H_
+
+#include "Vector.h"
+using namespace arda::Vector;
+
+#include <cmath>
+#include <string>
+
+namespace arda 
+   {
+   namespace Matrix
+      {
+      //////////////////////////////////////////////////////////////////////////
+      // Maxtrices
+      //
+      // Supported operations on matrices
+      //
+      // getstring*()	Returns a printable string representation of the matrix.
+      // assign*()	Assign the value of one matrix to another.
+      // add*()		Add two matrices
+      // subtract*()	A convenience function to avoid writing 
+      //			add(m1, scale(m2, -1))
+      // transpose*()	Returns the transpose of a matrix.
+      // scale*()	Multiplies a matrix by a scalar.
+      // multvec*()	Multiplies a matrix by a vector.
+      // multmat*()	Multiplies a matrix by another matrix.
+      // det*()		Calculates the determinant of a matrix.
+      // inverse*()	Returns the inverse of a non-singular matrix.
+      // setidentity*()	Sets a matrix to the identity matrix.
+      // getcol*()	Returns the vector that is column n of the matrix.
+      // setcol*()	Set the vector that is column n of the matrix.
+      //
+      // assign*(), add*(), and subtract*() are defined inline, though this might
+      // might not be the best thing.  Some real world testing is needed to decide
+      // whether or not htye should be remain inline.
+      
+      typedef int Matrix2i[4];
+      typedef float Matrix2f[4];
+      typedef double Matrix2d[4];
+      
+      typedef int Matrix3i[9];
+      typedef float Matrix3f[9];
+      typedef double Matrix3d[9];
+      
+      typedef int Matrix4i[16];
+      typedef float Matrix4f[16];
+      typedef double Matrix4d[16];
+      
+      //////////////////////////////////////////////////////////////////////////
+      // getstring*()
+      std::string & getstring2m(Matrix2i m, std::string & s);
+      std::string & getstring2m(Matrix2f m, std::string & s);
+      std::string & getstring2m(Matrix2d m, std::string & s);
+      
+      std::string & getstring3m(Matrix3i m, std::string & s);
+      std::string & getstring3m(Matrix3f m, std::string & s);
+      std::string & getstring3m(Matrix3d m, std::string & s);
+
+      std::string & getstring4m(Matrix4i m, std::string & s);
+      std::string & getstring4m(Matrix4f m, std::string & s);
+      std::string & getstring4m(Matrix4d m, std::string & s);
+
+      //////////////////////////////////////////////////////////////////////////
+      // assign*()
+
+      inline int * assign2m(Matrix2i m1, Matrix2i m2)
+ 	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; 
+	 m1[2] = m2[2]; m1[3] = m2[3]; 
+	 return m1; 
+	 }
+      inline float * assign2m(Matrix2f m1, Matrix2f m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; 
+	 m1[2] = m2[2]; m1[3] = m2[3]; 
+	 return m1; 
+	 }
+      inline double * assign2m(Matrix2d m1, Matrix2d m2)
+	 {
+	 m1[0] = m2[0]; m1[1] = m2[1]; 
+	 m1[2] = m2[2]; m1[3] = m2[3]; 
+	 return m1; 
+	 }
+
+      inline int * assign3m(Matrix3i m1, Matrix3i m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2];
+	 m1[3] = m2[3]; m1[4] = m2[4]; m1[5] = m2[5]; 	
+	 m1[6] = m2[6]; m1[7] = m2[7]; m1[8] = m2[8]; 
+	 return m1; 
+	 }
+      inline float * assign3m(Matrix3f m1, Matrix3f m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2];
+	 m1[3] = m2[3]; m1[4] = m2[4]; m1[5] = m2[5]; 	
+	 m1[6] = m2[6]; m1[7] = m2[7]; m1[8] = m2[8]; 
+	 return m1; 
+	 }
+      inline double * assign3m(Matrix3d m1, Matrix3d m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2];
+	 m1[3] = m2[3]; m1[4] = m2[4]; m1[5] = m2[5]; 	
+	 m1[6] = m2[6]; m1[7] = m2[7]; m1[8] = m2[8]; 
+	 return m1; 
+	 }
+
+      inline int * assign4m(Matrix4i m1, Matrix4i m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2]; m1[3] = m2[3];
+	 m1[4] = m2[4]; m1[5] = m2[5]; m1[6] = m2[6]; m1[7] = m2[7];
+	 m1[8] = m2[8]; m1[9] = m2[9]; m1[10] = m2[10]; m1[11] = m2[11];
+	 m1[12] = m2[12]; m1[13] = m2[13]; m1[14] = m2[14]; m1[15] = m2[15];
+	 return m1; 
+	 }
+      inline float * assign4m(Matrix4f m1, Matrix4f m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2]; m1[3] = m2[3];
+	 m1[4] = m2[4]; m1[5] = m2[5]; m1[6] = m2[6]; m1[7] = m2[7];
+	 m1[8] = m2[8]; m1[9] = m2[9]; m1[10] = m2[10]; m1[11] = m2[11];
+	 m1[12] = m2[12]; m1[13] = m2[13]; m1[14] = m2[14]; m1[15] = m2[15];
+	 return m1; 
+	 }
+      inline double * assign4m(Matrix4d m1, Matrix4d m2)
+	 { 
+	 m1[0] = m2[0]; m1[1] = m2[1]; m1[2] = m2[2]; m1[3] = m2[3];
+	 m1[4] = m2[4]; m1[5] = m2[5]; m1[6] = m2[6]; m1[7] = m2[7];
+	 m1[8] = m2[8]; m1[9] = m2[9]; m1[10] = m2[10]; m1[11] = m2[11];
+	 m1[12] = m2[12]; m1[13] = m2[13]; m1[14] = m2[14]; m1[15] = m2[15];
+	 return m1; 
+	 }
+
+      //////////////////////////////////////////////////////////////////////////
+      // add*()
+
+      inline int * add2m(Matrix2i m1, Matrix2i m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1];
+	 m1[2] += m2[2]; m1[3] += m2[3];
+	 return m1;
+	 }
+      inline float * add2m(Matrix2f m1, Matrix2f m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1];
+	 m1[2] += m2[2]; m1[3] += m2[3];
+	 return m1;
+	 }
+      inline double * add2m(Matrix2d m1, Matrix2d m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1];
+	 m1[2] += m2[2]; m1[3] += m2[3];
+	 return m1;
+	 }
+
+      inline int * add3m(Matrix3i m1, Matrix3i m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; 
+	 m1[3] += m2[3]; m1[4] += m2[4]; m1[5] += m2[5];
+	 m1[6] += m2[6]; m1[7] += m2[7]; m1[8] += m2[8];
+	 return m1;
+	 }
+      inline float * add3m(Matrix3f m1, Matrix3f m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; 
+	 m1[3] += m2[3]; m1[4] += m2[4]; m1[5] += m2[5];
+	 m1[6] += m2[6]; m1[7] += m2[7]; m1[8] += m2[8];
+	 return m1;
+	 }
+      inline double * add3m(Matrix3d m1, Matrix3d m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; 
+	 m1[3] += m2[3]; m1[4] += m2[4]; m1[5] += m2[5];
+	 m1[6] += m2[6]; m1[7] += m2[7]; m1[8] += m2[8];
+	 return m1;
+	 }
+
+      inline int * add4m(Matrix4i m1, Matrix4i m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; m1[3] += m2[3]; 
+	 m1[4] += m2[4]; m1[5] += m2[5]; m1[6] += m2[6]; m1[7] += m2[7]; 
+	 m1[8] += m2[8]; m1[9] += m2[9]; m1[10] += m2[10]; m1[11] += m2[11];
+	 m1[12] += m2[12]; m1[13] += m2[13]; m1[14] += m2[14]; m1[15] += m2[15];
+	 return m1;
+	 }
+      inline float * add4m(Matrix4f m1, Matrix4f m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; m1[3] += m2[3]; 
+	 m1[4] += m2[4]; m1[5] += m2[5]; m1[6] += m2[6]; m1[7] += m2[7]; 
+	 m1[8] += m2[8]; m1[9] += m2[9]; m1[10] += m2[10]; m1[11] += m2[11];
+	 m1[12] += m2[12]; m1[13] += m2[13]; m1[14] += m2[14]; m1[15] += m2[15];
+	 return m1;
+	 }
+      inline double * add4m(Matrix4d m1, Matrix4d m2)
+	 {
+	 m1[0] += m2[0]; m1[1] += m2[1]; m1[2] += m2[2]; m1[3] += m2[3]; 
+	 m1[4] += m2[4]; m1[5] += m2[5]; m1[6] += m2[6]; m1[7] += m2[7]; 
+	 m1[8] += m2[8]; m1[9] += m2[9]; m1[10] += m2[10]; m1[11] += m2[11];
+	 m1[12] += m2[12]; m1[13] += m2[13]; m1[14] += m2[14]; m1[15] += m2[15];
+	 return m1;
+	 }
+
+      //////////////////////////////////////////////////////////////////////////
+      // subtract*()
+
+      inline int * subtract2m(Matrix2i m1, Matrix2i m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1];
+	 m1[2] -= m2[2]; m1[3] -= m2[3];
+	 return m1;
+	 }
+      inline float * subtract2m(Matrix2f m1, Matrix2f m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1];
+	 m1[2] -= m2[2]; m1[3] -= m2[3];
+	 return m1;
+	 }
+      inline double * subtract2m(Matrix2d m1, Matrix2d m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1];
+	 m1[2] -= m2[2]; m1[3] -= m2[3];
+	 return m1;
+	 }
+
+      inline int * subtract3m(Matrix3i m1, Matrix3i m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; 
+	 m1[3] -= m2[3]; m1[4] -= m2[4]; m1[5] -= m2[5];
+	 m1[6] -= m2[6]; m1[7] -= m2[7]; m1[8] -= m2[8];
+	 return m1;
+	 }
+      inline float * subtract3m(Matrix3f m1, Matrix3f m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; 
+	 m1[3] -= m2[3]; m1[4] -= m2[4]; m1[5] -= m2[5];
+	 m1[6] -= m2[6]; m1[7] -= m2[7]; m1[8] -= m2[8];
+	 return m1;
+	 }
+      inline double * subtract3m(Matrix3d m1, Matrix3d m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; 
+	 m1[3] -= m2[3]; m1[4] -= m2[4]; m1[5] -= m2[5];
+	 m1[6] -= m2[6]; m1[7] -= m2[7]; m1[8] -= m2[8];
+	 return m1;
+	 }
+
+      inline int * subtract4m(Matrix4i m1, Matrix4i m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; m1[3] -= m2[3]; 
+	 m1[4] -= m2[4]; m1[5] -= m2[5]; m1[6] -= m2[6]; m1[7] -= m2[7]; 
+	 m1[8] -= m2[8]; m1[9] -= m2[9]; m1[10] -= m2[10]; m1[11] -= m2[11];
+	 m1[12] -= m2[12]; m1[13] -= m2[13]; m1[14] -= m2[14]; m1[15] -= m2[15];
+	 return m1;
+	 }
+      inline float * subtract4m(Matrix4f m1, Matrix4f m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; m1[3] -= m2[3]; 
+	 m1[4] -= m2[4]; m1[5] -= m2[5]; m1[6] -= m2[6]; m1[7] -= m2[7]; 
+	 m1[8] -= m2[8]; m1[9] -= m2[9]; m1[10] -= m2[10]; m1[11] -= m2[11];
+	 m1[12] -= m2[12]; m1[13] -= m2[13]; m1[14] -= m2[14]; m1[15] -= m2[15];
+	 return m1;
+	 }
+      inline double * subtract4m(Matrix4d m1, Matrix4d m2)
+	 {
+	 m1[0] -= m2[0]; m1[1] -= m2[1]; m1[2] -= m2[2]; m1[3] -= m2[3]; 
+	 m1[4] -= m2[4]; m1[5] -= m2[5]; m1[6] -= m2[6]; m1[7] -= m2[7]; 
+	 m1[8] -= m2[8]; m1[9] -= m2[9]; m1[10] -= m2[10]; m1[11] -= m2[11];
+	 m1[12] -= m2[12]; m1[13] -= m2[13]; m1[14] -= m2[14]; m1[15] -= m2[15];
+	 return m1;
+	 }
+
+      //////////////////////////////////////////////////////////////////////////
+      // scale*()
+
+      inline int * scale2m(Matrix2i m, int s)
+	 { 
+	 m[0] *= s; m[1] *= s; 
+	 m[2] *= s; m[3] *= s; 
+	 return m;
+	 }
+      inline float * scale2m(Matrix2f m, float s)
+	 { 
+	 m[0] *= s; m[1] *= s; 
+	 m[2] *= s; m[3] *= s; 
+	 return m;
+	 }
+      inline double * scale2m(Matrix2d m, double s)
+	 { 
+	 m[0] *= s; m[1] *= s; 
+	 m[2] *= s; m[3] *= s; 
+	 return m;
+	 }
+
+      inline int * scale3m(Matrix3i m, int s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s;
+	 m[3] *= s; m[4] *= s; m[5] *= s;
+	 m[4] *= s; m[5] *= s; m[6] *= s;
+	 return m;
+	 }
+      inline float * scale3m(Matrix3f m, float s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s;
+	 m[3] *= s; m[4] *= s; m[5] *= s;
+	 m[4] *= s; m[5] *= s; m[6] *= s;
+	 return m;
+	 }
+      inline double * scale3m(Matrix3d m, double s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s;
+	 m[3] *= s; m[4] *= s; m[5] *= s;
+	 m[4] *= s; m[5] *= s; m[6] *= s;
+	 return m;
+	 }
+
+      inline int * scale4m(Matrix4i m, int s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s; m[3] *= s; 
+	 m[4] *= s; m[5] *= s; m[4] *= s; m[5] *= s;
+	 m[6] *= s; m[7] *= s; m[8] *= s; m[9] *= s;
+	 m[10] *= s; m[11] *= s; m[12] *= s; m[13] *= s;
+	 return m;
+	 }
+      inline float * scale4m(Matrix4f m, float s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s; m[3] *= s; 
+	 m[4] *= s; m[5] *= s; m[4] *= s; m[5] *= s;
+	 m[6] *= s; m[7] *= s; m[8] *= s; m[9] *= s;
+	 m[10] *= s; m[11] *= s; m[12] *= s; m[13] *= s;
+	 return m;
+	 }
+      inline double * scale4m(Matrix4d m, double s)
+	 { 
+	 m[0] *= s; m[1] *= s; m[2] *= s; m[3] *= s; 
+	 m[4] *= s; m[5] *= s; m[6] *= s; m[7] *= s;
+	 m[8] *= s; m[9] *= s; m[10] *= s; m[11] *= s;
+	 m[12] *= s; m[13] *= s; m[14] *= s; m[15] *= s;
+	 return m;
+	 }
+
+      //////////////////////////////////////////////////////////////////////////
+      // multvec*mv(Matrix m, Vector v, Vector vres)
+      //	Right multiplies a matrix m by a (column) vector v and stores 
+      //	the result in vres.
+      //	Returns vres.
+      //
+      // multvec*vm(Vector v, Matrix m, Vector vres)
+      //	Left multiplies a (row) vector v by a matrix m and stores the
+      //	result in vres.
+      //	Returns vres.
+
+      int * multvec2mv(Matrix2i m, Vector2i v, Vector2i vres);
+      float * multvec2mv(Matrix2f m, Vector2f v, Vector2f vres);
+      double * multvec2mv(Matrix2d m, Vector2d v, Vector2d vres);
+
+      int * multvec3mv(Matrix3i m, Vector3i v, Vector3i vres);
+      float * multvec3mv(Matrix3f m, Vector3f v, Vector3f vres);
+      double * multvec3mv(Matrix3d m, Vector3d v, Vector3d vres);
+
+      int * multvec4mv(Matrix4i m, Vector4i v, Vector4i vres);
+      float * multvec4mv(Matrix4f m, Vector4f v, Vector4f vres);
+      double * multvec4mv(Matrix4d m, Vector4d v, Vector4d vres);
+
+      int * multvec2vm(Vector2i v, Matrix2i m, Vector2i vres);
+      float * multvec2vm(Vector2f v, Matrix2f m, Vector2f vres);
+      double * multvec2vm(Vector2d v, Matrix2d m, Vector2d vres);
+
+      int * multvec3vm(Vector3i v, Matrix3i m, Vector3i vres);
+      float * multvec3vm(Vector3f v, Matrix3f m, Vector3f vres);
+      double * multvec3vm(Vector3d v, Matrix3d m, Vector3d vres);
+
+      int * multvec4vm(Vector4i v, Matrix4i m, Vector4i vres);
+      float * multvec4vm(Vector4f v, Matrix4f m, Vector4f vres);
+      double * multvec4vm(Vector4d v, Matrix4d m, Vector4d vres);
+
+      //////////////////////////////////////////////////////////////////////////
+      // multmat*(Matrix m1, Matrix m2, Matrix mres)
+      //	Multiply two matrices together and return the result in mres.
+      //	returns mres.
+
+      int * multmat2(Matrix2i m1, Matrix2i m2, Matrix2i mres);
+      float * multmat2(Matrix2f m1, Matrix2f m2, Matrix2f mres);
+      double * multmat2(Matrix2d m1, Matrix2d m2, Matrix2d mres);
+
+      int * multmat3(Matrix3i m1, Matrix3i m2, Matrix3i mres);
+      float * multmat3(Matrix3f m1, Matrix3f m2, Matrix3f mres);
+      double * multmat3(Matrix3d m1, Matrix3d m2, Matrix3d mres);
+
+      int * multmat4(Matrix4i m1, Matrix4i m2, Matrix4i mres);
+      float * multmat4(Matrix4f m1, Matrix4f m2, Matrix4f mres);
+      double * multmat4(Matrix4d m1, Matrix4d m2, Matrix4d mres);
+
+      //////////////////////////////////////////////////////////////////////////
+      // transpose*(Matrix m, Matrix mres)
+      //	Compute the transpose of a matrix and store the result in mres.
+      //	Returns mres.
+
+      int * transpose2(Matrix2i m, Matrix2i mres);
+      float * transpose2(Matrix2f m, Matrix2f mres);
+      double * transpose2(Matrix2d m, Matrix2d mres);
+
+      int * transpose3(Matrix3i m, Matrix3i mres);
+      float * transpose3(Matrix3f m, Matrix3f mres);
+      double * transpose3(Matrix3d m, Matrix3d mres);
+
+      int * transpose4(Matrix4i m, Matrix4i mres);
+      float * transpose4(Matrix4f m, Matrix4f mres);
+      double * transpose4(Matrix4d m, Matrix4d mres);
+
+      } // namespace Matrix
+
+   } // namespace arda
+#endif // MATRIX_H_
diff -r f6e6f3c8f7eb -r 11e216148d1c include/Vector.h
--- a/include/Vector.h	Sun Sep 04 12:46:08 2011 -0500
+++ b/include/Vector.h	Tue Sep 06 11:26:39 2011 -0500
@@ -6,29 +6,28 @@
 
 namespace arda 
    {
-
-   ////////////////////////////////////////////////////////////////////////////////
-   // Vectors
-   //
-   // Supported operations on vectors
-   //
-   // assign*()		Assign the value of one vector to another.
-   // add*()		Add two vectors
-   // cross()		Cross product of two 3-space vectors.
-   // dot*()		Dot product.
-   // get_angle*()	Calculate angle between two vectors.
-   // length*()		Length of a vector (the Euclidean norm).
-   // normalize*()	Normalize the vector.
-   // proj*()		Calculate the projection of one vector onto another.
-   // scale*()		Multiply vector by a scalar.
-   // subtract*()	A convenience function to avoid writing add(v1, scale(v2, -1))
-   // vectostr*()	Returns a printable string representation of the vector.
-   //
-   // All functions except vectostr are defined inline.
-
    namespace Vector 
       {
 
+      ////////////////////////////////////////////////////////////////////////////////
+      // Vectors
+      //
+      // Supported operations on vectors
+      //
+      // assign*()	Assign the value of one vector to another.
+      // add*()		Add two vectors
+      // cross()	Cross product of two 3-space vectors.
+      // dot*()		Dot product.
+      // get_angle*()	Calculate angle between two vectors.
+      // length*()	Length of a vector (the Euclidean norm).
+      // normalize*()	Normalize the vector.
+      // proj*()	Calculate the projection of one vector onto another.
+      // scale*()	Multiply vector by a scalar.
+      // subtract*()	A convenience function to avoid writing add(v1, scale(v2, -1))
+      // getstring*()	Returns a printable string representation of the vector.
+      //
+      // All functions except getstring are defined inline.
+
       typedef int Vector2i[2];
       typedef float Vector2f[2];
       typedef double Vector2d[2];
@@ -42,89 +41,89 @@
       typedef double Vector4d[4];
 
       //////////////////////////////////////////////////////////////////////////
-      // vectostr*(Vector v)
+      // getstring*(Vector v)
       //	Returns a C string representation of the vector.
       //	This is one of the few non-inline Vector functions.
 
-      std::string & vectostr2(Vector2i v, std::string & s);
-      std::string & vectostr2(Vector2f v, std::string & s);
-      std::string & vectostr2(Vector2d v, std::string & s);
-      std::string & vectostr3(Vector3i v, std::string & s);
-      std::string & vectostr3(Vector3f v, std::string & s);
-      std::string & vectostr3(Vector3d v, std::string & s);
-      std::string & vectostr4(Vector4i v, std::string & s);
-      std::string & vectostr4(Vector4f v, std::string & s);
-      std::string & vectostr4(Vector4d v, std::string & s);
+      std::string & getstring2v(Vector2i v, std::string & s);
+      std::string & getstring2v(Vector2f v, std::string & s);
+      std::string & getstring2v(Vector2d v, std::string & s);
+      std::string & getstring3v(Vector3i v, std::string & s);
+      std::string & getstring3v(Vector3f v, std::string & s);
+      std::string & getstring3v(Vector3d v, std::string & s);
+      std::string & getstring4v(Vector4i v, std::string & s);
+      std::string & getstring4v(Vector4f v, std::string & s);
+      std::string & getstring4v(Vector4d v, std::string & s);
       
       //////////////////////////////////////////////////////////////////////////
       // scale*(Vector v, scalar s)
       //	Multiply a vector by a scalar.
       //	Returns v so that the result may be passed to other functions.
 
-      inline int * scale2(Vector2i v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline int * scale2(Vector2i v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline int * scale2(Vector2i v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
+      inline int * scale2v(Vector2i v, const int s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline int * scale2v(Vector2i v, const float s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline int * scale2v(Vector2i v, const double s)
+	 { v[0] *= s ; v[1] *= s; return v; }
 
-      inline float * scale2(Vector2f v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline float * scale2(Vector2f v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline float * scale2(Vector2f v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
+      inline float * scale2v(Vector2f v, const int s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline float * scale2v(Vector2f v, const float s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline float * scale2v(Vector2f v, const double s)
+	 { v[0] *= s ; v[1] *= s; return v; }
 
-      inline double * scale2(Vector2d v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline double * scale2(Vector2d v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
-      inline double * scale2(Vector2d v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; return v; }
+      inline double * scale2v(Vector2d v, const int s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline double * scale2v(Vector2d v, const float s)
+	 { v[0] *= s ; v[1] *= s; return v; }
+      inline double * scale2v(Vector2d v, const double s)
+	 { v[0] *= s ; v[1] *= s; return v; }
 
 
-      inline int * scale3(Vector3i v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline int * scale3(Vector3i v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline int * scale3(Vector3i v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
+      inline int * scale3v(Vector3i v, const int s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline int * scale3v(Vector3i v, const float s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline int * scale3v(Vector3i v, const double s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
  
-      inline float * scale3(Vector3f v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline float * scale3(Vector3f v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline float * scale3(Vector3f v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
+      inline float * scale3v(Vector3f v, const int s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline float * scale3v(Vector3f v, const float s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline float * scale3v(Vector3f v, const double s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
 
-      inline double * scale3(Vector3d v, const int s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline double *  scale3(Vector3d v, const float s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
-      inline double * scale3(Vector3d v, const double s)
-	 { v[0] = v[0] * s ; v[1] = v[1] * s; v[2] = v[2] * s; return v; }
+      inline double * scale3v(Vector3d v, const int s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline double *  scale3v(Vector3d v, const float s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
+      inline double * scale3v(Vector3d v, const double s)
+	 { v[0] *= s ; v[1] *= s; v[2] *= s; return v; }
 
 
-      inline int * scale4(Vector4i v, const int s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline int * scale4(Vector4i v, const float s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline int * scale4(Vector4i v, const double s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
+      inline int * scale4v(Vector4i v, const int s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline int * scale4v(Vector4i v, const float s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline int * scale4v(Vector4i v, const double s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
 
-      inline float * scale4(Vector4f v, const int s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline float * scale4(Vector4f v, const float s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline float * scale4(Vector4f v, const double s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
+      inline float * scale4v(Vector4f v, const int s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline float * scale4v(Vector4f v, const float s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline float * scale4v(Vector4f v, const double s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
 
-      inline double * scale4(Vector4d v, const int s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline double * scale4(Vector4d v, const float s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
-      inline double * scale4(Vector4d v, const double s) 
-	 { v[0] = v[0] * s ; v[1] = v[1] * s;  v[2] = v[2] * s;  v[3] = v[3] * s; return v; }
+      inline double * scale4v(Vector4d v, const int s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline double * scale4v(Vector4d v, const float s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
+      inline double * scale4v(Vector4d v, const double s) 
+	 { v[0] *= s ; v[1] *= s;  v[2] *= s;  v[3] *= s; return v; }
 
 
       //////////////////////////////////////////////////////////////////////////
@@ -132,25 +131,25 @@
       //	The dot product of two vectors, which must be the same type
       //	Returns a long for int vectors, a double for floating point vectors.
 
-      inline long dot2(Vector2i v1, Vector2i v2)
+      inline long dot2v(Vector2i v1, Vector2i v2)
 	 { return (long) v1[0]*v2[0] + (long) v1[1]*v2[1]; }
-      inline double dot2(Vector2f v1, Vector2f v2)
+      inline double dot2v(Vector2f v1, Vector2f v2)
 	 { return (double) v1[0]*v2[0] + (double) v1[1]*v2[1]; }
-      inline double dot2(Vector2d v1, Vector2d v2)
+      inline double dot2v(Vector2d v1, Vector2d v2)
 	 { return v1[0]*v2[0] + v1[1]*v2[1]; }
 
-      inline long dot3(Vector3i v1, Vector3i v2)
+      inline long dot3v(Vector3i v1, Vector3i v2)
 	 { return (long) v1[0]*v2[0] + (long) v1[1]*v2[1] + (long) v1[2]*v2[2]; }
-      inline double dot3(Vector3f v1, Vector3f v2)
+      inline double dot3v(Vector3f v1, Vector3f v2)
 	 { return (double) v1[0]*v2[0] + (double) v1[1]*v2[1] + (double) v1[2]*v2[2]; }
-      inline double dot3(Vector3d v1, Vector3d v2)
+      inline double dot3v(Vector3d v1, Vector3d v2)
 	 { return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; }
 
-      inline long dot4(Vector4i v1, Vector4i v2)
+      inline long dot4v(Vector4i v1, Vector4i v2)
 	 { return (long) v1[0]*v2[0] + (long) v1[1]*v2[1] + (long) v1[2]*v2[2] + (long) v1[3]*v2[3]; }
-      inline double dot4(Vector4f v1, Vector4f v2)
+      inline double dot4v(Vector4f v1, Vector4f v2)
 	 { return (double) v1[0]*v2[0] + (double) v1[1]*v2[1] + (double) v1[2]*v2[2] + (double) v1[3]*v2[3]; }
-      inline double dot4(Vector4d v1, Vector4d v2)
+      inline double dot4v(Vector4d v1, Vector4d v2)
 	 { return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3]; }
 
 
@@ -160,74 +159,74 @@
       //	Returns a float for float for float vectors, and a double for 
       //	int and double vectors.
 
-      inline double length2(Vector2i v)
-	 { return sqrt((dot2(v, v))); }
-      inline double length2(Vector2f v)
-	 { return sqrtf((dot2(v, v))); }
-      inline double length2(Vector2d v)
-	 { return sqrt((dot2(v, v))); }
+      inline double length2v(Vector2i v)
+	 { return sqrt((dot2v(v, v))); }
+      inline double length2v(Vector2f v)
+	 { return sqrtf((dot2v(v, v))); }
+      inline double length2v(Vector2d v)
+	 { return sqrt((dot2v(v, v))); }
 
-      inline double length3(Vector3i v)
-	 { return sqrt((dot3(v, v))); }
-      inline double length3(Vector3f v)
-	 { return sqrtf((dot3(v, v))); }
-      inline double length3(Vector3d v)
-	 { return sqrt((dot3(v, v))); }
+      inline double length3v(Vector3i v)
+	 { return sqrt((dot3v(v, v))); }
+      inline double length3v(Vector3f v)
+	 { return sqrtf((dot3v(v, v))); }
+      inline double length3v(Vector3d v)
+	 { return sqrt((dot3v(v, v))); }
 
-      inline double length4(Vector4i v)
-	 { return sqrt((dot4(v, v))); }
-      inline double length4(Vector4f v)
-	 { return sqrtf((dot4(v, v))); }
-      inline double length4(Vector4d v)
-	 { return sqrt((dot4(v, v))); }
+      inline double length4v(Vector4i v)
+	 { return sqrt((dot4v(v, v))); }
+      inline double length4v(Vector4f v)
+	 { return sqrtf((dot4v(v, v))); }
+      inline double length4v(Vector4d v)
+	 { return sqrt((dot4v(v, v))); }
 
 
       //////////////////////////////////////////////////////////////////////////
-      // add*(Vector v1, Vector v2, Vector vres)
-      //	Add two vectors together.
-      //        Returns vres so that the result may be used in further 
+      // add*(Vector v1, Vector v2)
+      //	Add v2 to v1.
+      //        Returns v1 so that the result may be used in further 
       //	calculations.
 
-      inline int * add2(Vector2i v1, Vector2i v2, Vector2i vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; return vres; }
-      inline float * add2(Vector2f v1, Vector2f v2, Vector2f vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; return vres; }
-      inline double * add2(Vector2d v1, Vector2d v2, Vector2d vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; return vres; }
+      inline int * add2v(Vector2i v1, Vector2i v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; return v1; }
+      inline float * add2v(Vector2f v1, Vector2f v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; return v1; }
+      inline double * add2v(Vector2d v1, Vector2d v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; return v1; }
 
-      inline int * add3(Vector3i v1, Vector3i v2, Vector3i vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2]; return vres; }
-      inline float * add3(Vector3f v1, Vector3f v2, Vector3f vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2];  return vres; }
-      inline double * add3(Vector3d v1, Vector3d v2, Vector3d vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2];  return vres; }
+      inline int * add3v(Vector3i v1, Vector3i v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2]; return v1; }
+      inline float * add3v(Vector3f v1, Vector3f v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2];  return v1; }
+      inline double * add3v(Vector3d v1, Vector3d v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2];  return v1; }
 
-      inline int * add4(Vector4i v1, Vector4i v2, Vector4i vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2]; vres[3] = v1[3] + v2[3]; return vres;}
-      inline float * add4(Vector4f v1, Vector4f v2, Vector4f vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2]; vres[3] = v1[3] + v2[3]; return vres; }
-      inline double * add4(Vector4d v1, Vector4d v2, Vector4d vres)
-	 { vres[0] = v1[0] + v2[0]; vres[1] = v1[1] + v2[1]; vres[2] = v1[2] + v2[2]; vres[3] = v1[3] + v2[3]; return vres; }
+      inline int * add4v(Vector4i v1, Vector4i v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2]; v1[3] += v2[3]; return v1; }
+      inline float * add4v(Vector4f v1, Vector4f v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2]; v1[3] += v2[3]; return v1; }
+      inline double * add4v(Vector4d v1, Vector4d v2)
+	 { v1[0] += v2[0]; v1[1] += v2[1]; v1[2] += v2[2]; v1[3] += v2[3]; return v1; }
 
       //////////////////////////////////////////////////////////////////////////
       // normalize*(Vector v)
       //	Normalizes the vecor.  Only defined for floating point vectors.
       //	Returns the vector.
 
-      inline float * normalize2(Vector2f v)
-	 { double l = length2(v); v[0] /= l; v[1] /= l; return v; }
-      inline double * normalize2(Vector2d v)
-	 { double l = length2(v); v[0] /= l; v[1] /= l; return v; }
+      inline float * normalize2v(Vector2f v)
+	 { double l = length2v(v); v[0] /= l; v[1] /= l; return v; }
+      inline double * normalize2v(Vector2d v)
+	 { double l = length2v(v); v[0] /= l; v[1] /= l; return v; }
 
-      inline float * normalize3(Vector3f v)
-	 { double l = length3(v); v[0] /= l; v[1] /= l; v[2] /= l; return v; }
-      inline double * normalize3(Vector3d v)
-	 { double l = length3(v); v[0] /= l; v[1] /= l; v[2] /= l; return v; }
+      inline float * normalize3v(Vector3f v)
+	 { double l = length3v(v); v[0] /= l; v[1] /= l; v[2] /= l; return v; }
+      inline double * normalize3v(Vector3d v)
+	 { double l = length3v(v); v[0] /= l; v[1] /= l; v[2] /= l; return v; }
 
-      inline float * normalize4(Vector4f v)
-	 { double l = length4(v); v[0] /= l; v[1] /= l; v[2] /= l; v[3] /= l; return v; }
-      inline double * normalize4(Vector4d v)
-	 { double l = length4(v); v[0] /= l; v[1] /= l; v[2] /= l; v[3] /= l; return v; }
+      inline float * normalize4v(Vector4f v)
+	 { double l = length4v(v); v[0] /= l; v[1] /= l; v[2] /= l; v[3] /= l; return v; }
+      inline double * normalize4v(Vector4d v)
+	 { double l = length4v(v); v[0] /= l; v[1] /= l; v[2] /= l; v[3] /= l; return v; }
 
       //////////////////////////////////////////////////////////////////////////
       // cross(Vector3 v)
@@ -261,10 +260,10 @@
       //
       // I started out with two version of get_angle3:
       // get_angle3_1 implemented as
-      //  return acos(dot3(v1, v2)/(length3(v1)*length3(v2)));
+      //  return acos(dot3v(v1, v2)/(length3v(v1)*length3v(v2)));
       // get_angle3_2 implemented as
-      //  double tmp = dot3(v1, v2); 
-      //  return acos(sqrtf(tmp*tmp/(dot3(v1,v1)*dot3(v2,v2))));
+      //  double tmp = dot3v(v1, v2); 
+      //  return acos(sqrtf(tmp*tmp/(dot3v(v1,v1)*dot3v(v2,v2))));
       // 
       // I ran timing tests, since I expected the second one - the one I chose to 
       // use, in the end - to be faster since it only resulted in one sqrt call
@@ -277,92 +276,92 @@
       // to worry about?  Eh.  I might as well use it.
 
       inline double get_angle2n(Vector2f v1, Vector2f v2)
-	 { return acos(dot2(v1, v2)); }
+	 { return acos(dot2v(v1, v2)); }
       inline double get_angle2n(Vector2d v1, Vector2d v2)
-	 { return acos(dot2(v1, v2)); }
-      inline double get_angle2(Vector2i v1, Vector2i v2)
+	 { return acos(dot2v(v1, v2)); }
+      inline double get_angle2v(Vector2i v1, Vector2i v2)
 	 { 
-	 double tmp = dot2(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot2(v1,v1)*dot2(v2,v2)))); 
+	 double tmp = dot2v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot2v(v1,v1)*dot2v(v2,v2)))); 
 	 }
-      inline double get_angle2(Vector2f v1, Vector2f v2)
+      inline double get_angle2v(Vector2f v1, Vector2f v2)
 	 { 
-	 double tmp = dot2(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot2(v1,v1)*dot2(v2,v2)))); 
+	 double tmp = dot2v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot2v(v1,v1)*dot2v(v2,v2)))); 
 	 }
-      inline double get_angle2(Vector2d v1, Vector2d v2)
+      inline double get_angle2v(Vector2d v1, Vector2d v2)
 	 { 
-	 double tmp = dot2(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot2(v1,v1)*dot2(v2,v2)))); 
+	 double tmp = dot2v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot2v(v1,v1)*dot2v(v2,v2)))); 
 	 }
 
       inline double get_angle3n(Vector3f v1, Vector3f v2)
-	 { return acos(dot3(v1, v2)); }
+	 { return acos(dot3v(v1, v2)); }
       inline double get_angle3n(Vector3d v1, Vector3d v2)
-	 { return acos(dot3(v1, v2)); }
-      inline double get_angle3(Vector3i v1, Vector3i v2)
+	 { return acos(dot3v(v1, v2)); }
+      inline double get_angle3v(Vector3i v1, Vector3i v2)
 	 { 
-	 double tmp = dot3(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot3(v1,v1)*dot3(v2,v2)))); 
+	 double tmp = dot3v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot3v(v1,v1)*dot3v(v2,v2)))); 
 	 }
-      inline double get_angle3(Vector3f v1, Vector3f v2)
+      inline double get_angle3v(Vector3f v1, Vector3f v2)
 	 { 
-	 double tmp = dot3(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot3(v1,v1)*dot3(v2,v2)))); 
+	 double tmp = dot3v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot3v(v1,v1)*dot3v(v2,v2)))); 
 	 }
-      inline double get_angle3(Vector3d v1, Vector3d v2)
+      inline double get_angle3v(Vector3d v1, Vector3d v2)
 	 { 
-	 double tmp = dot3(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot3(v1,v1)*dot3(v2,v2)))); 
+	 double tmp = dot3v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot3v(v1,v1)*dot3v(v2,v2)))); 
 	 }
 
       inline double get_angle4n(Vector4f v1, Vector4f v2)
-	 { return acos(dot4(v1, v2)); }
+	 { return acos(dot4v(v1, v2)); }
       inline double get_angle4n(Vector4d v1, Vector4d v2)
-	 { return acos(dot4(v1, v2)); }
-      inline double get_angle4(Vector4i v1, Vector4i v2)
+	 { return acos(dot4v(v1, v2)); }
+      inline double get_angle4v(Vector4i v1, Vector4i v2)
 	 { 
-	 double tmp = dot4(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot4(v1,v1)*dot4(v2,v2)))); 
+	 double tmp = dot4v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot4v(v1,v1)*dot4v(v2,v2)))); 
 	 }
-      inline double get_angle4(Vector4f v1, Vector4f v2)
+      inline double get_angle4v(Vector4f v1, Vector4f v2)
 	 { 
-	 double tmp = dot4(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot4(v1,v1)*dot4(v2,v2)))); 
+	 double tmp = dot4v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot4v(v1,v1)*dot4v(v2,v2)))); 
 	 }
-      inline double get_angle4(Vector4d v1, Vector4d v2)
+      inline double get_angle4v(Vector4d v1, Vector4d v2)
 	 { 
-	 double tmp = dot4(v1, v2); 
-	 return acos(sqrtf(tmp*tmp/(dot4(v1,v1)*dot4(v2,v2)))); 
+	 double tmp = dot4v(v1, v2); 
+	 return acos(sqrtf(tmp*tmp/(dot4v(v1,v1)*dot4v(v2,v2)))); 
 	 }
 
       //////////////////////////////////////////////////////////////////////////
-      // subtract*(Vector v1, Vector v2, Vector vres)
-      //	Subtract v2 from v1 and store result in vres.
-      //	Equivalent to add(v1, scale(v2, -1), vres).
-      //        Returns vres so that the result may be used in further 
+      // subtract*(Vector v1, Vector v2)
+      //	Subtract v2 from v1,
+      //	Equivalent to add(v1, scale(v2, -1)).
+      //        Returns v1 so that the result may be used in further 
       //	calculations.
 
-      inline int * subtract2(Vector2i v1, Vector2i v2, Vector2i vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; return vres; }
-      inline float * subtract2(Vector2f v1, Vector2f v2, Vector2f vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; return vres; }
-      inline double * subtract2(Vector2d v1, Vector2d v2, Vector2d vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; return vres; }
+      inline int * subtract2v(Vector2i v1, Vector2i v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; return v1; }
+      inline float * subtract2v(Vector2f v1, Vector2f v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; return v1; }
+      inline double * subtract2v(Vector2d v1, Vector2d v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; return v1; }
 
-      inline int * subtract3(Vector3i v1, Vector3i v2, Vector3i vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2]; return vres; }
-      inline float * subtract3(Vector3f v1, Vector3f v2, Vector3f vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2];  return vres; }
-      inline double * subtract3(Vector3d v1, Vector3d v2, Vector3d vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2];  return vres; }
+      inline int * subtract3v(Vector3i v1, Vector3i v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2]; return v1; }
+      inline float * subtract3v(Vector3f v1, Vector3f v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2];  return v1; }
+      inline double * subtract3v(Vector3d v1, Vector3d v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2];  return v1; }
 
-      inline int * subtract4(Vector4i v1, Vector4i v2, Vector4i vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2]; vres[3] = v1[3] - v2[3]; return vres;}
-      inline float * subtract4(Vector4f v1, Vector4f v2, Vector4f vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2]; vres[3] = v1[3] - v2[3]; return vres; }
-      inline double * subtract4(Vector4d v1, Vector4d v2, Vector4d vres)
-	 { vres[0] = v1[0] - v2[0]; vres[1] = v1[1] - v2[1]; vres[2] = v1[2] - v2[2]; vres[3] = v1[3] - v2[3]; return vres; }
+      inline int * subtract4v(Vector4i v1, Vector4i v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2]; v1[3] -= v2[3]; return v1;}
+      inline float * subtract4v(Vector4f v1, Vector4f v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2]; v1[3] -= v2[3]; return v1; }
+      inline double * subtract4v(Vector4d v1, Vector4d v2)
+	 { v1[0] -= v2[0]; v1[1] -= v2[1]; v1[2] -= v2[2]; v1[3] -= v2[3]; return v1; }
 
 
       //////////////////////////////////////////////////////////////////////////
@@ -370,25 +369,25 @@
       //	Copies the value of v2 into v1.
       //	Returns v1.
 
-      inline int * assign2(Vector2i v1, Vector2i v2)
+      inline int * assign2v(Vector2i v1, Vector2i v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; return v1; }
-      inline float * assign2(Vector2f v1, Vector2f v2)
+      inline float * assign2v(Vector2f v1, Vector2f v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; return v1; }
-      inline double * assign2(Vector2d v1, Vector2d v2)
+      inline double * assign2v(Vector2d v1, Vector2d v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; return v1; }
 
-      inline int * assign3(Vector3i v1, Vector3i v2)
+      inline int * assign3v(Vector3i v1, Vector3i v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; return v1; }
-      inline float * assign3(Vector3f v1, Vector3f v2)
+      inline float * assign3v(Vector3f v1, Vector3f v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; return v1; }
-      inline double * assign3(Vector3d v1, Vector3d v2)
+      inline double * assign3v(Vector3d v1, Vector3d v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; return v1; }
 
-      inline int * assign4(Vector4i v1, Vector4i v2)
+      inline int * assign4v(Vector4i v1, Vector4i v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; v1[3] = v2[3]; return v1; }
-      inline float * assign4(Vector4f v1, Vector4f v2)
+      inline float * assign4v(Vector4f v1, Vector4f v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; v1[3] = v2[3]; return v1; }
-      inline double * assign4(Vector4d v1, Vector4d v2)
+      inline double * assign4v(Vector4d v1, Vector4d v2)
 	 { v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2]; v1[3] = v2[3]; return v1; }
 
 
@@ -399,20 +398,20 @@
       //	Only defined for floating point vectors.
       //	Returns vres.
 
-      inline float * proj2(Vector2f v1, Vector2f v2, Vector2f vres)
-	 { assign2(vres, v2); scale2(vres, (dot2(v1, vres)/dot2(vres,vres))); return vres; }
-      inline double * proj2(Vector2d v1, Vector2d v2, Vector2d vres)
-	 { assign2(vres, v2); scale2(vres, (dot2(v1, vres)/dot2(vres,vres))); return vres; }
+      inline float * proj2v(Vector2f v1, Vector2f v2, Vector2f vres)
+	 { assign2v(vres, v2); scale2v(vres, (dot2v(v1, vres)/dot2v(vres,vres))); return vres; }
+      inline double * proj2v(Vector2d v1, Vector2d v2, Vector2d vres)
+	 { assign2v(vres, v2); scale2v(vres, (dot2v(v1, vres)/dot2v(vres,vres))); return vres; }
 
-      inline float * proj3(Vector3f v1, Vector3f v2, Vector3f vres)
-	 { assign3(vres, v2); scale3(vres, (dot3(v1, vres)/dot3(vres,vres))); return vres; }
-      inline double * proj3(Vector3d v1, Vector3d v2, Vector3d vres)
-	 { assign3(vres, v2); scale3(vres, (dot3(v1, vres)/dot3(vres,vres))); return vres; }
+      inline float * proj3v(Vector3f v1, Vector3f v2, Vector3f vres)
+	 { assign3v(vres, v2); scale3v(vres, (dot3v(v1, vres)/dot3v(vres,vres))); return vres; }
+      inline double * proj3v(Vector3d v1, Vector3d v2, Vector3d vres)
+	 { assign3v(vres, v2); scale3v(vres, (dot3v(v1, vres)/dot3v(vres,vres))); return vres; }
 
-      inline float * proj4(Vector4f v1, Vector4f v2, Vector4f vres)
-	 { assign4(vres, v2); scale4(vres, (dot4(v1, vres)/dot4(vres,vres))); return vres; }
-      inline double * proj4(Vector4d v1, Vector4d v2, Vector4d vres)
-	 { assign4(vres, v2); scale4(vres, (dot4(v1, vres)/dot4(vres,vres))); return vres; }
+      inline float * proj4v(Vector4f v1, Vector4f v2, Vector4f vres)
+	 { assign4v(vres, v2); scale4v(vres, (dot4v(v1, vres)/dot4v(vres,vres))); return vres; }
+      inline double * proj4v(Vector4d v1, Vector4d v2, Vector4d vres)
+	 { assign4v(vres, v2); scale4v(vres, (dot4v(v1, vres)/dot4v(vres,vres))); return vres; }
 
       } // Vector
 
diff -r f6e6f3c8f7eb -r 11e216148d1c notes
--- a/notes	Sun Sep 04 12:46:08 2011 -0500
+++ b/notes	Tue Sep 06 11:26:39 2011 -0500
@@ -11,8 +11,9 @@
 * normalize
 * cross
 * angle between two vectors
-assignment
-projection of one vector onto another.
+* subtraction
+* assignment
+* projection of one vector onto another.
 triple product i.e. dot(u, cross(v, w))   Is this really that common?
 
 Class with operator overloading?
@@ -20,13 +21,17 @@
 --------------------------------------------------------------------------------
 Matrices
 
+* getstring
+* assign
+* add
+* subtract
+* scalar multiplication
+* vector multiplication
+* matrix multiplication
+* transpose
 det
-transpose
 inverse
-add
-matrix multiplication
-vector multiplication
-scalar multiplication
 
 --------------------------------------------------------------------------------
 Quaternions
+
diff -r f6e6f3c8f7eb -r 11e216148d1c src/Math.cpp
--- a/src/Math.cpp	Sun Sep 04 12:46:08 2011 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,95 +0,0 @@
-#include "Math.h"
-
-#include <iostream>
-#include <string>
-#include <sstream>
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr2(Vector2i v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr2(Vector2f v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr2(Vector2d v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr3(Vector3i v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr3(Vector3f v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr3(Vector3d v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr4(Vector4i v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr4(Vector4f v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " }";
-   s = ss.str();
-   return s;   
-   }
-
-////////////////////////////////////////////////////////////////////////////////
-std::string & arda::Vector::vectostr4(Vector4d v, std::string & s)
-   {
-   s.clear();
-   std::stringstream ss;
-   ss << "{ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " }";
-   s = ss.str();
-   return s;   
-   }
diff -r f6e6f3c8f7eb -r 11e216148d1c src/Matrix.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Matrix.cpp	Tue Sep 06 11:26:39 2011 -0500
@@ -0,0 +1,515 @@
+#include "Matrix.h"
+
+#include <iostream>
+using namespace std;
+#include <string>
+#include <sstream>
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring2m(Matrix2i m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << " ], "
+	 "[ " << m[2] << ", " << m[3] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring2m(Matrix2f m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << " ], "
+	 "[ " << m[2] << ", " << m[3] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring2m(Matrix2d m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << " ], "
+	 "[ " << m[2] << ", " << m[3] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring3m(Matrix3i m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
+	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
+	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring3m(Matrix3f m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
+	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
+	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring3m(Matrix3d m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
+	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
+	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring4m(Matrix4i m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
+	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
+	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
+	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring4m(Matrix4f m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
+	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
+	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
+	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Matrix::getstring4m(Matrix4d m, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
+	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
+	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
+	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec2mv(Matrix2i m, Vector2i v, Vector2i vres)
+   {
+   vres[0] = m[0]*v[0] + m[2]*v[1];
+   vres[1] = m[1]*v[0] + m[3]*v[1];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec2mv(Matrix2f m, Vector2f v, Vector2f vres)
+   {
+   vres[0] = m[0]*v[0] + m[2]*v[1];
+   vres[1] = m[1]*v[0] + m[3]*v[1];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec2mv(Matrix2d m, Vector2d v, Vector2d vres)
+   {
+   vres[0] = m[0]*v[0] + m[2]*v[1];
+   vres[1] = m[1]*v[0] + m[3]*v[1];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec3mv(Matrix3i m, Vector3i v, Vector3i vres)
+   {
+   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
+   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
+   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec3mv(Matrix3f m, Vector3f v, Vector3f vres)
+   {
+   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
+   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
+   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec3mv(Matrix3d m, Vector3d v, Vector3d vres)
+   {
+   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
+   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
+   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec4mv(Matrix4i m, Vector4i v, Vector4i vres)
+   {
+   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
+   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
+   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
+   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec4mv(Matrix4f m, Vector4f v, Vector4f vres)
+   {
+   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
+   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
+   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
+   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
+   return vres;
+  
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec4mv(Matrix4d m, Vector4d v, Vector4d vres)
+   {
+   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
+   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
+   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
+   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
+   return vres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec2vm(Vector2i v, Matrix2i m, Vector2i vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1];
+   vres[1] = v[0]*m[2] + v[1]*m[3];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec2vm(Vector2f v, Matrix2f m, Vector2f vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1];
+   vres[1] = v[0]*m[2] + v[1]*m[3];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec2vm(Vector2d v, Matrix2d m, Vector2d vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1];
+   vres[1] = v[0]*m[2] + v[1]*m[3];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec3vm(Vector3i v, Matrix3i m, Vector3i vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
+   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
+   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec3vm(Vector3f v, Matrix3f m, Vector3f vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
+   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
+   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec3vm(Vector3d v, Matrix3d m, Vector3d vres)
+   {
+   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
+   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
+   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multvec4vm(Vector4i v, Matrix4i m, Vector4i vres)
+   {
+   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
+   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
+   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
+   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multvec4vm(Vector4f v, Matrix4f m, Vector4f vres)
+   {
+   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
+   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
+   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
+   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
+   return vres;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multvec4vm(Vector4d v, Matrix4d m, Vector4d vres)
+   {
+   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
+   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
+   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
+   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
+   return vres;   
+   }
+
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multmat2(Matrix2i m1, Matrix2i m2, Matrix2i mres)
+   {
+   const int size=2;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multmat2(Matrix2f m1, Matrix2f m2, Matrix2f mres)
+   {
+   const int size=2;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multmat2(Matrix2d m1, Matrix2d m2, Matrix2d mres)
+   {
+   const int size=2;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multmat3(Matrix3i m1, Matrix3i m2, Matrix3i mres)
+   {
+   const int size=3;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multmat3(Matrix3f m1, Matrix3f m2, Matrix3f mres)
+   {
+   const int size=3;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multmat3(Matrix3d m1, Matrix3d m2, Matrix3d mres)
+   {
+   const int size=3;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::multmat4(Matrix4i m1, Matrix4i m2, Matrix4i mres)
+   {
+   const int size=4;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::multmat4(Matrix4f m1, Matrix4f m2, Matrix4f mres)
+   {
+   const int size=4;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::multmat4(Matrix4d m1, Matrix4d m2, Matrix4d mres)
+   {
+   cout << "here" << endl;
+   const int size=4;
+   int i, j, k;
+   for (i=0; i<size*size; ++i)
+      mres[i] = 0;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; ++j)
+	 for (k=0; k<size; ++k)
+	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::transpose2(Matrix2i m, Matrix2i mres)
+   {
+   const int size=2;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::transpose2(Matrix2f m, Matrix2f mres)
+   {
+   const int size=2;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::transpose2(Matrix2d m, Matrix2d mres)
+   {
+   const int size=2;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::transpose3(Matrix3i m, Matrix3i mres)
+   {
+   const int size=3;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::transpose3(Matrix3f m, Matrix3f mres)
+   {
+   const int size=3;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::transpose3(Matrix3d m, Matrix3d mres)
+   {
+   const int size=3;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int * arda::Matrix::transpose4(Matrix4i m, Matrix4i mres)
+   {
+   const int size=4;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+float * arda::Matrix::transpose4(Matrix4f m, Matrix4f mres)
+   {
+   const int size=4;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+double * arda::Matrix::transpose4(Matrix4d m, Matrix4d mres)
+   {
+   const int size=4;
+   int i, j;
+   for (i=0; i<size; ++i)
+      for (j=0; j<size; j++)
+	 mres[size*i+j] = m[size*j+i];
+   return mres;
+   }
+
+
diff -r f6e6f3c8f7eb -r 11e216148d1c src/Vector.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Vector.cpp	Tue Sep 06 11:26:39 2011 -0500
@@ -0,0 +1,95 @@
+#include "Vector.h"
+
+#include <iostream>
+#include <string>
+#include <sstream>
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring2v(Vector2i v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring2v(Vector2f v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring2v(Vector2d v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring3v(Vector3i v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring3v(Vector3f v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring3v(Vector3d v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring4v(Vector4i v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring4v(Vector4f v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " ]";
+   s = ss.str();
+   return s;   
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+std::string & arda::Vector::getstring4v(Vector4d v, std::string & s)
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << " ]";
+   s = ss.str();
+   return s;   
+   }
diff -r f6e6f3c8f7eb -r 11e216148d1c src/main.cpp
--- a/src/main.cpp	Sun Sep 04 12:46:08 2011 -0500
+++ b/src/main.cpp	Tue Sep 06 11:26:39 2011 -0500
@@ -1,12 +1,15 @@
 #include <iostream>
 #include <iomanip>
+#include <limits>
 #include <string>
 using namespace std;
 
 #include "Math.h"
 using namespace arda::Vector;
+using namespace arda::Matrix;
 
-int main (int argc, char * argv[]) 
+////////////////////////////////////////////////////////////////////////////////
+void test_vector(void)
    {
    string s1, s2;
 
@@ -16,28 +19,28 @@
 	 "Testing Vector::scale*()" << endl;
 
    Vector2i v2i = { 1, 2 };
-   cout << setw(40) << "v2i: " << vectostr2(v2i, s1) << endl;
-   scale2(v2i, 2);
-   cout << setw(40) << "scale2(v2i, 2): " << vectostr2(v2i, s1) << endl;
-   scale2(v2i, 2.4);
-   cout << setw(40) << "scale2(vi, 2.4): " << vectostr2(v2i, s1) << endl;
-   scale2(scale2(v2i, 2), 3);
-   cout << setw(40) << "scale2(scale2(v2i, 2), 3): " << vectostr2(v2i, s1) << endl;
+   cout << setw(40) << "v2i: " << getstring2v(v2i, s1) << endl;
+   scale2v(v2i, 2);
+   cout << setw(40) << "scale2v(v2i, 2): " << getstring2v(v2i, s1) << endl;
+   scale2v(v2i, 2.4);
+   cout << setw(40) << "scale2v(v2i, 2.4): " << getstring2v(v2i, s1) << endl;
+   scale2v(scale2v(v2i, 2), 3);
+   cout << setw(40) << "scale2v(scale2v(v2i, 2), 3): " << getstring2v(v2i, s1) << endl;
 
 
    Vector3f v3f = { 3.0, 4.1, 5.2 };
-   cout << setw(40) << "v3f: " << vectostr3(v3f, s1) << endl;
-   scale3(v3f, 2.2f);
-   cout << setw(40) << "scale3(v3f,  2.2f): " << vectostr3(v3f, s1) << endl;
-   scale3(v3f, 2);
-   cout << setw(40) << "scale3(v2f, 2): " << vectostr3(v3f, s1) << endl;
+   cout << setw(40) << "v3f: " << getstring3v(v3f, s1) << endl;
+   scale3v(v3f, 2.2f);
+   cout << setw(40) << "scale3v(v3f,  2.2f): " << getstring3v(v3f, s1) << endl;
+   scale3v(v3f, 2);
+   cout << setw(40) << "scale3v(v2f, 2): " << getstring3v(v3f, s1) << endl;
 
    Vector4d v4d = { 6.0L, 7.1L, 8.2L, 9.3L };
-   cout << setw(40) << "vd4: " << vectostr4(v4d, s1) << endl;
-   scale4(v4d, 2.0);
-   cout << setw(40) << "scale4(v4d, 2.0): " << vectostr4(v4d, s1) << endl;
-   scale4(v4d, 2);
-   cout << setw(40) << "scale4(v4d, 2): " << vectostr4(v4d, s1) << endl;
+   cout << setw(40) << "vd4: " << getstring4v(v4d, s1) << endl;
+   scale4v(v4d, 2.0);
+   cout << setw(40) << "scale4v(v4d, 2.0): " << getstring4v(v4d, s1) << endl;
+   scale4v(v4d, 2);
+   cout << setw(40) << "scale4v(v4d, 2): " << getstring4v(v4d, s1) << endl;
    }
 
 
@@ -47,19 +50,19 @@
 	 "Testing Vector::dot*()" << endl;
 
    Vector2i v2i[2] = { { 1, 2 }, { 3, 4 } };
-   cout << setw(40) << "v2i[0]: " << vectostr2(v2i[0], s1) << endl;
-   cout << setw(40) << "v2i[1]: " << vectostr2(v2i[1], s1) << endl;
-   cout << setw(40) << "dot2(v2i[0], v2i[1]): " << dot2(v2i[0], v2i[1]) << endl;
+   cout << setw(40) << "v2i[0]: " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "v2i[1]: " << getstring2v(v2i[1], s1) << endl;
+   cout << setw(40) << "dot2v(v2i[0], v2i[1]): " << dot2v(v2i[0], v2i[1]) << endl;
    
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   cout << setw(40) << "dot3(v3f[0], v3f[1]): " << dot3(v3f[0], v3f[1]) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   cout << setw(40) << "dot3v(v3f[0], v3f[1]): " << dot3v(v3f[0], v3f[1]) << endl;
    
    Vector4d v4d[3] = { { 1.1, 2.2, 3.3, 4.4 }, { 5.5, 6.6, 7.7, 8.8 } };
-   cout << setw(40) << "v4d[0]: " << vectostr4(v4d[0], s1) << endl;
-   cout << setw(40) << "v4d[1]: " << vectostr4(v4d[1], s1) <<  endl;
-   cout << setw(40) << "dot4(v4d[0], v4d[1]): " << dot4(v4d[0], v4d[1]) << endl;
+   cout << setw(40) << "v4d[0]: " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[1]: " << getstring4v(v4d[1], s1) <<  endl;
+   cout << setw(40) << "dot4v(v4d[0], v4d[1]): " << dot4v(v4d[0], v4d[1]) << endl;
    
    }
 
@@ -69,16 +72,16 @@
 	 "Testing Vector::length*()" << endl;
 
    Vector2i v2i = { 1, 2 };
-   cout << setw(40) << "v2i: " << vectostr2(v2i, s1) << endl;
-   cout << setw(40) << "length2(v2i): " << length2(v2i) << endl;
+   cout << setw(40) << "v2i: " << getstring2v(v2i, s1) << endl;
+   cout << setw(40) << "length2v(v2i): " << length2v(v2i) << endl;
    
    Vector3f v3f = { 1.1f, 2.2f, 3.3f };
-   cout << setw(40) << "v3f: " << vectostr3(v3f, s1) << endl;
-   cout << setw(40) << "length3(v3f): " << length3(v3f) << endl;
+   cout << setw(40) << "v3f: " << getstring3v(v3f, s1) << endl;
+   cout << setw(40) << "length3v(v3f): " << length3v(v3f) << endl;
    
    Vector4d v4d = { 1.1, 2.2, 3.3, 4.4 };
-   cout << setw(40) << "v4d: " << vectostr4(v4d, s1) << endl;
-   cout << setw(40) << "length4(v4d): " << length4(v4d) << endl;
+   cout << setw(40) << "v4d: " << getstring4v(v4d, s1) << endl;
+   cout << setw(40) << "length4v(v4d): " << length4v(v4d) << endl;
    
    }
 
@@ -88,28 +91,25 @@
 	 "Testing Vector::add*()" << endl;
 
    Vector2i v2i[3] = { { 1, 2 }, { 3, 4 } };
-   Vector2i v2ires;
-   cout << setw(40) << "v2i[0]: " << vectostr2(v2i[0], s1) << endl;
-   cout << setw(40) << "v2i[1]: " << vectostr2(v2i[1], s1) << endl;
-   add2(v2i[0], v2i[1], v2ires);
-   cout << setw(40) << "add2(v2i[0], v2i[1], v2ires): " << vectostr2(v2ires, s1) << endl;
-   cout << setw(40) << "length2(add2(v2i[0], v2i[1], v2ires)): " << length2(add2(v2i[0], v2i[1], v2ires)) << endl;
+   cout << setw(40) << "v2i[0]: " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "v2i[1]: " << getstring2v(v2i[1], s1) << endl;
+   add2v(v2i[0], v2i[1]);
+   cout << setw(40) << "add2v(v2i[0], v2i[1]): " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "length2v(add2v(v2i[0], v2i[1])): " << length2v(add2v(v2i[0], v2i[1])) << endl;
 
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
-   Vector3f v3fres;
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   add3(v3f[0], v3f[1], v3fres);
-   cout << setw(40) << "add3(v3f[0], v3f[1], v3fres): " << vectostr3(v3fres, s1) << endl;
-   cout << setw(40) << "length3(add3(v3f[0], v3f[1], v3fres)): " << length3(add3(v3f[0], v3f[1], v3fres)) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   add3v(v3f[0], v3f[1]);
+   cout << setw(40) << "add3v(v3f[0], v3f[1]): " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "length3v(add3v(v3f[0], v3f[1])): " << length3v(add3v(v3f[0], v3f[1])) << endl;
 
    Vector4d v4d[3] = { { 1.1, 2.2, 3.3, 4.4 }, { 5.5, 6.6, 7.7, 8.8 } };
-   Vector4d v4dres;
-   cout << setw(40) << "v4d[0]: " << vectostr4(v4d[0], s1) << endl;
-   cout << setw(40) << "v4d[1]: " << vectostr4(v4d[1], s1) << endl;
-   add4(v4d[0], v4d[1], v4dres);
-   cout << setw(40) << "add4(v4d[0], v4d[1], v4dres): " << vectostr4(v4dres, s1) << endl;
-   cout << setw(40) << "length4(add4(v4d[0], v4d[1], v4dres)): " << length4(add4(v4d[0], v4d[1], v4dres)) << endl;
+   cout << setw(40) << "v4d[0]: " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[1]: " << getstring4v(v4d[1], s1) << endl;
+   add4v(v4d[0], v4d[1]);
+   cout << setw(40) << "add4v(v4d[0], v4d[1]): " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "length4v(add4v(v4d[0], v4d[1])): " << length4v(add4v(v4d[0], v4d[1])) << endl;
 
    }
 
@@ -119,25 +119,25 @@
 	 "Testing Vector::normalize*()" << endl;
 
    Vector2f v2f = { 1.1, 2.2 };
-   cout << setw(40) << "v2f: " << vectostr2(v2f, s1) << endl;
-   normalize2(v2f);
-   cout << setw(40) << "normalize2(v2f): " << vectostr2(v2f, s1) << endl;
-   scale2(normalize2(v2f), 2.0);
-   cout << setw(40) << "scale(normalize2(v2f), 2.0): " << vectostr2(v2f, s1) << endl;
+   cout << setw(40) << "v2f: " << getstring2v(v2f, s1) << endl;
+   normalize2v(v2f);
+   cout << setw(40) << "normalize2v(v2f): " << getstring2v(v2f, s1) << endl;
+   scale2v(normalize2v(v2f), 2.0);
+   cout << setw(40) << "scale(normalize2v(v2f), 2.0): " << getstring2v(v2f, s1) << endl;
    
    Vector3f v3f = { 1.1f, 2.2f, 3.3f };
-   cout << setw(40) << "v3f: " << vectostr3(v3f, s1) << endl;
-   normalize3(v3f);
-   cout << setw(40) << "normalize3(v3f): " << vectostr3(v3f, s1) << endl;
-   scale3(normalize3(v3f), 2.0);
-   cout << setw(40) << "scale(normalize3(v3f), 2.0): " << vectostr3(v3f, s1) << endl;
+   cout << setw(40) << "v3f: " << getstring3v(v3f, s1) << endl;
+   normalize3v(v3f);
+   cout << setw(40) << "normalize3v(v3f): " << getstring3v(v3f, s1) << endl;
+   scale3v(normalize3v(v3f), 2.0);
+   cout << setw(40) << "scale(normalize3v(v3f), 2.0): " << getstring3v(v3f, s1) << endl;
    
    Vector4d v4d = { 1.1, 2.2, 3.3, 4.4 };
-   cout << setw(40) << "v4d: " << vectostr4(v4d, s1) << endl;
-   normalize4(v4d);
-   cout << setw(40) << "normalize4(v4d): " << vectostr4(v4d, s1) << endl;
-   scale4(normalize4(v4d), 2.0);
-   cout << setw(40) << "scale4(normalize4(v4d), 2.0): " << vectostr4(v4d, s1) << endl;
+   cout << setw(40) << "v4d: " << getstring4v(v4d, s1) << endl;
+   normalize4v(v4d);
+   cout << setw(40) << "normalize4v(v4d): " << getstring4v(v4d, s1) << endl;
+   scale4v(normalize4v(v4d), 2.0);
+   cout << setw(40) << "scale4v(normalize4v(v4d), 2.0): " << getstring4v(v4d, s1) << endl;
    
    }
 
@@ -147,16 +147,16 @@
 	 "Testing Vector::cross()" << endl;
 
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f }, { 0.0f, 0.0f, 0.0f } };
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
    cross(v3f[0], v3f[1], v3f[2]);
-   cout << setw(40) << "cross(v3f[0], v3f[1], v3f[2]): " << vectostr3(v3f[2], s1) << endl;
+   cout << setw(40) << "cross(v3f[0], v3f[1], v3f[2]): " << getstring3v(v3f[2], s1) << endl;
 
    Vector3d v3d[3] = { { 1.1, 2.2, 3.3 }, { 4.4, 5.5, 6.6 }, { 0.0, 0.0, 0.0 } };
-   cout << setw(40) << "v3d[0]: " << vectostr3(v3d[0], s1) << endl;
-   cout << setw(40) << "v3d[1]: " << vectostr3(v3d[1], s1) << endl;
+   cout << setw(40) << "v3d[0]: " << getstring3v(v3d[0], s1) << endl;
+   cout << setw(40) << "v3d[1]: " << getstring3v(v3d[1], s1) << endl;
    cross(v3d[0], v3d[1], v3d[2]);
-   cout << setw(40) << "cross(v3d[0], v3d[1], v3d[2]): " << vectostr3(v3d[2], s1) << endl;
+   cout << setw(40) << "cross(v3d[0], v3d[1], v3d[2]): " << getstring3v(v3d[2], s1) << endl;
 
    }
 
@@ -167,28 +167,28 @@
 	 "Testing Vector::get_angle*()" << endl;
 
    Vector3i v3i[2] = { { 1, 2, 3 }, { 4, 5, 6 } };
-   cout << setw(40) << "v3i[0]: " << vectostr3(v3i[0], s1) << endl;
-   cout << setw(40) << "v3i[1]: " << vectostr3(v3i[1], s1) << endl;
-   cout << setw(40) << "get_angle3(v3i[0], v3i[1]): " << get_angle3(v3i[0], v3i[1]) << endl;
+   cout << setw(40) << "v3i[0]: " << getstring3v(v3i[0], s1) << endl;
+   cout << setw(40) << "v3i[1]: " << getstring3v(v3i[1], s1) << endl;
+   cout << setw(40) << "get_angle3v(v3i[0], v3i[1]): " << get_angle3v(v3i[0], v3i[1]) << endl;
 
    Vector3f v3f[2] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   cout << setw(40) << "get_angle3(v3f[0], v3f[1]): " << get_angle3(v3f[0], v3f[1]) << endl;
-   normalize3(v3f[0]);
-   normalize3(v3f[1]);
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   cout << setw(40) << "get_angle3v(v3f[0], v3f[1]): " << get_angle3v(v3f[0], v3f[1]) << endl;
+   normalize3v(v3f[0]);
+   normalize3v(v3f[1]);
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
    cout << setw(40) << "get_angle3n(v3f[0], v3f[1]): " << get_angle3n(v3f[0], v3f[1]) << endl;
 
    Vector3d v3d[2] = { { 1.1, 2.2, 3.3 }, { 4.4, 5.5, 6.6 } };
-   cout << setw(40) << "v3d[0]: " << vectostr3(v3d[0], s1) << endl;
-   cout << setw(40) << "v3d[1]: " << vectostr3(v3d[1], s1) << endl;
-   cout << setw(40) << "get_angle3(v3d[0], v3d[1]): " << get_angle3(v3d[0], v3d[1]) << endl;
-   normalize3(v3d[0]);
-   normalize3(v3d[1]);
-   cout << setw(40) << "v3d[0]: " << vectostr3(v3d[0], s1) << endl;
-   cout << setw(40) << "v3d[1]: " << vectostr3(v3d[1], s1) << endl;
+   cout << setw(40) << "v3d[0]: " << getstring3v(v3d[0], s1) << endl;
+   cout << setw(40) << "v3d[1]: " << getstring3v(v3d[1], s1) << endl;
+   cout << setw(40) << "get_angle3v(v3d[0], v3d[1]): " << get_angle3v(v3d[0], v3d[1]) << endl;
+   normalize3v(v3d[0]);
+   normalize3v(v3d[1]);
+   cout << setw(40) << "v3d[0]: " << getstring3v(v3d[0], s1) << endl;
+   cout << setw(40) << "v3d[1]: " << getstring3v(v3d[1], s1) << endl;
    cout << setw(40) << "get_angle3n(v3d[0], v3d[1]): " << get_angle3n(v3d[0], v3d[1]) << endl;
    }
 
@@ -198,28 +198,25 @@
 	 "Testing Vector::subtract*()" << endl;
 
    Vector2i v2i[3] = { { 1, 2 }, { 3, 4 } };
-   Vector2i v2ires;
-   cout << setw(40) << "v2i[0]: " << vectostr2(v2i[0], s1) << endl;
-   cout << setw(40) << "v2i[1]: " << vectostr2(v2i[1], s1) << endl;
-   subtract2(v2i[0], v2i[1], v2ires);
-   cout << setw(40) << "subtract2(v2i[0], v2i[1], v2ires): " << vectostr2(v2ires, s1) << endl;
-   cout << setw(40) << "length2(subtract2(v2i[0], v2i[1], v2ires)): " << length2(subtract2(v2i[0], v2i[1], v2ires)) << endl;
+   cout << setw(40) << "v2i[0]: " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "v2i[1]: " << getstring2v(v2i[1], s1) << endl;
+   subtract2v(v2i[0], v2i[1]);
+   cout << setw(40) << "subtract2v(v2i[0], v2i[1]): " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "length2v(subtract2v(v2i[0], v2i[1])): " << length2v(subtract2v(v2i[0], v2i[1])) << endl;
 
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
-   Vector3f v3fres;
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   subtract3(v3f[0], v3f[1], v3fres);
-   cout << setw(40) << "subtract3(v3f[0], v3f[1], v3fres): " << vectostr3(v3fres, s1) << endl;
-   cout << setw(40) << "length3(subtract3(v3f[0], v3f[1], v3fres)): " << length3(subtract3(v3f[0], v3f[1], v3fres)) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   subtract3v(v3f[0], v3f[1]);
+   cout << setw(40) << "subtract3v(v3f[0], v3f[1]): " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "length3v(subtract3v(v3f[0], v3f[1])): " << length3v(subtract3v(v3f[0], v3f[1])) << endl;
 
    Vector4d v4d[3] = { { 1.1, 2.2, 3.3, 4.4 }, { 5.5, 6.6, 7.7, 8.8 } };
-   Vector4d v4dres;
-   cout << setw(40) << "v4d[0]: " << vectostr4(v4d[0], s1) << endl;
-   cout << setw(40) << "v4d[1]: " << vectostr4(v4d[1], s1) << endl;
-   subtract4(v4d[0], v4d[1], v4dres);
-   cout << setw(40) << "subtract4(v4d[0], v4d[1], v4dres): " << vectostr4(v4dres, s1) << endl;
-   cout << setw(40) << "length4(subtract4(v4d[0], v4d[1], v4dres)): " << length4(subtract4(v4d[0], v4d[1], v4dres)) << endl;
+   cout << setw(40) << "v4d[0]: " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[1]: " << getstring4v(v4d[1], s1) << endl;
+   subtract4v(v4d[0], v4d[1]);
+   cout << setw(40) << "subtract4v(v4d[0], v4d[1]): " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "length4v(subtract4v(v4d[0], v4d[1])): " << length4v(subtract4v(v4d[0], v4d[1])) << endl;
 
    }
 
@@ -229,22 +226,22 @@
 	 "Testing Vector::assign*()" << endl;
 
    Vector2i v2i[3] = { { 1, 2 }, { 3, 4 } };
-   cout << setw(40) << "v2i[0]: " << vectostr2(v2i[0], s1) << endl;
-   cout << setw(40) << "v2i[1]: " << vectostr2(v2i[1], s1) << endl;
-   assign2(v2i[0], v2i[1]);
-   cout << setw(40) << "assign2(v2i[0], v2i[1]): " << vectostr2(v2i[0], s1) << endl;
+   cout << setw(40) << "v2i[0]: " << getstring2v(v2i[0], s1) << endl;
+   cout << setw(40) << "v2i[1]: " << getstring2v(v2i[1], s1) << endl;
+   assign2v(v2i[0], v2i[1]);
+   cout << setw(40) << "assign2v(v2i[0], v2i[1]): " << getstring2v(v2i[0], s1) << endl;
 
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   assign3(v3f[0], v3f[1]);
-   cout << setw(40) << "assign3(v3f[0], v3f[1]): " << vectostr3(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   assign3v(v3f[0], v3f[1]);
+   cout << setw(40) << "assign3v(v3f[0], v3f[1]): " << getstring3v(v3f[0], s1) << endl;
 
    Vector4d v4d[3] = { { 1.1, 2.2, 3.3, 4.4 }, { 5.5, 6.6, 7.7, 8.8 } };
-   cout << setw(40) << "v4d[0]: " << vectostr4(v4d[0], s1) << endl;
-   cout << setw(40) << "v4d[1]: " << vectostr4(v4d[1], s1) << endl;
-   assign4(v4d[0], v4d[1]);
-   cout << setw(40) << "assign4(v4d[0], v4d[1]): " << vectostr4(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[0]: " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[1]: " << getstring4v(v4d[1], s1) << endl;
+   assign4v(v4d[0], v4d[1]);
+   cout << setw(40) << "assign4v(v4d[0], v4d[1]): " << getstring4v(v4d[0], s1) << endl;
 
    }
 
@@ -255,24 +252,347 @@
 
    Vector2f v2f[3] = { { 1.1f, 2.2f }, { 3.3f, 4.4f } };
    Vector2f v2fres;
-   cout << setw(40) << "v2f[0]: " << vectostr2(v2f[0], s1) << endl;
-   cout << setw(40) << "v2f[1]: " << vectostr2(v2f[1], s1) << endl;
-   proj3(v2f[0], v2f[1], v2fres);
-   cout << setw(40) << "proj3(v2f[0], v2f[1], v2fres): " << vectostr2(v2fres, s1) << endl;
+   cout << setw(40) << "v2f[0]: " << getstring2v(v2f[0], s1) << endl;
+   cout << setw(40) << "v2f[1]: " << getstring2v(v2f[1], s1) << endl;
+   proj3v(v2f[0], v2f[1], v2fres);
+   cout << setw(40) << "proj3v(v2f[0], v2f[1], v2fres): " << getstring2v(v2fres, s1) << endl;
 
    Vector3f v3f[3] = { { 1.1f, 2.2f, 3.3f }, { 4.4f, 5.5f, 6.6f } };
    Vector3f v3fres;
-   cout << setw(40) << "v3f[0]: " << vectostr3(v3f[0], s1) << endl;
-   cout << setw(40) << "v3f[1]: " << vectostr3(v3f[1], s1) << endl;
-   proj3(v3f[0], v3f[1], v3fres);
-   cout << setw(40) << "proj3(v3f[0], v3f[1], v3fres): " << vectostr3(v3fres, s1) << endl;
+   cout << setw(40) << "v3f[0]: " << getstring3v(v3f[0], s1) << endl;
+   cout << setw(40) << "v3f[1]: " << getstring3v(v3f[1], s1) << endl;
+   proj3v(v3f[0], v3f[1], v3fres);
+   cout << setw(40) << "proj3v(v3f[0], v3f[1], v3fres): " << getstring3v(v3fres, s1) << endl;
 
    Vector4d v4d[3] = { { 1.1, 2.2, 3.3, 4.4 }, { 5.5, 6.6, 7.7, 8.8 } };
    Vector4d v4dres;
-   cout << setw(40) << "v4d[0]: " << vectostr4(v4d[0], s1) << endl;
-   cout << setw(40) << "v4d[1]: " << vectostr4(v4d[1], s1) << endl;
-   proj4(v4d[0], v4d[1], v4dres);
-   cout << setw(40) << "proj4(v4d[0], v4d[1], v4dres): " << vectostr4(v4dres, s1) << endl;
+   cout << setw(40) << "v4d[0]: " << getstring4v(v4d[0], s1) << endl;
+   cout << setw(40) << "v4d[1]: " << getstring4v(v4d[1], s1) << endl;
+   proj4v(v4d[0], v4d[1], v4dres);
+   cout << setw(40) << "proj4v(v4d[0], v4d[1], v4dres): " << getstring4v(v4dres, s1) << endl;
 
    }
+
    }
+
+////////////////////////////////////////////////////////////////////////////////
+void test_matrix(void)
+   {
+   std::string s;
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::assign*()" << endl;
+   Matrix2i m2i[2] = { { 0, 0, 
+			 0, 0 }, 
+		       { 1, 2, 
+			 3, 4 } };
+   Matrix3f m3f[2] = { { 0.0f, 0.0f, 0.0f, 
+			 0.0f, 0.0f, 0.0f, 
+			 0.0f, 0.0f, 0.0f }, 
+		       { 1.1f, 2.2f, 3.3f, 
+			 4.4f, 5.5f, 6.6f, 
+			 7.7f, 8.8f, 9.9f } };
+   Matrix4d m4d[2] = { { 0.0, 0.0, 0.0, 0.0, 
+			 0.0, 0.0, 0.0, 0.0, 
+			 0.0, 0.0, 0.0, 0.0, 
+			 0.0, 0.0, 0.0, 0.0 }, 
+		       { 1.1, 2.2, 3.3, 4.4, 
+			 5.5, 6.6, 7.7, 8.8, 
+			 9.9, 10.10, 11.11, 12.12,
+			 13.13, 14.14, 15.15, 16.16 } };
+   cout << setw(40) << "m2i[0]" << getstring2m(m2i[0], s) << endl;
+   cout << setw(40) << "m2i[1]" << getstring2m(m2i[1], s) << endl;
+   assign2m(m2i[0], m2i[1]);
+   cout << setw(40) << "assign2m(m2i[0], m2i[1])" << getstring2m(m2i[0], s) << endl;
+
+   cout << setw(40) << "m3f[0]" << getstring3m(m3f[0], s) << endl;
+   cout << setw(40) << "m3f[1]" << getstring3m(m3f[1], s) << endl;
+   assign3m(m3f[0], m3f[1]);
+   cout << setw(40) << "assign3m(m3f[0], m3f[1])" << getstring3m(m3f[0], s) << endl;
+
+   cout << setw(40) << "m4d[0]" << getstring4m(m4d[0], s) << endl;
+   cout << setw(40) << "m4d[1]" << getstring4m(m4d[1], s) << endl;
+   assign4m(m4d[0], m4d[1]);
+   cout << setw(40) << "assign4m(m3d[0], m3d[1])" << getstring4m(m4d[0], s) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::add*()" << endl;
+   Matrix2i m2i[2] = { { 1, 2, 
+			 3, 4 },
+		       { 5, 6, 
+			 7, 8 } };
+   Matrix3f m3f[2] = { { 1.1f, 2.2f, 3.3f, 
+			 4.4f, 5.5f, 6.6f, 
+			 7.7f, 8.8f, 9.9f },
+		       { 10.10f, 11.11f, 12.12f, 
+			 13.13f, 14.14f, 15.15f, 
+			 16.16f, 17.17f, 18.18f } };
+   Matrix4d m4d[2] = { { 1.1, 2.2, 3.3, 4.4, 
+			 5.5, 6.6, 7.7, 8.8, 
+			 9.9, 10.10, 11.11, 12.12,
+			 13.13, 14.14, 15.15, 16.16 },
+		       { 17.17, 18.18, 19.19, 20.20, 
+			 21.21, 22.22, 23.23, 24.24, 
+			 25.25, 26.26, 27.27, 28.28,
+			 29.29, 30.30, 31.31, 32.32 } };
+   cout << setw(40) << "m2i[0]" << getstring2m(m2i[0], s) << endl;
+   cout << setw(40) << "m2i[1]" << getstring2m(m2i[1], s) << endl;
+   add2m(m2i[0], m2i[1]);
+   cout << setw(40) << "add2m(m2i[0], m2i[1])" << getstring2m(m2i[0], s) << endl;
+
+   cout << setw(40) << "m3f[0]" << getstring3m(m3f[0], s) << endl;
+   cout << setw(40) << "m3f[1]" << getstring3m(m3f[1], s) << endl;
+   add3m(m3f[0], m3f[1]);
+   cout << setw(40) << "add3m(m3f[0], m3f[1])" << getstring3m(m3f[0], s) << endl;
+
+   cout << setw(40) << "m4d[0]" << getstring4m(m4d[0], s) << endl;
+   cout << setw(40) << "m4d[1]" << getstring4m(m4d[1], s) << endl;
+   add4m(m4d[0], m4d[1]);
+   cout << setw(40) << "add4m(m3d[0], m3d[1])" << getstring4m(m4d[0], s) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::subtract*()" << endl;
+   Matrix2i m2i[2] = { { 1, 2, 
+			 3, 4 },
+		       { 5, 6, 
+			 7, 8 } };
+   Matrix3f m3f[2] = { { 1.1f, 2.2f, 3.3f, 
+			 4.4f, 5.5f, 6.6f, 
+			 7.7f, 8.8f, 9.9f },
+		       { 10.10f, 11.11f, 12.12f, 
+			 13.13f, 14.14f, 15.15f, 
+			 16.16f, 17.17f, 18.18f } };
+   Matrix4d m4d[2] = { { 1.1, 2.2, 3.3, 4.4, 
+			 5.5, 6.6, 7.7, 8.8, 
+			 9.9, 10.10, 11.11, 12.12,
+			 13.13, 14.14, 15.15, 16.16 },
+		       { 17.17, 18.18, 19.19, 20.20, 
+			 21.21, 22.22, 23.23, 24.24, 
+			 25.25, 26.26, 27.27, 28.28,
+			 29.29, 30.30, 31.31, 32.32 } };
+   cout << setw(40) << "m2i[0]" << getstring2m(m2i[0], s) << endl;
+   cout << setw(40) << "m2i[1]" << getstring2m(m2i[1], s) << endl;
+   subtract2m(m2i[0], m2i[1]);
+   cout << setw(40) << "subtract2m(m2i[0], m2i[1])" << getstring2m(m2i[0], s) << endl;
+
+   cout << setw(40) << "m3f[0]" << getstring3m(m3f[0], s) << endl;
+   cout << setw(40) << "m3f[1]" << getstring3m(m3f[1], s) << endl;
+   subtract3m(m3f[0], m3f[1]);
+   cout << setw(40) << "subtract3m(m3f[0], m3f[1])" << getstring3m(m3f[0], s) << endl;
+
+   cout << setw(40) << "m4d[0]" << getstring4m(m4d[0], s) << endl;
+   cout << setw(40) << "m4d[1]" << getstring4m(m4d[1], s) << endl;
+   subtract4m(m4d[0], m4d[1]);
+   cout << setw(40) << "subtract4m(m3d[0], m3d[1])" << getstring4m(m4d[0], s) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::scale*()" << endl;
+   Matrix2i m2i = { 1, 2, 
+		    3, 4 };
+   Matrix3f m3f = { 1.1f, 2.2f, 3.3f, 
+		    4.4f, 5.5f, 6.6f, 
+		    7.7f, 8.8f, 9.9f };
+   Matrix4d m4d = { 1.1, 2.2, 3.3, 4.4, 
+		    5.5, 6.6, 7.7, 8.8, 
+		    9.9, 10.10, 11.11, 12.12,
+		    13.13, 14.14, 15.15, 16.16 };
+   cout << setw(40) << "m2i" << getstring2m(m2i, s) << endl;
+   cout << setw(40) << "scale2m(m2i, 2)" << getstring2m(scale2m(m2i, 2), s) << endl;
+
+   cout << setw(40) << "m3f" << getstring3m(m3f, s) << endl;
+   cout << setw(40) << "scale3m(m3f, 2.2f)" << getstring3m(scale3m(m3f, 2.2f), s) << endl;
+
+   cout << setw(40) << "m4d" << getstring4m(m4d, s) << endl;
+   cout << setw(40) << "scale4m(m3d, 2.2)" << getstring4m(scale4m(m4d, 2.2), s) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::multvec*()" << endl;
+   Vector2i v2ires;
+   Vector3f v3fres;
+   Vector4d v4dres;
+   Vector2i v2i = { 1, 2 };
+   Vector3f v3f = { 1.1f, 2.2f, 3.3f };
+   Vector4d v4d = { 1.1, 2.2, 3.3, 4.4 };
+   Matrix2i m2i = { 1, 2, 
+		    3, 4 };
+   Matrix3f m3f = { 1.1f, 2.2f, 3.3f, 
+		    4.4f, 5.5f, 6.6f, 
+		    7.7f, 8.8f, 9.9f };
+   Matrix4d m4d = { 1.1, 2.2, 3.3, 4.4, 
+		    5.5, 6.6, 7.7, 8.8, 
+		    9.9, 10.10, 11.11, 12.12,
+		    13.13, 14.14, 15.15, 16.16 };
+   cout << setw(40) << "v2i" << getstring2v(v2i, s) << endl;
+   cout << setw(40) << "m2i" << getstring2m(m2i, s) << endl;
+   cout << setw(40) << "multvec2mv(m2i, v2i, v2ires)" << 
+	 getstring2v(multvec2mv(m2i, v2i, v2ires), s) << endl;
+   cout << setw(40) << "multvec2mv(m2i, v2i, v2ires)" << 
+	 getstring2v(multvec2vm(v2i, m2i, v2ires), s) << endl;
+
+   cout << setw(40) << "v3f" << getstring3v(v3f, s) << endl;
+   cout << setw(40) << "m3f" << getstring3m(m3f, s) << endl;
+   cout << setw(40) << "multvec3mv(m3f, v3f, v3fres)" << 
+	 getstring3v(multvec3mv(m3f, v3f, v3fres), s) << endl;
+   cout << setw(40) << "multvec3mv(m3f, v3f, v3fres)" << 
+	 getstring3v(multvec3vm(v3f, m3f, v3fres), s) << endl;
+
+   cout << setw(40) << "v4d" << getstring4v(v4d, s) << endl;
+   cout << setw(40) << "m4d" << getstring4m(m4d, s) << endl;
+   cout << setw(40) << "multvec4mv(m4d, v4d, v4dres)" << 
+	 getstring4v(multvec4mv(m4d, v4d, v4dres), s) << endl;
+   cout << setw(40) << "multvec4mv(m4d, v4d, v4dres)" << 
+	 getstring4v(multvec4vm(v4d, m4d, v4dres), s) << endl;
+
+   }
+
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::multmat*()" << endl;
+   Matrix2i m2ires;
+   Matrix3f m3fres;
+   Matrix4d m4dres;
+   Matrix2i m2i[2] = { { 1, 2, 
+			 3, 4 },
+		       { 5, 6, 
+			 7, 8 } };
+   Matrix3f m3f[2] = { { 1.1f, 2.2f, 3.3f, 
+			 4.4f, 5.5f, 6.6f, 
+			 7.7f, 8.8f, 9.9f },
+		       { 10.10f, 11.11f, 12.12f, 
+			 13.13f, 14.14f, 15.15f, 
+			 16.16f, 17.17f, 18.18f } };
+   Matrix4d m4d[2] = { { 1.1, 2.2, 3.3, 4.4, 
+			 5.5, 6.6, 7.7, 8.8, 
+			 9.9, 10.10, 11.11, 12.12,
+			 13.13, 14.14, 15.15, 16.16 },
+		       { 17.17, 18.18, 19.19, 20.20, 
+			 21.21, 22.22, 23.23, 24.24, 
+			 25.25, 26.26, 27.27, 28.28,
+			 29.29, 30.30, 31.31, 32.32 } };
+   cout << setw(40) << "m2i[0]" << getstring2m(m2i[0], s) << endl;
+   cout << setw(40) << "m2i[1]" << getstring2m(m2i[1], s) << endl;
+   multmat2(m2i[0], m2i[1], m2ires);
+   cout << setw(40) << "multmat2(m2i[0], m2i[1], m21res)" << 
+	 getstring2m(m2ires, s) << endl;
+
+   cout << setw(40) << "m3f[0]" << getstring3m(m3f[0], s) << endl;
+   cout << setw(40) << "m3f[1]" << getstring3m(m3f[1], s) << endl;
+   multmat3(m3f[0], m3f[1], m3fres);
+   cout << setw(40) << "multmat3(m3f[0], m3f[1], m3fres)" << 
+	 getstring3m(m3fres, s) << endl;
+
+   cout << setw(40) << "m4d[0]" << getstring4m(m4d[0], s) << endl;
+   cout << setw(40) << "m4d[1]" << getstring4m(m4d[1], s) << endl;
+   multmat4(m4d[0], m4d[1], m4dres);
+   cout << setw(40) << "multmat4(m4d[0], m4d[1], m4dres)" << 
+	 getstring4m(m4dres, s) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix::transpose*()" << endl;
+   Matrix2i m2ires;
+   Matrix3f m3fres;
+   Matrix4d m4dres;
+   Matrix2i m2i = { 1, 2, 
+		    3, 4 };
+   Matrix3f m3f = { 1.1f, 2.2f, 3.3f, 
+		    4.4f, 5.5f, 6.6f, 
+		    7.7f, 8.8f, 9.9f };
+   Matrix4d m4d = { 1.1, 2.2, 3.3, 4.4, 
+		    5.5, 6.6, 7.7, 8.8, 
+		    9.9, 10.10, 11.11, 12.12,
+		    13.13, 14.14, 15.15, 16.16 };
+   cout << setw(40) << "m2i" << getstring2m(m2i, s) << endl;
+   cout << setw(40) << "transpose2(m2i, m2ires)" << 
+	 getstring2m(transpose2(m2i, m2ires), s) << endl;
+
+   cout << setw(40) << "m3f" << getstring3m(m3f, s) << endl;
+   cout << setw(40) << "transpose3(m3f, m3fres)" << 
+	 getstring3m(transpose3(m3f, m3fres), s) << endl;
+
+   cout << setw(40) << "m4d" << getstring4m(m4d, s) << endl;
+   cout << setw(40) << "transpose4(m4d, m4dres)" << 
+	 getstring4m(transpose4(m4d, m4dres), s) << endl;
+
+   }
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+void test_all(void)
+   {
+   test_vector();
+   test_matrix();
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+void show_menu(void)
+   {
+   cout << endl << endl
+	<< "Test what?" << endl
+	<< "0)      All" << endl << endl
+	<< "1)      Vector" << endl
+	<< "2)      Matrix" << endl
+	<< endl
+	<< "99      Quit" << endl;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+int main (int argc, char * argv[]) 
+   {
+   int retval = 0;
+
+   try
+      {
+      int choice = -1;
+      
+      while (choice != 99)
+	 {
+	 show_menu();
+	 cin >> choice;
+	 if(!cin)
+	    {
+	    cin.clear();
+	    cin.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
+	    }
+	 if (choice != 99)
+	    {
+	    switch (choice)
+	       {
+	       case 0:
+		  test_all();
+		  break;
+	       case 1:
+		  test_vector();
+		  break;
+	       case 2:
+		  test_matrix();
+		  break;
+	       default:
+		  cout << "Unrecognized choice.  Please try again." << endl;
+	       }
+	    choice = -1;
+	    }
+	 }
+      }
+   catch (const std::exception & error)
+      {
+      string e = "Caught exception: ";
+      e += error.what();
+      cerr << e << endl;
+      retval = 1;
+      }
+
+   return retval;;
+
+   }