+ if(flags==IDENTITY)
+ return other;
+ else if(other.flags==IDENTITY)
+ return *this;
+ else if(flags==TRANSLATE && !(other.flags&ARBITARY))
+ {
+ Matrix result = other;
+ result.matrix[12] += matrix[12];
+ result.matrix[13] += matrix[13];
+ result.matrix[14] += matrix[14];
+ result.flags |= flags;
+ return result;
+ }
+ else if(!(flags&ARBITARY) && other.flags==TRANSLATE)
+ {
+ Matrix result = *this;
+ const double *m = other.matrix;
+ result.matrix[12] += matrix[0]*m[12]+matrix[4]*m[13]+matrix[8]*m[14];
+ result.matrix[13] += matrix[1]*m[12]+matrix[5]*m[13]+matrix[9]*m[14];
+ result.matrix[14] += matrix[2]*m[12]+matrix[6]*m[13]+matrix[10]*m[14];
+ result.flags |= other.flags;
+ return result;
+ }
+ else
+ {
+ Matrix result;
+ fill(result.matrix, result.matrix+16, 0.0);
+ for(unsigned i=0; i<4; ++i)
+ for(unsigned j=0; j<4; ++j)
+ for(unsigned k=0; k<4; ++k)
+ result.matrix[i+j*4] += matrix[i+k*4]*other.matrix[k+j*4];
+ result.flags = flags|other.flags;
+ return result;
+ }
+}
+
+Matrix &Matrix::operator*=(const Matrix &other)
+{
+ multiply(other);
+ return *this;
+}
+
+Vector4 Matrix::operator*(const Vector4 &vec) const
+{
+ if(flags==IDENTITY)
+ return vec;
+ else if(flags==TRANSLATE)
+ return Vector4(vec.x+vec.w*matrix[12], vec.y+vec.w*matrix[13], vec.z+vec.w*matrix[14], vec.w);
+ else if(flags==SCALE)
+ return Vector4(vec.x*matrix[0], vec.y*matrix[5], vec.z*matrix[10], vec.w);
+ else
+ {
+ Vector4 result;
+ result.x = vec.x*matrix[0]+vec.y*matrix[4]+vec.z*matrix[8]+vec.w*matrix[12];
+ result.y = vec.x*matrix[1]+vec.y*matrix[5]+vec.z*matrix[9]+vec.w*matrix[13];
+ result.z = vec.x*matrix[2]+vec.y*matrix[6]+vec.z*matrix[10]+vec.w*matrix[14];
+ result.w = vec.x*matrix[3]+vec.y*matrix[7]+vec.z*matrix[11]+vec.w*matrix[15];
+ return result;
+ }
+}
+
+double Matrix::operator[](unsigned i) const
+{
+ if(i>=16)
+ throw out_of_range("Matrix::operator[]");
+ return matrix[i];
+}
+
+Matrix Matrix::translation(double x, double y, double z)
+{
+ Matrix result;
+ result.matrix[12] = x;
+ result.matrix[13] = y;
+ result.matrix[14] = z;
+ result.flags |= TRANSLATE;
+ return result;
+}
+
+Matrix Matrix::rotation(double a, double x, double y, double z)
+{
+ double l = sqrt(x*x+y*y+z*z);
+ x /= l;
+ y /= l;
+ z /= l;
+ double c = cos(a);
+ double s = sin(a);
+
+ // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_given_an_axis_and_an_angle
+ Matrix result;
+ result.matrix[0] = c+x*x*(1-c);
+ result.matrix[1] = y*x*(1-c)+z*s;
+ result.matrix[2] = z*x*(1-c)-y*s;
+ result.matrix[4] = x*y*(1-c)-z*s;
+ result.matrix[5] = c+y*y*(1-c);
+ result.matrix[6] = z*y*(1-c)+x*s;
+ result.matrix[8] = x*z*(1-c)+y*s;
+ result.matrix[9] = y*z*(1-c)-x*s;
+ result.matrix[10] = c+z*z*(1-c);
+ result.flags |= ROTATE;
+ return result;
+}
+
+Matrix Matrix::rotation_deg(double a, double x, double y, double z)
+{
+ return rotation(a*M_PI/180, x, y, z);
+}
+
+Matrix Matrix::scaling(double s)
+{
+ Matrix result;
+ result.matrix[0] = s;
+ result.matrix[5] = s;
+ result.matrix[10] = s;
+ result.flags |= SCALE;
+ return result;
+}
+
+Matrix Matrix::scaling(double x, double y, double z)
+{
+ Matrix result;
+ result.matrix[0] = x;
+ result.matrix[5] = y;
+ result.matrix[10] = z;
+ result.flags |= SCALE|ARBITARY;
+ return result;
+}
+
+Matrix Matrix::ortho(double l, double r, double b, double t, double n, double f)
+{
+ if(l==r || b==t || n==f)
+ throw invalid_argument("Matrix::ortho");
+
+ Matrix result;
+ result.matrix[0] = 2/(r-l);
+ result.matrix[5] = 2/(t-b);
+ result.matrix[10] = -2/(f-n);
+ result.matrix[12] = -(r+l)/(r-l);
+ result.matrix[13] = -(t+b)/(t-b);
+ result.matrix[14] = -(f+n)/(f-n);
+ result.flags = TRANSLATE|SCALE|ARBITARY;
+ return result;
+}
+
+Matrix Matrix::ortho_centered(double w, double h)
+{
+ return ortho(-w/2, w/2, -h/2, h/2, -1, 1);
+}
+
+Matrix Matrix::ortho_bottomleft(double w, double h)
+{
+ return ortho(0, w, 0, h, -1, 1);
+}
+
+Matrix Matrix::ortho_topleft(double w, double h)
+{
+ return ortho(0, w, h, 0, -1, 1);
+}
+
+Matrix Matrix::frustum(double l, double r, double b, double t, double n, double f)
+{
+ if(l==r || b==t || n<=0 || f<=n)
+ throw invalid_argument("Matrix::frustum");
+
+ Matrix result;
+ result.matrix[0] = 2*n/(r-l);
+ result.matrix[5] = 2*n/(t-b);
+ result.matrix[8] = (r+l)/(r-l);
+ result.matrix[9] = (t+b)/(t-b);
+ result.matrix[10] = -(f+n)/(f-n);
+ result.matrix[11] = -1;
+ result.matrix[14] = -2*f*n/(f-n);
+ result.matrix[15] = 0;
+ result.flags = ARBITARY;
+ return result;
+}
+
+Matrix Matrix::frustum_centered(double w, double h, double n, double f)
+{
+ return frustum(-w/2, w/2, -h/2, h/2, n, f);
+}
+
+Matrix Matrix::perspective(double h, double a, double n, double f)
+{
+ double hh = tan(h/2)*n;
+ return frustum(-hh*a, hh*a, -hh, hh, n, f);
+}
+
+
+GLenum MatrixStack::current_mode = GL_MODELVIEW;
+
+MatrixStack::MatrixStack(GLenum m):
+ mode(m)
+{
+ matrices.reserve(mode==GL_MODELVIEW ? 32 : 4);
+ matrices.push_back(Matrix());
+}
+
+MatrixStack::MatrixStack():
+ mode(0)
+{
+ matrices.reserve(32);
+ matrices.push_back(Matrix());
+}
+
+const Matrix &MatrixStack::top() const
+{
+ return matrices.back();
+}
+
+void MatrixStack::load(const Matrix &m)
+{
+ matrices.back() = m;
+ update();
+}
+
+void MatrixStack::multiply(const Matrix &m)
+{
+ matrices.back() *= m;
+ update();
+}
+
+void MatrixStack::push()
+{
+ matrices.push_back(top());
+}
+
+void MatrixStack::pop()
+{
+ if(matrices.size()==1)
+ throw stack_underflow("MatrixStack::pop()");
+
+ matrices.pop_back();
+ update();
+}
+
+void MatrixStack::update()
+{
+ if(!mode)
+ return;
+
+ if(mode!=current_mode)
+ {
+ glMatrixMode(mode);
+ current_mode = mode;
+ }
+
+ glLoadMatrixd(matrices.back().data());