unpredictable floating point errors from matrix decomposition

Question

I'm trying to decompose perspective matrix's near and far distances by using this formula:

near = m32 / (m22 - 1);
far  = m32 / (m22 + 1);

Here the perspective matrix test parameters:

  aspect   = 0.782f;
  fovy     = glm_rad(49.984f);
  nearDist = 0.1550385f;
  farDist  = 6000.340975f;

  glm_perspective(fovy, aspect, nearDist, farDist, proj);

Here what I'm doing to get near and far values (proj is column-major matrix):

far  = proj[3][2] / (proj[2][2] + 1.0f);
near = proj[3][2] / (proj[2][2] - 1.0f)

Results:

near = 0.155039 
far  = 5993.506348 

Near seems acceptable but far is not :/ If I use small value for far then I get more accurate results (right values are decomposed values):

farDist = 600.340975 (near, far): 0.155039 600.319885 
farDist = 60.340975f (near, far): 0.155039 60.340946 

Is there something wrong with math? What options do I have (without using double to store matrix)?

You can see perspective matrix formula here: https://www.khronos.org/registry/OpenGL-Refpages/gl2.1/xhtml/gluPerspective.xml

m22 = (near + far) / (near - far)
m32 = 2 * near * far / (near - far)

and implementation (line number may change by time): https://github.com/recp/cglm/blob/master/include/cglm/cam.h#L211

Answer 1

The problem is that the larger the far/near ratio gets, it requires more significant digits to extract far from the perspective matrix.

When the far/near ratio increases, m22 = (near+far)/(near-far) gets closer to 1.

For example, using double with near=0.155, and far=6,000, we get m22 = 1.0000516680014233. When this is stored as float, it is truncated to 1.0000516.

The important part of the result is the fraction. Even if all other opeartions are done with perfect accuracy, at this point you are already left with just 3 significant digits. This is very similar to catastrophic cancellation.

Essentially, you lose a significant digit every time far/near is multiplied by 10. When far is 6,000,000, the value of m22 would be truncated to 1.0 when stored as float, losing all information.

I tried to demonstrate it in a Jupyter Notebook.

But the real problem is not only that is is impossible to extract far without losing precision, but that the perspecive matrix itself is not accurate.

If you take a vector at z=6,000, apply the perspective matrix, you will not get z=1.0. Instead, applying the perspective matrix to a vector with the incorrect value of far, z=5993.506348 will give you z=1.0. The matrix itself is already wrong, so no method of extracting far can help.

TL;DR: If you want to extract near and far from a perspective matrix with reasonable accuracy, you have to use double.

Edited: added an explanation for the real problem, the original answer regarding catastrophic cancellation is just a second order effect.