tweaks to erff() threshold values

Steve Kargl sgk at troutmask.apl.washington.edu
Fri Aug 23 03:26:35 UTC 2013


On Thu, Aug 22, 2013 at 12:33:29PM -0700, Steve Kargl wrote:
> 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().
> 

The thresholds in erfcf() are also not well-chosen.
Here's a new diff fixing both erff() and erfcf().

-- 
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);
@@ -160,7 +160,7 @@
 	}
 
 	if(ix < 0x3f580000) {		/* |x|<0.84375 */
-	    if(ix < 0x23800000)  	/* |x|<2**-56 */
+	    if(ix < 0x33800000)  	/* |x|<2**-24*/
 		return one-x;
 	    z = x*x;
 	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
@@ -184,7 +184,7 @@
 		z = erx+P/Q; return one+z;
 	    }
 	}
-	if (ix < 0x41e00000) {		/* |x|<28 */
+	if (ix < 0x41131d83) {		/* |x|<9.194705 */
 	    x = fabsf(x);
  	    s = one/(x*x);
 	    if(ix< 0x4036DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
@@ -193,7 +193,7 @@
 	        S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
 				sa5+s*(sa6+s*(sa7+s*sa8)))))));
 	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
-		if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
+		if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
 	        R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
 				rb5+s*rb6)))));
 	        S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(


More information about the freebsd-numerics mailing list