Invert 4x4 matrix - Numerical most stable solution needed

Question

I want to invert a 4x4 matrix. My numbers are stored in fixed-point format (1.15.16 to be exact).

With floating-point arithmetic I usually just build the adjoint matrix and divide by the determinant (e.g. brute force the solution). That worked for me so far, but when dealing with fixed point numbers I get an unacceptable precision loss due to all of the multiplications used.

Note: In fixed point arithmetic I always throw away some of the least significant bits of immediate results.

So - What's the most numerical stable way to invert a matrix? I don't mind much about the performance, but simply going to floating-point would be to slow on my target architecture.

Answer 1

Meta-answer: Is it really a general 4x4 matrix? If your matrix has a special form, then there are direct formulas for inverting that would be fast and keep your operation count down.

For example, if it's a standard homogenous coordinate transform from graphics, like:

[ux vx wx tx]
[uy vy wy ty]
[uz vz wz tz]
[ 0  0  0  1]

(assuming a composition of rotation, scale, translation matrices)

then there's an easily-derivable direct formula, which is

[ux uy uz -dot(u,t)]
[vx vy vz -dot(v,t)]
[wx wy wz -dot(w,t)]
[ 0  0  0     1    ]

(ASCII matrices stolen from the linked page.)

You probably can't beat that for loss of precision in fixed point.

If your matrix comes from some domain where you know it has more structure, then there's likely to be an easy answer.