Re: libm trigonometric functions and i387

[Sverre Froyen <sverre%viewmark.com@localhost>]

On Tuesday 24 February 2009 09:59:24 Taylor R Campbell wrote:
> The i387 trigonometric instructions are broken by design: they return
> wildly incorrect values for arguments outside a very small range,
> because they reduce arguments based on an approximation of pi that is
> too coarse.  For example, on 2^62, the FSIN instruction gives about
> -0.70713 whereas the correct answer is closer to -0.70292.  Helpfully,
> when the arguments exceed about 2^63 in absolute value, the i387
> instructions set a processor flag that means `sorry, this hardware is
> too broken to give any remotely sensible answer'.  What NetBSD's libm
> does with this flag doesn't help the situation, though -- it proceeds
> to use other i387 instructions to reduce arguments in exactly the same
> way, with exactly the same approximation of pi, that the i387 said
> wouldn't give a remotely sensible answer.  For the sine of 1e22, this
> yields about 0.8740281, whereas the correct answer is closer to
> -0.8522008 (!).
>
> Looking at the source to NetBSD's libm, I see that there is one
> directory src/lib/libm/arch/i387/ for i387 miscomputations, and
> another directory src/lib/libm/src/ with C code to compute what
> appears to be better argument reduction based on 1584 bits of 2/pi,
> rather than the 66 or so bits of pi that the i387 uses.  However,
> /etc/ld.so.conf persuades the linker to favour libm387 over libm.
> This, I think, is a mistake, first of all: programmers who want their
> programs to use the i387's fast and loose trigonometric routines can
> always compile their programs with -ffast-math and -mfancy-math-387,
> but right now the default is to give wildly bogus answers.  Setting
> aside, though, the issue of whether NetBSD's libm should be fast or
> correct, how do I tell the linker that I want libm, not libm387,
> without editing /etc/ld.so.conf?

How do you accurately pass an argument of 2^63 to libm's trig functions?  A 
quick look suggests your only choices are float or double.  For an argument 
value of 2^63, the uncertainty exceeds 2*pi by several orders of magnitute.

Sverre