From owner-freebsd-numerics@FreeBSD.ORG Fri Aug 23 16:57:46 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 52EC113D for ; Fri, 23 Aug 2013 16:57:46 +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 EAACA240E for ; Fri, 23 Aug 2013 16:57:45 +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 r7NGvj5l047467; Fri, 23 Aug 2013 09:57:45 -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 r7NGviDs047466; Fri, 23 Aug 2013 09:57:44 -0700 (PDT) (envelope-from sgk) Date: Fri, 23 Aug 2013 09:57:44 -0700 From: Steve Kargl To: Bruce Evans Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130823165744.GA47369@troutmask.apl.washington.edu> References: <20130822213315.GA6708@troutmask.apl.washington.edu> <20130823202257.Q1593@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130823202257.Q1593@besplex.bde.org> 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: Fri, 23 Aug 2013 16:57:46 -0000 On Fri, Aug 23, 2013 at 09:12:33PM +1000, Bruce Evans wrote: > On Thu, 22 Aug 2013, Steve Kargl wrote: > > > (Having problems with freebsd-numerics list. Hopefully, not > > a duplicate.) > > It was a duplicate :-). Yep. The original finally showed up. hub.freebsd.org must have ben swamped. > > 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(). > > ... > > 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 */ > > This increases the number of incorrectly rounded cases slightly on i386, > but 2**-14 decreases it slightly. This is because i386 has extra > precision for floats which the nonlinear return expression can actually > use for approx. 2**-14 < |x| < 2**-13. I suppose I should update my test programs to count the number of results that fall within ULP bins. I typically only want to know do I make the worse case ULP better or worse. 2**-14 is certainly better than 2**-28, which is the double precision threshold. > > + if (ix < 0x04000000) /* |x|<0x1p-119 */ > > /*avoid underflow */ > > - return (float)0.125*((float)8.0*x+efx8*x); > > + return (8*x+efx8*x)/8; > > Also put the comment back on the same line as the return statement now > that it fits again (see the double version). OK. > Also add "spurious" to the comment. Underflow still occurs for most denormal > x. This is non-spurious. We just avoid it for |x| larger than about 8 > times the largest denormal, by multiplying efx by 8 so that it is larger > than 1 (efx8 = 8(efx)). OK. > Some modified lines: > > @ if(ix < 0x3f580000) { /* |x|<0.84375 */ > @ if(ix < 0x38800000) { /* |x|<2**-14 */ > @ if (ix < 0x04000000) /* |x|<0x1p-119 */ > @ return (8*x+efx8*x)/8; /*avoid spurious underflow */ > @ } > @ return x + efx*x; > @ } OK. I'll update my patch to use the above. Did you see my suggested changes to erfcf()? > 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. -- Steve