/* 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);
}