optimization/9654: extra-precision fp comparisons are less accurate

optimization/9654: extra-precision fp comparisons are less accurate

[richard]

>Number:         9654
>Category:       optimization
>Synopsis:       extra-precision fp comparisons are less accurate
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    unassigned
>State:          open
>Class:          wrong-code
>Submitter-Id:   net
>Arrival-Date:   Tue Feb 11 04:06:01 UTC 2003
>Closed-Date:
>Last-Modified:
>Originator:     richard@wetafx.co.nz
>Release:        g++3 (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7)
>Organization:
>Environment:
Linux 2.4.18-xfssmp #1 SMP i686
>Description:
//////////bug.cc:
#include <math.h>

double s(double t) throw()
{
  double d = 0.75+t;
  double o = d < 1.0 ? acos(d) : 1.0/67108864.0;
  return sin(o*0.5)/sin(o);
}

int main(void)
{
  return isnan(s(0.25-1.0/18014398509481984.0));
}
//////////GNUmakefile:
bug: bug.cc; g++3 -O1 -o bug bug.cc -lm 
>How-To-Repeat:
With the two specified files, use gmake to build.

With optimization turned on, the comparison between
d and 1.0 on line 6 succeeds when it should fail.

Essentially, the value that is computed for d is
still in an fp register (with extra precision)
when it is compared with 1.0.  That value is less
then 1.0, but by the time it is actually converted
to a double, it gets rounded to 1.0.  Thus, the
situation that the comparison is trying to trap
does not get trapped.

Corresponding problems also appear with floats.

This problem also shows up with:

gcc version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
gcc version 2.96 20000731 (Red Hat Linux 7.2 2.96-112.7.2)
>Fix:
There are various ways to work around this.  Obviously,
the code can be changed slightly so that the intended
effect is correct.  However, other similar problems may
be far more subtle than generating a NaN.

The problem can be avoid using the -ffloat-store option,
but that has detrimental effects on performance.

It would be great if floating point comparisons (float
or double) were guaranteed to be on values at the
precision specified by the source code, and not on
values with extra precision.

If this cannot be the default, perhaps there should
at least be an option like -ffloat-compare so that
there is a mid-ground between the performance hit of
using -ffloat-store and the reliability hit of not
using -ffloat-store.
>Release-Note:
>Audit-Trail:
>Unformatted:

Re: optimization/9654: extra-precision fp comparisons are less accurate

[ebotcazou]

Synopsis: extra-precision fp comparisons are less accurate

State-Changed-From-To: open->closed
State-Changed-By: ebotcazou
State-Changed-When: Sun Feb 16 15:30:11 2003
State-Changed-Why:
    Not a bug.

http://gcc.gnu.org/cgi-bin/gnatsweb.pl?cmd=view%20audit-trail&database=gcc&pr=9654

Re: optimization/9654: extra-precision fp comparisons are less accurate

[Richard Earnshaw]

The following reply was made to PR optimization/9654; it has been noted by GNATS.

From: Richard Earnshaw <rearnsha@arm.com>
To: Richard Addison-Wood <richard@wetafx.co.nz>
Cc: rearnsha@gcc.gnu.org, gcc-gnats@gcc.gnu.org, gcc-bugs@gcc.gnu.org
Subject: Re: optimization/9654: extra-precision fp comparisons are less 
 accurate
