# HG changeset patch
# User Eris Caffee <discordia@eldalin.com>
# Date 1318003718 18000
# Node ID 81d2aa42a8607a817ae75fc9669ee44c7fe380b5
# Parent  5377337c7de098bcfb875c3ff631d12cd4c0b463
Matrix22 working, but a few methods not yet written.  I started thinking about refactoring again t omake things even cleaner in design, with less repetition.

diff -r 5377337c7de0 -r 81d2aa42a860 CMakeLists.txt
--- a/CMakeLists.txt	Wed Sep 28 12:53:57 2011 -0500
+++ b/CMakeLists.txt	Fri Oct 07 11:08:38 2011 -0500
@@ -148,6 +148,7 @@
 set (SRCS_${App_Name}
   include/Math.h
   include/Vector.h
+  include/Matrix.h
 )
 
 include_directories (
diff -r 5377337c7de0 -r 81d2aa42a860 include/Math.h
--- a/include/Math.h	Wed Sep 28 12:53:57 2011 -0500
+++ b/include/Math.h	Fri Oct 07 11:08:38 2011 -0500
@@ -2,7 +2,7 @@
 #define MATH_H_
 
 #include "Vector.h"
-// #include "Matrix.h"
+#include "Matrix.h"
 
 namespace arda 
    {
@@ -29,6 +29,35 @@
       const double PI = 3.14159265358979323846264338327950288;
       }
 
+   ////////////////////////////////////////////////////////////////////////////////
+   // Mixed Vector / Matrix operations.
+   // * *=
+   //     Defined for Vector * Matrix and Matrix * Vector, but these
+   //     operators are not part of the class.  They are defined independently.
+
+   // Vector * Matrix
+   template <typename T> inline Vector2<T>& operator*=(Vector2<T>& v, const Matrix22<T>& m)
+      {
+      Vector2<T> vres(0);
+      vres[0] = v[0]*m[0] + v[1]*m[1];
+      vres[1] = v[0]*m[2] + v[1]*m[3];
+      v = vres;
+      return v;
+      }
+   template <typename T> inline Vector2<T> operator*(const Vector2<T>& v, const Matrix22<T>& m)
+      { return Vector2<T>(v) *= m; }
+
+   // Matrix * Vector
+   // Note: can't define Matrix *= Vector since the result is a Vector.
+   template <typename T> inline Vector2<T> operator*(const Matrix22<T>& m, const Vector2<T>& v)
+      { 
+      Vector2<T> vres(0);
+      vres[0] = m[0]*v[0] + m[2]*v[1];
+      vres[1] = m[1]*v[0] + m[3]*v[1];
+      return vres;
+      }
+
+
    } // namespace arda
 
 #endif
diff -r 5377337c7de0 -r 81d2aa42a860 include/Matrix.h
--- a/include/Matrix.h	Wed Sep 28 12:53:57 2011 -0500
+++ b/include/Matrix.h	Fri Oct 07 11:08:38 2011 -0500
@@ -2,520 +2,339 @@
 #define MATRIX_H_
 
 #include "Vector.h"
-using namespace arda::Vector;
 
 #include <cmath>
+#include <cstring>
 #include <string>
