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