From owner-freebsd-numerics@FreeBSD.ORG Sat Aug 24 20:21:08 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id BF18A218 for ; Sat, 24 Aug 2013 20:21:08 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id 83A132165 for ; Sat, 24 Aug 2013 20:21:08 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.7/8.14.7) with ESMTP id r7OKL2Yj055014; Sat, 24 Aug 2013 13:21:02 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.7/8.14.7/Submit) id r7OKL2Nd055013; Sat, 24 Aug 2013 13:21:02 -0700 (PDT) (envelope-from sgk) Date: Sat, 24 Aug 2013 13:21:02 -0700 From: Steve Kargl To: Bruce Evans Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130824202102.GA54981@troutmask.apl.washington.edu> References: <20130822213315.GA6708@troutmask.apl.washington.edu> <20130823202257.Q1593@besplex.bde.org> <20130823165744.GA47369@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130823165744.GA47369@troutmask.apl.washington.edu> User-Agent: Mutt/1.5.21 (2010-09-15) Cc: freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 24 Aug 2013 20:21:08 -0000 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; }