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

Richard Earnshaw rearnsha@arm.com
Fri Feb 14 12:22:00 GMT 2003


> 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.



More information about the Gcc-bugs mailing list