+#include <sstream>
 
 namespace arda 
    {
-   namespace Matrix
+   //////////////////////////////////////////////////////////////////////////
+   // Maxtrices
+   //
+   // Supported operations on matrices
+   //
+   // []
+   //     If you have a Matrix named m you can access it's member elements as
+   //     m[0], m[1], m[2], etc.  The indexing is COLUMN MAJOR.  So for a 2x2
+   //     matrix, m[0] and m[1] are the elements of the first column, not the
+   //     first row.
+   // == != + += - -=
+   //     Defined for operations on two matrices of the same type.
+   // * *= / /=
+   //     Defined for scalar multiplication/division with int, float, and double.
+   //     * is defined as a standalone template operator for the scalar * Matrix form.
+   // * *=
+   //     Also defined for Matrix * Matrix
+   //
+   // assign()          Assign the value of one matrix to another.
+   // det()             Calculates the determinant of a matrix.
+   // getstring()       Returns a printable string representation of the matrix.
+   // transpose()       Returns the transpose of a matrix.
+   // setidentity()     Sets a matrix to the identity matrix.
+
+   // inverse()         Returns the inverse of a non-singular matrix.
+   // getcol()          Returns the vector that is column n of the matrix.
+   // setcol()          Set the vector that is column n of the matrix.
+      
+   //////////////////////////////////////////////////////////////////////////
+   template <typename T> class Matrix22
       {
-      //////////////////////////////////////////////////////////////////////////
-      // 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);
+      public:
+      T m[4];
 
-      std::string & getstring4m(Matrix4i m, std::string & s);
-      std::string & getstring4m(Matrix4f m, std::string & s);
-      std::string & getstring4m(Matrix4d m, std::string & s);
+      // Constructors
+      Matrix22() {}
+      Matrix22(T a) { m[0] = m[1] = m[2] = m[3] = a; }
+      Matrix22(T a0, T a1, T a2, T a3) 
+	 { 
+	 m[0] = a0; m[2] = a2;
+	 m[1] = a1; m[3] = a3; 
+	 } 
 
-      //////////////////////////////////////////////////////////////////////////
-      // assign*()
+      // Array indexing
+      // Remember: column major order is used.
+      inline T& operator[](unsigned int i) 
+	 { assert (i<4); return m[i]; }
+      inline const T& operator[](unsigned int i) const
+	 { assert (i<4); return m[i]; }
 
-      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; 
+      // Assignment
+      inline Matrix22<T>& operator=(const Matrix22<T>& m2)
+	 { 
+	 m[0] = m2[0]; m[2] = m2[2]; 
+	 m[1] = m2[1]; m[3] = m2[3]; 
+	 return *this;
 	 }
-      inline float * assign2m(Matrix2f m1, Matrix2f m2)
+      inline Matrix22<T>& assign(T a0, T a1, T a2, T a3)
 	 { 
-	 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; 
+	 m[0] = a0; m[2] = a2;
+	 m[1] = a1; m[3] = a3; 
+	 return *this;
 	 }
 
-      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; 
+      // Comparison
+      inline bool operator==(Matrix22<T>& m2) const
+	 {
+	 return  
+	       m[0] == m2[0] && m[2] == m2[2] &&
+	       m[1] == m2[1] && m[3] == m2[3]; 
 	 }
-      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 bool operator!=(Matrix22<T>& m2) const
+	 { return ! (*this == m2); } 
+
+      // Matrix addition
+      inline Matrix22<T>& operator+=(const Matrix22<T>& m2)
+	 {
+	 m[0] += m2[0]; m[2] += m2[2]; 
+	 m[1] += m2[1]; m[3] += m2[3];
+	 return *this;
 	 }
-      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 Matrix22<T> operator+(const Matrix22<T>& m2) const
+	 { return Matrix22<T>(*this) += m2; }
+
+      // Matrix subtraction
+      inline Matrix22<T>& operator-=(const Matrix22<T>& m2)
+	 {
+	 m[0] -= m2[0]; m[2] -= m2[2]; 
+	 m[1] -= m2[1]; m[3] -= m2[3];
+	 return *this;
+	 }
+      inline Matrix22<T> operator-(const Matrix22<T>& m2) const
+	 { return Matrix22<T>(*this) -= m2; }
+      
+      // Scalar multiplication
+      inline Matrix22<T>& operator*=(const int a)
+	 {
+	 m[0] *= a; m[2] *= a;
+	 m[1] *= a; m[3] *= a;
+	 return *this;
+	 }
+      inline Matrix22<T>& operator*=(const float a)
+	 {
+	 m[0] *= a; m[2] *= a;
+	 m[1] *= a; m[3] *= a;
+	 return *this;
+	 }
+      inline Matrix22<T>& operator*=(const double a)
+	 {
+	 m[0] *= a; m[2] *= a;
+	 m[1] *= a; m[3] *= a;
+	 return *this;
 	 }
 
-      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; 
+      // Scalar division
+      inline Matrix22<T>& operator/=(const int a)
+	 {
+	 assert(a!=0);
+	 m[0] /= a; m[2] /= a;
+	 m[1] /= a; m[3] /= a;
+	 return *this;
 	 }
-      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 Matrix22<T>& operator/=(const float a)
+	 {
+	 assert(a!=0);
+	 m[0] /= a; m[2] /= a;
+	 m[1] /= a; m[3] /= a;
+	 return *this;
 	 }
-      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; 
+      inline Matrix22<T>& operator/=(const double a)
+	 {
+	 assert(a!=0);
+	 m[0] /= a; m[2] /= a;
+	 m[1] /= a; m[3] /= a;
+	 return *this;
+	 }
+      inline Matrix22<T> operator/(const int a)
+	 { return Matrix22<T>(*this) /= a; }
+      inline Matrix22<T> operator/(const float a)
+	 { return Matrix22<T>(*this) /= a; }
+      inline Matrix22<T> operator/(const double a)
+	 { return Matrix22<T>(*this) /= a; }
+
+      // Matrix multiplication
+      // Not sure if this should be inlined at all.  Sure the compiler will
+      // probably ignoe the inline request here, but maybe it won't and maybe
+      // that would be bad.  Needs real testing.
+      inline Matrix22<T>& operator*=(Matrix22<T> m2)
+	 {
+	 const int size=2;
+	 Matrix22<T> mres(0);
+	 int i, j, k;
+	 for (i=0; i<size; ++i)
+	    for (j=0; j<size; ++j)
+	       for (k=0; k<size; ++k)
+		  mres[size*i+j] += (*this)[size*k+j] * m2[size*i+k];
+	 *this = mres;
+	 return *this;
+	 }
+      inline Matrix22<T> operator*(Matrix22<T> m2) const
+	 { return Matrix22<T>(*this) *= m2; }
+
+      // methods
+      inline double det() const
+	 { return (double) (m[0] * m[3]) - (double) (m[1] * m[2]); }
+
+      std::string & getstring(std::string & s) const;
+
+      inline Matrix22<T>& setidentity()
+	 {
+	 memset(m, 0, sizeof(m));
+	 m[0] = m[3] = 1;
+	 return *this;
 	 }
 
-      //////////////////////////////////////////////////////////////////////////
-      // add*()
+      inline Matrix22<T> transpose() const
+	 {
+	 const int size=2;
+	 Matrix22<T> mres(0);
+	 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;
+	 }
+      };
 
-      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;
-	 }
+   // Scalar multiplication continued
+   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const int a)
+      { return Matrix22<T>(m) *= a; }
+   template <typename T> inline Matrix22<T> operator*(const int a, const Matrix22<T>& m)
+      { return Matrix22<T>(m) *= a; }
+   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const float a)
+      { return Matrix22<T>(m) *= a; }
+   template <typename T> inline Matrix22<T> operator*(const float a, const Matrix22<T>& m)
+      { return Matrix22<T>(m) *= a; }
+   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const double a)
+      { return Matrix22<T>(m) *= a; }
+   template <typename T> inline Matrix22<T> operator*(const double a, const Matrix22<T>& m)
+      { return Matrix22<T>(m) *= a; }
 
-      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;
-	 }
+   //////////////////////////////////////////////////////////////////////////
+   typedef Matrix22<int> Matrix22i;
+   typedef Matrix22<float> Matrix22f;
+   typedef Matrix22<double> Matrix22d;
 
-      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;
-	 }
+#ifdef NOTHING   
+   typedef Matrix33<int> Matrix33i;
+   typedef Matrix33<float> Matrix33f;
+   typedef Matrix33<double> Matrix33d;
+      
+   typedef Matrix44<int> Matrix44i;
+   typedef Matrix44<float> Matrix44f;
+   typedef Matrix44<double> Matrix44d;
+      
+#endif
+   } // namespace arda
 
-      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;
-	 }
+////////////////////////////////////////////////////////////////////////////////
+template <typename T> std::string & arda::Matrix22<T>::getstring(std::string & s) const
+   {
+   s.clear();
+   std::stringstream ss;
+   ss << "[ "
+	 "[ " << m[0] << ", " << m[1] << " ], "
+	 "[ " << m[2] << ", " << m[3] << " ] ]";
+   s = ss.str();
+   return s;   
+   }
 
