cvs commit: src/sys/i386/include _types.h

[Bruce Evans]

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