typedef GLfloat vec4[4];
inline void vec4_add(vec4 a, vec4 b)
{
    int i;
    for(i=0; i<4; ++i)
        a[i] += b[i];
}
inline void vec4_sub(vec4 a, vec4 b)
{
    int i;
    for(i=0; i<4; ++i)
        a[i] -= b[i];
}
inline void vec4_scale(vec4 v, GLfloat s)
{
    int i;
    for(i=0; i<4; ++i)
        v[i] *= s;
}
inline GLfloat vec4_inner_product(vec4 a, vec4 b)
{
    GLfloat p = 0.;
    int i;
    for(i=0; i<4; ++i)
        p += b[i]*a[i];
    return p;
}
inline GLfloat vec4_length(vec4 v)
{
    return sqrtf(vec4_inner_product(v,v));
}
inline void vec4_normalize(vec4 v)
{
    GLfloat k = 1.0 / vec4_length(v);
    vec4_scale(v, k);
}
inline void vec4_cross(vec4 a, vec4 b)
{
    vec4 c;
    c[0] = a[1]*b[2] - a[2]*b[1];
    c[1] = a[2]*b[0] - a[0]*b[2];
    c[2] = a[0]*b[1] - a[1]*b[0];
    c[3] = 0.;
    memcpy(a, c, sizeof(a));
}

typedef vec4 mat4x4[4];
inline void mat4x4_identity(mat4x4 M)
{
    int i, j;
    M[0][0] = 1; M[1][0] = 0; M[2][0] = 0; M[3][0] = 0;
    M[0][1] = 0; M[1][1] = 1; M[2][1] = 0; M[3][1] = 0;
    M[0][2] = 0; M[1][2] = 0; M[2][2] = 1; M[3][2] = 0;
    M[0][3] = 0; M[1][3] = 0; M[2][3] = 0; M[3][3] = 1;
    /*for(j=0; j<4; ++j)
        for(i=0; i<4; ++i) {
            M[i][j] = i==j ? 1 : 0;
    }*/
}
inline void mat4x4_cpy(mat4x4 M, mat4x4 N)
{
    int i, j;
    for(j=0; j<4; ++j) {
        for(i=0; i<4; ++i) {
            M[i][j] = N[i][j];
        }
    }
}
inline void mat4x4_mul(mat4x4 M, mat4x4 b)
{
    mat4x4 a;
    int i, j, k;
    memcpy(a, M, sizeof(a));
    for(j=0; j<4; ++j) {
        for(i=0; i<4; ++i) {
            M[i][j] = 0;
            for(k=0; k<4; ++k) {
                M[i][j] += a[i][k]*b[k][j];
            }
        }
    }
}
inline void mat4x4_trans(mat4x4 M, GLfloat x, GLfloat y, GLfloat z)
{
    mat4x4 T; mat4x4_identity(T);
    T[3][0] = x;
    T[3][1] = y;
    T[3][2] = z;
    mat4x4_mul(M, T);
}
inline void mat4x4_rot(mat4x4 M, GLfloat x, GLfloat y, GLfloat z, GLfloat angle)
{
}
inline void mat4x4_rot_X(mat4x4 M, GLfloat angle)
{
    GLfloat s = sinf(angle);
    GLfloat c = cosf(angle);
    mat4x4 R = {
        {1, 0, 0, 0},
        {0, c, s, 0},
        {0,-s, c, 0},
        {0, 0, 0, 1}
    };
    mat4x4_mul(M, R);
}
inline void mat4x4_rot_Y(mat4x4 M, GLfloat angle)
{
    GLfloat s = sinf(angle);
    GLfloat c = cosf(angle);
    mat4x4 R = {
        {c, 0, s, 0},
        {0, 1, 0, 0},
        {-s, 0, c, 0},
        {0, 0, 0, 1}
    };
    mat4x4_mul(M, R);
}
inline void mat4x4_rot_Z(mat4x4 M, GLfloat angle)
{
    GLfloat s = sinf(angle);
    GLfloat c = cosf(angle);
    mat4x4 R = {
        {c, s, 0, 0},
        {-s, c, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1}
    };
    mat4x4_mul(M, R);
}
inline void mat4x4_row(vec4 r, mat4x4 M, int i)
{
    int k;
    for(k=0; k<4; ++k)
        r[k] = M[k][i];
}
inline void mat4x4_col(vec4 r, mat4x4 M, int i)
{
    int k;
    for(k=0; k<4; ++k)
        r[k] = M[i][k];
}
inline void mat4x4_cpy_T(mat4x4 M, mat4x4 N)
{
    int i, j;
    for(j=0; j<4; ++j) {
        for(i=0; i<4; ++i) {
            M[i][j] = N[j][i];
        }
    }
}