-      //////////////////////////////////////////////////////////////////////////
-      // 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.
+#ifdef NOTHING
+   //////////////////////////////////////////////////////////////////////////
+   // det*(Matrix m)
+   //	Computer the determinant of a matrix.
+   //	Returns the determinant.
 
-      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);
-
-      //////////////////////////////////////////////////////////////////////////
-      // det*(Matrix m)
-      //	Computer the determinant of a matrix.
-      //	Returns the determinant.
-
-      inline double det2(Matrix2i m)
-	 { return m[0] * m[3] - m[1] * m[2]; }
-      inline double det2(Matrix2f m)
-	 { return m[0] * m[3] - m[1] * m[2]; }
-      inline double det2(Matrix2d m)
-	 { return m[0] * m[3] - m[1] * m[2]; }
-
-      inline double det3(Matrix3i m)
-	 { return m[0] * (  m[4] * m[8]  
-			  - m[5] * m[7])  
-	       -  m[3] * (  m[1] * m[8]
-			  - m[2] * m[7]) 
-	       +  m[6] * (  m[1] * m[5] 
-			  - m[2] * m[4]); }
-      inline double det3(Matrix3f m)
-	 { return m[0] * (  m[4] * m[8] 
-			  - m[5] * m[7])  
-	       -  m[3] * (  m[1] * m[8] 
-                          - m[2] * m[7]) 
-	       +  m[6] * (  m[1] * m[5] 
-			    - m[2] * m[4]); }
-      inline double det3(Matrix3d m)
-	 { return m[0] * (  m[4] * m[8] 
-			  - m[5] * m[7])
-	       - m[3] * (  m[1] * m[8] 
-			 - m[2] * m[7])
-	       + m[6] * (m[1] * m[5] 
+   inline double det3(Matrix3i m)
+      { return m[0] * (  m[4] * m[8]  
+			 - m[5] * m[7])  
+	    -  m[3] * (  m[1] * m[8]
+			 - m[2] * m[7]) 
+	    +  m[6] * (  m[1] * m[5] 
 			 - m[2] * m[4]); }
 
-      inline double det4(Matrix4i m)
-	 { return m[0] * (m[5] * (m[10] * m[15] - 
-				  m[11] * m[14]) - 
-			  m[9] * (m[6] * m[15] - 
-				  m[7] * m[14]) +
-			  m[13] * (m[6] * m[11] - 
-				   m[7] * m[10])) -
-	       m[4] * (m[1] * (m[10] * m[15] - 
+   inline double det4(Matrix4i m)
+      { return m[0] * (m[5] * (m[10] * m[15] - 
 			       m[11] * m[14]) - 
-		       m[9] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[11] - 
-				m[3] * m[10])) +
-	       m[8] * (m[1] * (m[6] * m[15] - 
-			       m[7] * m[14]) - 
-		       m[5] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[7] - 
-				m[3] * m[6])) +
-	       m[12] * (m[1] * (m[6] * m[11] - 
-				m[7] * m[10]) - 
-			m[5] * (m[2] * m[11] - 
-				m[3] * m[10]) +
-			m[9] * (m[2] * m[7] - 
-				m[3] * m[6])); }
-      inline double det4(Matrix4f m)
-	 { return m[0] * (m[5] * (m[10] * m[15] - 
-				  m[11] * m[14]) - 
-			  m[9] * (m[6] * m[15] - 
-				  m[7] * m[14]) +
-			  m[13] * (m[6] * m[11] - 
-				   m[7] * m[10])) -
-	       m[4] * (m[1] * (m[10] * m[15] - 
-			       m[11] * m[14]) - 
-		       m[9] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[11] - 
-				m[3] * m[10])) +
-	       m[8] * (m[1] * (m[6] * m[15] - 
-			       m[7] * m[14]) - 
-		       m[5] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[7] - 
-				m[3] * m[6])) +
-	       m[12] * (m[1] * (m[6] * m[11] - 
-				m[7] * m[10]) - 
-			m[5] * (m[2] * m[11] - 
-				m[3] * m[10]) +
-			m[9] * (m[2] * m[7] - 
-				m[3] * m[6])); }
-      inline double det4(Matrix4d m)
-	 { return m[0] * (m[5] * (m[10] * m[15] - 
-				  m[11] * m[14]) - 
-			  m[9] * (m[6] * m[15] - 
-				  m[7] * m[14]) +
-			  m[13] * (m[6] * m[11] - 
-				   m[7] * m[10])) -
-	       m[4] * (m[1] * (m[10] * m[15] - 
-			       m[11] * m[14]) - 
-		       m[9] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[11] - 
-				m[3] * m[10])) +
-	       m[8] * (m[1] * (m[6] * m[15] - 
-			       m[7] * m[14]) - 
-		       m[5] * (m[2] * m[15] - 
-			       m[3] * m[14]) +
-		       m[13] * (m[2] * m[7] - 
-				m[3] * m[6])) +
-	       m[12] * (m[1] * (m[6] * m[11] - 
-				m[7] * m[10]) - 
-			m[5] * (m[2] * m[11] - 
-				m[3] * m[10]) +
-			m[9] * (m[2] * m[7] - 
-				m[3] * m[6])); }
-      } // namespace Matrix
+		       m[9] * (m[6] * m[15] - 
+			       m[7] * m[14]) +
+		       m[13] * (m[6] * m[11] - 
+				m[7] * m[10])) -
+	    m[4] * (m[1] * (m[10] * m[15] - 
+			    m[11] * m[14]) - 
+		    m[9] * (m[2] * m[15] - 
+			    m[3] * m[14]) +
+		    m[13] * (m[2] * m[11] - 
+			     m[3] * m[10])) +
+	    m[8] * (m[1] * (m[6] * m[15] - 
+			    m[7] * m[14]) - 
+		    m[5] * (m[2] * m[15] - 
+			    m[3] * m[14]) +
+		    m[13] * (m[2] * m[7] - 
+			     m[3] * m[6])) +
+	    m[12] * (m[1] * (m[6] * m[11] - 
+			     m[7] * m[10]) - 
+		     m[5] * (m[2] * m[11] - 
+			     m[3] * m[10]) +
+		     m[9] * (m[2] * m[7] - 
+			     m[3] * m[6])); }
 
