typedef float real;
class Matrix;
class Vector {
public:
	enum {length=4};
	Vector();
	Vector(real x,real y,real z);
	Vector(const Vector& rhs);
	real x() const { return data[0];}
	real y() const { return data[1];}
	real z() const { return data[2];}
	real& x() { return data[0];}
	real& y() { return data[1];}
	real& z() { return data[2];}
	Vector mul(const Vector &rhs);
	Vector mul(const Matrix& rhs);
	Vector mul(const real scale);
	real dot(const Vector &rhs);
	Vector cross(const Vector& rhs);
	real mag();
	Vector unit();
	real get(int i) const { return data[i];}
private:
	real data[length];
};

class Matrix {
	public:
		enum {length=4};
		Matrix mul(const Matrix& rhs);
		Matrix();
		Matrix(const Matrix& rhs);
		real get(int i,int j) const { return data[i][j];}
		void set(int i,int j,real val)  { data[i][j]=val;}
	private:
		real data[length][length];
};
#include <cmath>
#include <cstring>
Matrix::Matrix() {
	for(int i=0;i<length;++i)
		for(int j=0;j<length;++j)
			data[i][j]=i==j ? 1:0;
}
Matrix::Matrix(const Matrix&rhs) {
	memcpy(data,rhs.data,sizeof(data));
}

Matrix Matrix::mul(const Matrix &rhs) {
		Matrix ret;
		for (int i = 0; i < length; ++i) {
			for (int j = 0; j < length; ++j) {
				real cell = 0;
				for (int r = 0; r < length; ++r) {
					cell += this->get(i, r) * rhs.get(r, j);
				}
				ret.set(i, j, cell);
			}
		}
		return ret;
	}
Vector Vector::mul(const Matrix& mat) {

		 int m = 1;
		int n = length;
		int p = n;
		Vector ret;

			for (int j = 0; j < p; ++j) {
				int cell = 0;
				for (int r = 0; r < n; ++r) {
					cell += this->get(r) * mat.get(r, j);
				}
				ret.data[j]= cell;
			}

		return ret;
	}

Vector::Vector() {
	for(int i=0;i<length-1;++i)
		data[i]=0;
	data[length-1]=1;
}
Vector::Vector(real x,real y,real z) {
	data[0]=x;
	data[1]=y;
	data[2]=z;
	data[3]=1;
}
Vector Vector::mul(real scale) {
	Vector v;
	for(int i=0;i<length;++i) {
		v.data[i] = this->data[i] * scale;
	}
	return v;
}
Vector Vector::cross(const Vector& rhs) {
	return Vector(y()*rhs.z() - rhs.y()*z() , z()*rhs.x()-rhs.z()*x(), x()*rhs.y()-rhs.x()*y());
}
Vector::Vector(const Vector &rhs) {
	for(int i=0;i<length;++i)
		data[i]=rhs.data[i];
}
real Vector::mag() {
	real tot=0;
	for(int i=0;i<length;++i)
		tot+= data[i]*data[i];
	return sqrt(tot);
}
Vector Vector::unit() {
	real l=this->mag();
	real s=1/l;
	return this->mul(s);
}
real Vector::dot(const Vector& rhs) {
	float d=0;
	for(int i=0;i<length;++i)
		d += data[i]*rhs.data[i];
	return d;
}
Matrix frustrum(real l,real r,real b,real t,real n,real f) {
	Matrix m;
	m.set(0,0, (2.0f*n)/(r-l));
	m.set(0,2, (r+l)/(r-l));
	m.set(1,1, (2.0f*n)/(t-b));
	m.set(1,2, (t+b)/(t-b));
	m.set(2,2,(-(f+n))/(f-n));
	m.set(2,3,(-(2.0f*f*n))/(f-n));
	m.set(3,2,-1);
	return m;
}

#include <cstdio>

void printMat(const Matrix& m) {
	int l=m.length;
	for(int i=0;i<l;++i) {
		printf("[");
		for(int j=0;j<l;++j) {
			printf(" %5.02f ",m.get(i,j));
		}
		printf("]\n");
	}
}
void printVec(const Vector& v) {
	for(int i=0;i<v.length;++i) {
		printf("[%5.02f]\n",v.get(i));
	}
}
void printVec2(const Vector&v) {
	for(int i=0;i<v.length-1;++i) {
		printf("[%5.02f]\n",v.get(i)/v.get(v.length-1));
	}
}
int main(void) {
			 // simple projection test
	Vector X(1,0,0);
	Vector Y(0,1,0);
	Vector G(5,1,191);
	printf("X.Y = %.2f\n",X.dot(Y));
	printf("X.G = %.2f\n",X.dot(G));
	printf("Y.G = %.2f\n",Y.dot(G));
	real view_height=4;
	Matrix f = frustrum(.5f, -.5f, -1.f , 1.f , 1.f, 90.f);
	printMat(f);
	printVec(G);
	Vector s=G.mul(f);
	printVec2(s);

return 0;
}