On Thu, 6 Mar 2008, Colin Percival wrote:
Bruce Evans wrote:
On Wed, 5 Mar 2008, Colin Percival wrote:
You could have stopped this sentence here -- for all practical purposes,
correctly rounded trigonometric functions are not feasible.
Nah, it's quite feasible, especially for transcendental functions since
transcendental numbers in general can be approximated more closely by
rationals than algebraic numbers.
In fact, one of the common ways to prove that a real number is transcendental
is to show that it can be approximated more closely by rationals than any
algebraic function can be.
I've only studied a few pages related to this.
The evaluation just needs to be
done in extra precision
How much extra precision do you need? If the exact value of a transcendental
function is extremely close to halfway between two consecutive floating-point
values, how do you decide which way to round it?
Probably not much more than the same number of extra bits as there are in
the original precision. Values of trig functions mod 2^N are apparently
sort of uniformly distributed. Thus each extra bit of precision gives
correct rounding for about half of the remaining cases and N extra bits
gives correct rounding for about 2^N cases, which is all cases if there
were only 2^N cases to start with. This works out in practice for float
precision -- with about 24 bits of extra precision, all except about 10
cases are perfectly rounded, and it is easy to evaluate those cases in
extra precision to see that only a few more bits are needed.
If the exact value is very close to half-way then you may need many
more bits but would be unlucky to need many more. AFAIK (not far),
no one knows exactly how many bits are needed for the worst case for
even double precision. The domain is just too large to search (almost
2^64 bits). However, the large size of the domain means that very-near-
half-way cases are very unlikely to occur in practice, so calling a
slow multi-precision library to handle such cases would be a good way
to handle them and sending mail about them would be a not so good way.
But detecting such cases is too costly to do routinely, at least without
hardware support.
Computing transcendental functions to within (0.5+x)ulp for small positive
values of x is certainly feasible, but computing provably correctly rounded
transcendental functions is very very hard.
I think it is only very hard globally and for general functions.
Elementary transcendental functions like exp() have rational coefficients
and well-known properties which make them easier to bound.
Fortunately, library functions aren't required to have any particular
error bounds.
We require error bounds of < 1 ulp and aim for < 0.5[0]1 ulps.
What standard states that those bounds are required (or can be relied upon
when proving that code is correct)?
Just the principle of least error.
Intel FPUs over the years have been
variously documented to compute transcendental functions to within 1.0, 2.5,
or 3.5 ulp -- are you saying that the 8087 doesn't comply with the standard
which it inspired?
I didn't say that before, but in fact the ix87 has always had multi-gulp
errors where its docs indicate errors of < 1 ulp, mainly for trig
functions. Its trig functions can't possibly have errors of < 1 ulp,
since its approximation to pi has only 66-68 bits, but an approximation
with thousands of bits is needed (see libm). Multi-gulp errors occur
as early as for x near pi/2 in 64-bit precision, since subtracting a
64-bit x near pi/2 from a 66 or 68 bit approx to pi/2 gives about 2
or 4 bits of accuracy or about 62 or 60 bits of error. 2^60 is about
1 exa-ulp (10^9 gulps). The errors are much larger for larger x, and
give multi-gulp errors even for float precision for not very large x.
Backwards compatibility has apparently resulted in these errors being
perfectly preserved in all implementations of the ix87.
Other ix87 errors are relatively minor AFAIK. E.g., exp(x) is not
implemented directly in hardware. It requires several imprecise
calculations like x*log2(e), so the hardware is good enough for double
precision but not for long double precision where I guess it has an
error of a few ulps. It has to switch the mode to long double precision
to be good enough for double precision. log(x) is implemented more
directly in hardware and I guess it has an error of about 1 ulp in long
double precision.
Bruce
_______________________________________________
cvs-all@freebsd.org mailing list
http://lists.freebsd.org/mailman/listinfo/cvs-all
To unsubscribe, send any mail to "[EMAIL PROTECTED]"