-   } // namespace arda
+////////////////////////////////////////////////////////////////////////////////
+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;
+   }
+
+
+////////////////////////////////////////////////////////////////////////////////
+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;
+   }
+
+////////////////////////////////////////////////////////////////////////////////
+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::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;   
+   }
+#endif
+
+
+
+
 #endif // MATRIX_H_
diff -r 5377337c7de0 -r 81d2aa42a860 include/Vector.h
--- a/include/Vector.h	Wed Sep 28 12:53:57 2011 -0500
+++ b/include/Vector.h	Fri Oct 07 11:08:38 2011 -0500
@@ -32,9 +32,7 @@
    //     Defined for operations on two vectors of the same type.
    // * *= / /=
    //     Defined for scalar multiplication/division with int, float, and double.
-   //     *= and /= are defined as member functions, but * and / are defined
-   //     as standalone template operators so that I can support v*a as well as
-   //     a*v order of multiplication.
+   //     * is defined as a standalone template operator for the scalar * Vector form.
    // assign(T a, T b [, T b [, T d]])
    //     Directly assign values to the vector.
    // cross(Vector3& v2)
@@ -83,15 +81,16 @@
    //
    // Scalar multiplication and division are overloaded.  Why?
    // I'm not sure I need this.  The idea is to support all likely scalar
-   // types without having to rely on implicit or explicity type conversion,
+   // types without having to rely on implicit or explicit type conversion,
    // especially not a conversion that might result in loss of precision.
-   // Multiplying a float Vector by a double scala, for example.
-   // Is there a better way to do this?
+   // Multiplying a float Vector by a double scalar, for example.
+   //
+   // Is there a better way to do this?  Am I worrying about nothing?
    //
    //
    //
    // I suppose I could make this into a single Vector template on two parameters
-   // template <class T, unsigned int N> class Vector
+   // template <typename T, unsigned int N> class Vector
    //    {
    //    T v[N];
    //
@@ -102,13 +101,13 @@
 
 
    /////////////////////////////////////////////////////////////////////////////
-   template <class T> class Vector2
+   template <typename T> class Vector2
       {
       public:
       T x, y;
 
       // Constructors
-      Vector2() : x (0), y(0) {}
+      Vector2() {}
       Vector2(T a) : x (a), y(a) {}
       Vector2(T a, T b) : x (a), y(b) {}
 
@@ -128,19 +127,19 @@
       inline bool operator==(Vector2<T>& v2) const
 	 { return ((x == v2.x) && (y == v2.y)); }
       inline bool operator!=(Vector2<T>& v2) const
-	 { return ! ((x == v2.x) && (y == v2.y)); }
+	 { return ! (*this == v2); }
 
       // Vector addition
-      inline Vector2<T> operator+(const Vector2<T>& v2) const
-	 { Vector2<T> vres; vres.x = x + v2.x; vres.y = y + v2.y; return vres; }
       inline Vector2<T>& operator+=(const Vector2<T>& v2)
 	 { x += v2.x; y += v2.y; return *this; }
+      inline Vector2<T> operator+(const Vector2<T>& v2) const
+	 { return Vector2<T>(*this) += v2; }
 
       // Vector subtraction
-      inline Vector2<T> operator-(const Vector2<T>& v2) const
-	 { Vector2<T> vres; vres.x = x - v2.x; vres.y = y - v2.y; return vres; }
       inline Vector2<T>& operator-=(const Vector2<T>& v2)
 	 { x -= v2.x; y -= v2.y; return *this; }
+      inline Vector2<T> operator-(const Vector2<T>& v2) const
+	 { return Vector2<T>(*this) -= v2; }
 
       // Scalar multiplication
       inline Vector2<T>& operator*=(const int a)
@@ -149,6 +148,12 @@
 	 { x *= a; y *= a; return *this; }
       inline Vector2<T>& operator*=(const double a)
 	 { x *= a; y *= a; return *this; }
+      inline Vector2<T> operator*(const int a)
+	 { return Vector2<T>(*this) *= a;}
+      inline Vector2<T> operator*(const float a)
+	 { return Vector2<T>(*this) *= a;}
+      inline Vector2<T> operator*(const double a)
+	 { return Vector2<T>(*this) *= a;}
 
       // Scalar division
       inline Vector2<T>& operator/=(const int a)
@@ -157,6 +162,12 @@
 	 { assert(a!=0); x /= a; y /= a; return *this; }
       inline Vector2<T>& operator/=(const double a)
 	 { assert(a!=0); x /= a; y /= a; return *this; }
+      inline Vector2<T> operator/(const int a)
+	 { return Vector2<T>(*this) /= a;}
+      inline Vector2<T> operator/(const float a)
+	 { return Vector2<T>(*this) /= a;}
+      inline Vector2<T> operator/(const double a)
+	 { return Vector2<T>(*this) /= a;}
 
 
       // methods
@@ -188,41 +199,22 @@
       };
 
    // Scalar multiplication, continued
-   template <class T> inline Vector2<T> operator*(const Vector2<T>& v, const int a)
-         { return Vector2<T>(v) *= a;}
-   template <class T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
-         { return Vector2<T>(v) *= a;}
-   template <class T> inline Vector2<T> operator*(const Vector2<T> v, const float a)
+   template <typename T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
       { return Vector2<T>(v) *= a;}
-   template <class T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
+   template <typename T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
       { return Vector2<T>(v) *= a;}
-   template <class T> inline Vector2<T> operator*(const Vector2<T>& v, const double a)
-      { return Vector2<T>(v) *= a;}
-   template <class T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
+   template <typename T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
       { return Vector2<T>(v) *= a;}
 
