Date: Sat, 24 Aug 2013 13:21:02 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@FreeBSD.org Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130824202102.GA54981@troutmask.apl.washington.edu> In-Reply-To: <20130823165744.GA47369@troutmask.apl.washington.edu> References: <20130822213315.GA6708@troutmask.apl.washington.edu> <20130823202257.Q1593@besplex.bde.org> <20130823165744.GA47369@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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;
}
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130824202102.GA54981>