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>