[Bug math/22906] new version of exp doesn't always round correctly

[woodard at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=22906

--- Comment #3 from Ben Woodard <woodard at redhat dot com> ---
Glibc on x86_64 managed 1/2ULP which is practically correctly rounded. 
In my testing of the new exp function, I found 22 cases where it was not
correctly rounded.

I made a mode in my little test program to filter out all the cases where the
answer is within 1ULP. All of the answers, were within the 1ULP

$ ./compare_vals -t -c 200000 input-decks/exp-inputs-hard
/home/ben/Work/glibc/build/glibc/math/libm.so 
Between: 0x1.0000015823067p+0 - 0x1.0000015884ae7p+0    1 bad values    max
error: 0x1p-51
Between: 0x1.ffffffffcf2ep-1 - 0x1.0000000030dap+0      1 bad values    max
error: 0x1p-51
Between: 0x1.0fffffffcf2cp+6 - 0x1.100000004f42ep+6     1 bad values    max
error: 0x1p+46
Between: 0x1.cbffffffcf2cp+6 - 0x1.cc00000058d1ep+6     1 bad values    max
error: 0x1p+113
Between: 0x1.62001f9860ef3p+9 - 0x1.62001f98e88bep+9    1 bad values    max
error: 0x1p+969
Between: 0x1.89ffffffcf2cp+7 - 0x1.8a00000044405p+7     2 bad values    max
error: 0x1p+232
Between: 0x1.ebffffffcf2cp+6 - 0x1.ec00000035fbep+6     1 bad values    max
error: 0x1p+125
Between: 0x1.0dffffffca266p+7 - 0x1.0e00000030d4p+7     1 bad values    max
error: 0x1p+142
Between: 0x1.1bffffffb4856p+8 - 0x1.1c00000037bffp+8    2 bad values    max
error: 0x1p+357
Between: 0x1.62002fffd6992p+9 - 0x1.620030019d156p+9    1 bad values    max
error: 0x1p+969
Between: 0x1.ffffffffa6609p-1 - 0x1.000000038abdp+0     4 bad values    max
error: 0x1p-51
Between: 0x1.6200315aa0262p+9 - 0x1.6200315fa7813p+9    3 bad values    max
error: 0x1p+969
Between: 0x1.ffffffffcf2cp-54 - 0x1.0000000039273p-53   1 bad values    max
error: 0x1p-52
Between: 0x1.ffffffffcf7cp-54 - 0x1.000000003562p-53    1 bad values    max
error: 0x1p-52
Between: 0x1.00000003df24cp+0 - 0x1.0000000729649p+0    4 bad values    max
error: 0x1p-51
-------
exp(x) = glibc, /home/ben/Work/glibc/build/glibc/math/libm.so
exp(0x1.00000002502f2p+0) = 0x1.5bf0a8b46a533p+1, 0x1.5bf0a8b46a532p+1
difference: 0x1p-51
exp(0x1.0000000266056p+0) = 0x1.5bf0a8b488011p+1, 0x1.5bf0a8b488012p+1
difference: 0x1p-51
exp(0x1.000000027bdbap+0) = 0x1.5bf0a8b4a5af1p+1, 0x1.5bf0a8b4a5afp+1
difference: 0x1p-51
exp(0x1.000000041cf9ep+0) = 0x1.5bf0a8b6dc9b3p+1, 0x1.5bf0a8b6dc9b2p+1
difference: 0x1p-51
exp(0x1.00000004552a6p+0) = 0x1.5bf0a8b728f9bp+1, 0x1.5bf0a8b728f9cp+1
difference: 0x1p-51
exp(0x1.00000004f2923p+0) = 0x1.5bf0a8b7fee97p+1, 0x1.5bf0a8b7fee96p+1
difference: 0x1p-51
exp(0x1.00000006c9b25p+0) = 0x1.5bf0a8ba7f3cdp+1, 0x1.5bf0a8ba7f3cep+1
difference: 0x1p-51
exp(0x1.0000015853da7p+0) = 0x1.5bf0aa8542d27p+1, 0x1.5bf0aa8542d26p+1
difference: 0x1p-51
exp(0x1.620030003123fp+9) = 0x1.58f0ae7a25891p+1021, 0x1.58f0ae7a2589p+1021
difference: 0x1p+969
exp(0x1.6200315b10d8dp+9) = 0x1.58f455448a857p+1021, 0x1.58f455448a858p+1021
difference: 0x1p+969
exp(0x1.6200315bd5b13p+9) = 0x1.58f4575707c8bp+1021, 0x1.58f4575707c8ap+1021
difference: 0x1p+969
exp(0x1.6200315ce81efp+9) = 0x1.58f45a3a9a9c3p+1021, 0x1.58f45a3a9a9c4p+1021
difference: 0x1p+969
exp(0x1p-53) = 0x1.0000000000001p+0, 0x1p+0 difference: 0x1p-52
exp(0x1.ffffffffd7349p-1) = 0x1.5bf0a8b129bd9p+1, 0x1.5bf0a8b129bd8p+1
difference: 0x1p-51
exp(0x1.100000001e6eep+6) = 0x1.12fec76184ad3p+98, 0x1.12fec76184ad2p+98
difference: 0x1p+46
exp(0x1.62001f98b7b7ep+9) = 0x1.58c47c6081583p+1021, 0x1.58c47c6081582p+1021
difference: 0x1p+969
exp(0x1.cc00000027fdep+6) = 0x1.e1032031f2bafp+165, 0x1.e1032031f2baep+165
difference: 0x1p+113
exp(0x1.8a000000136c5p+7) = 0x1.284d3bbd648f7p+284, 0x1.284d3bbd648f6p+284
difference: 0x1p+232
exp(0x1.ec0000000527ep+6) = 0x1.5e1152fb6b48fp+177, 0x1.5e1152fb6b49p+177
difference: 0x1p+125
exp(0x1.0dfffffffafa6p+7) = 0x1.b2afab7b3c48bp+194, 0x1.b2afab7b3c48ap+194
difference: 0x1p+142
exp(0x1.1bffffffe5596p+8) = 0x1.a741e7506e9b1p+409, 0x1.a741e7506e9bp+409
difference: 0x1p+357
exp(0x1.1c00000006ebfp+8) = 0x1.a741e787f06b9p+409, 0x1.a741e787f06b8p+409
difference: 0x1p+357
22/908 failed

filtering for 1ULP:

$ ./compare_vals -t1 -c 200000 input-decks/exp-inputs-hard
/home/ben/Work/glibc/build/glibc/math/libm.so 
-------
exp(x) = glibc, /home/ben/Work/glibc/build/glibc/math/libm.so
0/908 failed

If we look at intel's math library it is even worse:
$ ./compare_vals -t input-decks/exp-inputs-hard
/usr/tce/packages/intel/intel-18.0.1/lib/intel64_lin/libimf.so
<snip - see attached file if you care>
104/908 failed

But once again if I filter out all the values which are within 1ULP rather than
1/2ULP then there are no errors:

$ ./compare_vals -t1 input-decks/exp-inputs-hard
/usr/tce/packages/intel/intel-18.0.1/lib/intel64_lin/libimf.so
-------
exp(x) = glibc, /usr/tce/packages/intel/intel-18.0.1/lib/intel64_lin/libimf.so
0/908 failed

As for performance:

The original libc:
$ ./time_math2 input-decks/exp-inputs-hard 

10 Ranges: 
1-99 c:99       m:1547  r:1437:1607 
100-116 c:17    m:2003  r:1990:2097 
117-214 c:98    m:2710  r:2433:2832 
215-357 c:143   m:3373  r:3141:3757 
361-546 c:186   m:301475        r:230357:406579 
547-558 c:12    m:466542        r:425866:501855 
559-565 c:7     m:539577        r:538734:540095 
566-891 c:326   m:2425683       r:2198307:3031765 
892-905 c:14    m:3636118       r:3225192:3953467 
906-907 c:1     m:4838777       r:4838777:4841767 

Intel's math library:

$ ./time_math2 -A
/usr/tce/packages/intel/intel-18.0.1/lib/intel64_lin/libimf.so -a exp
input-decks/exp-inputs-hard 

12 Ranges: 
1-12 c:12       m:349   r:347:357 
16-100 c:85     m:563   r:537:639 
101-132 c:32    m:820   r:809:881 
133-137 c:5     m:981   r:977:994 
139-141 c:3     m:1258  r:1222:1278 
142-657 c:516   m:1567  r:1448:1774 
658-815 c:158   m:1952  r:1898:2078 
816-818 c:3     m:2368  r:2297:2410 
819-821 c:3     m:2587  r:2548:2613 
822-825 c:4     m:3280  r:3204:3355 
826-905 c:80    m:3804  r:3553:3900 
906-907 c:1     m:7765  r:7765:9851 

crlibm

$ ./time_math2 -A ../crlibm/.libs/libcrlibm.so.0.0.0 -a exp_rn
input-decks/exp-inputs-hard 

9 Ranges: 
2-99 c:98       m:791   r:681:955 
100-119 c:20    m:2001  r:1997:2013 
120-260 c:141   m:3694  r:3243:4327 
261-263 c:3     m:12988         r:12974:13012 
265-360 c:96    m:17260         r:15246:20282 
369-372 c:4     m:119446        r:114979:121839 
373-376 c:4     m:138689        r:133139:146041 
381-394 c:14    m:229740        r:205505:246519 
395-907 c:512   m:353069        r:263795:568649 

and the new libm version of exp:
$ ./time_math2 -A ../glibc/build/math/libm.so.6 -a exp
input-decks/exp-inputs-hard

5 Ranges: 
1-102 c:102     m:1257  r:1231:1315 
103-105 c:3     m:1413  r:1394:1430 
107-160 c:54    m:1669  r:1595:1954 
161-905 c:745   m:2631  r:2101:3078 
906-907 c:1     m:3467  r:3467:4316

-- 
You are receiving this mail because:
You are on the CC list for the bug.