From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 31237 invoked by alias); 14 Feb 2003 12:26:00 -0000 Mailing-List: contact gcc-prs-help@gcc.gnu.org; run by ezmlm Precedence: bulk List-Archive: List-Post: List-Help: Sender: gcc-prs-owner@gcc.gnu.org Received: (qmail 31223 invoked by uid 71); 14 Feb 2003 12:26:00 -0000 Date: Fri, 14 Feb 2003 12:26:00 -0000 Message-ID: <20030214122600.31222.qmail@sources.redhat.com> To: nobody@gcc.gnu.org Cc: gcc-prs@gcc.gnu.org, From: Richard Earnshaw Subject: Re: optimization/9654: extra-precision fp comparisons are less accurate Reply-To: Richard Earnshaw X-SW-Source: 2003-02/txt/msg00610.txt.bz2 List-Id: The following reply was made to PR optimization/9654; it has been noted by GNATS. From: Richard Earnshaw To: Richard Addison-Wood 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.