//*************************************************************************////  File Name: Matrix.h //  Author   : Justin Legakis (Original Author)//			   Ali BaderEddin (Applied some changes + added useful comments)//  //  Description: A class that represents a matrix and supports many matrix//               operations.//  //       0 1 2 3 //    0 [1 0 0 5]//    1 [0 1 0 0]//    2 [0 0 1 0]  ==> matrix [0][3] is equal to 5 					//    3 [0 0 0 1]      (that's how matrix is represened)////*************************************************************************#ifndef _MATRIX_H_#define _MATRIX_H_#include <math.h>#include <assert.h>#include "vectors.h"class Matrix {	private:	  // REPRESENTATION	  float	data[4][4];	public:	  // CONSTRUCTORS & DESTRUCTOR	  Matrix() { Clear(); }	  Matrix(const Matrix& A) 	  {		for (int y = 0; y < 4; y++) 		  for (int x = 0; x < 4; x++) 			data[y][x] = A.data[y][x]; 	  }	  ~Matrix() {}  	  // ACCESSOR	  float Get (int x, int y) const 	  { 		assert (x >= 0 && x < 4);		assert (y >= 0 && y < 4);		return data[y][x]; 	  }  	  // MODIFIERS	  void Set(int x, int y, float v) 	  {		assert (x >= 0 && x < 4);		assert (y >= 0 && y < 4);		data[y][x] = v; 	  }	  void SetToIdentity() 	  {		for (int y = 0; y < 4; y++)		  for (int x = 0; x < 4; x++)			data[y][x] = (x == y); 	  }	  void Clear() 	  {		for (int y = 0; y < 4; y++)		  for (int x = 0; x < 4; x++)			data[y][x] = 0; 	  }  	  // OVERLOADED OPERATORS	  Matrix& operator=(const Matrix& A) 	  {		for (int y=0; y<4; y++) 		  for (int x=0; x<4; x++) 			data[y][x] = A.data[y][x]; 		return (*this); 	  }	  int operator==(const Matrix& A) const 	  {		for (int y=0; y<4; y++) 		  for (int x=0; x<4; x++) 			if (this->data[y][x] != A.data[y][x]) 				return 0; 		return 1; 	  }	  int operator!=(const Matrix &A) const { return !(*this==A); }  	  friend Matrix operator+(const Matrix &A, const Matrix &B);	  friend Matrix operator-(const Matrix &A, const Matrix &B);	  friend Matrix operator*(const Matrix &A, const Matrix &B);	  friend Matrix operator*(const Matrix &A, float f);	  friend Matrix operator*(float f, const Matrix &A) { return A * f; }  	  Matrix& operator+=(const Matrix& A) 	  {		*this = *this + A;		return *this; 	  }	  Matrix& operator-=(const Matrix& A) 	  {		*this = *this - A;		return *this; 	  }	  Matrix& operator*=(const float d) 	  {		*this = *this * d;		return *this; 	  }	  Matrix& operator*=(const Matrix& A) 	  {		*this = *this * A;		return *this; 	  }	  // TRANSFORMATIONS	  void Translate(float x, float y, float z);	  void Scale(float x, float y, float z);	  void Scale(float s) { Scale(s, s, s); }	  void XRotate(float theta);	  void YRotate(float theta);	  void ZRotate(float theta);	  void Transform(Vec4f &v);	  void Transform(Vec3f &v) 	  {		Vec4f v2 = Vec4f(v.x(),v.y(),v.z(),1);		Transform(v2);		v.Set(v2.x(),v2.y(),v2.z()); 	  }	  void Transform(Vec2f &v) 	  {		Vec4f v2 = Vec4f(v.x(),v.y(),1,1);		Transform(v2);		v.Set(v2.x(),v2.y()); 	  }	  // INPUT / OUTPUT (Edited by Ali BaderEddin)	  void Write(FILE *F = stdout) const 	  {		assert (F != NULL);		for (int y = 0; y < 4; y++) 		{		  for (int x = 0; x < 4; x++) 			fprintf (F, "%.2f ", data[y][x]); 		  fprintf (F,"\n"); 		} 		fprintf (F,"\n"); 	  }	  void Write3x3(FILE *F = stdout) const 	  {		assert (F != NULL);		for (int y = 0; y < 4; y++) 		{		  if (y == 2) continue;		  for (int x = 0; x < 4; x++) 		  {			if (x == 2) continue;			fprintf (F, "%.2f ", data[y][x]); 		  }		  fprintf (F,"\n"); 		} 		fprintf (F,"\n"); 	  }	  void Read(FILE *F) 	  {		assert (F != NULL);		for (int y = 0; y < 4; y++) 		  for (int x = 0; x < 4; x++) 		  {			int scanned = fscanf (F,"%f",&data[y][x]);			assert (scanned == 1); 		  } 	  }	  void Read3x3(FILE *F) 	  {		assert (F != NULL);		Clear();		for (int y = 0; y < 4; y++) 		{		  if (y == 2) continue;		  for (int x = 0; x < 4; x++) 		  {			if (x == 2) continue;			int scanned = fscanf (F,"%f",&data[y][x]);			assert (scanned == 1); 		  } 		} 	  }	  	  /***** Stuff Added By Ali BaderEddin *****/	  //  Loading OpenGL matrix into our Matrix	  void LoadOpenGLMatrix (float *mat)	  {		  for (int i = 0; i < 16; i++)			  this->data[i%4][i/4] = mat;	  }	  //  Saving our matrix into an OpenGL Matrix	  float * GetOpenGLMatrix ()	  {		  float * mat;		  		  mat = (float *) malloc (16 * sizeof (float));		  for (int i = 0; i < 16; i++)			  mat = this->data[i%4][i/4];		  return mat;	  }};#endif

