help requested with long double issue on 32-bit FreeBSD

[Bruce Evans]

On Fri, 9 Aug 2013, Zimmermann Paul wrote:

> [this mail was sent on the MPFR list, but Steve Kargl suggested to post it
> also here]
>
> we need help to investigate an error we get with the development version of
> MPFR on i686-freebsd. The error occurs on the hydra.nixos.org continuous
> integration platform. Unfortunately we don't have direct access to the
> corresponding computer, thus isolating the issue is not easy.
>
> The error we get can be seen at http://hydra.nixos.org/build/5666171/log/raw.
> Starting from the long double value in 3, and subtracting e=0.25, we get the
> long double value in 4, which is obviously wrong in the last bits:
>
> 3: d=2.632048527199790448306230411512629530e-01
> 4:  d=1.320485271997904469509776959057489876e-02
>                        ^^^^^^^^^^^^^^^^^^^^
>
> We believe this is a compiler bug (GCC 4.2.1 was used in this test).
>
> If you have access to a 32-bit computer under FreeBSD, please could you try
> to compile the development version and run the "tset_ld" test, after doing
> export GMP_CHECK_RANDOMIZE=1376689137?
>
> If you can reproduce the error, please tell Vincent Lefèvre and myself
> (Vincent.Lefevre at ens-lyon.fr and Paul.Zimmermann at inria.fr); we'll help
> you to isolate the problem (e.g., trying different optimization levels).
>
> If you can't reproduce (on a 32-bit FreeBSD machine), please tell us too;
> this might indicate the bug was fixed in a later gcc version.

I didn't run the tests, but it was easy to see that this is just the
result of FreeBSD's default rounding precision in 32-bit mode being
53 bits and not switching to 64 bits in the tests.  Long doubles having
64 mantissa bits obviously cannot work right with a rounding precision
of 53, so something must switch to 64 bits in any program or part of
a program that uses long doubles and actually needs them.  Also, gcc
supports the default rounding precision to a fault.  It rounds even
literal constants to 53 bits at compile time, so it is difficult to
even initialize d to the above values.

The following programs demonstrates several variations of the problem.

% #include <ieeefp.h>
% 
% unsigned char cd[10] = {
%     0x05, 0x7d, 0x5f, 0x29, 0x55, 0xc9, 0xc2, 0x86, 0xfd, 0x3f,
% };
% 
% main()
% {
% 	long double d = 2.632048527199790448306230411512629530e-01L;
% 	long double e = 1.320485271997904469509776959057489876e-02L;
% 
% #ifdef FIX_RUNTIME_ROUNDING_PRECISION
% 	fpsetprec(FP_PE);
% #endif
% #ifdef FIX_COMPILE_TIME_ROUNDING_PRECISION
% 	d = *(long double *)&cd;
% #endif
% 	printf("%La\n", d);
% 	printf("%.37Le\n", d);
% 	printf("%La\n", e);
% 	printf("%.37Le\n", e);
% 	printf("%La\n", d - 0.25);
% 	printf("%.37Le\n", d - 0.25);
% }

This is easiest to test using CC=clang, since clang is not bug-for-bug
compatible with gcc4-2.1 and doesn't round even long double constant
expressions to 53 bits at compile time.  So it gives your 3:d directly.
gcc needs the magic FIX_COMPILE_TIME_ROUNDING_PRECISION to give same
3:d.  (Even hex constants are rounded to 53 bits by gcc.)  Then the
result of d - 0.25 is the same as your 4:d in the default rounding
precision, but FIX_RUNTIME_ROUNDING_PRECISION fixes this.

The bug is also easy to see by printing the values in hex on a system
without the bug.  On FreeBSD-amd64:

% 0x1.0d8592aa52befa0ap-2
% 2.6320485271997904483062304115126295301e-01
% 0x1.b0b2554a57df4p-7
% 1.3204852719979044695097769590574898757e-02
% 0x1.b0b2554a57df414p-7
% 1.3204852719979044830623041151262953008e-02

The intermediate result 4:d is the final result rounded to 53 bits.
Note that this was printed by a program compiled by clang.  clang
presumably rounding 4:d to 64 bits, but it ends up with only 53 nonzero
bits since it was produced by a program that rounded it to that many
and the decimal format has more than enough precision to recover the
53-bit value when rounded to 64 bits.

It is interesting that the decimal digits for d and e line up but the
binary digits don't.  I prefer a different %a format that lines up the
binary digits more often, but this is unfortunately not supported by
printf().

This is hard to fix.  Programs using long double precision should switch
the rounding precision at runtime using fpsetprec() as above.  But
naive programmers won't know to do this, and apparently even people
at INRIA don't know this :-).  Newer long double functions in FreeBSD
libm do the switch at runtime using macros that hide most of the
complications including avoiding switching the rounding precision if
it is already 64 bits (switching the precision is very slow, but testing
to avoid it when unnecessary only takes large code when it is in fact
unnecessary).

Compile-time rounding of constant expressions is even harder to fix.
FreeBSD libm used to mostly use the method of writing long double
constants as the sum of 2 double constants, with a volatile hack to
force the evaluation at runtime (in 64-bit precision iff someone
remembered to switch the rounding precision).  Newer parts mostly
use macro that hides some of the details of initializing a union
containing the necessary values in bits.

Bruce