[cfarm-users] cfarm26 FP inaccuracies

[Bruno Haible]

Jacob Bachmeyer wrote:
> > and the paranoia.c stringent IEEE floating point test reports one flaw:
> >
> > Checking rounding on multiply, divide and add/subtract.
> > * is neither chopped nor correctly rounded.
> > / is neither chopped nor correctly rounded.
> > Addition/Subtraction neither rounds nor chops.
> > Sticky bit used incorrectly or not at all.
> > FLAW:  lack(s) of guard digits or failure(s) to correctly round or chop
> >
> > On a machine where the test succeeds, paranoia.c reports no flaws.
> 
> To make this even weirder, cfarm26 and cfarm27 are different VMs on the 
> same physical box, yet they give different results from paranoia.c.

Similarly, GNU clisp's floatparam.c [1] reports
  #define double_rounds_correctly 0
on cfarm26 (x86 with glibc), but
  #define double_rounds_correctly 1
on cfarm27 (x86 with musl) and on non-x86 platforms.

And when you compile it with gcc options or "-ffloat-store" or
"-msse2 -mfpmath=sse", even on cfarm26, it reports
  #define double_rounds_correctly 1

> Perhaps there is something weird about Debian's libm on cfarm26? Or the 
> compiler defaults are different?
> 
> Are long doubles being used internally somewhere doubles are expected to 
> be used throughout?  Paranoia's indication of about 11 extra bits would 
> be consistent with the x87 80-bit format vs. standard 64-bit double...

Yes. x86 CPUs have two floating-point units, that are available for
'float' and 'double' computations: the 387 and the MMX/SSE unit.

While the 'float' computations are done the same way in both units,
for 'double' computations there is a difference: In the 387 unit
they are performed with excess precision (as 80-bit numbers instead
of 64-bit numbers), if a certain bit in the FPU control word permits it.

musl libc disables this bit in the FPU control word by default, thus
making it impossible to do computations with 80-bit numbers. Thus, on musl:
  - 'long double' is the same as 'double', not longer,
  - but 'double' computations produce the same reproducible results as on
    most platforms.

glibc, on the other hand, enables this bit. Thus computations with
'long double' are possible. But 'double' computations have excess
precision. There are several ways to avoid this excess precision:
  - Use 'volatile double' instead of 'double' everywhere. This comes
    at a performance cost.
  - Use gcc option '-ffloat-store'. This comes at the same performance cost.
  - Instruct gcc to use the MMX/SSE unit instead; see above.

Btw, on x86_64 you have this problem as well, if you use the gcc option
"-mfpmath=387". Fortunately this option is not the default, which is why
this problem with the 387 is not seen in practice on x86_64 systems.

Bruno

[1] https://gitlab.com/gnu-clisp/clisp/-/raw/master/src/floatparam.c?ref_type=heads