Date: Thu, 22 Aug 2013 14:33:15 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: (2nd time) tweaks to erff() threshold values Message-ID: <20130822213315.GA6708@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
(Having problems with freebsd-numerics list. Hopefully, not a duplicate.) I would like to apply the patch, which follows my .sig, to msun/src/s_erff.c. It tweaks the threshold values to those appropriate for erff(). For both the old values and the new values suggested by me, exhaustive testing in the interval [0,0x1p-11) shows % ./testf -m 0 -M 0x1p-11 -E double Max ULP: 0.77011 x Max ULP: 4.00666904e-04, 0x1.a42134p-12 Here, the assumed accurate solution is libm's erf(). For 1000000 values, uniformly distributed in the interval, speed tests show % ./testf -m 0 -M 0x1p-11 -s old: 0.024801 secs, 0.02480 usecs/call new: 0.021259 secs, 0.02126 usecs/call So, we see a slight speed up with the new values and no change in the ULP. The biggest changes comes from changing the interval 6 <= |x| < inf to 4 <= |x| < inf. Accuracy tests in the interval [2,6) give % ./testf -m 2 -M 6 -E double Max ULP: 0.51126 x Max ULP: 2.00804758e+00, 0x1.0107b4p+1 for both the old value and the new value. The improvement comes in the speed testing of 1000000 values in [2,6). % ./testf -m 2 -M 6 -s old: 0.105391 secs, 0.10539 usecs/call new: 0.061907 secs, 0.06191 usecs/call Any objections to committing the change? -- Steve Index: src/s_erff.c =================================================================== --- src/s_erff.c (revision 1358) +++ src/s_erff.c (working copy) @@ -107,10 +107,10 @@ } if(ix < 0x3f580000) { /* |x|<0.84375 */ - if(ix < 0x31800000) { /* |x|<2**-28 */ - if (ix < 0x04000000) + if(ix < 0x39000000) { /* |x|<2**-13 */ + if (ix < 0x04000000) /* |x|<0x1p-119 */ /*avoid underflow */ - return (float)0.125*((float)8.0*x+efx8*x); + return (8*x+efx8*x)/8; return x + efx*x; } z = x*x; @@ -125,7 +125,7 @@ Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); if(hx>=0) return erx + P/Q; else return -erx - P/Q; } - if (ix >= 0x40c00000) { /* inf>|x|>=6 */ + if (ix >= 0x40800000) { /* inf>|x|>=4 */ if(hx>=0) return one-tiny; else return tiny-one; } x = fabsf(x); -- Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130822213315.GA6708>