//*************************************************************************////  File Name: Matrix.cpp //  Author   : Justin Legakis (Original Author)//			   Ali BaderEddin (Applied some changes + added useful comments)//  //  Description: A class that represents a matrix and supports many matrix//               operations.//  //       0 1 2 3 //    0 [1 0 0 5]//    1 [0 1 0 0]//    2 [0 0 1 0]  ==> matrix [0][3] is equal to 5 					//    3 [0 0 0 1]      (that's how matrix is represened)////*************************************************************************/#include <stdlib.h>#include <stdio.h>#include <math.h>#include <string.h>#include "matrix.h"#include "vectors.h"// OVERLOADED OPERATORSMatrix operator+(const Matrix& A, const Matrix& B) {  Matrix answer;  for (int y=0; y<4; y++)     for (int x=0; x<4; x++)       answer.data[y][x] = A.data[y][x] + B.data[y][x];  return answer; }Matrix operator-(const Matrix& A, const Matrix& B) {  Matrix answer;  for (int y=0; y<4; y++)     for (int x=0; x<4; x++)       answer.data[y][x] = A.data[y][x] - B.data[y][x];  return answer; }Matrix operator*(const Matrix& A, const Matrix& B) {  Matrix answer;  for (int y=0; y<4; y++)     for (int x=0; x<4; x++)       for (int i=0; i<4; i++) 		answer.data[y][x] += A.data[y] * B.data[x];  return answer;}Matrix operator*(const Matrix& A, float f) {  Matrix answer;  for (int y=0; y<4; y++)     for (int x=0; x<4; x++)       answer.data[y][x] = A.data[y][x] * f;  return answer;}// TRANSFORMATIONS/*  Translation Matrix       0 1 2  3     0 [1 0 0 t_x]    1 [0 1 0 t_y]    2 [0 0 1 t_z]  					    3 [0 0 0   1]*/void Matrix::Translate(float t_x, float t_y, float t_z) {  Matrix t;  t.SetToIdentity();  t.data[0][3] = t_x;  t.data[1][3] = t_y;  t.data[2][3] = t_z;  t *= *this;  *this = t;}/*  Scaling Matrix       0   1   2   3     0 [s_z 0   0   0]    1 [0   s_y 0   0]    2 [0   0   s_z 0]  					    3 [0   0   0   1]*/void Matrix::Scale(float s_x, float s_y, float s_z) {  Matrix s;   s.SetToIdentity();  s.data[0][0] = s_x;  s.data[1][1] = s_y;  s.data[2][2] = s_z;  s.data[3][3] = 1;  s *= *this;  *this = s;}/*  X Rotation Matrix       0 1   2    3     0 [1 0   0    0]    1 [0 cos -sin 0]    2 [0 sin cos  0]  					    3 [0 0   0    1]*/void Matrix::XRotate(float theta) {  Matrix rx;  rx.SetToIdentity();  rx.data[1][1]= (float)cos((float)theta);  rx.data[1][2]=-(float)sin((float)theta);  rx.data[2][1]= (float)sin((float)theta);  rx.data[2][2]= (float)cos((float)theta);  rx *= *this;  *this = rx;}/*  Y Rotation Matrix       0    1   2    3     0 [cos  0   sin  0]    1 [0    1   0    0]    2 [-sin 0   cos  0]  					    3 [0    0   0    1]*/void Matrix::YRotate(float theta) {  Matrix ry;  ry.SetToIdentity();  ry.data[0][0]= (float)cos((float)theta);  ry.data[0][2]= (float)sin((float)theta);  ry.data[2][0]=-(float)sin((float)theta);  ry.data[2][2]= (float)cos((float)theta);  ry *= *this;  *this = ry;}/*  Z Rotation Matrix       0     1  2 3     0 [cos -sin 0 0]    1 [sin cos  0 0]    2 [0   0    1 0]  					    3 [0   0    0 1]*/void Matrix::ZRotate(float theta) {  Matrix rz;  rz.SetToIdentity();  rz.data[0][0]= (float)cos((float)theta);  rz.data[0][1]=-(float)sin((float)theta);  rz.data[1][0]= (float)sin((float)theta);  rz.data[1][1]= (float)cos((float)theta);  rz *= *this;  *this = rz;}void Matrix::Transform(Vec4f &v) {  Vec4f answer;  for (int y=0; y<4; y++)   {    answer.data[y] = 0;    for (int i=0; i<4; i++)       answer.data[y] += data[y] * v;  }  v = answer;}