Date: Fri, 14 Feb 2003 12:21:27 +0000

 > Of course, this discourse into acos(), sin(),
 > and sqrt() is beside the point.  With code
 > like this:
 > 
 >   if (d < 1.0)
 >   {
 >     arglessthanone(d);
 >   }
 > 
 > it should be the case that the arglessthanone()
 > function should not see an argument value that
 > is greater than or equal to 1.0.
 
 If we were talking about the pure mathematical operations then I would 
 agree with you.  But we aren't.  We are talking about a finite 
 representation of the operations.
 
 d can be less than one, but to such a limited extent that by the time we 
 have called into arglessthanone() it's representation has been reduced to 
 the same as one.  If your algorithm is sensitive to this then it should 
 really be considered unstable.  Going back to your original code, you had:
 
 double s(double t) throw()
 {
   double d = 0.75+t;
   double o = d < 1.0 ? acos(d) : 1.0/67108864.0;
   return sin(o*0.5)/sin(o);
 }
 
 
 If we look at the Taylor expansion for sin(x) we can express the above as
 
 	o/2 - ((o/2)^3)/3! + ....
         ------------------------
         o - (o^3)/3! + ...
 
 Now clearly this can be simplified to
 
        1/2 - (o^2)/48 + ...
        --------------------
        1 - (o^2)/6 + ...
 
 From this it should be clear that if o is less than DBL_EPSILON the answer 
 is going to be 0.5 (since o^2 is much smaller than can be represented as a 
 difference from 0.5, let alone from 1), so we can now rewrite your 
 function as
 
 double s(double t) throw
 {
   double d = 0.75+t;
   double o = d < 1.0 ? acos(d) : 1.0/67108864.0;
   if (o < DBL_EPSILON)
     return 0.5;		// Limit value for small o
   return sin(o*0.5) / sin(o);
 }    
 
 So what does all this show? 
 
 1) That we need to be careful if there is a possibility that we may be 
 dealing with a discontinuity in an intermediate expression
 2) Careful analysis shows us that often we don't need to go anywhere near 
 the discontinuity in order to get accurate answers; we just need to take a 
 bit more care about how we calculate them.
 
 R.
 

Re: optimization/9654: extra-precision fp comparisons are less accurate

[Richard Addison-Wood]

The following reply was made to PR optimization/9654; it has been noted by GNATS.

From: Richard Addison-Wood <richard@wetafx.co.nz>
To: rearnsha@gcc.gnu.org
Cc: gcc-gnats@gcc.gnu.org, gcc-bugs@gcc.gnu.org
Subject: Re: optimization/9654: extra-precision fp comparisons are less accurate
Date: Fri, 14 Feb 2003 16:48:38 +1300

 It's quite reasonable to assume that
 
   acos(1-0.5*DBL_EPSILON)
 
 will not be 0.  For a value of d that is
 nearly 1:
 
   acos(d) == sqrt(1-d*d)
 
 is a good approximation (certainly within
 representational limits ieee754 double
 precision).  The point to notice here is that
 if d is distinguishable from 1, d*d will be
 even more distinguishable from 1.  All told,
 it very reasonable to expect:
 
   acos(1-0.5*DBL_EPSILON) == sqrt(DBL_EPSILON)
 
 which is a value that is not particularly
 close to 0.  However, it is still close
 enough to 0 so that:
 
   sin(sqrt(DBL_EPSILON)) == sqrt(DBL_EPSILON)
 
 is a good approximation.
 
 Of course, this discourse into acos(), sin(),
 and sqrt() is beside the point.  With code
 like this:
 
   if (d < 1.0)
   {
     arglessthanone(d);
   }
 
 it should be the case that the arglessthanone()
 function should not see an argument value that
 is greater than or equal to 1.0.

Re: optimization/9654: extra-precision fp comparisons are less accurate

[rearnsha]

The following reply was made to PR optimization/9654; it has been noted by GNATS.

From: <rearnsha@gcc.gnu.org>
To: richard@wetafx.co.nz, gcc-gnats@gcc.gnu.org, gcc-bugs@gcc.gnu.org
Cc:  
Subject: Re: optimization/9654: extra-precision fp comparisons are less accurate
Date: Wed, 12 Feb 2003 17:49:17 +0000

 http://gcc.gnu.org/cgi-bin/gnatsweb.pl?cmd=view%20audit-trail&database=gcc&pr=9654
 
 I think it's very dangerous to assume that
 
 acos(something-nearly-equal-to-1)
 
 won't return 0, or, more importantly
 
 sin (acos(something-nearly-equal-to-1)
 
 won't.
 
 A simple test after the assignment of o would suffice to avoid this 
 problem, eg:
 
 #include <math.h>
 
 double s(double t) throw()
 {
    double d = 0.75+t;
    double o = d < 1.0 ? acos(d) : 1.0/67108864.0;
    if (o == 0.0) o = 1.0 / 67108864.0;
    return sin(o*0.5)/sin(o);
 }
 
 int main(void)
 {
    return isnan(s(0.25-1.0/18014398509481984.0));
 }
 
 If you are really paranoid, then a range check on small o and applying 
 L'Hoptial's Rule might not be a bad idea.
 
 I really don't think you'd like the alternatives to -ffloat-store any 
 more than you like -ffloat-store.