[Bug target/46080] New: incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[Bug target/46080] New: incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[vincent at vinc17 dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

           Summary: incorrect precision of sqrtf builtin for x87
                    arithmetic (-mfpmath=387)
           Product: gcc
           Version: 4.4.5
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: target
        AssignedTo: unassigned@gcc.gnu.org
        ReportedBy: vincent@vinc17.org


With -mfpmath=387 (tested on a x86_64 platform), the first sqrtf is computed in
double precision instead of single. On

#include <stdio.h>
#include <math.h>

float x = (float) M_PI;

int main(void)
{
  printf ("%.60f\n", sqrtf(x));
  printf ("%.60f\n", sqrtf(x));
  printf ("%.60f\n", sqrtf(x));
  return 0;
}

I get with various gcc versions (including gcc version 4.6.0 20101009
(experimental) [trunk revision 165234] (Debian 20101009-1)):

$ gcc -Wall -mfpmath=387 bug.c -o bug -lm -O0; ./bug
1.772453875567026715387441981874871999025344848632812500000000
1.772453904151916503906250000000000000000000000000000000000000
1.772453904151916503906250000000000000000000000000000000000000

The bug is also present with -O1 when the sqrtf calls are grouped in a single
printf and disappears if the builtin is disabled with -fno-builtin-sqrtf.

For the first sqrtf, the asm code shows a fsqrt followed by a "call sqrtf"
under some condition, but the condition is not satisfied. The other occurrences
just have a "call sqrtf", which is correct.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[rguenth at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

Richard Guenther <rguenth at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Target|                            |i?86-*-*
             Status|UNCONFIRMED                 |NEW
      Known to work|                            |4.3.4
           Keywords|                            |wrong-code
   Last reconfirmed|                            |2010.10.19 10:49:12
     Ever Confirmed|0                           |1
            Summary|incorrect precision of      |[4.4/4.5/4.6 Regression]
                   |sqrtf builtin for x87       |incorrect precision of
                   |arithmetic (-mfpmath=387)   |sqrtf builtin for x87
                   |                            |arithmetic (-mfpmath=387)
   Target Milestone|---                         |4.4.6
      Known to fail|                            |4.4.3, 4.5.1, 4.6.0

--- Comment #1 from Richard Guenther <rguenth at gcc dot gnu.org> 2010-10-19 10:49:12 UTC ---
Confirmed.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[vincent at vinc17 dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #2 from Vincent Lefèvre <vincent at vinc17 dot org> 2010-10-20 01:51:56 UTC ---
Created attachment 22089
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=22089
sh script to test sqrtf

Similar problems can also be found with:

  printf ("%.60f\n%.60f\n%.60f\n", sqrtf(x), sqrtf(x), sqrtf(x));

I've found that every GCC version I could test was showing some incorrect
behavior (but GCC 4.2.4 was the most consistent one). With the attached script,
I get:

           -DSEP  -O0   -O1   -O2
GCC 3.4.6   SSS   SSS   SDD   SDD
GCC 4.1.3   SSS   SSS   DSS   DDS
GCC 4.2.4   SSS   SSS   DDD   DDD   (x86)
GCC 4.3.5   SSS   SSS   DSS   DDD   (ditto with GCC 4.3.2 on x86)
GCC 4.4.5   DSS   SSD   DSS   DDD

where S means that one gets the result in single precision (as expected) and D
means that one gets the result in double precision.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[ubizjak at gmail dot com]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #3 from Uros Bizjak <ubizjak at gmail dot com> 2010-10-20 08:28:57 UTC ---
(In reply to comment #2)
> Created attachment 22089 [details]
> sh script to test sqrtf
> 
> Similar problems can also be found with:
> 
>   printf ("%.60f\n%.60f\n%.60f\n", sqrtf(x), sqrtf(x), sqrtf(x));
> 
> I've found that every GCC version I could test was showing some incorrect
> behavior (but GCC 4.2.4 was the most consistent one). With the attached script,
> I get:
> 
>            -DSEP  -O0   -O1   -O2
> GCC 3.4.6   SSS   SSS   SDD   SDD
> GCC 4.1.3   SSS   SSS   DSS   DDS
> GCC 4.2.4   SSS   SSS   DDD   DDD   (x86)
> GCC 4.3.5   SSS   SSS   DSS   DDD   (ditto with GCC 4.3.2 on x86)
> GCC 4.4.5   DSS   SSD   DSS   DDD
> 
> where S means that one gets the result in single precision (as expected) and D
> means that one gets the result in double precision.

You should use -ffloat-store to remove excess precision.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jakub at gcc dot gnu.org

--- Comment #4 from Jakub Jelinek <jakub at gcc dot gnu.org> 2010-10-20 09:39:29 UTC ---
Or in 4.6 better yet -std=c99 or -fexcess-precision=standard

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[zimmerma+gcc at loria dot fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #5 from Paul Zimmermann <zimmerma+gcc at loria dot fr> 2010-10-20 09:54:01 UTC ---
(In reply to comment #3)
> You should use -ffloat-store to remove excess precision.

the main issue is not the excess of precision, but the fact that identical
calls
to printf with identical arguments give different values in the same program!

The fact that different versions of GCC behave differently vs precision is a
different issue.

Paul

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #6 from Jakub Jelinek <jakub at gcc dot gnu.org> 2010-10-20 10:46:14 UTC ---
That depends on the register allocation, whether the result of fsqrt needs to
be flushed into memory or can stay in the i387 register stack.  In the former
case it gets rounded from long double to float, in the latter case it doesn't.

OT, I wonder about:
  /* Test the result; if it is NaN, set errno=EDOM because
     the argument was not in the domain.  */
  do_compare_rtx_and_jump (target, target, EQ, 0, GET_MODE (target),
                           NULL_RTX, NULL_RTX, lab,
                           /* The jump is very likely.  */
                           REG_BR_PROB_BASE - (REG_BR_PROB_BASE / 2000 - 1));
Why do we use (EQ target target) here instead of (ORDERED target target)?
At least on i?87 the latter is cheaper than the former, we eventually simplify
it, but it would be IMHO better to avoid generating useless rtx when the
comparison often needs to be split into ORDERED && UNEQ.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[vincent at vinc17 dot org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #7 from Vincent Lefèvre <vincent at vinc17 dot org> 2010-10-20 23:43:33 UTC ---
But there's something strange in the generated code: sometimes the fsqrt
instruction is used, sometimes "call sqrtf" is used (for the same sqrtf() call
in the C source). This is not consistent.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[pinskia at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #8 from Andrew Pinski <pinskia at gcc dot gnu.org> 2010-10-20 23:49:18 UTC ---
%.60f

You really should use hex float to see the diferences.  I bet it is just the
final digit of the hex float that is different and only by one.  This is
actually ok IIRC.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[zimmerma+gcc at loria dot fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #9 from Paul Zimmermann <zimmerma+gcc at loria dot fr> 2010-10-21 19:26:11 UTC ---
(In reply to comment #8)

> You really should use hex float to see the diferences.  I bet it is just the
> final digit of the hex float that is different and only by one.  This is
> actually ok IIRC.

see comment 5. (Moreover the sqrt function must return correct rounding
according
to IEEE 754, thus even a difference of 1 ulp is *not* ok.)

Paul Zimmermann

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #10 from Jakub Jelinek <jakub at gcc dot gnu.org> 2010-10-21 20:55:58 UTC ---
fsqrt insn is always used, and if the result is NaN, it calls library sqrtf
function so that errno is set correctly.  The (conditional) call causes (some
of) the values to be forced into stack and thus rounded to IEEE single
precision, if they aren't forced into stack, they will have long double
precision.
You can use -fno-errno-math if you don't want errno to be set, then there will
be no calls to sqrtf and all 3 calls should at least when optimizing evaluate
in extended precision.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[zimmerma+gcc at loria dot fr]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #11 from Paul Zimmermann <zimmerma+gcc at loria dot fr> 2010-10-22 06:56:20 UTC ---
(In reply to comment #10)
> You can use -fno-errno-math if you don't want errno to be set, then there will
> be no calls to sqrtf and all 3 calls should at least when optimizing evaluate
> in extended precision.

indeed with -fno-math-errno I get three identical results on a 64-bit Core 2
with gcc 4.4.4:


tarte% cat bug.c
#include <stdio.h>
#include <math.h>

float x=(float) 2.0;
int main() {
  printf ("%.10f\n%.10f\n%.10f\n",sqrtf(x),sqrtf(x),sqrtf(x));
  return 0;
}
tarte% gcc -mfpmath=387 bug.c -lm
tarte% ./a.out 
1.4142135382
1.4142135382
1.4142135624

tarte% gcc -mfpmath=387 -fno-math-errno bug.c -lm
tarte% ./a.out 
1.4142135624
1.4142135624
1.4142135624

However setting errno should not have side effects on the results.

Paul

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

--- Comment #12 from Jakub Jelinek <jakub at gcc dot gnu.org> 2010-10-26 13:41:10 UTC ---
Well, the setting of errno by calling another function has similar effects like
calling any other function in between, if you do
  float a = sqrtf (x);
  foo ();
  float b = sqrtf (x);
  bar ();
  float c = sqrtf (x);
  printf ("%f %f %f\n", a, b, c);
then there is similar problem - every call in between on i?86/x86_64 will
clobber the i387 register stack in the ABI, so the values need to be flushed
into memory and thus with the fast extended precision mode rounded to floating
point precision, but the c value is probably kept in reg-stack and saved as a
double instead of float, because of varargs in printf.  Thus it can have
different value.

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[rguenth at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

Richard Guenther <rguenth at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Priority|P3                          |P2

[Bug target/46080] [4.4/4.5/4.6 Regression] incorrect precision of sqrtf builtin for x87 arithmetic (-mfpmath=387)

[jakub at gcc dot gnu.org]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=46080

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |INVALID

--- Comment #13 from Jakub Jelinek <jakub at gcc dot gnu.org> 2010-11-16 10:54:36 UTC ---
Even if we add code to hardcode EDOM value for Linux targets and teach gcc that
errno is *__errno_location (), the first sqrtf could still be rounded to float
precision (because there is an optional __errno_location () call after it and
according to the ABI it might clobber %st(0)) while the remaining two could
very well be kept in extended precision and just rounded down to double
precision for varargs passing.

So I'm afraid there is nothing we can do about this in GCC and there is no GCC
bug, but simply the unpredictability of excess precision arithmetics.

Use -ffloat-store/-fexcess-precision=standard/volatile float to ensure it works
as expected if you care about it.