The state of glibc libm

["Joseph S. Myers" <joseph at codesourcery dot com>]

I've reviewed many (not yet all) of glibc's open "math" component bugs.  I 
hope some actual summary information on what the current state of libm 
looks like may be of interest to the people involved in the various past 
discussions of better libm for GCC or glibc - and those interested in 
fixing things, whether through patches to existing code, new 
implementations or both.

I would say the actual main issues with the present glibc libm are the 
following (if there are others, they are not well-reflected in the open 
bugs).  When I refer to a function by name I'm generally referred to all 
of the float, double and various long double implementations (x86 
extended, IEEE binary128, IBM double+double) - and fixes will often need 
to fix multiple versions of a function in similar ways.

(a) Most libm functions are not correctly rounded - and do not make an 
attempt at being correctly rounded.

A full fix would likely require new (automatically generated and tuned) 
implementations such as proposed at 
<http://gcc.gnu.org/ml/gcc/2012-02/msg00298.html>.  As I understand it, 
correct rounding (proved correct) is generally feasible for functions of 
one binary32, binary64 or x86 extended argument.  For functions of two 
arguments, or one binary128 argument, it may not be feasible to search for 
worst cases for correct rounding, although it may be possible to produce 
implementations that are "probably" correctly rounding.  For functions of 
complex arguments or IBM long double, correct rounding may be less 
feasible (even complex multiplication and division, and all of +-*/ on IBM 
long double, are not correctly rounding, and correct rounding isn't so 
well-defined for IBM long double with its possibility of discontiguous 
mantissa bits).

Known inaccuracies in functions are indicated by the libm-test-ulps files 
in glibc.  Given my patch 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00770.html> to reduce some 
old figures that look like having been left behind after tests or the 
library were fixed, all listed errors are 24ulp or below (9ulp or below 
for x86 and x86_64; I haven't tested the larger errors on other 
architectures to see if they are actually still applicable).

(b) Where functions do make attempts at being correctly rounded 
(especially the IBM Accurate Mathematical Library functions), they tend to 
be sufficiently slow that the slowness attracts bug reports.  Again, this 
would likely be addressed by new implementations that use careful error 
bounds and information about worst cases to reduce the cost of being 
correctly rounding.

(c) Various functions do not set errno correctly (many cases) or raise the 
proper floating-point exceptions (a smaller number of cases - both 
spurious exceptions where not permitted by ISO C, and failing to raise 
required overflow/underflow exceptions).  In general this is a separate 
bug for each function (filed as many separate bugs in glibc Bugzilla) and 
can be fixed by a separate local patch for each function (adding a 
testcase, of course - note that glibc's main libm-test.inc presently only 
tests invalid and divide-by-zero exceptions, so if working on these error 
handling issues it might be useful to extend it to cover other exceptions 
as well as errno values).

(d) There are some specific bugs filed for functions such as nexttoward 
whose results are precisely specified by ISO C; in general these should be 
fixable by local patches.

(e) Various functions, mainly IBM Accurate Mathematical Library ones, 
produce wildly wrong results or crash in non-default rounding modes.  I 
have a patch for exp pending review at 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00748.html> and I expect 
others can be fixed similarly.

(f) Various trigonometrical functions are inaccurate on x86 and x86_64 
(see glibc bug 13658 and recent discussions).

(g) Bessel function implementations handle large inputs in different ways 
to other functions, as I discuss at 
<http://sourceware.org/ml/libc-alpha/2012-02/msg00512.html>.

(h) Various complex functions have problems such as inaccuracy or wrong 
branch cuts.

(i) Some real functions have particular issues (which should be fixable by 
local changes short of new correctly rounded implementations):

  - erfc (my patch 
    <http://sourceware.org/ml/libc-alpha/2012-02/msg00640.html> is pending 
    review).

  - pow (bugs 369, 706, 2678, 3866).  The assembly implementations may 
    complicate fixing these issues, though it's probably possible to fix 
    only some bugs (in all relevant implementations, with plenty of 
    testcases) rather than a patch needing to fix both all issues and all 
    implementations at once.

  - lgamma, in cases where the result is close to 0 and there is a lot of 
    cancellation in the present calculations.

  - tgamma, in cases of results of large magnitude (where the approach of 
    using exp (lgamma) leads to large errors).

-- 
Joseph S. Myers
joseph@codesourcery.com