changeset 11:398894c52b69

Changed mind. Sticking with current design. Matrix22<T> is done.
author Eris Caffee <discordia@eldalin.com>
date Fri, 07 Oct 2011 13:33:48 -0500
parents 139d189290e0
children 1d2b25d4517f
files include/Math.h include/Matrix.h include/Vector.h src/Matrix.cpp src/main.cpp
diffstat 5 files changed, 733 insertions(+), 1122 deletions(-) [+]
line diff
     1.1 --- a/include/Math.h	Fri Oct 07 11:09:15 2011 -0500
     1.2 +++ b/include/Math.h	Fri Oct 07 13:33:48 2011 -0500
     1.3 @@ -22,41 +22,42 @@
     1.4     ////////////////////////////////////////////////////////////////////////////////
     1.5     // Angles
     1.6  
     1.7 -   namespace Angles
     1.8 +   namespace Math
     1.9        {
    1.10        // This many decimals is enough to cover even IEEE quad precision (128 bit)
    1.11        // reals.  Run on PA-RISC lately?  :)
    1.12        const double PI = 3.14159265358979323846264338327950288;
    1.13 -      }
    1.14  
    1.15 -   ////////////////////////////////////////////////////////////////////////////////
    1.16 -   // Mixed Vector / Matrix operations.
    1.17 -   // * *=
    1.18 -   //     Defined for Vector * Matrix and Matrix * Vector, but these
    1.19 -   //     operators are not part of the class.  They are defined independently.
    1.20 +      //////////////////////////////////////////////////////////////////////////
    1.21 +      // Mixed Vector / Matrix operations.
    1.22 +      // * *=
    1.23 +      //     Defined for Vector * Matrix and Matrix * Vector, but these
    1.24 +      //     operators are not part of the class.  They are defined independently.
    1.25 +      //
    1.26  
    1.27 -   // Vector * Matrix
    1.28 -   template <typename T> inline Vector2<T>& operator*=(Vector2<T>& v, const Matrix22<T>& m)
    1.29 -      {
    1.30 -      Vector2<T> vres(0);
    1.31 -      vres[0] = v[0]*m[0] + v[1]*m[1];
    1.32 -      vres[1] = v[0]*m[2] + v[1]*m[3];
    1.33 -      v = vres;
    1.34 -      return v;
    1.35 -      }
    1.36 -   template <typename T> inline Vector2<T> operator*(const Vector2<T>& v, const Matrix22<T>& m)
    1.37 -      { return Vector2<T>(v) *= m; }
    1.38 +      // Vector * Matrix
    1.39 +      template <typename T> inline Vector2<T>& operator*=(Vector2<T>& v, const Matrix22<T>& m)
    1.40 +	 {
    1.41 +	 Vector2<T> vres(0);
    1.42 +	 vres[0] = v[0]*m[0] + v[1]*m[1];
    1.43 +	 vres[1] = v[0]*m[2] + v[1]*m[3];
    1.44 +	 v = vres;
    1.45 +	 return v;
    1.46 +	 }
    1.47 +      template <typename T> inline Vector2<T> operator*(const Vector2<T>& v, const Matrix22<T>& m)
    1.48 +	 { return Vector2<T>(v) *= m; }
    1.49  
    1.50 -   // Matrix * Vector
    1.51 -   // Note: can't define Matrix *= Vector since the result is a Vector.
    1.52 -   template <typename T> inline Vector2<T> operator*(const Matrix22<T>& m, const Vector2<T>& v)
    1.53 -      { 
    1.54 -      Vector2<T> vres(0);
    1.55 -      vres[0] = m[0]*v[0] + m[2]*v[1];
    1.56 -      vres[1] = m[1]*v[0] + m[3]*v[1];
    1.57 -      return vres;
    1.58 -      }
    1.59 +      // Matrix * Vector
    1.60 +      // Note: can't define Matrix *= Vector since the result is a Vector.
    1.61 +      template <typename T> inline Vector2<T> operator*(const Matrix22<T>& m, const Vector2<T>& v)
    1.62 +	 { 
    1.63 +	 Vector2<T> vres(0);
    1.64 +	 vres[0] = m[0]*v[0] + m[2]*v[1];
    1.65 +	 vres[1] = m[1]*v[0] + m[3]*v[1];
    1.66 +	 return vres;
    1.67 +	 }
    1.68  
    1.69 +      } // namespace Math
    1.70  
    1.71     } // namespace arda
    1.72  
     2.1 --- a/include/Matrix.h	Fri Oct 07 11:09:15 2011 -0500
     2.2 +++ b/include/Matrix.h	Fri Oct 07 13:33:48 2011 -0500
     2.3 @@ -10,182 +10,283 @@
     2.4  
     2.5  namespace arda 
     2.6     {
     2.7 -   //////////////////////////////////////////////////////////////////////////
     2.8 -   // Maxtrices
     2.9 -   //
    2.10 -   // Supported operations on matrices
    2.11 -   //
    2.12 -   // []
    2.13 -   //     If you have a Matrix named m you can access it's member elements as
    2.14 -   //     m[0], m[1], m[2], etc.  The indexing is COLUMN MAJOR.  So for a 2x2
    2.15 -   //     matrix, m[0] and m[1] are the elements of the first column, not the
    2.16 -   //     first row.
    2.17 -   // == != + += - -=
    2.18 -   //     Defined for operations on two matrices of the same type.
    2.19 -   // * *= / /=
    2.20 -   //     Defined for scalar multiplication/division with int, float, and double.
    2.21 -   //     * is defined as a standalone template operator for the scalar * Matrix form.
    2.22 -   // * *=
    2.23 -   //     Also defined for Matrix * Matrix
    2.24 -   //
    2.25 -   // assign()          Assign the value of one matrix to another.
    2.26 -   // det()             Calculates the determinant of a matrix.
    2.27 -   // getstring()       Returns a printable string representation of the matrix.
    2.28 -   // transpose()       Returns the transpose of a matrix.
    2.29 -   // setidentity()     Sets a matrix to the identity matrix.
    2.30 +   namespace Math
    2.31 +      {
    2.32 +      //////////////////////////////////////////////////////////////////////////
    2.33 +      // Maxtrices
    2.34 +      //
    2.35 +      // Supported operations on matrices
    2.36 +      //
    2.37 +      // []
    2.38 +      //     If you have a Matrix named m you can access it's member elements as
    2.39 +      //     m[0], m[1], m[2], etc.  The indexing is COLUMN MAJOR.  So for a 2x2
    2.40 +      //     matrix, m[0] and m[1] are the elements of the first column, not the
    2.41 +      //     first row.
    2.42 +      // == != + += - -=
    2.43 +      //     Defined for operations on two matrices of the same type.
    2.44 +      // * *= / /=
    2.45 +      //     Defined for scalar multiplication/division with int, float, and double.
    2.46 +      //     * is defined as a standalone template operator for the scalar * Matrix
    2.47 +      //     form.
    2.48 +      // * *=
    2.49 +      //     Also defined for Matrix * Matrix
    2.50 +      //
    2.51 +      // assign()          Assign the value of one matrix to another.
    2.52 +      // getcol()          Returns the vector that is column n of the matrix.
    2.53 +      // getstring()       Returns a printable string representation of the matrix.
    2.54 +      // setcol()          Set the vector that is column n of the matrix.
    2.55 +      //                   Overloaded so that you can set from a vector or from
    2.56 +      //                   specific values.
    2.57 +      // setidentity()     Sets a matrix to the identity matrix.
    2.58 +      //
    2.59 +      //
    2.60 +      //
    2.61 +      // The following are functions, not class methods.  This is because I think
    2.62 +      // that the notation det(m) and transpose(m), for example, is more natural
    2.63 +      // to write and read than m.det() and m.transpose().  Using functional
    2.64 +      // notation just more naturally mimics the standard mathematical notation.
    2.65 +      //
    2.66 +      // det()             Calculates the determinant of a matrix.
    2.67 +      // transpose()       Returns the transpose of a matrix.
    2.68 +      
    2.69 +      //////////////////////////////////////////////////////////////////////////
    2.70 +      template <typename T> class Matrix22
    2.71 +	 {
    2.72 +	 public:
    2.73 +	 T m[4];
    2.74  
    2.75 -   // inverse()         Returns the inverse of a non-singular matrix.
    2.76 -   // getcol()          Returns the vector that is column n of the matrix.
    2.77 -   // setcol()          Set the vector that is column n of the matrix.
    2.78 +	 // Constructors
    2.79 +	 Matrix22() {}
    2.80 +	 explicit Matrix22(T a) { m[0] = m[1] = m[2] = m[3] = a; }
    2.81 +	 Matrix22(T a0, T a1, T a2, T a3) 
    2.82 +	    { 
    2.83 +	    m[0] = a0; m[2] = a2;
    2.84 +	    m[1] = a1; m[3] = a3; 
    2.85 +	    } 
    2.86 +
    2.87 +	 // Array indexing
    2.88 +	 // Remember: column major order is used.
    2.89 +	 inline T& operator[](unsigned int i) 
    2.90 +	    { assert (i<4); return m[i]; }
    2.91 +	 inline T operator[](unsigned int i) const
    2.92 +	    { assert (i<4); return m[i]; }
    2.93 +
    2.94 +	 // Assignment
    2.95 +	 inline Matrix22<T>& operator=(const Matrix22<T>& m2)
    2.96 +	    { 
    2.97 +	    m[0] = m2[0]; m[2] = m2[2]; 
    2.98 +	    m[1] = m2[1]; m[3] = m2[3]; 
    2.99 +	    return *this;
   2.100 +	    }
   2.101 +	 inline Matrix22<T>& assign(T a0, T a1, T a2, T a3)
   2.102 +	    { 
   2.103 +	    m[0] = a0; m[2] = a2;
   2.104 +	    m[1] = a1; m[3] = a3; 
   2.105 +	    return *this;
   2.106 +	    }
   2.107 +
   2.108 +	 // Comparison
   2.109 +	 inline bool operator==(Matrix22<T>& m2) const
   2.110 +	    {
   2.111 +	    return  
   2.112 +		  m[0] == m2[0] && m[2] == m2[2] &&
   2.113 +		  m[1] == m2[1] && m[3] == m2[3]; 
   2.114 +	    }
   2.115 +	 inline bool operator!=(Matrix22<T>& m2) const
   2.116 +	    { return ! (*this == m2); } 
   2.117 +
   2.118 +	 // Matrix addition
   2.119 +	 inline Matrix22<T>& operator+=(const Matrix22<T>& m2)
   2.120 +	    {
   2.121 +	    m[0] += m2[0]; m[2] += m2[2]; 
   2.122 +	    m[1] += m2[1]; m[3] += m2[3];
   2.123 +	    return *this;
   2.124 +	    }
   2.125 +	 inline Matrix22<T> operator+(const Matrix22<T>& m2) const
   2.126 +	    { return Matrix22<T>(*this) += m2; }
   2.127 +
   2.128 +	 // Matrix subtraction
   2.129 +	 inline Matrix22<T>& operator-=(const Matrix22<T>& m2)
   2.130 +	    {
   2.131 +	    m[0] -= m2[0]; m[2] -= m2[2]; 
   2.132 +	    m[1] -= m2[1]; m[3] -= m2[3];
   2.133 +	    return *this;
   2.134 +	    }
   2.135 +	 inline Matrix22<T> operator-(const Matrix22<T>& m2) const
   2.136 +	    { return Matrix22<T>(*this) -= m2; }
   2.137        
   2.138 -   //////////////////////////////////////////////////////////////////////////
   2.139 -   template <typename T> class Matrix22
   2.140 -      {
   2.141 -      public:
   2.142 -      T m[4];
   2.143 +	 // Scalar multiplication
   2.144 +	 inline Matrix22<T>& operator*=(const int a)
   2.145 +	    {
   2.146 +	    m[0] *= a; m[2] *= a;
   2.147 +	    m[1] *= a; m[3] *= a;
   2.148 +	    return *this;
   2.149 +	    }
   2.150 +	 inline Matrix22<T>& operator*=(const float a)
   2.151 +	    {
   2.152 +	    m[0] *= a; m[2] *= a;
   2.153 +	    m[1] *= a; m[3] *= a;
   2.154 +	    return *this;
   2.155 +	    }
   2.156 +	 inline Matrix22<T>& operator*=(const double a)
   2.157 +	    {
   2.158 +	    m[0] *= a; m[2] *= a;
   2.159 +	    m[1] *= a; m[3] *= a;
   2.160 +	    return *this;
   2.161 +	    }
   2.162  
   2.163 -      // Constructors
   2.164 -      Matrix22() {}
   2.165 -      Matrix22(T a) { m[0] = m[1] = m[2] = m[3] = a; }
   2.166 -      Matrix22(T a0, T a1, T a2, T a3) 
   2.167 -	 { 
   2.168 -	 m[0] = a0; m[2] = a2;
   2.169 -	 m[1] = a1; m[3] = a3; 
   2.170 -	 } 
   2.171 +	 // Scalar division
   2.172 +	 inline Matrix22<T>& operator/=(const int a)
   2.173 +	    {
   2.174 +	    assert(a!=0);
   2.175 +	    m[0] /= a; m[2] /= a;
   2.176 +	    m[1] /= a; m[3] /= a;
   2.177 +	    return *this;
   2.178 +	    }
   2.179 +	 inline Matrix22<T>& operator/=(const float a)
   2.180 +	    {
   2.181 +	    assert(a!=0);
   2.182 +	    m[0] /= a; m[2] /= a;
   2.183 +	    m[1] /= a; m[3] /= a;
   2.184 +	    return *this;
   2.185 +	    }
   2.186 +	 inline Matrix22<T>& operator/=(const double a)
   2.187 +	    {
   2.188 +	    assert(a!=0);
   2.189 +	    m[0] /= a; m[2] /= a;
   2.190 +	    m[1] /= a; m[3] /= a;
   2.191 +	    return *this;
   2.192 +	    }
   2.193 +	 inline Matrix22<T> operator/(const int a)
   2.194 +	    { return Matrix22<T>(*this) /= a; }
   2.195 +	 inline Matrix22<T> operator/(const float a)
   2.196 +	    { return Matrix22<T>(*this) /= a; }
   2.197 +	 inline Matrix22<T> operator/(const double a)
   2.198 +	    { return Matrix22<T>(*this) /= a; }
   2.199  
   2.200 -      // Array indexing
   2.201 -      // Remember: column major order is used.
   2.202 -      inline T& operator[](unsigned int i) 
   2.203 -	 { assert (i<4); return m[i]; }
   2.204 -      inline const T& operator[](unsigned int i) const
   2.205 -	 { assert (i<4); return m[i]; }
   2.206 +	 // Matrix multiplication
   2.207 +	 // Not sure if this should be inlined at all.  Sure the compiler will
   2.208 +	 // probably ignore the inline request here, but maybe it won't and maybe
   2.209 +	 // that would be bad.  Needs real testing.
   2.210 +	 inline Matrix22<T>& operator*=(Matrix22<T> m2)
   2.211 +	    {
   2.212 +	    const int size=2;
   2.213 +	    Matrix22<T> mres(0);
   2.214 +	    int i, j, k;
   2.215 +	    for (i=0; i<size; ++i)
   2.216 +	       for (j=0; j<size; ++j)
   2.217 +		  for (k=0; k<size; ++k)
   2.218 +		     mres[size*i+j] += (*this)[size*k+j] * m2[size*i+k];
   2.219 +	    *this = mres;
   2.220 +	    return *this;
   2.221 +	    }
   2.222 +	 inline Matrix22<T> operator*(Matrix22<T> m2) const
   2.223 +	    { return Matrix22<T>(*this) *= m2; }
   2.224  
   2.225 -      // Assignment
   2.226 -      inline Matrix22<T>& operator=(const Matrix22<T>& m2)
   2.227 -	 { 
   2.228 -	 m[0] = m2[0]; m[2] = m2[2]; 
   2.229 -	 m[1] = m2[1]; m[3] = m2[3]; 
   2.230 -	 return *this;
   2.231 -	 }
   2.232 -      inline Matrix22<T>& assign(T a0, T a1, T a2, T a3)
   2.233 -	 { 
   2.234 -	 m[0] = a0; m[2] = a2;
   2.235 -	 m[1] = a1; m[3] = a3; 
   2.236 -	 return *this;
   2.237 -	 }
   2.238  
   2.239 -      // Comparison
   2.240 -      inline bool operator==(Matrix22<T>& m2) const
   2.241 -	 {
   2.242 -	 return  
   2.243 -	       m[0] == m2[0] && m[2] == m2[2] &&
   2.244 -	       m[1] == m2[1] && m[3] == m2[3]; 
   2.245 -	 }
   2.246 -      inline bool operator!=(Matrix22<T>& m2) const
   2.247 -	 { return ! (*this == m2); } 
   2.248 +	 inline Vector2<T> getcol(const unsigned int i) const
   2.249 +	    {
   2.250 +	    assert(i<2);
   2.251 +	    const int size=2;
   2.252 +	    return Vector2<T>(m[i*size], m[i*size+1]);
   2.253 +	    }
   2.254  
   2.255 -      // Matrix addition
   2.256 -      inline Matrix22<T>& operator+=(const Matrix22<T>& m2)
   2.257 -	 {
   2.258 -	 m[0] += m2[0]; m[2] += m2[2]; 
   2.259 -	 m[1] += m2[1]; m[3] += m2[3];
   2.260 -	 return *this;
   2.261 -	 }
   2.262 -      inline Matrix22<T> operator+(const Matrix22<T>& m2) const
   2.263 -	 { return Matrix22<T>(*this) += m2; }
   2.264 +	 inline Vector2<T> getrow(const unsigned int i) const
   2.265 +	    {
   2.266 +	    assert(i<2);
   2.267 +	    const int size=2;
   2.268 +	    return Vector2<T>(m[i], m[i+size]);
   2.269 +	    }
   2.270 +	 
   2.271 +	 inline Matrix22<T>& setcol(const unsigned int i, const Vector2<T>& v)
   2.272 +	    {
   2.273 +	    assert(i<2);
   2.274 +	    const int size=2;
   2.275 +	    m[i*size] = v[0];
   2.276 +	    m[i*size+1] = v[1];
   2.277 +	    return *this;
   2.278 +	    }
   2.279  
   2.280 -      // Matrix subtraction
   2.281 -      inline Matrix22<T>& operator-=(const Matrix22<T>& m2)
   2.282 -	 {
   2.283 -	 m[0] -= m2[0]; m[2] -= m2[2]; 
   2.284 -	 m[1] -= m2[1]; m[3] -= m2[3];
   2.285 -	 return *this;
   2.286 -	 }
   2.287 -      inline Matrix22<T> operator-(const Matrix22<T>& m2) const
   2.288 -	 { return Matrix22<T>(*this) -= m2; }
   2.289 +	 inline Matrix22<T>& setcol(const unsigned int i, const T a, const T b)
   2.290 +	    {
   2.291 +	    assert(i<2);
   2.292 +	    const int size=2;
   2.293 +	    m[i*size] = a;
   2.294 +	    m[i*size+1] = b;
   2.295 +	    return *this;
   2.296 +	    }
   2.297 +
   2.298 +	 inline Matrix22<T>& setrow(const unsigned int i, const Vector2<T>& v)
   2.299 +	    {
   2.300 +	    assert(i<2);
   2.301 +	    const int size=2;
   2.302 +	    m[i] = v[0];
   2.303 +	    m[i+size] = v[1];
   2.304 +	    return *this;
   2.305 +	    }
   2.306 +
   2.307 +	 inline Matrix22<T>& setrow(const unsigned int i, const T a, const T b)
   2.308 +	    {
   2.309 +	    assert(i<2);
   2.310 +	    const int size=2;
   2.311 +	    m[i] = a;
   2.312 +	    m[i+size] = b;
   2.313 +	    return *this;
   2.314 +	    }
   2.315 +
   2.316 +	 std::string & getstring(std::string & s) const;
   2.317 +
   2.318 +	 inline Matrix22<T>& setidentity()
   2.319 +	    {
   2.320 +	    memset(m, 0, sizeof(m));
   2.321 +	    m[0] = m[3] = 1;
   2.322 +	    return *this;
   2.323 +	    }
   2.324 +
   2.325 +	 };
   2.326 +
   2.327 +      // Scalar multiplication continued
   2.328 +      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const int a)
   2.329 +	 { return Matrix22<T>(m) *= a; }
   2.330 +      template <typename T> inline Matrix22<T> operator*(const int a, const Matrix22<T>& m)
   2.331 +	 { return Matrix22<T>(m) *= a; }
   2.332 +      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const float a)
   2.333 +	 { return Matrix22<T>(m) *= a; }
   2.334 +      template <typename T> inline Matrix22<T> operator*(const float a, const Matrix22<T>& m)
   2.335 +	 { return Matrix22<T>(m) *= a; }
   2.336 +      template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const double a)
   2.337 +	 { return Matrix22<T>(m) *= a; }
   2.338 +      template <typename T> inline Matrix22<T> operator*(const double a, const Matrix22<T>& m)
   2.339 +	 { return Matrix22<T>(m) *= a; }
   2.340 +
   2.341 +
   2.342 +
   2.343 +
   2.344 +      //////////////////////////////////////////////////////////////////////////
   2.345 +      typedef Matrix22<int> Matrix22i;
   2.346 +      typedef Matrix22<float> Matrix22f;
   2.347 +      typedef Matrix22<double> Matrix22d;
   2.348 +
   2.349 +#ifdef NOTHING   
   2.350 +      typedef Matrix33<int> Matrix33i;
   2.351 +      typedef Matrix33<float> Matrix33f;
   2.352 +      typedef Matrix33<double> Matrix33d;
   2.353        
   2.354 -      // Scalar multiplication
   2.355 -      inline Matrix22<T>& operator*=(const int a)
   2.356 -	 {
   2.357 -	 m[0] *= a; m[2] *= a;
   2.358 -	 m[1] *= a; m[3] *= a;
   2.359 -	 return *this;
   2.360 -	 }
   2.361 -      inline Matrix22<T>& operator*=(const float a)
   2.362 -	 {
   2.363 -	 m[0] *= a; m[2] *= a;
   2.364 -	 m[1] *= a; m[3] *= a;
   2.365 -	 return *this;
   2.366 -	 }
   2.367 -      inline Matrix22<T>& operator*=(const double a)
   2.368 -	 {
   2.369 -	 m[0] *= a; m[2] *= a;
   2.370 -	 m[1] *= a; m[3] *= a;
   2.371 -	 return *this;
   2.372 -	 }
   2.373 +      typedef Matrix44<int> Matrix44i;
   2.374 +      typedef Matrix44<float> Matrix44f;
   2.375 +      typedef Matrix44<double> Matrix44d;
   2.376 +      
   2.377 +#endif
   2.378  
   2.379 -      // Scalar division
   2.380 -      inline Matrix22<T>& operator/=(const int a)
   2.381 -	 {
   2.382 -	 assert(a!=0);
   2.383 -	 m[0] /= a; m[2] /= a;
   2.384 -	 m[1] /= a; m[3] /= a;
   2.385 -	 return *this;
   2.386 -	 }
   2.387 -      inline Matrix22<T>& operator/=(const float a)
   2.388 -	 {
   2.389 -	 assert(a!=0);
   2.390 -	 m[0] /= a; m[2] /= a;
   2.391 -	 m[1] /= a; m[3] /= a;
   2.392 -	 return *this;
   2.393 -	 }
   2.394 -      inline Matrix22<T>& operator/=(const double a)
   2.395 -	 {
   2.396 -	 assert(a!=0);
   2.397 -	 m[0] /= a; m[2] /= a;
   2.398 -	 m[1] /= a; m[3] /= a;
   2.399 -	 return *this;
   2.400 -	 }
   2.401 -      inline Matrix22<T> operator/(const int a)
   2.402 -	 { return Matrix22<T>(*this) /= a; }
   2.403 -      inline Matrix22<T> operator/(const float a)
   2.404 -	 { return Matrix22<T>(*this) /= a; }
   2.405 -      inline Matrix22<T> operator/(const double a)
   2.406 -	 { return Matrix22<T>(*this) /= a; }
   2.407  
   2.408 -      // Matrix multiplication
   2.409 -      // Not sure if this should be inlined at all.  Sure the compiler will
   2.410 -      // probably ignoe the inline request here, but maybe it won't and maybe
   2.411 -      // that would be bad.  Needs real testing.
   2.412 -      inline Matrix22<T>& operator*=(Matrix22<T> m2)
   2.413 -	 {
   2.414 -	 const int size=2;
   2.415 -	 Matrix22<T> mres(0);
   2.416 -	 int i, j, k;
   2.417 -	 for (i=0; i<size; ++i)
   2.418 -	    for (j=0; j<size; ++j)
   2.419 -	       for (k=0; k<size; ++k)
   2.420 -		  mres[size*i+j] += (*this)[size*k+j] * m2[size*i+k];
   2.421 -	 *this = mres;
   2.422 -	 return *this;
   2.423 -	 }
   2.424 -      inline Matrix22<T> operator*(Matrix22<T> m2) const
   2.425 -	 { return Matrix22<T>(*this) *= m2; }
   2.426  
   2.427 -      // methods
   2.428 -      inline double det() const
   2.429 +      //////////////////////////////////////////////////////////////////////////
   2.430 +      // Matrix functions
   2.431 +
   2.432 +      template <typename T> inline double det(const Matrix22<T> m)
   2.433  	 { return (double) (m[0] * m[3]) - (double) (m[1] * m[2]); }
   2.434  
   2.435 -      std::string & getstring(std::string & s) const;
   2.436 -
   2.437 -      inline Matrix22<T>& setidentity()
   2.438 -	 {
   2.439 -	 memset(m, 0, sizeof(m));
   2.440 -	 m[0] = m[3] = 1;
   2.441 -	 return *this;
   2.442 -	 }
   2.443 -
   2.444 -      inline Matrix22<T> transpose() const
   2.445 +      template <typename T> inline Matrix22<T> transpose(const Matrix22<T>& m)
   2.446  	 {
   2.447  	 const int size=2;
   2.448  	 Matrix22<T> mres(0);
   2.449 @@ -195,44 +296,15 @@
   2.450  	       mres[size*i+j] = m[size*j+i];
   2.451  	 return mres;
   2.452  	 }
   2.453 -      };
   2.454  
   2.455 -   // Scalar multiplication continued
   2.456 -   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const int a)
   2.457 -      { return Matrix22<T>(m) *= a; }
   2.458 -   template <typename T> inline Matrix22<T> operator*(const int a, const Matrix22<T>& m)
   2.459 -      { return Matrix22<T>(m) *= a; }
   2.460 -   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const float a)
   2.461 -      { return Matrix22<T>(m) *= a; }
   2.462 -   template <typename T> inline Matrix22<T> operator*(const float a, const Matrix22<T>& m)
   2.463 -      { return Matrix22<T>(m) *= a; }
   2.464 -   template <typename T> inline Matrix22<T> operator*(const Matrix22<T>& m, const double a)
   2.465 -      { return Matrix22<T>(m) *= a; }
   2.466 -   template <typename T> inline Matrix22<T> operator*(const double a, const Matrix22<T>& m)
   2.467 -      { return Matrix22<T>(m) *= a; }
   2.468 -
   2.469 -
   2.470 -
   2.471 -
   2.472 -   //////////////////////////////////////////////////////////////////////////
   2.473 -   typedef Matrix22<int> Matrix22i;
   2.474 -   typedef Matrix22<float> Matrix22f;
   2.475 -   typedef Matrix22<double> Matrix22d;
   2.476 -
   2.477 -#ifdef NOTHING   
   2.478 -   typedef Matrix33<int> Matrix33i;
   2.479 -   typedef Matrix33<float> Matrix33f;
   2.480 -   typedef Matrix33<double> Matrix33d;
   2.481 -      
   2.482 -   typedef Matrix44<int> Matrix44i;
   2.483 -   typedef Matrix44<float> Matrix44f;
   2.484 -   typedef Matrix44<double> Matrix44d;
   2.485 -      
   2.486 -#endif
   2.487 +      } // namespace Math
   2.488     } // namespace arda
   2.489  
   2.490 +
   2.491  ////////////////////////////////////////////////////////////////////////////////
   2.492 -template <typename T> std::string & arda::Matrix22<T>::getstring(std::string & s) const
   2.493 +// Matrix methods
   2.494 +
   2.495 +template <typename T> std::string & arda::Math::Matrix22<T>::getstring(std::string & s) const
   2.496     {
   2.497     s.clear();
   2.498     std::stringstream ss;
     3.1 --- a/include/Vector.h	Fri Oct 07 11:09:15 2011 -0500
     3.2 +++ b/include/Vector.h	Fri Oct 07 13:33:48 2011 -0500
     3.3 @@ -10,449 +10,453 @@
     3.4  
     3.5  namespace arda 
     3.6     {
     3.7 -   ////////////////////////////////////////////////////////////////////////////////
     3.8 -   // Vectors
     3.9 -   //
    3.10 -   // Templated classes for 2, 3, and 4 dimensional vectors of any type,
    3.11 -   //     although these types are really intended to be int, float, and double 
    3.12 -   //     only.  (i.e. if you roll your own base type, it had better act like
    3.13 -   //     a numeric scalar.)
    3.14 -   //
    3.15 -   // Supported operations on vectors
    3.16 -   //
    3.17 -   // Construction
    3.18 -   //     Vector v;             // Default zero vector.
    3.19 -   //     Vector v2(v);         // Construct from other vector.
    3.20 -   //     T a, b;
    3.21 -   //     Vector v(a, b)        // Construct from values.
    3.22 -   // []
    3.23 -   //     If you have a Vector named v you can access it's member elements as
    3.24 -   //     v[0], v[1], v[2], and v[3]. in addtion to v.x, v.y, v.z, and v.w
    3.25 -   // == != + += - -=
    3.26 -   //     Defined for operations on two vectors of the same type.
    3.27 -   // * *= / /=
    3.28 -   //     Defined for scalar multiplication/division with int, float, and double.
    3.29 -   //     * is defined as a standalone template operator for the scalar * Vector form.
    3.30 -   // assign(T a, T b [, T b [, T d]])
    3.31 -   //     Directly assign values to the vector.
    3.32 -   // cross(Vector3& v2)
    3.33 -   // cross(Vector3& v2, Vector3& vres)
    3.34 -   //     Cross product of the vector with v2.  Only defined for Vector3 types.
    3.35 -   //     The version that takes two arguments stores the result in the second
    3.36 -   //     argument.  This may be faster since no temporary object is created
    3.37 -   //     during the operation.
    3.38 -   // dot(Vector& v2)
    3.39 -   //     The dot product of the vector with v2.
    3.40 -   // get_angle(Vector& v2)
    3.41 -   //     Returns the angle (in radians) between the vector and v2.
    3.42 -   //     Not meaningful for int vectors.
    3.43 -   // get_anglen(Vector& v2)
    3.44 -   //     Like get_angle(), but the vectors must already be normalized.
    3.45 -   //     Not meaningful for int vectors.
    3.46 -   // getstring(std::string & s)
    3.47 -   //     Returns a C string representation of the vector.
    3.48 -   //     This is one of the few non-inline Vector functions.
    3.49 -   // length()
    3.50 -   //     Get the length of a vector.
    3.51 -   //     Returns a double for all vector types.
    3.52 -   // normalize()
    3.53 -   //     Normalizes the vecor.
    3.54 -   //     Not meaningful for int vectors.
    3.55 -   // proj(Vector& v2, Vector& vres)
    3.56 -   //     Calculates the projection of the vector onto v2 and returns the result 
    3.57 -   //     in vres. 
    3.58 -   //     Not meaningful for int vectors.
    3.59 -   //     Returns vres.
    3.60 -   //
    3.61 -   //
    3.62 -   // Note that the following typedefs are defined for convenience:
    3.63 -   //
    3.64 -   // typedef Vector2<int> Vector2i;
    3.65 -   // typedef Vector2<float> Vector2f;
    3.66 -   // typedef Vector2<double> Vector2d;
    3.67 -   // typedef Vector3<int> Vector3i;
    3.68 -   // typedef Vector3<float> Vector3f;
    3.69 -   // typedef Vector3<double> Vector3d;
    3.70 -   // typedef Vector4<int> Vector4i;
    3.71 -   // typedef Vector4<float> Vector4f;
    3.72 -   // typedef Vector4<double> Vector4d;
    3.73 -   //
    3.74 -   //
    3.75 -   //
    3.76 -   // Scalar multiplication and division are overloaded.  Why?
    3.77 -   // I'm not sure I need this.  The idea is to support all likely scalar
    3.78 -   // types without having to rely on implicit or explicit type conversion,
    3.79 -   // especially not a conversion that might result in loss of precision.
    3.80 -   // Multiplying a float Vector by a double scalar, for example.
    3.81 -   //
    3.82 -   // Is there a better way to do this?  Am I worrying about nothing?
    3.83 -   //
    3.84 -   //
    3.85 -   //
    3.86 -   // I suppose I could make this into a single Vector template on two parameters
    3.87 -   // template <typename T, unsigned int N> class Vector
    3.88 -   //    {
    3.89 -   //    T v[N];
    3.90 -   //
    3.91 -   //    etc...
    3.92 -   //
    3.93 -   // This might buy perfect generality, but it would come at the expense of
    3.94 -   // loops in all of the code, so I don't think it would really be worth it.
    3.95 +   namespace Math
    3.96 +      {
    3.97 +      ////////////////////////////////////////////////////////////////////////////////
    3.98 +      // Vectors
    3.99 +      //
   3.100 +      // Templated classes for 2, 3, and 4 dimensional vectors of any type,
   3.101 +      //     although these types are really intended to be int, float, and double 
   3.102 +      //     only.  (i.e. if you roll your own base type, it had better act like
   3.103 +      //     a numeric scalar.)
   3.104 +      //
   3.105 +      // Supported operations on vectors
   3.106 +      //
   3.107 +      // Construction
   3.108 +      //     Vector v;             // Default zero vector.
   3.109 +      //     Vector v2(v);         // Construct from other vector.
   3.110 +      //     T a, b;
   3.111 +      //     Vector v(a, b)        // Construct from values.
   3.112 +      // []
   3.113 +      //     If you have a Vector named v you can access it's member elements as
   3.114 +      //     v[0], v[1], v[2], and v[3]. in addtion to v.x, v.y, v.z, and v.w
   3.115 +      // == != + += - -=
   3.116 +      //     Defined for operations on two vectors of the same type.
   3.117 +      // * *= / /=
   3.118 +      //     Defined for scalar multiplication/division with int, float, and double.
   3.119 +      //     * is defined as a standalone template operator for the scalar * Vector form.
   3.120 +      // assign(T a, T b [, T b [, T d]])
   3.121 +      //     Directly assign values to the vector.
   3.122 +      // cross(Vector3& v2)
   3.123 +      // cross(Vector3& v2, Vector3& vres)
   3.124 +      //     Cross product of the vector with v2.  Only defined for Vector3 types.
   3.125 +      //     The version that takes two arguments stores the result in the second
   3.126 +      //     argument.  This may be faster since no temporary object is created
   3.127 +      //     during the operation.
   3.128 +      // dot(Vector& v2)
   3.129 +      //     The dot product of the vector with v2.
   3.130 +      // get_angle(Vector& v2)
   3.131 +      //     Returns the angle (in radians) between the vector and v2.
   3.132 +      //     Not meaningful for int vectors.
   3.133 +      // get_anglen(Vector& v2)
   3.134 +      //     Like get_angle(), but the vectors must already be normalized.
   3.135 +      //     Not meaningful for int vectors.
   3.136 +      // getstring(std::string & s)
   3.137 +      //     Returns a C string representation of the vector.
   3.138 +      //     This is one of the few non-inline Vector functions.
   3.139 +      // length()
   3.140 +      //     Get the length of a vector.
   3.141 +      //     Returns a double for all vector types.
   3.142 +      // normalize()
   3.143 +      //     Normalizes the vecor.
   3.144 +      //     Not meaningful for int vectors.
   3.145 +      // proj(Vector& v2, Vector& vres)
   3.146 +      //     Calculates the projection of the vector onto v2 and returns the result 
   3.147 +      //     in vres. 
   3.148 +      //     Not meaningful for int vectors.
   3.149 +      //     Returns vres.
   3.150 +      //
   3.151 +      //
   3.152 +      // Note that the following typedefs are defined for convenience:
   3.153 +      //
   3.154 +      // typedef Vector2<int> Vector2i;
   3.155 +      // typedef Vector2<float> Vector2f;
   3.156 +      // typedef Vector2<double> Vector2d;
   3.157 +      // typedef Vector3<int> Vector3i;
   3.158 +      // typedef Vector3<float> Vector3f;
   3.159 +      // typedef Vector3<double> Vector3d;
   3.160 +      // typedef Vector4<int> Vector4i;
   3.161 +      // typedef Vector4<float> Vector4f;
   3.162 +      // typedef Vector4<double> Vector4d;
   3.163 +      //
   3.164 +      //
   3.165 +      //
   3.166 +      // Scalar multiplication and division are overloaded.  Why?
   3.167 +      // I'm not sure I need this.  The idea is to support all likely scalar
   3.168 +      // types without having to rely on implicit or explicit type conversion,
   3.169 +      // especially not a conversion that might result in loss of precision.
   3.170 +      // Multiplying a float Vector by a double scalar, for example.
   3.171 +      //
   3.172 +      // Is there a better way to do this?  Am I worrying about nothing?
   3.173 +      //
   3.174 +      //
   3.175 +      //
   3.176 +      // I suppose I could make this into a single Vector template on two parameters
   3.177 +      // template <typename T, unsigned int N> class Vector
   3.178 +      //    {
   3.179 +      //    T v[N];
   3.180 +      //
   3.181 +      //    etc...
   3.182 +      //
   3.183 +      // This might buy perfect generality, but it would come at the expense of
   3.184 +      // loops in all of the code, so I don't think it would really be worth it.
   3.185 +      
   3.186 +      
   3.187 +      /////////////////////////////////////////////////////////////////////////////
   3.188 +      template <typename T> class Vector2
   3.189 +	 {
   3.190 +	 public:
   3.191 +	 T x, y;
   3.192 +	 
   3.193 +	 // Constructors
   3.194 +	 Vector2() {}
   3.195 +	 explicit Vector2(const T a) : x (a), y(a) {}
   3.196 +	 Vector2(const T a, const T b) : x (a), y(b) {}
   3.197  
   3.198 +	 // Array indexing
   3.199 +	 inline T& operator[](unsigned int i)
   3.200 +	    { assert (i<2); if (i==0) return x; return y; }
   3.201 +	 inline T operator[](unsigned int i) const
   3.202 +	    { assert (i<2); if (i==0) return x; return y; }
   3.203  
   3.204 -   /////////////////////////////////////////////////////////////////////////////
   3.205 -   template <typename T> class Vector2
   3.206 -      {
   3.207 -      public:
   3.208 -      T x, y;
   3.209 +	 // Assignment
   3.210 +	 inline Vector2<T>& operator=(const Vector2<T>& v2)
   3.211 +	    { x = v2.x; y = v2.y; return *this; }
   3.212 +	 inline Vector2<T>& assign(const T a, const T b)
   3.213 +	    { x = a; y = b; return *this; }
   3.214  
   3.215 -      // Constructors
   3.216 -      Vector2() {}
   3.217 -      Vector2(T a) : x (a), y(a) {}
   3.218 -      Vector2(T a, T b) : x (a), y(b) {}
   3.219 +	 // Comparison
   3.220 +	 inline bool operator==(Vector2<T>& v2) const
   3.221 +	    { return ((x == v2.x) && (y == v2.y)); }
   3.222 +	 inline bool operator!=(Vector2<T>& v2) const
   3.223 +	    { return ! (*this == v2); }
   3.224  
   3.225 -      // Array indexing
   3.226 -      inline T& operator[](unsigned int i)
   3.227 -	 { assert (i<2); if (i==0) return x; return y; }
   3.228 -      inline const T& operator[](unsigned int i) const
   3.229 -	 { assert (i<2); if (i==0) return x; return y; }
   3.230 +	 // Vector addition
   3.231 +	 inline Vector2<T>& operator+=(const Vector2<T>& v2)
   3.232 +	    { x += v2.x; y += v2.y; return *this; }
   3.233 +	 inline Vector2<T> operator+(const Vector2<T>& v2) const
   3.234 +	    { return Vector2<T>(*this) += v2; }
   3.235  
   3.236 -      // Assignment
   3.237 -      inline Vector2<T>& operator=(const Vector2<T>& v2)
   3.238 -	 { x = v2.x; y = v2.y; return *this; }
   3.239 -      inline Vector2<T>& assign(const T a, const T b)
   3.240 -	 { x = a; y = b; return *this; }
   3.241 +	 // Vector subtraction
   3.242 +	 inline Vector2<T>& operator-=(const Vector2<T>& v2)
   3.243 +	    { x -= v2.x; y -= v2.y; return *this; }
   3.244 +	 inline Vector2<T> operator-(const Vector2<T>& v2) const
   3.245 +	    { return Vector2<T>(*this) -= v2; }
   3.246  
   3.247 -      // Comparison
   3.248 -      inline bool operator==(Vector2<T>& v2) const
   3.249 -	 { return ((x == v2.x) && (y == v2.y)); }
   3.250 -      inline bool operator!=(Vector2<T>& v2) const
   3.251 -	 { return ! (*this == v2); }
   3.252 +	 // Scalar multiplication
   3.253 +	 inline Vector2<T>& operator*=(const int a)
   3.254 +	    { x *= a; y *= a; return *this; }
   3.255 +	 inline Vector2<T>& operator*=(const float a)
   3.256 +	    { x *= a; y *= a; return *this; }
   3.257 +	 inline Vector2<T>& operator*=(const double a)
   3.258 +	    { x *= a; y *= a; return *this; }
   3.259 +	 inline Vector2<T> operator*(const int a)
   3.260 +	    { return Vector2<T>(*this) *= a;}
   3.261 +	 inline Vector2<T> operator*(const float a)
   3.262 +	    { return Vector2<T>(*this) *= a;}
   3.263 +	 inline Vector2<T> operator*(const double a)
   3.264 +	    { return Vector2<T>(*this) *= a;}
   3.265  
   3.266 -      // Vector addition
   3.267 -      inline Vector2<T>& operator+=(const Vector2<T>& v2)
   3.268 -	 { x += v2.x; y += v2.y; return *this; }
   3.269 -      inline Vector2<T> operator+(const Vector2<T>& v2) const
   3.270 -	 { return Vector2<T>(*this) += v2; }
   3.271 +	 // Scalar division
   3.272 +	 inline Vector2<T>& operator/=(const int a)
   3.273 +	    { assert(a!=0); x /= a; y /= a; return *this; }
   3.274 +	 inline Vector2<T>& operator/=(const float a)
   3.275 +	    { assert(a!=0); x /= a; y /= a; return *this; }
   3.276 +	 inline Vector2<T>& operator/=(const double a)
   3.277 +	    { assert(a!=0); x /= a; y /= a; return *this; }
   3.278 +	 inline Vector2<T> operator/(const int a)
   3.279 +	    { return Vector2<T>(*this) /= a;}
   3.280 +	 inline Vector2<T> operator/(const float a)
   3.281 +	    { return Vector2<T>(*this) /= a;}
   3.282 +	 inline Vector2<T> operator/(const double a)
   3.283 +	    { return Vector2<T>(*this) /= a;}
   3.284  
   3.285 -      // Vector subtraction
   3.286 -      inline Vector2<T>& operator-=(const Vector2<T>& v2)
   3.287 -	 { x -= v2.x; y -= v2.y; return *this; }
   3.288 -      inline Vector2<T> operator-(const Vector2<T>& v2) const
   3.289 -	 { return Vector2<T>(*this) -= v2; }
   3.290  
   3.291 -      // Scalar multiplication
   3.292 -      inline Vector2<T>& operator*=(const int a)
   3.293 -	 { x *= a; y *= a; return *this; }
   3.294 -      inline Vector2<T>& operator*=(const float a)
   3.295 -	 { x *= a; y *= a; return *this; }
   3.296 -      inline Vector2<T>& operator*=(const double a)
   3.297 -	 { x *= a; y *= a; return *this; }
   3.298 -      inline Vector2<T> operator*(const int a)
   3.299 -	 { return Vector2<T>(*this) *= a;}
   3.300 -      inline Vector2<T> operator*(const float a)
   3.301 -	 { return Vector2<T>(*this) *= a;}
   3.302 -      inline Vector2<T> operator*(const double a)
   3.303 -	 { return Vector2<T>(*this) *= a;}
   3.304 +	 // methods
   3.305  
   3.306 -      // Scalar division
   3.307 -      inline Vector2<T>& operator/=(const int a)
   3.308 -	 { assert(a!=0); x /= a; y /= a; return *this; }
   3.309 -      inline Vector2<T>& operator/=(const float a)
   3.310 -	 { assert(a!=0); x /= a; y /= a; return *this; }
   3.311 -      inline Vector2<T>& operator/=(const double a)
   3.312 -	 { assert(a!=0); x /= a; y /= a; return *this; }
   3.313 -      inline Vector2<T> operator/(const int a)
   3.314 -	 { return Vector2<T>(*this) /= a;}
   3.315 -      inline Vector2<T> operator/(const float a)
   3.316 -	 { return Vector2<T>(*this) /= a;}
   3.317 -      inline Vector2<T> operator/(const double a)
   3.318 -	 { return Vector2<T>(*this) /= a;}
   3.319 +	 inline T dot(const Vector2<T>& v2) const
   3.320 +	    { return x*v2.x + y*v2.y; }
   3.321  
   3.322 +	 inline double get_angle(Vector2<T>& v2)
   3.323 +	    { 
   3.324 +	    double tmp = dot(v2); 
   3.325 +	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.326 +	    }
   3.327  
   3.328 -      // methods
   3.329 +	 inline double get_anglen(Vector2<T>& v2)
   3.330 +	    { return acos((double) dot(v2)); }
   3.331  
   3.332 -      inline T dot(const Vector2<T>& v2) const
   3.333 -	 { return x*v2.x + y*v2.y; }
   3.334 +	 std::string & getstring(std::string & s) const;
   3.335  
   3.336 -      inline double get_angle(Vector2<T>& v2)
   3.337 -	 { 
   3.338 -	 double tmp = dot(v2); 
   3.339 -	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.340 -	 }
   3.341 +	 inline double length(void) const
   3.342 +	    { return sqrt((dot(*this))); }
   3.343  
   3.344 -      inline double get_anglen(Vector2<T>& v2)
   3.345 -	 { return acos((double) dot(v2)); }
   3.346 +	 inline Vector2<T>& normalize()
   3.347 +	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; return *this; }
   3.348  
   3.349 -      std::string & getstring(std::string & s) const;
   3.350 +	 inline Vector2<T> proj(Vector2<T>& v2)
   3.351 +	    { Vector2<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.352 +	 inline Vector2<T>& proj(Vector2<T>& v2, Vector2<T>& vres)
   3.353 +	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.354 +	 };
   3.355  
   3.356 -      inline double length(void) const
   3.357 -	 { return sqrt((dot(*this))); }
   3.358 +      // Scalar multiplication, continued
   3.359 +      template <typename T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
   3.360 +	 { return Vector2<T>(v) *= a;}
   3.361 +      template <typename T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
   3.362 +	 { return Vector2<T>(v) *= a;}
   3.363 +      template <typename T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
   3.364 +	 { return Vector2<T>(v) *= a;}
   3.365  
   3.366 -      inline Vector2<T>& normalize()
   3.367 -	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; return *this; }
   3.368  
   3.369 -      inline Vector2<T> proj(Vector2<T>& v2)
   3.370 -	 { Vector2<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.371 -      inline Vector2<T>& proj(Vector2<T>& v2, Vector2<T>& vres)
   3.372 -	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.373 -      };
   3.374 +      /////////////////////////////////////////////////////////////////////////////
   3.375 +      template <typename T> class Vector3
   3.376 +	 {
   3.377 +	 public:
   3.378 +	 T x, y, z;
   3.379  
   3.380 -   // Scalar multiplication, continued
   3.381 -   template <typename T> inline Vector2<T> operator*(const int a, const Vector2<T>& v)
   3.382 -      { return Vector2<T>(v) *= a;}
   3.383 -   template <typename T> inline Vector2<T> operator*(const float a, const Vector2<T> v)
   3.384 -      { return Vector2<T>(v) *= a;}
   3.385 -   template <typename T> inline Vector2<T> operator*(const double a, const Vector2<T>& v)
   3.386 -      { return Vector2<T>(v) *= a;}
   3.387 +	 // Constructors
   3.388 +	 Vector3() {}
   3.389 +	 explicit Vector3(const T a) : x (a), y(a), z(a) {}
   3.390 +	 Vector3(const T a, const T b, const T c) : x (a), y(b), z(c) {}
   3.391  
   3.392 +	 // Array indexing
   3.393 +	 inline T& operator[](unsigned int i)
   3.394 +	    { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
   3.395 +	 inline T operator[](unsigned int i) const
   3.396 +	    { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
   3.397  
   3.398 -   /////////////////////////////////////////////////////////////////////////////
   3.399 -   template <typename T> class Vector3
   3.400 -      {
   3.401 -      public:
   3.402 -      T x, y, z;
   3.403 +	 // Assignment
   3.404 +	 inline Vector3<T>& operator=(const Vector3<T>& v2)
   3.405 +	    { x = v2.x; y = v2.y; z = v2.z; return *this; }
   3.406 +	 inline Vector3<T>& assign(const T a, const T b, const T c)
   3.407 +	    { x = a; y = b; z = c; return *this; }
   3.408  
   3.409 -      // Constructors
   3.410 -      Vector3() {}
   3.411 -      Vector3(T a) : x (a), y(a), z(a) {}
   3.412 -      Vector3(T a, T b, T c) : x (a), y(b), z(c) {}
   3.413 +	 // Comparison
   3.414 +	 inline bool operator==(Vector3<T>& v2) const
   3.415 +	    { return ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
   3.416 +	 inline bool operator!=(Vector3<T>& v2) const
   3.417 +	    { return ! (*this == v2); }
   3.418  
   3.419 -      // Array indexing
   3.420 -      inline T& operator[](unsigned int i)
   3.421 -	 { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
   3.422 -      inline const T& operator[](unsigned int i) const
   3.423 -	 { assert (i<3); if (i==0) return x; if (i==1) return y; return z; }
   3.424 +	 // Vector addition
   3.425 +	 inline Vector3<T>& operator+=(const Vector3<T>& v2)
   3.426 +	    { x += v2.x; y += v2.y; z += v2.z; return *this; }
   3.427 +	 inline Vector3<T> operator+(const Vector3<T>& v2) const
   3.428 +	    { return Vector3<T>(*this) += v2; }
   3.429  
   3.430 -      // Assignment
   3.431 -      inline Vector3<T>& operator=(const Vector3<T>& v2)
   3.432 -	 { x = v2.x; y = v2.y; z = v2.z; return *this; }
   3.433 -      inline Vector3<T>& assign(const T a, const T b, const T c)
   3.434 -	 { x = a; y = b; z = c; return *this; }
   3.435 +	 // Vector subtraction
   3.436 +	 inline Vector3<T>& operator-=(const Vector3<T>& v2)
   3.437 +	    { x -= v2.x; y -= v2.y; z -= v2.z; return *this; }
   3.438 +	 inline Vector3<T> operator-(const Vector3<T>& v2) const
   3.439 +	    { return Vector3<T>(*this) -= v2; }
   3.440  
   3.441 -      // Comparison
   3.442 -      inline bool operator==(Vector3<T>& v2) const
   3.443 -	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z)); }
   3.444 -      inline bool operator!=(Vector3<T>& v2) const
   3.445 -	 { return ! (*this == v2); }
   3.446 +	 // Scalar multiplication
   3.447 +	 inline Vector3<T>& operator*=(const int a)
   3.448 +	    { x *= a; y *= a; z *= a; return *this; }
   3.449 +	 inline Vector3<T>& operator*=(const float a)
   3.450 +	    { x *= a; y *= a; z *= a; return *this; }
   3.451 +	 inline Vector3<T>& operator*=(const double a)
   3.452 +	    { x *= a; y *= a; z *= a; return *this; }
   3.453 +	 inline Vector3<T> operator*(const int a)
   3.454 +	    { return Vector3<T>(*this) *= a;}
   3.455 +	 inline Vector3<T> operator*(const float a)
   3.456 +	    { return Vector3<T>(*this) *= a;}
   3.457 +	 inline Vector3<T> operator*(const double a)
   3.458 +	    { return Vector3<T>(*this) *= a;}
   3.459  
   3.460 -      // Vector addition
   3.461 -      inline Vector3<T>& operator+=(const Vector3<T>& v2)
   3.462 -	 { x += v2.x; y += v2.y; z += v2.z; return *this; }
   3.463 -      inline Vector3<T> operator+(const Vector3<T>& v2) const
   3.464 -	 { return Vector3<T>(*this) += v2; }
   3.465 +	 // Scalar division
   3.466 +	 inline Vector3<T>& operator/=(const int a)
   3.467 +	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.468 +	 inline Vector3<T>& operator/=(const float a)
   3.469 +	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.470 +	 inline Vector3<T>& operator/=(const double a)
   3.471 +	    { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.472 +	 inline Vector3<T> operator/(const int a)
   3.473 +	    { return Vector3<T>(*this) /= a;}
   3.474 +	 inline Vector3<T> operator/(const float a)
   3.475 +	    { return Vector3<T>(*this) /= a;}
   3.476 +	 inline Vector3<T> operator/(const double a)
   3.477 +	    { return Vector3<T>(*this) /= a;}
   3.478  
   3.479 -      // Vector subtraction
   3.480 -      inline Vector3<T>& operator-=(const Vector3<T>& v2)
   3.481 -	 { x -= v2.x; y -= v2.y; z -= v2.z; return *this; }
   3.482 -      inline Vector3<T> operator-(const Vector3<T>& v2) const
   3.483 -	 { return Vector3<T>(*this) -= v2; }
   3.484  
   3.485 -      // Scalar multiplication
   3.486 -      inline Vector3<T>& operator*=(const int a)
   3.487 -	 { x *= a; y *= a; z *= a; return *this; }
   3.488 -      inline Vector3<T>& operator*=(const float a)
   3.489 -	 { x *= a; y *= a; z *= a; return *this; }
   3.490 -      inline Vector3<T>& operator*=(const double a)
   3.491 -	 { x *= a; y *= a; z *= a; return *this; }
   3.492 -      inline Vector3<T> operator*(const int a)
   3.493 -	 { return Vector3<T>(*this) *= a;}
   3.494 -      inline Vector3<T> operator*(const float a)
   3.495 -	 { return Vector3<T>(*this) *= a;}
   3.496 -      inline Vector3<T> operator*(const double a)
   3.497 -	 { return Vector3<T>(*this) *= a;}
   3.498 +	 // methods
   3.499  
   3.500 -      // Scalar division
   3.501 -      inline Vector3<T>& operator/=(const int a)
   3.502 -	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.503 -      inline Vector3<T>& operator/=(const float a)
   3.504 -	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.505 -      inline Vector3<T>& operator/=(const double a)
   3.506 -	 { assert(a!=0); x /= a; y /= a; z /= a; return *this; }
   3.507 -      inline Vector3<T> operator/(const int a)
   3.508 -	 { return Vector3<T>(*this) /= a;}
   3.509 -      inline Vector3<T> operator/(const float a)
   3.510 -	 { return Vector3<T>(*this) /= a;}
   3.511 -      inline Vector3<T> operator/(const double a)
   3.512 -	 { return Vector3<T>(*this) /= a;}
   3.513 +	 inline Vector3<T>& cross(Vector3<T>& v2, Vector3<T>& vres)
   3.514 +	    { 
   3.515 +	    vres.x =  y*v2.z - v2.y*z;
   3.516 +	    vres.y = -x*v2.z + v2.x*z;
   3.517 +	    vres.z =  x*v2.y - v2.x*y;
   3.518 +	    return vres;
   3.519 +	    }
   3.520 +	 inline Vector3<T> cross(Vector3<T>& v2)
   3.521 +	    { 
   3.522 +	    Vector3<T> vres;
   3.523 +	    vres.x =  y*v2.z - v2.y*z;
   3.524 +	    vres.y = -x*v2.z + v2.x*z;
   3.525 +	    vres.z =  x*v2.y - v2.x*y;
   3.526 +	    return vres;
   3.527 +	    }
   3.528  
   3.529 +	 inline T dot(const Vector3<T>& v2) const
   3.530 +	    { return x*v2.x + y*v2.y + z*v2.z; }
   3.531  
   3.532 -      // methods
   3.533 +	 inline double get_angle(Vector3<T>& v2)
   3.534 +	    { 
   3.535 +	    double tmp = dot(v2); 
   3.536 +	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.537 +	    }
   3.538  
   3.539 -      inline Vector3<T>& cross(Vector3<T>& v2, Vector3<T>& vres)
   3.540 -	 { 
   3.541 -	 vres.x =  y*v2.z - v2.y*z;
   3.542 -	 vres.y = -x*v2.z + v2.x*z;
   3.543 -	 vres.z =  x*v2.y - v2.x*y;
   3.544 -	 return vres;
   3.545 -	 }
   3.546 -      inline Vector3<T> cross(Vector3<T>& v2)
   3.547 -	 { 
   3.548 -	 Vector3<T> vres;
   3.549 -	 vres.x =  y*v2.z - v2.y*z;
   3.550 -	 vres.y = -x*v2.z + v2.x*z;
   3.551 -	 vres.z =  x*v2.y - v2.x*y;
   3.552 -	 return vres;
   3.553 -	 }
   3.554 +	 inline double get_anglen(Vector3<T>& v2)
   3.555 +	    { return acos((double) dot(v2)); }
   3.556  
   3.557 -      inline T dot(const Vector3<T>& v2) const
   3.558 -	 { return x*v2.x + y*v2.y + z*v2.z; }
   3.559 +	 std::string & getstring(std::string & s) const;
   3.560  
   3.561 -      inline double get_angle(Vector3<T>& v2)
   3.562 -	 { 
   3.563 -	 double tmp = dot(v2); 
   3.564 -	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.565 -	 }
   3.566 +	 inline double length(void) const
   3.567 +	    { return sqrt((dot(*this))); }
   3.568  
   3.569 -      inline double get_anglen(Vector3<T>& v2)
   3.570 -	 { return acos((double) dot(v2)); }
   3.571 +	 inline Vector3<T>& normalize()
   3.572 +	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; return *this; }
   3.573  
   3.574 -      std::string & getstring(std::string & s) const;
   3.575 +	 inline Vector3<T> proj(Vector3<T>& v2)
   3.576 +	    { Vector3<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.577 +	 inline Vector3<T>& proj(Vector3<T>& v2, Vector3<T>& vres)
   3.578 +	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.579 +	 };
   3.580  
   3.581 -      inline double length(void) const
   3.582 -	 { return sqrt((dot(*this))); }
   3.583 +      // Scalar multiplication, continued
   3.584 +      template <typename T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
   3.585 +	 { return Vector3<T>(v) *= a;}
   3.586 +      template <typename T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
   3.587 +	 { return Vector3<T>(v) *= a;}
   3.588 +      template <typename T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
   3.589 +	 { return Vector3<T>(v) *= a;}
   3.590  
   3.591 -      inline Vector3<T>& normalize()
   3.592 -	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; return *this; }
   3.593 +      /////////////////////////////////////////////////////////////////////////////
   3.594 +      template <typename T> class Vector4
   3.595 +	 {
   3.596 +	 public:
   3.597 +	 T x, y, z, w;
   3.598  
   3.599 -      inline Vector3<T> proj(Vector3<T>& v2)
   3.600 -	 { Vector3<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.601 -      inline Vector3<T>& proj(Vector3<T>& v2, Vector3<T>& vres)
   3.602 -	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.603 -      };
   3.604 +	 // Constructors
   3.605 +	 Vector4() {}
   3.606 +	 explicit Vector4(const T a) : x (a), y(a), z(a), w(a) {}
   3.607 +	 Vector4(const T a, const T b, const T c, const T d) : x (a), y(b), z(c), w(d) {}
   3.608  
   3.609 -   // Scalar multiplication, continued
   3.610 -   template <typename T> inline Vector3<T> operator*(const int a, const Vector3<T>& v)
   3.611 -      { return Vector3<T>(v) *= a;}
   3.612 -   template <typename T> inline Vector3<T> operator*(const float a, const Vector3<T> v)
   3.613 -      { return Vector3<T>(v) *= a;}
   3.614 -   template <typename T> inline Vector3<T> operator*(const double a, const Vector3<T>& v)
   3.615 -      { return Vector3<T>(v) *= a;}
   3.616 +	 // Array indexing
   3.617 +	 inline T& operator[](unsigned int i)
   3.618 +	    { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
   3.619 +	 inline T operator[](unsigned int i) const
   3.620 +	    { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
   3.621  
   3.622 -   /////////////////////////////////////////////////////////////////////////////
   3.623 -   template <typename T> class Vector4
   3.624 -      {
   3.625 -      public:
   3.626 -      T x, y, z, w;
   3.627 +	 // Assignment
   3.628 +	 inline Vector4<T>& operator=(const Vector4<T>& v2)
   3.629 +	    { x = v2.x; y = v2.y; z = v2.z; w = v2.w; return *this; }
   3.630 +	 inline Vector4<T>& assign(const T a, const T b, const T c, const T d)
   3.631 +	    { x = a; y = b; z = c; w = d; return *this; }
   3.632  
   3.633 -      // Constructors
   3.634 -      Vector4() {}
   3.635 -      Vector4(T a) : x (a), y(a), z(a), w(a) {}
   3.636 -      Vector4(T a, T b, T c, T d) : x (a), y(b), z(c), w(d) {}
   3.637 +	 // Comparison
   3.638 +	 inline bool operator==(Vector4<T>& v2) const
   3.639 +	    { return ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
   3.640 +	 inline bool operator!=(Vector4<T>& v2) const
   3.641 +	    { return ! (*this == v2); }
   3.642  
   3.643 -      // Array indexing
   3.644 -      inline T& operator[](unsigned int i)
   3.645 -	 { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
   3.646 -      inline const T& operator[](unsigned int i) const
   3.647 -	 { assert (i<4); if (i==0) return x; if (i==1) return y; if (i==2) return z; return w; }
   3.648 +	 // Vector addition
   3.649 +	 inline Vector4<T>& operator+=(const Vector4<T>& v2)
   3.650 +	    { x += v2.x; y += v2.y; z += v2.z; w += v2.w; return *this; }
   3.651 +	 inline Vector4<T> operator+(const Vector4<T>& v2) const
   3.652 +	    { return Vector4<T>(*this) += v2; }
   3.653  
   3.654 -      // Assignment
   3.655 -      inline Vector4<T>& operator=(const Vector4<T>& v2)
   3.656 -	 { x = v2.x; y = v2.y; z = v2.z; w = v2.w; return *this; }
   3.657 -      inline Vector4<T>& assign(const T a, const T b, const T c, const T d)
   3.658 -	 { x = a; y = b; z = c; w = d; return *this; }
   3.659 +	 // Vector subtraction
   3.660 +	 inline Vector4<T>& operator-=(const Vector4<T>& v2)
   3.661 +	    { x -= v2.x; y -= v2.y; z -= v2.z; w -= v2.w; return *this; }
   3.662 +	 inline Vector4<T> operator-(const Vector4<T>& v2) const
   3.663 +	    { return Vector4<T>(*this) -= v2; }
   3.664  
   3.665 -      // Comparison
   3.666 -      inline bool operator==(Vector4<T>& v2) const
   3.667 -	 { return ((x == v2.x) && (y == v2.y) && (z == v2.z) && (w == v2.w)); }
   3.668 -      inline bool operator!=(Vector4<T>& v2) const
   3.669 -	 { return ! (*this == v2); }
   3.670 +	 // Scalar multiplication
   3.671 +	 inline Vector4<T>& operator*=(const int a)
   3.672 +	    { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.673 +	 inline Vector4<T>& operator*=(const float a)
   3.674 +	    { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.675 +	 inline Vector4<T>& operator*=(const double a)
   3.676 +	    { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.677 +	 inline Vector4<T> operator*(const int a)
   3.678 +	    { return Vector4<T>(*this) *= a;}
   3.679 +	 inline Vector4<T> operator*(const float a)
   3.680 +	    { return Vector4<T>(*this) *= a;}
   3.681 +	 inline Vector4<T> operator*(const double a)
   3.682 +	    { return Vector4<T>(*this) *= a;}
   3.683  
   3.684 -      // Vector addition
   3.685 -      inline Vector4<T>& operator+=(const Vector4<T>& v2)
   3.686 -	 { x += v2.x; y += v2.y; z += v2.z; w += v2.w; return *this; }
   3.687 -      inline Vector4<T> operator+(const Vector4<T>& v2) const
   3.688 -	 { return Vector4<T>(*this) += v2; }
   3.689 +	 // Scalar division
   3.690 +	 inline Vector4<T>& operator/=(const int a)
   3.691 +	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.692 +	 inline Vector4<T>& operator/=(const float a)
   3.693 +	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.694 +	 inline Vector4<T>& operator/=(const double a)
   3.695 +	    { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.696 +	 inline Vector4<T> operator/(const int a)
   3.697 +	    { return Vector4<T>(*this) /= a;}
   3.698 +	 inline Vector4<T> operator/(const float a)
   3.699 +	    { return Vector4<T>(*this) /= a;}
   3.700 +	 inline Vector4<T> operator/(const double a)
   3.701 +	    { return Vector4<T>(*this) /= a;}
   3.702  
   3.703 -      // Vector subtraction
   3.704 -      inline Vector4<T>& operator-=(const Vector4<T>& v2)
   3.705 -	 { x -= v2.x; y -= v2.y; z -= v2.z; w -= v2.w; return *this; }
   3.706 -      inline Vector4<T> operator-(const Vector4<T>& v2) const
   3.707 -	 { return Vector4<T>(*this) -= v2; }
   3.708  
   3.709 -      // Scalar multiplication
   3.710 -      inline Vector4<T>& operator*=(const int a)
   3.711 -	 { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.712 -      inline Vector4<T>& operator*=(const float a)
   3.713 -	 { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.714 -      inline Vector4<T>& operator*=(const double a)
   3.715 -	 { x *= a; y *= a; z *= a; w *= a; return *this; }
   3.716 -      inline Vector4<T> operator*(const int a)
   3.717 -	 { return Vector4<T>(*this) *= a;}
   3.718 -      inline Vector4<T> operator*(const float a)
   3.719 -	 { return Vector4<T>(*this) *= a;}
   3.720 -      inline Vector4<T> operator*(const double a)
   3.721 -	 { return Vector4<T>(*this) *= a;}
   3.722 +	 // methods
   3.723  
   3.724 -      // Scalar division
   3.725 -      inline Vector4<T>& operator/=(const int a)
   3.726 -	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.727 -      inline Vector4<T>& operator/=(const float a)
   3.728 -	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.729 -      inline Vector4<T>& operator/=(const double a)
   3.730 -	 { assert(a!=0); x /= a; y /= a; z /= a; w /= a; return *this; }
   3.731 -      inline Vector4<T> operator/(const int a)
   3.732 -	 { return Vector4<T>(*this) /= a;}
   3.733 -      inline Vector4<T> operator/(const float a)
   3.734 -	 { return Vector4<T>(*this) /= a;}
   3.735 -      inline Vector4<T> operator/(const double a)
   3.736 -	 { return Vector4<T>(*this) /= a;}
   3.737 +	 inline T dot(const Vector4<T>& v2) const
   3.738 +	    { return x*v2.x + y*v2.y + z*v2.z + w*v2.w; }
   3.739  
   3.740 +	 inline double get_angle(Vector4<T>& v2)
   3.741 +	    { 
   3.742 +	    double tmp = dot(v2); 
   3.743 +	    return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.744 +	    }
   3.745  
   3.746 -      // methods
   3.747 +	 inline double get_anglen(Vector4<T>& v2)
   3.748 +	    { return acos((double) dot(v2)); }
   3.749  
   3.750 -      inline T dot(const Vector4<T>& v2) const
   3.751 -	 { return x*v2.x + y*v2.y + z*v2.z + w*v2.w; }
   3.752 +	 std::string & getstring(std::string & s) const;
   3.753  
   3.754 -      inline double get_angle(Vector4<T>& v2)
   3.755 -	 { 
   3.756 -	 double tmp = dot(v2); 
   3.757 -	 return acos(sqrt(tmp*tmp/(dot(*this)*v2.dot(v2))));
   3.758 -	 }
   3.759 +	 inline double length(void) const
   3.760 +	    { return sqrt((dot(*this))); }
   3.761  
   3.762 -      inline double get_anglen(Vector4<T>& v2)
   3.763 -	 { return acos((double) dot(v2)); }
   3.764 +	 inline Vector4<T>& normalize()
   3.765 +	    { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; w /= l; return *this; }
   3.766  
   3.767 -      std::string & getstring(std::string & s) const;
   3.768 +	 inline Vector4<T> proj(Vector4<T>& v2)
   3.769 +	    { Vector4<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.770 +	 inline Vector4<T>& proj(Vector4<T>& v2, Vector4<T>& vres)
   3.771 +	    { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.772 +	 };
   3.773  
   3.774 -      inline double length(void) const
   3.775 -	 { return sqrt((dot(*this))); }
   3.776 +      // Scalar multiplication, continued
   3.777 +      template <typename T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
   3.778 +	 { return Vector4<T>(v) *= a;}
   3.779 +      template <typename T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
   3.780 +	 { return Vector4<T>(v) *= a;}
   3.781 +      template <typename T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
   3.782 +	 { return Vector4<T>(v) *= a;}
   3.783  
   3.784 -      inline Vector4<T>& normalize()
   3.785 -	 { double l = length(); if (l == 0.0) return *this; x /= l; y /= l; z /= l; w /= l; return *this; }
   3.786 +      /////////////////////////////////////////////////////////////////////////////
   3.787  
   3.788 -      inline Vector4<T> proj(Vector4<T>& v2)
   3.789 -	 { Vector4<T> vres; vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.790 -      inline Vector4<T>& proj(Vector4<T>& v2, Vector4<T>& vres)
   3.791 -	 { vres = v2 * dot(v2)/v2.dot(v2); return vres; }
   3.792 -      };
   3.793 +      typedef Vector2<int> Vector2i;
   3.794 +      typedef Vector2<float> Vector2f;
   3.795 +      typedef Vector2<double> Vector2d;
   3.796  
   3.797 -   // Scalar multiplication, continued
   3.798 -   template <typename T> inline Vector4<T> operator*(const int a, const Vector4<T>& v)
   3.799 -      { return Vector4<T>(v) *= a;}
   3.800 -   template <typename T> inline Vector4<T> operator*(const float a, const Vector4<T> v)
   3.801 -      { return Vector4<T>(v) *= a;}
   3.802 -   template <typename T> inline Vector4<T> operator*(const double a, const Vector4<T>& v)
   3.803 -      { return Vector4<T>(v) *= a;}
   3.804 +      typedef Vector3<int> Vector3i;
   3.805 +      typedef Vector3<float> Vector3f;
   3.806 +      typedef Vector3<double> Vector3d;
   3.807  
   3.808 -   /////////////////////////////////////////////////////////////////////////////
   3.809 +      typedef Vector4<int> Vector4i;
   3.810 +      typedef Vector4<float> Vector4f;
   3.811 +      typedef Vector4<double> Vector4d;
   3.812  
   3.813 -   typedef Vector2<int> Vector2i;
   3.814 -   typedef Vector2<float> Vector2f;
   3.815 -   typedef Vector2<double> Vector2d;
   3.816 -
   3.817 -   typedef Vector3<int> Vector3i;
   3.818 -   typedef Vector3<float> Vector3f;
   3.819 -   typedef Vector3<double> Vector3d;
   3.820 -
   3.821 -   typedef Vector4<int> Vector4i;
   3.822 -   typedef Vector4<float> Vector4f;
   3.823 -   typedef Vector4<double> Vector4d;
   3.824 +      } // namespace Math
   3.825  
   3.826     } // namespace arda
   3.827  
   3.828  ////////////////////////////////////////////////////////////////////////////////
   3.829 -template <typename T> std::string & arda::Vector2<T>::getstring(std::string & s) const
   3.830 +template <typename T> std::string & arda::Math::Vector2<T>::getstring(std::string & s) const
   3.831     {
   3.832     s.clear();
   3.833     std::stringstream ss;
   3.834 @@ -462,7 +466,7 @@
   3.835     }
   3.836  
   3.837  ////////////////////////////////////////////////////////////////////////////////
   3.838 -template <typename T> std::string & arda::Vector3<T>::getstring(std::string & s) const
   3.839 +template <typename T> std::string & arda::Math::Vector3<T>::getstring(std::string & s) const
   3.840     {
   3.841     s.clear();
   3.842     std::stringstream ss;
   3.843 @@ -472,7 +476,7 @@
   3.844     }
   3.845  
   3.846  ////////////////////////////////////////////////////////////////////////////////
   3.847 -template <typename T> std::string & arda::Vector4<T>::getstring(std::string & s) const
   3.848 +template <typename T> std::string & arda::Math::Vector4<T>::getstring(std::string & s) const
   3.849     {
   3.850     s.clear();
   3.851     std::stringstream ss;
     4.1 --- a/src/Matrix.cpp	Fri Oct 07 11:09:15 2011 -0500
     4.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.3 @@ -1,515 +0,0 @@
     4.4 -#include "Matrix.h"
     4.5 -
     4.6 -#include <iostream>
     4.7 -using namespace std;
     4.8 -#include <string>
     4.9 -#include <sstream>
    4.10 -
    4.11 -////////////////////////////////////////////////////////////////////////////////
    4.12 -std::string & arda::Matrix::getstring2m(Matrix2i m, std::string & s)
    4.13 -   {
    4.14 -   s.clear();
    4.15 -   std::stringstream ss;
    4.16 -   ss << "[ "
    4.17 -	 "[ " << m[0] << ", " << m[1] << " ], "
    4.18 -	 "[ " << m[2] << ", " << m[3] << " ] ]";
    4.19 -   s = ss.str();
    4.20 -   return s;   
    4.21 -   }
    4.22 -
    4.23 -////////////////////////////////////////////////////////////////////////////////
    4.24 -std::string & arda::Matrix::getstring2m(Matrix2f m, std::string & s)
    4.25 -   {
    4.26 -   s.clear();
    4.27 -   std::stringstream ss;
    4.28 -   ss << "[ "
    4.29 -	 "[ " << m[0] << ", " << m[1] << " ], "
    4.30 -	 "[ " << m[2] << ", " << m[3] << " ] ]";
    4.31 -   s = ss.str();
    4.32 -   return s;   
    4.33 -   }
    4.34 -
    4.35 -////////////////////////////////////////////////////////////////////////////////
    4.36 -std::string & arda::Matrix::getstring2m(Matrix2d m, std::string & s)
    4.37 -   {
    4.38 -   s.clear();
    4.39 -   std::stringstream ss;
    4.40 -   ss << "[ "
    4.41 -	 "[ " << m[0] << ", " << m[1] << " ], "
    4.42 -	 "[ " << m[2] << ", " << m[3] << " ] ]";
    4.43 -   s = ss.str();
    4.44 -   return s;   
    4.45 -   }
    4.46 -
    4.47 -////////////////////////////////////////////////////////////////////////////////
    4.48 -std::string & arda::Matrix::getstring3m(Matrix3i m, std::string & s)
    4.49 -   {
    4.50 -   s.clear();
    4.51 -   std::stringstream ss;
    4.52 -   ss << "[ "
    4.53 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
    4.54 -	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
    4.55 -	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
    4.56 -   s = ss.str();
    4.57 -   return s;   
    4.58 -   }
    4.59 -
    4.60 -////////////////////////////////////////////////////////////////////////////////
    4.61 -std::string & arda::Matrix::getstring3m(Matrix3f m, std::string & s)
    4.62 -   {
    4.63 -   s.clear();
    4.64 -   std::stringstream ss;
    4.65 -   ss << "[ "
    4.66 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
    4.67 -	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
    4.68 -	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
    4.69 -   s = ss.str();
    4.70 -   return s;   
    4.71 -   }
    4.72 -
    4.73 -////////////////////////////////////////////////////////////////////////////////
    4.74 -std::string & arda::Matrix::getstring3m(Matrix3d m, std::string & s)
    4.75 -   {
    4.76 -   s.clear();
    4.77 -   std::stringstream ss;
    4.78 -   ss << "[ "
    4.79 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << " ], "
    4.80 -	 "[ " << m[3] << ", " << m[4] << ", " << m[5] << " ], "
    4.81 -	 "[ " << m[6] << ", " << m[7] << ", " << m[8] << " ] ]";
    4.82 -   s = ss.str();
    4.83 -   return s;   
    4.84 -   }
    4.85 -
    4.86 -////////////////////////////////////////////////////////////////////////////////
    4.87 -std::string & arda::Matrix::getstring4m(Matrix4i m, std::string & s)
    4.88 -   {
    4.89 -   s.clear();
    4.90 -   std::stringstream ss;
    4.91 -   ss << "[ "
    4.92 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
    4.93 -	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
    4.94 -	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
    4.95 -	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
    4.96 -   s = ss.str();
    4.97 -   return s;   
    4.98 -   }
    4.99 -
   4.100 -////////////////////////////////////////////////////////////////////////////////
   4.101 -std::string & arda::Matrix::getstring4m(Matrix4f m, std::string & s)
   4.102 -   {
   4.103 -   s.clear();
   4.104 -   std::stringstream ss;
   4.105 -   ss << "[ "
   4.106 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
   4.107 -	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
   4.108 -	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
   4.109 -	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
   4.110 -   s = ss.str();
   4.111 -   return s;   
   4.112 -   }
   4.113 -
   4.114 -////////////////////////////////////////////////////////////////////////////////
   4.115 -std::string & arda::Matrix::getstring4m(Matrix4d m, std::string & s)
   4.116 -   {
   4.117 -   s.clear();
   4.118 -   std::stringstream ss;
   4.119 -   ss << "[ "
   4.120 -	 "[ " << m[0] << ", " << m[1] << ", " << m[2] << ", " << m[3] << " ], "
   4.121 -	 "[ " << m[4] << ", " << m[5] << ", " << m[6] << ", " << m[7] <<  " ], "
   4.122 -	 "[ " << m[8] << ", " << m[9] << ", " << m[10] << ", " << m[11] <<  " ], "
   4.123 -	 "[ " << m[12] << ", " << m[13] << ", " << m[14] << ", " << m[15] <<  " ] ]";
   4.124 -   s = ss.str();
   4.125 -   return s;   
   4.126 -   }
   4.127 -
   4.128 -////////////////////////////////////////////////////////////////////////////////
   4.129 -int * arda::Matrix::multvec2mv(Matrix2i m, Vector2i v, Vector2i vres)
   4.130 -   {
   4.131 -   vres[0] = m[0]*v[0] + m[2]*v[1];
   4.132 -   vres[1] = m[1]*v[0] + m[3]*v[1];
   4.133 -   return vres;
   4.134 -   }
   4.135 -
   4.136 -////////////////////////////////////////////////////////////////////////////////
   4.137 -float * arda::Matrix::multvec2mv(Matrix2f m, Vector2f v, Vector2f vres)
   4.138 -   {
   4.139 -   vres[0] = m[0]*v[0] + m[2]*v[1];
   4.140 -   vres[1] = m[1]*v[0] + m[3]*v[1];
   4.141 -   return vres;   
   4.142 -   }
   4.143 -
   4.144 -////////////////////////////////////////////////////////////////////////////////
   4.145 -double * arda::Matrix::multvec2mv(Matrix2d m, Vector2d v, Vector2d vres)
   4.146 -   {
   4.147 -   vres[0] = m[0]*v[0] + m[2]*v[1];
   4.148 -   vres[1] = m[1]*v[0] + m[3]*v[1];
   4.149 -   return vres;
   4.150 -   }
   4.151 -
   4.152 -////////////////////////////////////////////////////////////////////////////////
   4.153 -int * arda::Matrix::multvec3mv(Matrix3i m, Vector3i v, Vector3i vres)
   4.154 -   {
   4.155 -   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
   4.156 -   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
   4.157 -   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
   4.158 -   return vres;
   4.159 -   }
   4.160 -
   4.161 -////////////////////////////////////////////////////////////////////////////////
   4.162 -float * arda::Matrix::multvec3mv(Matrix3f m, Vector3f v, Vector3f vres)
   4.163 -   {
   4.164 -   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
   4.165 -   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
   4.166 -   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
   4.167 -   return vres;
   4.168 -   }
   4.169 -
   4.170 -////////////////////////////////////////////////////////////////////////////////
   4.171 -double * arda::Matrix::multvec3mv(Matrix3d m, Vector3d v, Vector3d vres)
   4.172 -   {
   4.173 -   vres[0] = m[0]*v[0] + m[3]*v[1] + m[6]*v[2];
   4.174 -   vres[1] = m[1]*v[0] + m[4]*v[1] + m[7]*v[2];
   4.175 -   vres[2] = m[2]*v[0] + m[5]*v[1] + m[8]*v[2];
   4.176 -   return vres;
   4.177 -   }
   4.178 -
   4.179 -////////////////////////////////////////////////////////////////////////////////
   4.180 -int * arda::Matrix::multvec4mv(Matrix4i m, Vector4i v, Vector4i vres)
   4.181 -   {
   4.182 -   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
   4.183 -   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
   4.184 -   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
   4.185 -   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
   4.186 -   return vres;
   4.187 -   }
   4.188 -
   4.189 -////////////////////////////////////////////////////////////////////////////////
   4.190 -float * arda::Matrix::multvec4mv(Matrix4f m, Vector4f v, Vector4f vres)
   4.191 -   {
   4.192 -   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
   4.193 -   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
   4.194 -   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
   4.195 -   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
   4.196 -   return vres;
   4.197 -  
   4.198 -   }
   4.199 -
   4.200 -////////////////////////////////////////////////////////////////////////////////
   4.201 -double * arda::Matrix::multvec4mv(Matrix4d m, Vector4d v, Vector4d vres)
   4.202 -   {
   4.203 -   vres[0] = m[0]*v[0] + m[4]*v[1] + m[8] *v[2] + m[12]*v[3];
   4.204 -   vres[1] = m[1]*v[0] + m[5]*v[1] + m[9]*v[2]  + m[13]*v[3];
   4.205 -   vres[2] = m[2]*v[0] + m[6]*v[1] + m[10]*v[2] + m[14]*v[3];
   4.206 -   vres[3] = m[3]*v[0] + m[7]*v[1] + m[11]*v[2] + m[15]*v[3];
   4.207 -   return vres;
   4.208 -   }
   4.209 -
   4.210 -////////////////////////////////////////////////////////////////////////////////
   4.211 -int * arda::Matrix::multvec2vm(Vector2i v, Matrix2i m, Vector2i vres)
   4.212 -   {
   4.213 -   vres[0] = v[0]*m[0] + v[1]*m[1];
   4.214 -   vres[1] = v[0]*m[2] + v[1]*m[3];
   4.215 -   return vres;   
   4.216 -   }
   4.217 -
   4.218 -////////////////////////////////////////////////////////////////////////////////
   4.219 -float * arda::Matrix::multvec2vm(Vector2f v, Matrix2f m, Vector2f vres)
   4.220 -   {
   4.221 -   vres[0] = v[0]*m[0] + v[1]*m[1];
   4.222 -   vres[1] = v[0]*m[2] + v[1]*m[3];
   4.223 -   return vres;   
   4.224 -   }
   4.225 -
   4.226 -////////////////////////////////////////////////////////////////////////////////
   4.227 -double * arda::Matrix::multvec2vm(Vector2d v, Matrix2d m, Vector2d vres)
   4.228 -   {
   4.229 -   vres[0] = v[0]*m[0] + v[1]*m[1];
   4.230 -   vres[1] = v[0]*m[2] + v[1]*m[3];
   4.231 -   return vres;   
   4.232 -   }
   4.233 -
   4.234 -////////////////////////////////////////////////////////////////////////////////
   4.235 -int * arda::Matrix::multvec3vm(Vector3i v, Matrix3i m, Vector3i vres)
   4.236 -   {
   4.237 -   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
   4.238 -   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
   4.239 -   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
   4.240 -   return vres;   
   4.241 -   }
   4.242 -
   4.243 -////////////////////////////////////////////////////////////////////////////////
   4.244 -float * arda::Matrix::multvec3vm(Vector3f v, Matrix3f m, Vector3f vres)
   4.245 -   {
   4.246 -   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
   4.247 -   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
   4.248 -   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
   4.249 -   return vres;   
   4.250 -   }
   4.251 -
   4.252 -////////////////////////////////////////////////////////////////////////////////
   4.253 -double * arda::Matrix::multvec3vm(Vector3d v, Matrix3d m, Vector3d vres)
   4.254 -   {
   4.255 -   vres[0] = v[0]*m[0] + v[1]*m[1] + v[2]*m[2];
   4.256 -   vres[1] = v[0]*m[3] + v[1]*m[4] + v[2]*m[5];
   4.257 -   vres[2] = v[0]*m[6] + v[1]*m[7] + v[2]*m[8];
   4.258 -   return vres;   
   4.259 -   }
   4.260 -
   4.261 -////////////////////////////////////////////////////////////////////////////////
   4.262 -int * arda::Matrix::multvec4vm(Vector4i v, Matrix4i m, Vector4i vres)
   4.263 -   {
   4.264 -   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
   4.265 -   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
   4.266 -   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
   4.267 -   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
   4.268 -   return vres;   
   4.269 -   }
   4.270 -
   4.271 -////////////////////////////////////////////////////////////////////////////////
   4.272 -float * arda::Matrix::multvec4vm(Vector4f v, Matrix4f m, Vector4f vres)
   4.273 -   {
   4.274 -   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
   4.275 -   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
   4.276 -   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
   4.277 -   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
   4.278 -   return vres;   
   4.279 -   }
   4.280 -
   4.281 -////////////////////////////////////////////////////////////////////////////////
   4.282 -double * arda::Matrix::multvec4vm(Vector4d v, Matrix4d m, Vector4d vres)
   4.283 -   {
   4.284 -   vres[0] = v[0]*m[0]  + v[1]*m[1]  + v[2]*m[2]  + v[3]*m[3];
   4.285 -   vres[1] = v[0]*m[4]  + v[1]*m[5]  + v[2]*m[6]  + v[3]*m[7];
   4.286 -   vres[2] = v[0]*m[8]  + v[1]*m[9]  + v[2]*m[10] + v[3]*m[11];
   4.287 -   vres[3] = v[0]*m[12] + v[1]*m[13] + v[2]*m[14] + v[3]*m[15];
   4.288 -   return vres;   
   4.289 -   }
   4.290 -
   4.291 -
   4.292 -////////////////////////////////////////////////////////////////////////////////
   4.293 -int * arda::Matrix::multmat2(Matrix2i m1, Matrix2i m2, Matrix2i mres)
   4.294 -   {
   4.295 -   const int size=2;
   4.296 -   int i, j, k;
   4.297 -   for (i=0; i<size*size; ++i)
   4.298 -      mres[i] = 0;
   4.299 -   for (i=0; i<size; ++i)
   4.300 -      for (j=0; j<size; ++j)
   4.301 -	 for (k=0; k<size; ++k)
   4.302 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.303 -   return mres;
   4.304 -   }
   4.305 -
   4.306 -////////////////////////////////////////////////////////////////////////////////
   4.307 -float * arda::Matrix::multmat2(Matrix2f m1, Matrix2f m2, Matrix2f mres)
   4.308 -   {
   4.309 -   const int size=2;
   4.310 -   int i, j, k;
   4.311 -   for (i=0; i<size*size; ++i)
   4.312 -      mres[i] = 0;
   4.313 -   for (i=0; i<size; ++i)
   4.314 -      for (j=0; j<size; ++j)
   4.315 -	 for (k=0; k<size; ++k)
   4.316 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.317 -   return mres;
   4.318 -   }
   4.319 -
   4.320 -////////////////////////////////////////////////////////////////////////////////
   4.321 -double * arda::Matrix::multmat2(Matrix2d m1, Matrix2d m2, Matrix2d mres)
   4.322 -   {
   4.323 -   const int size=2;
   4.324 -   int i, j, k;
   4.325 -   for (i=0; i<size*size; ++i)
   4.326 -      mres[i] = 0;
   4.327 -   for (i=0; i<size; ++i)
   4.328 -      for (j=0; j<size; ++j)
   4.329 -	 for (k=0; k<size; ++k)
   4.330 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.331 -   return mres;
   4.332 -   }
   4.333 -
   4.334 -////////////////////////////////////////////////////////////////////////////////
   4.335 -int * arda::Matrix::multmat3(Matrix3i m1, Matrix3i m2, Matrix3i mres)
   4.336 -   {
   4.337 -   const int size=3;
   4.338 -   int i, j, k;
   4.339 -   for (i=0; i<size*size; ++i)
   4.340 -      mres[i] = 0;
   4.341 -   for (i=0; i<size; ++i)
   4.342 -      for (j=0; j<size; ++j)
   4.343 -	 for (k=0; k<size; ++k)
   4.344 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.345 -   return mres;
   4.346 -   }
   4.347 -
   4.348 -////////////////////////////////////////////////////////////////////////////////
   4.349 -float * arda::Matrix::multmat3(Matrix3f m1, Matrix3f m2, Matrix3f mres)
   4.350 -   {
   4.351 -   const int size=3;
   4.352 -   int i, j, k;
   4.353 -   for (i=0; i<size*size; ++i)
   4.354 -      mres[i] = 0;
   4.355 -   for (i=0; i<size; ++i)
   4.356 -      for (j=0; j<size; ++j)
   4.357 -	 for (k=0; k<size; ++k)
   4.358 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.359 -   return mres;
   4.360 -   }
   4.361 -
   4.362 -////////////////////////////////////////////////////////////////////////////////
   4.363 -double * arda::Matrix::multmat3(Matrix3d m1, Matrix3d m2, Matrix3d mres)
   4.364 -   {
   4.365 -   const int size=3;
   4.366 -   int i, j, k;
   4.367 -   for (i=0; i<size*size; ++i)
   4.368 -      mres[i] = 0;
   4.369 -   for (i=0; i<size; ++i)
   4.370 -      for (j=0; j<size; ++j)
   4.371 -	 for (k=0; k<size; ++k)
   4.372 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.373 -   return mres;
   4.374 -   }
   4.375 -
   4.376 -////////////////////////////////////////////////////////////////////////////////
   4.377 -int * arda::Matrix::multmat4(Matrix4i m1, Matrix4i m2, Matrix4i mres)
   4.378 -   {
   4.379 -   const int size=4;
   4.380 -   int i, j, k;
   4.381 -   for (i=0; i<size*size; ++i)
   4.382 -      mres[i] = 0;
   4.383 -   for (i=0; i<size; ++i)
   4.384 -      for (j=0; j<size; ++j)
   4.385 -	 for (k=0; k<size; ++k)
   4.386 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.387 -   return mres;
   4.388 -   }
   4.389 -
   4.390 -////////////////////////////////////////////////////////////////////////////////
   4.391 -float * arda::Matrix::multmat4(Matrix4f m1, Matrix4f m2, Matrix4f mres)
   4.392 -   {
   4.393 -   const int size=4;
   4.394 -   int i, j, k;
   4.395 -   for (i=0; i<size*size; ++i)
   4.396 -      mres[i] = 0;
   4.397 -   for (i=0; i<size; ++i)
   4.398 -      for (j=0; j<size; ++j)
   4.399 -	 for (k=0; k<size; ++k)
   4.400 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.401 -   return mres;
   4.402 -   }
   4.403 -
   4.404 -////////////////////////////////////////////////////////////////////////////////
   4.405 -double * arda::Matrix::multmat4(Matrix4d m1, Matrix4d m2, Matrix4d mres)
   4.406 -   {
   4.407 -   cout << "here" << endl;
   4.408 -   const int size=4;
   4.409 -   int i, j, k;
   4.410 -   for (i=0; i<size*size; ++i)
   4.411 -      mres[i] = 0;
   4.412 -   for (i=0; i<size; ++i)
   4.413 -      for (j=0; j<size; ++j)
   4.414 -	 for (k=0; k<size; ++k)
   4.415 -	    mres[size*i+j] += m1[size*k+j] * m2[size*i+k];
   4.416 -   return mres;
   4.417 -   }
   4.418 -
   4.419 -////////////////////////////////////////////////////////////////////////////////
   4.420 -int * arda::Matrix::transpose2(Matrix2i m, Matrix2i mres)
   4.421 -   {
   4.422 -   const int size=2;
   4.423 -   int i, j;
   4.424 -   for (i=0; i<size; ++i)
   4.425 -      for (j=0; j<size; j++)
   4.426 -	 mres[size*i+j] = m[size*j+i];
   4.427 -   return mres;
   4.428 -   }
   4.429 -
   4.430 -////////////////////////////////////////////////////////////////////////////////
   4.431 -float * arda::Matrix::transpose2(Matrix2f m, Matrix2f mres)
   4.432 -   {
   4.433 -   const int size=2;
   4.434 -   int i, j;
   4.435 -   for (i=0; i<size; ++i)
   4.436 -      for (j=0; j<size; j++)
   4.437 -	 mres[size*i+j] = m[size*j+i];
   4.438 -   return mres;
   4.439 -   }
   4.440 -
   4.441 -////////////////////////////////////////////////////////////////////////////////
   4.442 -double * arda::Matrix::transpose2(Matrix2d m, Matrix2d mres)
   4.443 -   {
   4.444 -   const int size=2;
   4.445 -   int i, j;
   4.446 -   for (i=0; i<size; ++i)
   4.447 -      for (j=0; j<size; j++)
   4.448 -	 mres[size*i+j] = m[size*j+i];
   4.449 -   return mres;
   4.450 -   }
   4.451 -
   4.452 -////////////////////////////////////////////////////////////////////////////////
   4.453 -int * arda::Matrix::transpose3(Matrix3i m, Matrix3i mres)
   4.454 -   {
   4.455 -   const int size=3;
   4.456 -   int i, j;
   4.457 -   for (i=0; i<size; ++i)
   4.458 -      for (j=0; j<size; j++)
   4.459 -	 mres[size*i+j] = m[size*j+i];
   4.460 -   return mres;
   4.461 -   }
   4.462 -
   4.463 -////////////////////////////////////////////////////////////////////////////////
   4.464 -float * arda::Matrix::transpose3(Matrix3f m, Matrix3f mres)
   4.465 -   {
   4.466 -   const int size=3;
   4.467 -   int i, j;
   4.468 -   for (i=0; i<size; ++i)
   4.469 -      for (j=0; j<size; j++)
   4.470 -	 mres[size*i+j] = m[size*j+i];
   4.471 -   return mres;
   4.472 -   }
   4.473 -
   4.474 -////////////////////////////////////////////////////////////////////////////////
   4.475 -double * arda::Matrix::transpose3(Matrix3d m, Matrix3d mres)
   4.476 -   {
   4.477 -   const int size=3;
   4.478 -   int i, j;
   4.479 -   for (i=0; i<size; ++i)
   4.480 -      for (j=0; j<size; j++)
   4.481 -	 mres[size*i+j] = m[size*j+i];
   4.482 -   return mres;
   4.483 -   }
   4.484 -
   4.485 -////////////////////////////////////////////////////////////////////////////////
   4.486 -int * arda::Matrix::transpose4(Matrix4i m, Matrix4i mres)
   4.487 -   {
   4.488 -   const int size=4;
   4.489 -   int i, j;
   4.490 -   for (i=0; i<size; ++i)
   4.491 -      for (j=0; j<size; j++)
   4.492 -	 mres[size*i+j] = m[size*j+i];
   4.493 -   return mres;
   4.494 -   }
   4.495 -
   4.496 -////////////////////////////////////////////////////////////////////////////////
   4.497 -float * arda::Matrix::transpose4(Matrix4f m, Matrix4f mres)
   4.498 -   {
   4.499 -   const int size=4;
   4.500 -   int i, j;
   4.501 -   for (i=0; i<size; ++i)
   4.502 -      for (j=0; j<size; j++)
   4.503 -	 mres[size*i+j] = m[size*j+i];
   4.504 -   return mres;
   4.505 -   }
   4.506 -
   4.507 -////////////////////////////////////////////////////////////////////////////////
   4.508 -double * arda::Matrix::transpose4(Matrix4d m, Matrix4d mres)
   4.509 -   {
   4.510 -   const int size=4;
   4.511 -   int i, j;
   4.512 -   for (i=0; i<size; ++i)
   4.513 -      for (j=0; j<size; j++)
   4.514 -	 mres[size*i+j] = m[size*j+i];
   4.515 -   return mres;
   4.516 -   }
   4.517 -
   4.518 -
     5.1 --- a/src/main.cpp	Fri Oct 07 11:09:15 2011 -0500
     5.2 +++ b/src/main.cpp	Fri Oct 07 13:33:48 2011 -0500
     5.3 @@ -6,7 +6,7 @@
     5.4  using namespace std;
     5.5  
     5.6  #include "Math.h"
     5.7 -using namespace arda;
     5.8 +using namespace arda::Math;
     5.9  
    5.10  ////////////////////////////////////////////////////////////////////////////////
    5.11  void test_vector(void)
    5.12 @@ -588,8 +588,8 @@
    5.13     double d;
    5.14  
    5.15     cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
    5.16 -   d = m22i_1.det();
    5.17 -   cout << setw(40) << "d = m22i_1.det(): " << d << endl;
    5.18 +   d = det(m22i_1);
    5.19 +   cout << setw(40) << "d = det(m22i_1): " << d << endl;
    5.20     }
    5.21  
    5.22     {
    5.23 @@ -599,8 +599,8 @@
    5.24     Matrix22i m22i_2;
    5.25  
    5.26     cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
    5.27 -   m22i_2 = m22i_1.transpose();
    5.28 -   cout << setw(40) << "m22i_2 = m22i_1.transpose(): " << m22i_2.getstring(s1) << endl;
    5.29 +   m22i_2 = transpose(m22i_1);
    5.30 +   cout << setw(40) << "m22i_2 = transpose(m22i_1): " << m22i_2.getstring(s1) << endl;
    5.31     }
    5.32  
    5.33     {
    5.34 @@ -612,6 +612,55 @@
    5.35     m22i_1.setidentity();
    5.36     cout << setw(40) << "m22i_1.setidentity(): " << m22i_1.getstring(s1) << endl;
    5.37     }
    5.38 +
    5.39 +   {
    5.40 +   cout << "===============================================================" << endl <<
    5.41 +	 "Testing Matrix getrow() and getcol()" << endl;
    5.42 +   Matrix22i m22i_1(1, 2, 3, 4);
    5.43 +   Vector2i v2i_1(0);
    5.44 +
    5.45 +   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
    5.46 +   v2i_1 = m22i_1.getrow(0);
    5.47 +   cout << setw(40) << "v2i_1 = m22i_1.getrow(0): " << v2i_1.getstring(s1) << endl;
    5.48 +   v2i_1 = m22i_1.getrow(1);
    5.49 +   cout << setw(40) << "v2i_1 = m22i_1.getrow(1): " << v2i_1.getstring(s1) << endl;
    5.50 +   v2i_1 = m22i_1.getcol(0);
    5.51 +   cout << setw(40) << "v2i_1 = m22i_1.getcol(0): " << v2i_1.getstring(s1) << endl;
    5.52 +   v2i_1 = m22i_1.getcol(1);
    5.53 +   cout << setw(40) << "v2i_1 = m22i_1.getcol(1): " << v2i_1.getstring(s1) << endl;
    5.54 +   }
    5.55 +
    5.56 +   {
    5.57 +   cout << "===============================================================" << endl <<
    5.58 +	 "Testing Matrix setrow() and setcol()" << endl;
    5.59 +   Matrix22i m22i_1;
    5.60 +   Vector2i v2i_1(2, 3);
    5.61 +
    5.62 +   m22i_1.setidentity();
    5.63 +   cout << setw(40) << "m22i_1: " << m22i_1.getstring(s1) << endl;
    5.64 +   m22i_1.setrow(0, v2i_1);
    5.65 +   cout << setw(40) << "m22i_1.setrow(0, v2i_1): " << m22i_1.getstring(s1) << endl;
    5.66 +
    5.67 +   m22i_1.setidentity();
    5.68 +   m22i_1.setrow(1, v2i_1);
    5.69 +   cout << setw(40) << "m22i_1.setrow(1, v2i_1): " << m22i_1.getstring(s1) << endl;
    5.70 +
    5.71 +   m22i_1.setidentity();
    5.72 +   m22i_1.setrow(1, 4, 5);
    5.73 +   cout << setw(40) << "m22i_1.setrow(1, 4, 5): " << m22i_1.getstring(s1) << endl;
    5.74 +
    5.75 +   m22i_1.setidentity();
    5.76 +   m22i_1.setcol(0, v2i_1);
    5.77 +   cout << setw(40) << "m22i_1.setcol(0, v2i_1): " << m22i_1.getstring(s1) << endl;
    5.78 +
    5.79 +   m22i_1.setidentity();
    5.80 +   m22i_1.setcol(1, v2i_1);
    5.81 +   cout << setw(40) << "m22i_1.setcol(1, v2i_1): " << m22i_1.getstring(s1) << endl;
    5.82 +
    5.83 +   m22i_1.setidentity();
    5.84 +   m22i_1.setcol(1, 4, 5);
    5.85 +   cout << setw(40) << "m22i_1.setcol(1, 4, 5): " << m22i_1.getstring(s1) << endl;
    5.86 +   }
    5.87     }
    5.88  
    5.89  ////////////////////////////////////////////////////////////////////////////////