-   // Scalar division, continued
-   template <class T> inline Vector2<T> operator/(const Vector2<T>& v, const int a)
-      { return Vector2<T>(v) /= a;}
-   template <class T> inline Vector2<T> operator/(const int a, const Vector2<T>& v)
-      { return Vector2<T>(v) /= a;}
-   template <class T> inline Vector2<T> operator/(const Vector2<T>& v, const float a)
-      { return Vector2<T>(v) /= a;}
-   template <class T> inline Vector2<T> operator/(const float a, const Vector2<T>& v)
-      { return Vector2<T>(v) /= a;}
-   template <class T> inline Vector2<T> operator/(const Vector2<T>& v, const double a)
-      { return Vector2<T>(v) /= a;}
-   template <class T> inline Vector2<T> operator/(const double a, const Vector2<T>& v)
-      { return Vector2<T>(v) /= a;}
 
    /////////////////////////////////////////////////////////////////////////////
-   template <class T> class Vector3
+   template <typename T> class Vector3
       {
       public:
       T x, y, z;
 
       // Constructors
-      Vector3() : x (0), y(0), z(0) {}
+      Vector3() {}
       Vector3(T a) : x (a), y(a), z(a) {}
       Vector3(T a, T b, T c) : x (a), y(b), z(c) {}
 
@@ -242,19 +234,19 @@
       inline bool operator==(Vector3<T>& v2) const
 	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
       inline bool operator!=(Vector3<T>& v2) const
-	 { return ! ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
+	 { return ! (*this == v2); }
 
       // Vector addition
-      inline Vector3<T> operator+(const Vector3<T>& v2) const
-	 { Vector3<T> vres; vres.x = x + v2.x; vres.y = y + v2.y; vres.z = z + v2.z; return vres; }
       inline Vector3<T>& operator+=(const Vector3<T>& v2)
 	 { x += v2.x; y += v2.y; z += v2.z; return *this; }
+      inline Vector3<T> operator+(const Vector3<T>& v2) const
+	 { return Vector3<T>(*this) += v2; }
 
       // Vector subtraction
-      inline Vector3<T> operator-(const Vector3<T>& v2) const
-	 { Vector3<T> vres; vres.x = x - v2.x; vres.y = y - v2.y; vres.z = z - v2.z; return vres; }
       inline Vector3<T>& operator-=(const Vector3<T>& v2)
 	 { x -= v2.x; y -= v2.y; z -= v2.z; return *this; }
+      inline Vector3<T> operator-(const Vector3<T>& v2) const
+	 { return Vector3<T>(*this) -= v2; }
 
       // Scalar multiplication
       inline Vector3<T>& operator*=(const int a)
@@ -263,6 +255,12 @@
 	 { x *= a; y *= a; z *= a; return *this; }
       inline Vector3<T>& operator*=(const double a)
 	 { x *= a; y *= a; z *= a; return *this; }
+      inline Vector3<T> operator*(const int a)
+	 { return Vector3<T>(*this) *= a;}
+      inline Vector3<T> operator*(const float a)
+	 { return Vector3<T>(*this) *= a;}
+      inline Vector3<T> operator*(const double a)
+	 { return Vector3<T>(*this) *= a;}
 
       // Scalar division
       inline Vector3<T>& operator/=(const int a)
@@ -271,6 +269,12 @@
 	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
       inline Vector3<T>& operator/=(const double a)
 	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
+      inline Vector3<T> operator/(const int a)
+	 { return Vector3<T>(*this) /= a;}
+      inline Vector3<T> operator/(const float a)
+	 { return Vector3<T>(*this) /= a;}
+      inline Vector3<T> operator/(const double a)
+	 { return Vector3<T>(*this) /= a;}
 
 
       // methods
@@ -318,41 +322,21 @@
       };
 
    // Scalar multiplication, continued
-   template <class T> inline Vector3<T> operator*(const Vector3<T>& v, const int a)
-         { return Vector3<T>(v) *= a;}
-   template <class T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
-         { return Vector3<T>(v) *= a;}
-   template <class T> inline Vector3<T> operator*(const Vector3<T> v, const float a)
+   template <typename T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
       { return Vector3<T>(v) *= a;}
-   template <class T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
+   template <typename T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
       { return Vector3<T>(v) *= a;}
-   template <class T> inline Vector3<T> operator*(const Vector3<T>& v, const double a)
-      { return Vector3<T>(v) *= a;}
-   template <class T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
+   template <typename T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
       { return Vector3<T>(v) *= a;}
 
-   // Scalar division, continued
-   template <class T> inline Vector3<T> operator/(const Vector3<T>& v, const int a)
-      { return Vector3<T>(v) /= a;}
-   template <class T> inline Vector3<T> operator/(const int a, const Vector3<T>& v)
-      { return Vector3<T>(v) /= a;}
-   template <class T> inline Vector3<T> operator/(const Vector3<T>& v, const float a)
-      { return Vector3<T>(v) /= a;}
-   template <class T> inline Vector3<T> operator/(const float a, const Vector3<T>& v)
-      { return Vector3<T>(v) /= a;}
-   template <class T> inline Vector3<T> operator/(const Vector3<T>& v, const double a)
-      { return Vector3<T>(v) /= a;}
-   template <class T> inline Vector3<T> operator/(const double a, const Vector3<T>& v)
-      { return Vector3<T>(v) /= a;}
-
    /////////////////////////////////////////////////////////////////////////////
-   template <class T> class Vector4
+   template <typename T> class Vector4
       {
       public:
       T x, y, z, w;
 
       // Constructors
-      Vector4() : x (0), y(0), z(0), w(0) {}
+      Vector4() {}
       Vector4(T a) : x (a), y(a), z(a), w(a) {}
       Vector4(T a, T b, T c, T d) : x (a), y(b), z(c), w(d) {}
 
@@ -372,19 +356,19 @@
       inline bool operator==(Vector4<T>& v2) const
 	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
       inline bool operator!=(Vector4<T>& v2) const
-	 { return ! ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
+	 { return ! (*this == v2); }
 
       // Vector addition
-      inline Vector4<T> operator+(const Vector4<T>& v2) const
-	 { Vector4<T> vres; vres.x = x + v2.x; vres.y = y + v2.y; vres.z = z + v2.z; vres.w = w + v2.w; return vres; }
       inline Vector4<T>& operator+=(const Vector4<T>& v2)
 	 { x += v2.x; y += v2.y; z += v2.z; w += v2.w; return *this; }
+      inline Vector4<T> operator+(const Vector4<T>& v2) const
+	 { return Vector4<T>(*this) += v2; }
 
       // Vector subtraction
-      inline Vector4<T> operator-(const Vector4<T>& v2) const
-	 { Vector4<T> vres; vres.x = x - v2.x; vres.y = y - v2.y; vres.z = z - v2.z; vres.w = w - v2.w; return vres; }
       inline Vector4<T>& operator-=(const Vector4<T>& v2)
 	 { x -= v2.x; y -= v2.y; z -= v2.z; w -= v2.w; return *this; }
+      inline Vector4<T> operator-(const Vector4<T>& v2) const
+	 { return Vector4<T>(*this) -= v2; }
 
       // Scalar multiplication
       inline Vector4<T>& operator*=(const int a)
@@ -393,6 +377,12 @@
 	 { x *= a; y *= a; z *= a; w *= a; return *this; }
       inline Vector4<T>& operator*=(const double a)
 	 { x *= a; y *= a; z *= a; w *= a; return *this; }
+      inline Vector4<T> operator*(const int a)
+	 { return Vector4<T>(*this) *= a;}
+      inline Vector4<T> operator*(const float a)
+	 { return Vector4<T>(*this) *= a;}
+      inline Vector4<T> operator*(const double a)
+	 { return Vector4<T>(*this) *= a;}
 
       // Scalar division
       inline Vector4<T>& operator/=(const int a)
@@ -401,6 +391,12 @@
 	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
       inline Vector4<T>& operator/=(const double a)
 	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
+      inline Vector4<T> operator/(const int a)
+	 { return Vector4<T>(*this) /= a;}
+      inline Vector4<T> operator/(const float a)
+	 { return Vector4<T>(*this) /= a;}
+      inline Vector4<T> operator/(const double a)
+	 { return Vector4<T>(*this) /= a;}
 
 
       // methods
@@ -432,32 +428,12 @@
       };
 
    // Scalar multiplication, continued
-   template <class T> inline Vector4<T> operator*(const Vector4<T>& v, const int a)
-         { return Vector4<T>(v) *= a;}
-   template <class T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
-         { return Vector4<T>(v) *= a;}
-   template <class T> inline Vector4<T> operator*(const Vector4<T> v, const float a)
+   template <typename T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
       { return Vector4<T>(v) *= a;}
-   template <class T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
+   template <typename T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
       { return Vector4<T>(v) *= a;}
-   template <class T> inline Vector4<T> operator*(const Vector4<T>& v, const double a)
+   template <typename T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
       { return Vector4<T>(v) *= a;}
-   template <class T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
-      { return Vector4<T>(v) *= a;}
-
-   // Scalar division, continued
-   template <class T> inline Vector4<T> operator/(const Vector4<T>& v, const int a)
-      { return Vector4<T>(v) /= a;}
-   template <class T> inline Vector4<T> operator/(const int a, const Vector4<T>& v)
-      { return Vector4<T>(v) /= a;}
-   template <class T> inline Vector4<T> operator/(const Vector4<T>& v, const float a)
-      { return Vector4<T>(v) /= a;}
-   template <class T> inline Vector4<T> operator/(const float a, const Vector4<T>& v)
-      { return Vector4<T>(v) /= a;}
-   template <class T> inline Vector4<T> operator/(const Vector4<T>& v, const double a)
-      { return Vector4<T>(v) /= a;}
-   template <class T> inline Vector4<T> operator/(const double a, const Vector4<T>& v)
-      { return Vector4<T>(v) /= a;}
 
    /////////////////////////////////////////////////////////////////////////////
 
@@ -476,7 +452,7 @@
    } // namespace arda
 
 ////////////////////////////////////////////////////////////////////////////////
-template <class T> std::string & arda::Vector2<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Vector2<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
@@ -486,7 +462,7 @@
    }
 
 ////////////////////////////////////////////////////////////////////////////////
