/* 3x3 matrices are easiest to invert analytically.  Compute the
 * adjoint matrix of cofactors, the determinant of the original
 * matrix, and multiply by 1/determinant.  Note the use of cyclic
 * indexing to avoid the "-1 in even cells" rules that most
 * treatments use. */
#define LOOP1Di for(i=0; i<3; i++)
#define LOOP2Dij for(i=0; i<3; i++) for(j=0; j<3; j++)
void invert(double m[3][3], double out[3][3])
{
#define M(i,j) (m[(i)%3][(j)%3])
    int i, j;
    double det = 0;
    LOOP2Dij out[j][i] = M(i+1, j+1) * M(i+2, j+2) - M(i+2, j+1) * M(i+1, j+2);
    LOOP1Di  det += m[0][i] * out[i][0];
    LOOP2Dij out[i][j] *= (1/det);
}