Date: Thu, 8 Dec 2005 14:50:55 +1100 (EST) From: Bruce Evans <bde@zeta.org.au> To: Mike Silbersack <silby@silby.com> Cc: src-committers@freebsd.org, Andre Oppermann <andre@freebsd.org>, cvs-src@freebsd.org, cvs-all@freebsd.org, Bruce Evans <bde@freebsd.org>, Steve Kargl <sgk@troutmask.apl.washington.edu>, Andrey Chernov <ache@freebsd.org> Subject: Re: cvs commit: src/lib/msun/src e_lgammaf_r.c Message-ID: <20051208141941.L63825@delplex.bde.org> In-Reply-To: <20051204231731.N43418@odysseus.silby.com> References: <200511280832.jAS8WGvs059057@repoman.freebsd.org> <438AD8FB.A8B96AB6@freebsd.org> <20051128172718.GA59929@troutmask.apl.washington.edu> <20051129110058.T33820@delplex.bde.org> <20051129012102.GA84108@nagual.pp.ru> <20051129184901.I34802@delplex.bde.org> <20051204231731.N43418@odysseus.silby.com>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051208141941.L63825>