From owner-freebsd-numerics@FreeBSD.ORG Thu Aug 22 21:33:16 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 638BBDAF for ; Thu, 22 Aug 2013 21:33:16 +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 4522724EC for ; Thu, 22 Aug 2013 21:33:16 +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 r7MLXFgx006874 for ; Thu, 22 Aug 2013 14:33:15 -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 r7MLXFHV006873 for freebsd-numerics@freebsd.org; Thu, 22 Aug 2013 14:33:15 -0700 (PDT) (envelope-from sgk) Date: Thu, 22 Aug 2013 14:33:15 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: (2nd time) tweaks to erff() threshold values Message-ID: <20130822213315.GA6708@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.5.21 (2010-09-15) 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: Thu, 22 Aug 2013 21:33:16 -0000 (Having problems with freebsd-numerics list. Hopefully, not a duplicate.) 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(). For both the old values and the new values suggested by me, exhaustive testing in the interval [0,0x1p-11) shows % ./testf -m 0 -M 0x1p-11 -E double Max ULP: 0.77011 x Max ULP: 4.00666904e-04, 0x1.a42134p-12 Here, the assumed accurate solution is libm's erf(). For 1000000 values, uniformly distributed in the interval, speed tests show % ./testf -m 0 -M 0x1p-11 -s old: 0.024801 secs, 0.02480 usecs/call new: 0.021259 secs, 0.02126 usecs/call So, we see a slight speed up with the new values and no change in the ULP. The biggest changes comes from changing the interval 6 <= |x| < inf to 4 <= |x| < inf. Accuracy tests in the interval [2,6) give % ./testf -m 2 -M 6 -E double Max ULP: 0.51126 x Max ULP: 2.00804758e+00, 0x1.0107b4p+1 for both the old value and the new value. The improvement comes in the speed testing of 1000000 values in [2,6). % ./testf -m 2 -M 6 -s old: 0.105391 secs, 0.10539 usecs/call new: 0.061907 secs, 0.06191 usecs/call Any objections to committing the change? -- 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); -- Steve