Re: Changing x87 precision to full 63bit as default

[Joerg Sonnenberger <joerg%britannica.bec.de@localhost>]

On Wed, Nov 13, 2013 at 07:40:31AM +0900, tsugutomo.enami%jp.sony.com@localhost 
wrote:
> Joerg Sonnenberger <joerg%britannica.bec.de@localhost> writes:
> 
> > What exactly have you tried?
> 
> This code:
> 
> #include <stdio.h>
> #include <ieeefp.h>
> 
> int
> main(int argc, char *argv[])
> {
>         volatile double x, y, r;
> 
>         if (argc > 1) {
>                 fp_prec_t oldprec;
>                 oldprec = fpsetprec(FP_PD);
>                 printf("prec %d -> %d\n", oldprec, FP_PD);
>         }
>         x = 3002399751580332.0;
>         y = 3002399751580331.0;
>         r = x / y;
>         printf("%.16f, %.16f, %.16f\n", x, y, r);
> 
>         return 0;
> }
> 
> And result was
> % ./t
> 3002399751580332.0000000000000000, 3002399751580331.0000000000000000, 
> 1.0000000000000004
> % ./t setprec
> prec 3 -> 2
> 3002399751580332.0000000000000000, 3002399751580331.0000000000000000, 
> 1.0000000000000002
> 
> Compiling with -ffloat-store and/or -fexcess-precision=standard was
> some result.

Ah, thanks. Took a moment to understand what is happening here. The
machine code produced is actually correct -- the problem is the division
itself. With the choosen values, the exact value of r is slightly
smaller than (1 + DBL_EPSILON / 2). In double precision mode, this gets
rounded down. In extended precision mode, this gets rounded *up* to (1+
DBL_EPSILON / 2) for the *internal* 80bit representation and then stored
to double, which results in another rounding up. So yes, I think for
this specific case wrapping the computations with fpsetprec if
FP_EVAL_METHOD == 2 and restoring the old value afterwards is the best
approach.

Joerg