(2nd time) tweaks to erff() threshold values

[Steve Kargl]

On Fri, Aug 23, 2013 at 09:57:44AM -0700, Steve Kargl wrote:
> On Fri, Aug 23, 2013 at 09:12:33PM +1000, Bruce Evans wrote:
> 
>> The whole erf implementation is probably not the best way for
>> float precision, but is good for understanding higher precisions.
> 
> I'm fairly certain that the implementation of erff() is not very
> efficient.  The polynomials, used in the rational approximations,
> are the same order as those used in the double precision approximation.
> I'll check the polys when update my P(x)/Q(x) remes algorithm for
> erfl and erfcl.

I seem to be right (although I haven't iterated on the P's and Q's).

erff() with the currently used pp and qq coefficients in [0,0.84375).
./testf -m 0 -M 0.84375 -E double
    Max ULP: 0.62437
  x Max ULP: 1.04175471e-38, 0x1.c5bf88p-127

./testf -m 0 -M 0.84375 -s
[0., 0.84375), 1000000 calls, 0.030670 secs, 0.03067 usecs/call



erff() with P and Q coefficients in [0,0.84375) from my the
initial iteration of my rational polynomial approximation code.

./testf -m 0 -M 0.84375 -E double
Interval tested: [0.0000,0.8438]
    Max ULP: 0.62437
  x Max ULP: 1.04175471e-38, 0x1.c5bf88p-127

./testf -m 0 -M 0.84375 -s
[0., 0.84375], 1000000 calls, 0.023354 secs, 0.02335 usecs/call

So, max ulp doesn't change and the code is faster.
The reported error estimate from solving the matrix is 

err = 7.028164541356344965e-11 0x1.351a191817c3cp-34

-- 
Steve

Index: src/s_erff.c
===================================================================
--- src/s_erff.c	(revision 1361)
+++ src/s_erff.c	(working copy)
@@ -41,6 +41,12 @@ qq2  =  6.5022252500e-02, /* 0x3d852a63 
 qq3  =  5.0813062117e-03, /* 0x3ba68116 */
 qq4  =  1.3249473704e-04, /* 0x390aee49 */
 qq5  = -3.9602282413e-06, /* 0xb684e21a */
+#define P0  1.28379166e-01F	/*  0x1.06eba8p-3 */
+#define P1 -3.36702317e-01F	/* -0x1.58c87ep-2 */
+#define P2 -1.00184719e-04F	/* -0x1.a43486p-14 */
+#define Q1  3.07090849e-01F	/*  0x1.3a7606p-2 */
+#define Q2  1.99971031e-02F	/*  0x1.47a1eep-6 */
+#define Q3 -2.07959116e-03F	/* -0x1.109380p-9 */
 /*
  * Coefficients for approximation to  erf  in [0.84375,1.25]
  */
@@ -114,8 +120,13 @@ erff(float x)
 		return x + efx*x;
 	    }
 	    z = x*x;
+#if 0
 	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
 	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
+#else
+	    r = P0 + z * (P1 + z * P2);
+	    s = 1 + z * (Q1 + z * (Q2 + z * Q3));
+#endif
 	    y = r/s;
 	    return x + x*y;
 	}