4
votes

I am trying to implement emulated double-precision in GLSL, and I observe a strange behaviour difference leading to subtle floating point errors in GLSL.

Consider the following fragment shader, writing to a 4-float texture to print the output.

layout (location = 0) out vec4 Output
uniform float s;
void main()
{
  float a = 0.1f;
  float b = s;

  const float split = 8193.0; // = 2^13 + 1

  float ca = split * a;
  float cb = split * b;

  float v1a = ca - (ca - a);
  float v1b = cb - (cb - b);

  Output = vec4(a,b,v1a,v1b);
}

This is the output I observe

GLSL output with uniform :

a = 0.1            0x3dcccccd
b = 2.86129e-06    0x36400497
v1a = 0.0999756    0x3dccc000
v1b = 2.86129e-06  0x36400497

Now, with the given values of b1 and b2 as inputs, the value of v2b does not have the expected result. Or at least it does not have the same result as on CPU (as seen here):

C++ output :

a = 0.100000     0x3dcccccd
b = 0.000003     0x36400497
v1a = 0.099976   0x3dccc000
v1b = 0.000003   0x36400000

Note the discrepancy for the value of v1b (0x36400497 vs 0x36400000).

So in an effort to figure out what was happening (and who was right) I attempted to redo the computation in GLSL, replacing the uniform value by a constant, using a slightly modified shader, where I replaced the uniform by it value.

layout (location = 0) out vec4 Output
void main()
{
  float a = 0.1f;
  float b = uintBitsToFloat(0x36400497u);

  const float split = 8193.0; // = 2^13 + 1

  float ca = split * a;
  float cb = split * b;

  float v1a = ca - (ca - a);
  float v1b = cb - (cb - b);

  Output = vec4(a,b,v1a,v1b);
}

This time, I get the same output as the C++ version of the same computation.

GLSL output with constants :

a = 0.1            0x3dcccccd
b = 2.86129e-06    0x36400497
v1a = 0.0999756    0x3dccc000
v1b = 2.86102e-06  0x36400000

My question is, what makes the floating point computation behave differently between a uniform variable and a constant ? Is this some kind of behind-the-scenes compiler optimization ?

Here are my OpenGL vendor strings from my laptop's intel GPU, but I also observed the same behaviour on a nVidia card as well.

Renderer : Intel(R) HD Graphics 520
Vendor   : Intel
OpenGL   : 4.5.0 - Build 23.20.16.4973
GLSL     : 4.50 - Build 23.20.16.4973
2
you are comparing different things. If I get it right the discrepancy is on CPU not on GLSL ... const float might get optimized by compiler and compute at compile time as double hence the difference between variable and constant. There is usually also difference between CPU and GPU results (not your case however) because GPU floats are not always standard floats they have sometimes different number of bits ...Spektre
If you change float v1b = cb - (cb - b); to float tempb = cb - b; float v1b = cb - tempb; (and similarly for v1a), does the behavior change?Eric Postpischil
Try adding the precise modifier to v1a and v1b as this looks like a case where floating-point re-association or FMA contraction are being applied as a performance optimization (allowed by default, I think), changing the result.njuffa
re: @Spektre : the problem is that the first version has the wrong result (as it does not conform to IEEE 754 standard the other two follow).Louen

2 Answers

3
votes

GPU do not necessarily have/use IEEE 754 some implementations have smaller number of bits so there is no brainer the result will be different. Its the same as you would compare float vs. double results on FPU. However you can try to enforce precision if your GLSL implementation allows it see:

In worse case use double and dvec if your GPU allows it but beware there are no 64bit interpolators yet (at least to my knowledge).

To rule out rounding due to passing results by texture see:

You can also check the number of mantissa bits on your GPU simply by printing

1.0+1.0/2.0
1.0+1.0/4.0
1.0+1.0/8.0
1.0+1.0/16.0
...
1.0+1.0/2.0^i

The i of last number not printed as 1.0 is the number of the mantissa bits. So you can check if it is 23 or not ...

3
votes

So, as mentioned by @njuffa, in the comments, the problem was solved by using the precise modifier on the values which depended on rigorous IEEE754 operations:

layout (location = 0) out vec4 Output
uniform float s;
void main()
{
  float a = 0.1f;
  float b = s;

  const float split = 8193.0; // = 2^13 + 1

  precise float ca = split * a;
  precise float cb = split * b;

  precise float v1a = ca - (ca - a);
  precise float v1b = cb - (cb - b);

  Output = vec4(a,b,v1a,v1b);
}

Output :

a = 0.1            0x3dcccccd
b = 2.86129e-06    0x36400497
v1a = 0.0999756    0x3dccc000
v1b = 2.86102e-06  0x36400000

Edit: it is higly probable that only the last precise are needed to constrain the operations leading to its computation to avoid unwanted optimizations.