Use of C99 extra long double math functions after r236148
Stephen Montgomery-Smith
stephen at missouri.edu
Sun Aug 12 23:11:04 UTC 2012
On 07/20/2012 04:19 AM, Bruce Evans wrote:
> On Fri, 20 Jul 2012, Bruce Evans wrote:
> The only way for x*x + y*y to be _very_ near 1 in infinite precision
> is for |x| or y| to be 1 (I think). Other cases are bounded away from
> 1, and if you are lucky the bound is fairly far from 1 so that sloppier
> approximations work OK. Mathematicians should determine the bound
> exactly using continued fractions or something like they do for
> approximations to N*Pi/2. This becomes especially interesting in high
> precisions where you can't hope to get near the worst case by random
> testing.
>
> % x0 = ldexp(floor(ldexp(x, 24)), -24);
> % x1 = x - x0;
> % y0 = ldexp(floor(ldexp(y, 24)), -24);
> % y1 = y - y0;
>
> This has a chance of working iff the bound away from 1 is something like
> 2**-24. Otherwise, multiplying by 2**24 and flooring a positive value
> will just produce 0. 2**-24 seems much too small a bound. My test
> coverage is not wide enough to hit many bad cases.
This is meant to cover a situation where x = cos(t) and y = sin(t) for
some t not a multiple of PI/2.
Now, hypot(x,y) will be 1, but only to within machine precision, i.e. an
error of about 1e-17.
So log(hypot(x,y)) will be about 1e-17. The true answer being 0, the
ULP will be infinite.
BUT (and this goes with Goldberg's paper when he considers the quadratic
formula when the quadratic equation has nearly equal roots), suppose
x = (double)(cos(t))
that is, x is not exactly cos(t), but it is the number "cos(t) written
in IEEE double precision". Similarly for y. That is, even though the
formula that produce x and y isn't exact, let's pretend that x and y are
exact.
Again log(hypot(x,y)) will be about 1e-17. But the true answer will
also be about 1e-17. But they won't be the same, and the ULP will be
about 1e17.
What my formula does is deal with the second case, but reduce the ULP to
about 1e8! That is, if x and y are exact numbers, and it so happens
that hypot(x,y) is very close to 1, my method will get you about 8 extra
digits of accuracy.
Now you have special formulas that handle the cases when z is close to
1, -1, I and -I. Earlier this morning I sent you a formula, which I
think might be slightly more accurate than yours, for when z is close to
1. I think similar formulas can be produced for when z is close to -1,
I and -I.
To get ULP of about 1 when x and y are exact, and it happens that
hypot(x,y) is close to 1, but z is not close to 1, -1, I or -I, would
require, I think, hypot(x,y)-1 being computed using double double
precision (i.e. a mantissa of 108 bits), and then feeding this into log1p.
Of course, you may argue that situations when x and y are exact, not
close to 1 or 0, and hypot(x,y) is close to 1, are so very rare that
extra consideration is not required.
But my algorithm produces better answers than the naive formula even
when the distance between 1 and hypot is about 1/10. The naive formula
has a ULP of about 10, and I get it down to less than 2. And when the
distance between hypot(x,y) and 1 is about 1e-5, the naive formula has a
ULP of about 1e5, and I still manage to get a ULP of about less than 2.
Stephen
More information about the freebsd-numerics
mailing list