-template <class T> std::string & arda::Vector3<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Vector3<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
@@ -496,7 +472,7 @@
    }
 
 ////////////////////////////////////////////////////////////////////////////////
-template <class T> std::string & arda::Vector4<T>::getstring(std::string & s) const
+template <typename T> std::string & arda::Vector4<T>::getstring(std::string & s) const
    {
    s.clear();
    std::stringstream ss;
diff -r 5377337c7de0 -r 81d2aa42a860 src/main.cpp
--- a/src/main.cpp	Wed Sep 28 12:53:57 2011 -0500
+++ b/src/main.cpp	Fri Oct 07 11:08:38 2011 -0500
@@ -2,6 +2,7 @@
 #include <iomanip>
 #include <limits>
 #include <string>
+
 using namespace std;
 
 #include "Math.h"
@@ -15,7 +16,7 @@
    {
    cout << "===============================================================" << endl <<
 	 "Testing Vector constructors" << endl;
-   Vector2i v2i_1;
+   Vector2i v2i_1(0);
    Vector2i v2i_2(1);
    Vector2i v2i_3(1, 2);
    Vector2i v2i_4(v2i_3);
@@ -25,7 +26,7 @@
    cout << setw(40) << "v2i_3: " << v2i_3.getstring(s1) << endl;
    cout << setw(40) << "v2i_4: " << v2i_4.getstring(s1) << endl;
 
-   Vector3f v3f_1;
+   Vector3f v3f_1(0);
    Vector3f v3f_2(1);
    Vector3f v3f_3(1, 2, 3);
    Vector3f v3f_4(v3f_3);
@@ -35,7 +36,7 @@
    cout << setw(40) << "v3f_3: " << v3f_3.getstring(s1) << endl;
    cout << setw(40) << "v3f_4: " << v3f_4.getstring(s1) << endl;
 
-   Vector4d v4d_1;
+   Vector4d v4d_1(0);
    Vector4d v4d_2(1);
    Vector4d v4d_3(1, 2, 3, 4);
    Vector4d v4d_4(v4d_3);
@@ -64,9 +65,9 @@
    Vector2i v2i(1,2);
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector4d v4d(1.1, 2.2, 3.3, 4.4);
-   Vector2i v2i_2;
-   Vector3f v3f_2;
-   Vector4d v4d_2;
+   Vector2i v2i_2(0);
+   Vector3f v3f_2(0);
+   Vector4d v4d_2(0);
 
    cout << "Before assignment" << endl;
    cout << setw(40) << "v2i: " << v2i.getstring(s1) << endl;
@@ -99,9 +100,12 @@
    Vector2i v2i(1,2);
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector4d v4d(1.1, 2.2, 3.3, 4.4);
-   Vector2i v2i_2;
-   Vector3f v3f_2;
-   Vector4d v4d_2;
+   Vector2i v2i_2(1,2);
+   Vector3f v3f_2(1.1f, 2.2f, 3.3f);
+   Vector4d v4d_2(1.1, 2.2, 3.3, 4.4);
+   Vector2i v2i_3(0);
+   Vector3f v3f_3(0);
+   Vector4d v4d_3(0);
 
    cout << setw(40) << "v2i: " << v2i.getstring(s1) << endl;
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
@@ -109,13 +113,22 @@
    cout << setw(40) << "v2i_2: " << v2i_2.getstring(s1) << endl;
    cout << setw(40) << "v3f_2: " << v3f_2.getstring(s1) << endl;
    cout << setw(40) << "v4d_2: " << v4d_2.getstring(s1) << endl;
+   cout << setw(40) << "v2i_3: " << v2i_3.getstring(s1) << endl;
+   cout << setw(40) << "v3f_3: " << v3f_3.getstring(s1) << endl;
+   cout << setw(40) << "v4d_3: " << v4d_3.getstring(s1) << endl;
    cout << boolalpha;
    cout << setw(40) << "v2i == v2i_2: " << (v2i == v2i_2) << endl;
+   cout << setw(40) << "v2i == v2i_3: " << (v2i == v2i_3) << endl;
+   cout << setw(40) << "v2i != v2i_2: " << (v2i != v2i_2) << endl;
+   cout << setw(40) << "v2i != v2i_3: " << (v2i != v2i_3) << endl;
    cout << setw(40) << "v3f == v3f_2: " << (v3f == v3f_2) << endl;
+   cout << setw(40) << "v3f == v3f_3: " << (v3f == v3f_3) << endl;
+   cout << setw(40) << "v3f != v3f_2: " << (v3f != v3f_2) << endl;
+   cout << setw(40) << "v3f != v3f_3: " << (v3f != v3f_3) << endl;
    cout << setw(40) << "v4d == v4d_2: " << (v4d == v4d_2) << endl;
-   cout << setw(40) << "v2i != v2i_2: " << (v2i != v2i_2) << endl;
-   cout << setw(40) << "v3f != v3f_2: " << (v3f != v3f_2) << endl;
+   cout << setw(40) << "v4d == v4d_3: " << (v4d == v4d_3) << endl;
    cout << setw(40) << "v4d != v4d_2: " << (v4d != v4d_2) << endl;
+   cout << setw(40) << "v4d != v4d_3: " << (v4d != v4d_3) << endl;
    }
 
    {
@@ -127,9 +140,9 @@
    Vector2i v2i_2(3,4);
    Vector3f v3f_2(4.4, 5.5, 6.6);
    Vector4d v4d_2(5.5, 6.6, 7.7, 8.8);
-   Vector2i v2i_3;
-   Vector3f v3f_3 ;
-   Vector4d v4d_3;
+   Vector2i v2i_3(0);
+   Vector3f v3f_3(0);
+   Vector4d v4d_3(0);
 
    cout << setw(40) << "v2i: " << v2i.getstring(s1) << endl;
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
@@ -160,9 +173,9 @@
    Vector2i v2i_2(3,4);
    Vector3f v3f_2(4.4, 5.5, 6.6);
    Vector4d v4d_2(5.5, 6.6, 7.7, 8.8);
-   Vector2i v2i_3;
-   Vector3f v3f_3;
-   Vector4d v4d_3;
+   Vector2i v2i_3(0);
+   Vector3f v3f_3(0);
+   Vector4d v4d_3(0);
 
    cout << setw(40) << "v2i: " << v2i.getstring(s1) << endl;
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
@@ -190,13 +203,16 @@
    Vector2i v2i(1,2);
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector4d v4d(1.1, 2.2, 3.3, 4.4);
-   Vector2i v2i_2;
-   Vector3f v3f_2;
-   Vector4d v4d_2;
+   Vector2i v2i_2(0);
+   Vector3f v3f_2(0);
+   Vector4d v4d_2(0);
    int i = 2;
    float f = 2.f;
    double d = 2.0;
 
+   cout << setw(40) << "i: " << i << endl;
+   cout << setw(40) << "f: " << f << endl;
+   cout << setw(40) << "d: " << d << endl;
    cout << setw(40) << "v2i: " << v2i.getstring(s1) << endl;
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
    cout << setw(40) << "v4d: " << v4d.getstring(s1) << endl;
@@ -223,9 +239,9 @@
    Vector2i v2i(1,2);
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector4d v4d(1.1, 2.2, 3.3, 4.4);
-   Vector2i v2i_2;
-   Vector3f v3f_2;
-   Vector4d v4d_2;
+   Vector2i v2i_2(0);
+   Vector3f v3f_2(0);
+   Vector4d v4d_2(0);
    int i = 2;
    float f = 2.f;
    double d = 2.0;
@@ -255,8 +271,8 @@
 	 "Testing Vector cross product" << endl;
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector3f v3f_2(4.4, 5.5, 6.6);
-   Vector3f v3f_3;
-   Vector3f v3f_4;
+   Vector3f v3f_3(0);
+   Vector3f v3f_4(0);
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
    cout << setw(40) << "v3f_2: " << v3f_2.getstring(s1) << endl;
    cout << setw(40) << "v3f_3: " << v3f_3.getstring(s1) << endl;
@@ -388,13 +404,13 @@
 	 "Testing Vector get_proj" << endl;
    Vector2f v2f(1.1, 2.2);
    Vector2f v2f_2(3.3, 4.4);
-   Vector2f v2f_3;
+   Vector2f v2f_3(0);
    Vector3f v3f(1.1f, 2.2f, 3.3f);
    Vector3f v3f_2(4.4, 5.5, 6.6);
-   Vector3f v3f_3;
+   Vector3f v3f_3(0);
    Vector4d v4d(1.1, 2.2, 3.3, 4.4);
    Vector4d v4d_2(5.5, 6.6, 7.7, 8.8);
-   Vector4d v4d_3;
+   Vector4d v4d_3(0);
    cout << setw(40) << "v2f: " << v2f.getstring(s1) << endl;
    cout << setw(40) << "v2f_2: " << v2f_2.getstring(s1) << endl;
    cout << setw(40) << "v3f: " << v3f.getstring(s1) << endl;
@@ -423,7 +439,179 @@
 ////////////////////////////////////////////////////////////////////////////////
 void test_matrix(void)
    {
-   string s;
+   string s1, s2;
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix constructors" << endl;
+   Matrix22i m22i_1(0);
+   Matrix22i m22i_2(1);
+   Matrix22i m22i_3(1, 2, 3, 4);
+   Matrix22i m22i_4(m22i_3);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   cout << setw(40) << "m22i_3: " << m22i_3.getstring(s1) << endl;
+   cout << setw(40) << "m22i_4: " << m22i_4.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix array indexing" << endl;
+   Matrix22i m22i(1, 2, 3, 4);
+
+   cout << setw(40) << "m22i: " << m22i[0] << " "  << m22i[1] << " " << m22i[2] << " " << m22i[3] << " "<< endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix array assignment" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(5, 6, 7, 8);
+
+   cout << "Before assignment" << endl;
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+
+   m22i_2 = m22i_1;
+      cout << "After assignment by = " << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+
+   m22i_2.assign(9, 10, 11, 12);
+   cout << "After assignment by assign" << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix comparison" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(1, 2, 3, 4);
+   Matrix22i m22i_3(0);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   cout << setw(40) << "m22i_3: " << m22i_3.getstring(s1) << endl;
+   cout << boolalpha;
+   cout << setw(40) << "m22i_1 == m22i_2: " << (m22i_1 == m22i_2) << endl;
+   cout << setw(40) << "m22i_1 == m22i_3: " << (m22i_2 == m22i_3) << endl;
+   cout << setw(40) << "m22i_1 != m22i_2: " << (m22i_1 != m22i_2) << endl;
+   cout << setw(40) << "m22i_1 != m22i_3: " << (m22i_2 != m22i_3) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix addition" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(5, 6, 7, 8);
+   Matrix22i m22i_3(0);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   cout << setw(40) << "m22i_3: " << m22i_3.getstring(s1) << endl;
+   m22i_3 = m22i_1 + m22i_2;
+   cout << setw(40) << "m22i_3 = m22i_1 + m22i_2: " << m22i_3.getstring(s1) << endl;
+
+   m22i_3 += m22i_1;
+   cout << setw(40) << "m22i_3 += m22i_1: " << m22i_3.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix subtraction" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(5, 6, 7, 8);
+   Matrix22i m22i_3(0);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   cout << setw(40) << "m22i_3: " << m22i_3.getstring(s1) << endl;
+   m22i_3 = m22i_1 - m22i_2;
+   cout << setw(40) << "m22i_3 = m22i_1 - m22i_2: " << m22i_3.getstring(s1) << endl;
+
+   m22i_3 -= m22i_1;
+   cout << setw(40) << "m22i_3 -= m22i_1: " << m22i_3.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix scalar multiplication" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(0);
+   int i = 2;
+
+   cout << setw(40) << "i: " << i << endl;
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   m22i_2 = m22i_1 *i;
+   cout << setw(40) << "m22i_2 = m22i_1 * i: " << m22i_2.getstring(s1) << endl;
+   m22i_2 = i * m22i_1;
+   cout << setw(40) << "m22i_2 = i * m22i_1: " << m22i_2.getstring(s1) << endl;
+   m22i_2 *= i;
+   cout << setw(40) << "m22i_2 *= i: " << m22i_2.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Matrix scalar division" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(0);
+   int i = 2;
+
+   cout << setw(40) << "i: " << i << endl;
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   m22i_2 = m22i_1 / i;
+   cout << setw(40) << "m22i_2 = m22i_1 / i: " << m22i_2.getstring(s1) << endl;
+   m22i_1 /= i;
+   cout << setw(40) << "m22i_1 /= i: " << m22i_2.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix multiplication" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2(5, 6, 7, 8);
+   Matrix22i m22i_3(0);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   cout << setw(40) << "m22i_2: " << m22i_2.getstring(s1) << endl;
+   cout << setw(40) << "m22i_3: " << m22i_3.getstring(s1) << endl;
+   m22i_3 = m22i_1 * m22i_2;
+   cout << setw(40) << "m22i_3 = m22i_1 * m22i_2: " << m22i_3.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix det" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   double d;
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   d = m22i_1.det();
+   cout << setw(40) << "d = m22i_1.det(): " << d << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix transpose" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+   Matrix22i m22i_2;
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   m22i_2 = m22i_1.transpose();
+   cout << setw(40) << "m22i_2 = m22i_1.transpose(): " << m22i_2.getstring(s1) << endl;
+   }
+
+   {
+   cout << "===============================================================" << endl <<
+	 "Testing Maxtrix setidentity" << endl;
+   Matrix22i m22i_1(1, 2, 3, 4);
+
+   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
+   m22i_1.setidentity();
+   cout << setw(40) << "m22i_1.setidentity(): " << m22i_1.getstring(s1) << endl;
+   }
    }
 
 ////////////////////////////////////////////////////////////////////////////////