Re: [ballen@uwm.edu] libc/2269: triginometric argument reductionerror in libm

[Stephen L Moshier <moshier at mediaone dot net>]

On 20 May 2001, Andreas Jaeger wrote:
> here's a report about problems with the trigonmetric range reduction
> in glibc.
> 
> Any comments or suggestions?



> The IEEE standard clearly states that math functions
> should return the closest IEEE float to the true correct value.

No, I don't think IEEE 754 requires anything at all about math functions
beyond arithmetic and the square root.


> /usr/lib/libc/sysdeps/ieee754/flt-32/e_rem_pio2f.c

This is not the function that is invoked on x86. The one actually used
is sysdeps/i386/fpu/s_sinf.S.  The fdlibm function e_rem_pio2f.c would
be used on sparcstations or other computers for which no hardware
function routine has been provided.

For Intel's hardware routine the typical absolute error in sin(x) increases
linearly with the argument, approximately as 4e-22 x.  When the result
is rounded to double precision, the function suffers partial loss of
precision (> 1 ulp absolute) for x >= 2^18, and total loss of precision (50%
of full scale) for x >= 2^65.

Intel's hardware behaves that way on purpose.  It is not a mistake.
The design reflects an engineering tradeoff choice among speed,
accuracy, memory requirement, and the size of your electricity bill.
The reason it makes sense arises from an analysis of the error
propagation.  In real life, function arguments are usually not exact
values but are inexact.  If the result of your input being off by half
a bit is that the function value changes 100%, that is not a good
situation.  Your obelisk, bridge, or skyscraper will fall down.  That
is a reason to go back to the drawing board and think of another way
to solve the problem.

As for the value judgments, I personally wouldn't classify the
Intel behavior as a bug.  I do think the user should be offered a
choice of routines, if that is possible.  I have previously advocated
that glibc should provide the user with a range of choices for math
libraries.  The response was that the maintainers do not physically
have the time and energy to offer more than the support that exists
now.  A possible compromise in this case might be to test the input
argument and, if larger than the partial-loss-of-precision threshold,
to call Sun's rem_pio2 routine (or perhaps the new "IBM Accurate
Mathematical Library" function that might recently have replaced Sun's
program).  But I haven't looked to see how much work that would
entail, and I don't think it is worth a high priority.  I think that a
better project would be to provide a complete, separate, fdlibm and/or
IBM library in a form that can be used on i386.  Whether such
libraries would be distributed with libc is a management or political
decision.  There is no technical reason why it would need to be part of libc.

Regarding s_sinf.S, the function is checking for Intel's total loss of
precision flag and it jumps to some range reduction code on that flag.
But the code looks like a futile gesture since the built-in value of
pi is no more accurate than the one that already failed.  I offer
no opinion on whether this is worth fixing.  There are numerous
design errors like this in libm, and on my private list I put
rewriting all the single precision fdlibm functions at a higher priority.