//*************************************************************************////  File Name: Vector.h //  Author   : Justin Legakis (Original Author)//			   Ali BaderEddin (Applied some changes)//  //  Description: A class that represents a 2d, 3d, and 4d vector //  			 and supports many vector operations.//  //*************************************************************************/#ifndef _VECTORS_H_#define _VECTORS_H_#include <stdlib.h>#include <stdio.h>#include <math.h>#include <assert.h>class Matrix;class Vec2f {	private:	  // REPRESENTATION	  float		data[2];	public:	  // CONSTRUCTORS & DESTRUCTOR	  Vec2f() { data[0] = data[1] = 0; }	  Vec2f(const Vec2f &V) 	  {		data[0] = V.data[0];		data[1] = V.data[1]; 	  }	  Vec2f(float d0, float d1) 	  {		data[0] = d0;		data[1] = d1; 	  }	  Vec2f(const Vec2f &V1, const Vec2f &V2) 	  {		data[0] = V1.data[0] - V2.data[0];		data[1] = V1.data[1] - V2.data[1]; 	  }	  ~Vec2f() { }	  // ACCESSORS	  void Get(float &d0, float &d1) const 	  {		d0 = data[0];		d1 = data[1]; 	  }	  float operator[](int i) const 	  { 		assert (i >= 0 && i < 2); 		return data; 	  }	  float x() const { return data[0]; }	  float y() const { return data[1]; }	  float Length() const 	  {		float l = (float)sqrt( data[0] * data[0] + data[1] * data[1] );		return l; 	  }	  // MODIFIERS	  void Set(float d0, float d1) 	  {		data[0] = d0;		data[1] = d1; 	  }	  void Scale(float d0, float d1) 	  {		data[0] *= d0;		data[1] *= d1; 	  }	  void Divide(float d0, float d1) 	  {		data[0] /= d0;		data[1] /= d1; 	  }	  void Negate() 	  {		data[0] = -data[0];		data[1] = -data[1]; 	  }	  // OVERLOADED OPERATORS	  Vec2f& operator=(const Vec2f &V) 	  {		data[0] = V.data[0];		data[1] = V.data[1];		return *this; 	  }	  int operator==(const Vec2f &V) const 	  {		return ((data[0] == V.data[0]) &&			(data[1] == V.data[1])); 	  }	  int operator!=(const Vec2f &V) 	  {		return ((data[0] != V.data[0]) ||			(data[1] != V.data[1])); 	  }	  Vec2f& operator+=(const Vec2f &V) 	  {		data[0] += V.data[0];		data[1] += V.data[1];		return *this; 	  }	  Vec2f& operator-=(const Vec2f &V) 	  {		data[0] -= V.data[0];		data[1] -= V.data[1];		return *this; 	  }	  Vec2f& operator*=(float f) 	  {		data[0] *= f;		data[1] *= f;		return *this; 	  }	  Vec2f& operator/=(float f) 	  {		data[0] /= f;		data[1] /= f;		return *this; 	  }  	  // OPERATIONS	  float Dot2(const Vec2f &V) const 	  {		return data[0] * V.data[0] + data[1] * V.data[1] ;	  }	  // STATIC OPERATIONS	  static void Add(Vec2f &a, const Vec2f &b, const Vec2f &c ) 	  {		a.data[0] = b.data[0] + c.data[0];		a.data[1] = b.data[1] + c.data[1]; 	  }	  static void Sub(Vec2f &a, const Vec2f &b, const Vec2f &c ) 	  {		a.data[0] = b.data[0] - c.data[0];		a.data[1] = b.data[1] - c.data[1]; 	  }	  static void CopyScale(Vec2f &a, const Vec2f &b, float c ) 	  {		a.data[0] = b.data[0] * c;		a.data[1] = b.data[1] * c; 	  }	  static void AddScale(Vec2f &a, const Vec2f &b, const Vec2f &c, float d ) 	  {		a.data[0] = b.data[0] + c.data[0] * d;		a.data[1] = b.data[1] + c.data[1] * d; 	  }	  static void Average(Vec2f &a, const Vec2f &b, const Vec2f &c ) 	  {		a.data[0] = (b.data[0] + c.data[0]) * 0.5f;		a.data[1] = (b.data[1] + c.data[1]) * 0.5f; 	  }	  static void WeightedSum(Vec2f &a, const Vec2f &b, float c, const Vec2f &d, float e ) 	  {		a.data[0] = b.data[0] * c + d.data[0] * e;		a.data[1] = b.data[1] * c + d.data[1] * e; 	  }	  // INPUT / OUTPUT	  void Write(FILE *F = stdout) 	  {		fprintf (F, "%f %f\n",data[0],data[1]); 	  }};// ====================================================================// ====================================================================class Vec3f {	private:	  // REPRESENTATION	  float		data[3];	public:	  // CONSTRUCTORS & DESTRUCTOR	  Vec3f() { data[0] = data[1] = data[2] = 0; }	  Vec3f(const Vec3f &V) 	  {		data[0] = V.data[0];		data[1] = V.data[1];		data[2] = V.data[2]; 	  }	  Vec3f(float d0, float d1, float d2) 	  {		data[0] = d0;		data[1] = d1;		data[2] = d2; 	  }	  Vec3f(const Vec3f &V1, const Vec3f &V2) 	  {		data[0] = V1.data[0] - V2.data[0];		data[1] = V1.data[1] - V2.data[1];		data[2] = V1.data[2] - V2.data[2]; 	  }	  ~Vec3f() { }	  // ACCESSORS	  void Get(float &d0, float &d1, float &d2) const 	  {		d0 = data[0];		d1 = data[1];		d2 = data[2]; 	  }	  float operator[](int i) const 	  { 		assert (i >= 0 && i < 3); 		return data; 	  }	  float x() const { return data[0]; }	  float y() const { return data[1]; }	  float z() const { return data[2]; }	  float r() const { return data[0]; }	  float g() const { return data[1]; }	  float b() const { return data[2]; }	  float Length() const 	  {		float l = (float)sqrt( data[0] * data[0] +				data[1] * data[1] +				data[2] * data[2] );		return l; 	  }	  // MODIFIERS	  void Set(float d0, float d1, float d2) 	  {		data[0] = d0;		data[1] = d1;		data[2] = d2; 	  }	  void Scale(float d0, float d1, float d2) 	  {		data[0] *= d0;		data[1] *= d1;		data[2] *= d2; 	  }	  void Divide(float d0, float d1, float d2) 	  {		data[0] /= d0;		data[1] /= d1;		data[2] /= d2; 	  }	  void Normalize() 	  {		float l = Length();		if (l > 0) 		{		  data[0] /= l;		  data[1] /= l;		  data[2] /= l; 		}	  }	  void Negate() 	  {		data[0] = -data[0];		data[1] = -data[1];		data[2] = -data[2]; 	  }	  // OVERLOADED OPERATORS	  Vec3f& operator=(const Vec3f &V) 	  {		data[0] = V.data[0];		data[1] = V.data[1];		data[2] = V.data[2];		return *this; 	  }	  int operator==(const Vec3f &V) 	  {		return ((data[0] == V.data[0]) &&			(data[1] == V.data[1]) &&			(data[2] == V.data[2])); 	  }	  int operator!=(const Vec3f &V) 	  {		return ((data[0] != V.data[0]) ||			(data[1] != V.data[1]) ||			(data[2] != V.data[2])); 	  }	  Vec3f& operator+=(const Vec3f &V) 	  {		data[0] += V.data[0];		data[1] += V.data[1];		data[2] += V.data[2];		return *this; 	  }	  Vec3f& operator-=(const Vec3f &V) 	  {		data[0] -= V.data[0];		data[1] -= V.data[1];		data[2] -= V.data[2];		return *this; 	  }	  Vec3f& operator*=(int i) 	  {		data[0] = float(data[0] * i);		data[1] = float(data[1] * i);		data[2] = float(data[2] * i);		return *this; 	  }	  Vec3f& operator*=(float f) 	  {		data[0] *= f;		data[1] *= f;		data[2] *= f;		return *this;	  }	  Vec3f& operator/=(int i) 	  {		data[0] = float(data[0] / i);		data[1] = float(data[1] / i);		data[2] = float(data[2] / i);		return *this; 	  }	  Vec3f& operator/=(float f) 	  {		data[0] /= f;		data[1] /= f;		data[2] /= f;		return *this; 	  }  	  // OPERATIONS	  float Dot3(const Vec3f &V) const 	  {		return data[0] * V.data[0] +		  data[1] * V.data[1] +		  data[2] * V.data[2] ; 	  }	  // STATIC OPERATIONS	  friend void Add(Vec3f &a, const Vec3f &b, const Vec3f &c ) 	  {		a.data[0] = b.data[0] + c.data[0];		a.data[1] = b.data[1] + c.data[1];		a.data[2] = b.data[2] + c.data[2]; 	  }	  friend void Sub(Vec3f &a, const Vec3f &b, const Vec3f &c ) 	  {		a.data[0] = b.data[0] - c.data[0];		a.data[1] = b.data[1] - c.data[1];		a.data[2] = b.data[2] - c.data[2]; 	  }	  friend void CopyScale(Vec3f &a, const Vec3f &b, float c ) 	  {		a.data[0] = b.data[0] * c;		a.data[1] = b.data[1] * c;		a.data[2] = b.data[2] * c; 	  }	  friend void AddScale(Vec3f &a, const Vec3f &b, const Vec3f &c, float d ) 	  {		a.data[0] = b.data[0] + c.data[0] * d;		a.data[1] = b.data[1] + c.data[1] * d;		a.data[2] = b.data[2] + c.data[2] * d; 	  }	  friend void Average(Vec3f &a, const Vec3f &b, const Vec3f &c ) 	  {		a.data[0] = (b.data[0] + c.data[0]) * 0.5f;		a.data[1] = (b.data[1] + c.data[1]) * 0.5f;		a.data[2] = (b.data[2] + c.data[2]) * 0.5f; 	  }	  friend void WeightedSum(Vec3f &a, const Vec3f &b, float c, const Vec3f &d, float e ) 	  {		a.data[0] = b.data[0] * c + d.data[0] * e;		a.data[1] = b.data[1] * c + d.data[1] * e;		a.data[2] = b.data[2] * c + d.data[2] * e; 	  }	  static void Cross3(Vec3f &c, const Vec3f &v1, const Vec3f &v2) 	  {		float x = v1.data[1]*v2.data[2] - v1.data[2]*v2.data[1];		float y = v1.data[2]*v2.data[0] - v1.data[0]*v2.data[2];		float z = v1.data[0]*v2.data[1] - v1.data[1]*v2.data[0];		c.data[0] = x; c.data[1] = y; c.data[2] = z; 	  }  	  // INPUT / OUTPUT	  void Write(FILE *F = stdout) 	  {		fprintf (F, "%f %f %f\n",data[0],data[1],data[2]); 	  }  };// ====================================================================// ====================================================================class Vec4f {private:  friend class Matrix;  // REPRESENTATION  float		data[4];public:  // CONSTRUCTORS & DESTRUCTOR  Vec4f() { data[0] = data[1] = data[2] = data[3] = 0; }  Vec4f(const Vec4f &V)   {    data[0] = V.data[0];    data[1] = V.data[1];    data[2] = V.data[2];    data[3] = V.data[3];   }  Vec4f(float d0, float d1, float d2, float d3)   {    data[0] = d0;    data[1] = d1;    data[2] = d2;    data[3] = d3;   }  Vec4f(const Vec3f &V, float w)   {    data[0] = V.x();    data[1] = V.y();    data[2] = V.z();    data[3] = w;   }  Vec4f(const Vec4f &V1, const Vec4f &V2)   {    data[0] = V1.data[0] - V2.data[0];    data[1] = V1.data[1] - V2.data[1];    data[2] = V1.data[2] - V2.data[2];    data[3] = V1.data[3] - V2.data[3];   }  ~Vec4f() { }  // ACCESSORS  void Get(float &d0, float &d1, float &d2, float &d3) const   {    d0 = data[0];    d1 = data[1];    d2 = data[2];    d3 = data[3];   }  float operator[](int i) const   {     assert (i >= 0 && i < 4);     return data;   }  float x() const { return data[0]; }  float y() const { return data[1]; }  float z() const { return data[2]; }  float w() const { return data[3]; }  float r() const { return data[0]; }  float g() const { return data[1]; }  float b() const { return data[2]; }  float a() const { return data[3]; }  float Length() const   {    float l = (float)sqrt( data[0] * data[0] +		     data[1] * data[1] +		     data[2] * data[2] +		     data[3] * data[3] );    return l;   }  // MODIFIERS  void Set(float d0, float d1, float d2, float d3)  {    data[0] = d0;    data[1] = d1;    data[2] = d2;    data[3] = d3;   }  void Scale(float d0, float d1, float d2, float d3)   {    data[0] *= d0;    data[1] *= d1;    data[2] *= d2;    data[3] *= d3;   }  void Divide(float d0, float d1, float d2, float d3)   {    data[0] /= d0;    data[1] /= d1;    data[2] /= d2;    data[3] /= d3;   }  void Negate()   {    data[0] = -data[0];    data[1] = -data[1];    data[2] = -data[2];    data[3] = -data[3];   }  void Normalize()   {    float l = Length();    if (l > 0) 	{      data[0] /= l;      data[1] /= l;      data[2] /= l; 	}  }  void DivideByW()   {    if (data[3] != 0) 	{      data[0] /= data[3];      data[1] /= data[3];      data[2] /= data[3];    } 	else 	{      data[0] = data[1] = data[2] = 0; 	}    data[3] = 1;  }  // OVERLOADED OPERATORS  Vec4f& operator=(const Vec4f &V)   {    data[0] = V.data[0];    data[1] = V.data[1];    data[2] = V.data[2];    data[3] = V.data[3];    return *this;   }  int operator==(const Vec4f &V) const   {    return ((data[0] == V.data[0]) &&	    (data[1] == V.data[1]) &&	    (data[2] == V.data[2]) &&	    (data[3] == V.data[3]));   }  int operator!=(const Vec4f &V) const   {    return ((data[0] != V.data[0]) ||	    (data[1] != V.data[1]) ||	    (data[2] != V.data[2]) ||	    (data[3] != V.data[3]));   }  Vec4f& operator+=(const Vec4f &V)   {    data[0] += V.data[0];    data[1] += V.data[1];    data[2] += V.data[2];    data[3] += V.data[3];    return *this;   }  Vec4f& operator-=(const Vec4f &V)   {    data[0] -= V.data[0];    data[1] -= V.data[1];    data[2] -= V.data[2];    data[3] -= V.data[3];    return *this;   }  Vec4f& operator*=(float f)   {    data[0] *= f;    data[1] *= f;    data[2] *= f;    data[3] *= f;    return *this;   }  Vec4f& operator/=(float f)   {    data[0] /= f;    data[1] /= f;    data[2] /= f;    data[3] /= f;    return *this;  }  // OPERATIONS  float Dot2(const Vec4f &V) const   {    return data[0] * V.data[0] +      data[1] * V.data[1];   }  float Dot3(const Vec4f &V) const   {    return data[0] * V.data[0] +      data[1] * V.data[1] +      data[2] * V.data[2];   }  float Dot4(const Vec4f &V) const   {    return data[0] * V.data[0] +      data[1] * V.data[1] +      data[2] * V.data[2] +      data[3] * V.data[3];   }    // STATIC OPERATIONS  static void Add(Vec4f &a, const Vec4f &b, const Vec4f &c )   {    a.data[0] = b.data[0] + c.data[0];    a.data[1] = b.data[1] + c.data[1];    a.data[2] = b.data[2] + c.data[2];    a.data[3] = b.data[3] + c.data[3];   }  static void Sub(Vec4f &a, const Vec4f &b, const Vec4f &c )   {    a.data[0] = b.data[0] - c.data[0];    a.data[1] = b.data[1] - c.data[1];    a.data[2] = b.data[2] - c.data[2];    a.data[3] = b.data[3] - c.data[3];   }  static void CopyScale(Vec4f &a, const Vec4f &b, float c )   {    a.data[0] = b.data[0] * c;    a.data[1] = b.data[1] * c;    a.data[2] = b.data[2] * c;    a.data[3] = b.data[3] * c;   }  static void AddScale(Vec4f &a, const Vec4f &b, const Vec4f &c, float d )   {    a.data[0] = b.data[0] + c.data[0] * d;    a.data[1] = b.data[1] + c.data[1] * d;    a.data[2] = b.data[2] + c.data[2] * d;    a.data[3] = b.data[3] + c.data[3] * d;   }  static void Average(Vec4f &a, const Vec4f &b, const Vec4f &c )   {    a.data[0] = (b.data[0] + c.data[0]) * 0.5f;    a.data[1] = (b.data[1] + c.data[1]) * 0.5f;    a.data[2] = (b.data[2] + c.data[2]) * 0.5f;    a.data[3] = (b.data[3] + c.data[3]) * 0.5f;   }  static void WeightedSum(Vec4f &a, const Vec4f &b, float c, const Vec4f &d, float e )   {    a.data[0] = b.data[0] * c + d.data[0] * e;    a.data[1] = b.data[1] * c + d.data[1] * e;    a.data[2] = b.data[2] * c + d.data[2] * e;    a.data[3] = b.data[3] * c + d.data[3] * e;   }  static void Cross3(Vec4f &c, const Vec4f &v1, const Vec4f &v2)   {    float x = v1.data[1]*v2.data[2] - v1.data[2]*v2.data[1];    float y = v1.data[2]*v2.data[0] - v1.data[0]*v2.data[2];    float z = v1.data[0]*v2.data[1] - v1.data[1]*v2.data[0];    c.data[0] = x; c.data[1] = y; c.data[2] = z;   }  // INPUT / OUTPUT  void Write(FILE *F = stdout)   {    fprintf (F, "%f %f %f %f\n",data[0],data[1],data[2],data[3]);   }};#endif

#include "vectors.h"#include "matrix.h"void main (){	Matrix A, B, C;	A.SetToIdentity ();	A.Scale (2);	B.SetToIdentity ();	B.Translate (1, 1, 1);		C = A + B;	A.Write ();	B.Write ();	C.Write ();}