cvs commit: src/lib/msun/src e_lgammaf_r.c

[Bruce Evans]

On Sun, 4 Dec 2005, Mike Silbersack wrote:

> On Tue, 29 Nov 2005, Bruce Evans wrote:
>
>> On Tue, 29 Nov 2005, Andrey Chernov wrote:
>> 
>>> On Tue, Nov 29, 2005 at 11:49:13AM +1100, Bruce Evans wrote:
>>>> I don't like writing papers, and rarely read them these days.
>>> 
>>> Not writting the paper about your libm changes will increase chances your
>>> changes will be simple lost at some step. Possible scenario: 1) One of
>> 
>> They won't be lost, becaue they are in cvs :-).
>
> Would it be possible to write regression tests which verify that accuracy is 
> not lost due to future changes to libm?  That would be the best way to ensure 
> that your work is not lost in the future.

Whoever makes the changes would write the regression tests :-).  Mine
are't sufficently general to commit.  In batch mode which takes about
10 hours on a 2GHz Athlon to check all cases for floats, they are
currently reporting the following errors on amd64:

% acos:  max_er = 0x1cbc92d0 0.8980, avg_er = 0.170, #>=1:0.5 = 0:5422147
% acosh: max_er = 0x40018f09 2.0002, avg_er = 0.080, #>=1:0.5 = 99094:244658831

All hyperbolic functions use algorithms that aren't quite good enough
for <1 ulp precision in either float or double precision, and are not
easy to change to get that precision.

% asin:  max_er = 0x1d9cc0b9 0.9254, avg_er = 0.013, #>=1:0.5 = 0:2357938
% asinh: max_er = 0x38e6127e 1.7781, avg_er = 0.176, #>=1:0.5 = 1122862:542122336
% atan:  max_er = 0x1b44773c 0.8521, avg_er = 0.186, #>=1:0.5 = 0:6406812
% atanh: max_er = 0x371d20ea 1.7223, avg_er = 0.017, #>=1:0.5 = 1203812:52062350
% cbrt:  max_er = 0x71e2e7fe 3.5589, avg_er = 0.506, #>=1:0.5 = 502518670:1799139486

cbrt was broken on translation to float precision.  Fixed locally.

% ceil:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% ceil:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% ceil:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% ceil:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0

All cases are not actually checked except for functions like ceil() which
are supposed to work in all rounding modes.  Cases not tested mostly
work OK, especially considering that we don't try hard to make nonstandard
rounding modes always work.

% cos:   max_er = 0x10075fce 0.5009, avg_er = 0.138, #>=1:0.5 = 0:647600
% cosh:  max_er = 0x5017fcd2 2.5029, avg_er = 0.021, #>=1:0.5 = 417078:23901640
% erf:   max_er = 0x1ef93235 0.9679, avg_er = 0.127, #>=1:0.5 = 0:124418716
% exp:   max_er = 0x1d1f5c1d 0.9101, avg_er = 0.034, #>=1:0.5 = 0:17928157
% exp2:  max_er = 0x11004a89 0.5313, avg_er = 0.033, #>=1:0.5 = 0:1113072
% expm1: max_er = 0x1a026f9e 0.8128, avg_er = 0.032, #>=1:0.5 = 0:12920601
% fabs:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% fabs:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% fabs:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% fabs:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% floor: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% floor: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% floor: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% floor: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% j0:    max_er = 0x6f8686c5c9f26 3654467.3863, avg_er = 0.418, #>=1:0.5 = 671562746:1332944948
% j1:    max_er = 0x449026e9c286f 2246675.4566, avg_er = 0.423, #>=1:0.5 = 674510414:1330993378
% lgamma:max_er = 0xe4cf242400000 7497618.0703, avg_er = 0.274, #>=1:0.5 = 70033450:508702225

All bessel functions, and lgamma, use algorithms that aren't nowhere near
good enough for <1 ulp precision in either float or double precision, and
are hard to fix to get that precision.

% log:   max_er = 0x1c6a0bec 0.8879, avg_er = 0.125, #>=1:0.5 = 0:13361409
% log10: max_er = 0x42e34742 2.0902, avg_er = 0.126, #>=1:0.5 = 1129782:30059364

logs are easy to do but this one is apparently done too easily on args
near 1.

% log1p: max_er = 0x1ad7d5e0 0.8388, avg_er = 0.089, #>=1:0.5 = 0:11534114
% logb:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% logb:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% logb:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% logb:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% rint:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% rint:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% rint:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% rint:  max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% round: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% round: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% round: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% round: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% significand:max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% significand:max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% significand:max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% significand:max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% sin:   max_er = 0x10075fb5 0.5009, avg_er = 0.137, #>=1:0.5 = 0:625118
% sinh:  max_er = 0x5017fcd2 2.5029, avg_er = 0.024, #>=1:0.5 = 2182602:72386346
% sqrt:  max_er =  0xffffffc 0.5000, avg_er = 0.125, #>=1:0.5 = 0:0
% sqrt:  max_er = 0x1fffffef 1.0000, avg_er = 0.249, #>=1:0.5 = 0:1069630165
% sqrt:  max_er = 0x1fffff84 1.0000, avg_er = 0.249, #>=1:0.5 = 0:1069202731
% sqrt:  max_er = 0x1fffffef 1.0000, avg_er = 0.249, #>=1:0.5 = 0:1069630165

sqrt is supposed to be prefectly rounded in all rounding modes, and
apparently is.  This batch only tested whatever is in libm, so it
tested hardware sqrt here.  IIRC, fdlibm sqrt is prefectly rounded in
all modes too.  But the regression test only understands round-to-nearest
mode and perfect rounding can only be guessed from its output.  (Perfect
rounding for round-to-nearest is 0 in the "#>=0.5" column.).

% tan:   max_er = 0x19986d82 0.7999, avg_er = 0.150, #>=1:0.5 = 0:303818258
% tanh:  max_er = 0x4608a402 2.1886, avg_er = 0.035, #>=1:0.5 = 37533748:118674314

tgamma is omitted since tgammaf doesn't exist.

% trunc: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% trunc: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% trunc: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0
% trunc: max_er =          0 0.0000, avg_er = 0.000, #>=1:0.5 = 0:0

Testing double and long precision is harder since exhaustive testing is
impossible.  I don't have any ideas better than testing on every value
representable as a float, and on a few nearby and special values.  testing
on every 0x10000'th float indicates that sparse testing works quite well --
output like the above can be generated in closer to 1 second than overnight,
and it shows the same problem cases as above, including ones that ucbtest
doesn't find in tests that take a few minutes.

Bruce