Skip site navigation (1)Skip section navigation (2)
Date:      Thu, 22 Aug 2013 12:33:29 -0700
From:      Steve Kargl <sgk@troutmask.apl.washington.edu>
To:        freebsd-numerics@freebsd.org
Subject:   tweaks to erff() threshold values
Message-ID:  <20130822193329.GA81453@troutmask.apl.washington.edu>

Next in thread | Raw E-Mail | Index | Archive | Help
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);



Want to link to this message? Use this URL: <http://docs.FreeBSD.org/cgi/mid.cgi?20130822193329.GA81453>