From owner-freebsd-numerics@FreeBSD.ORG Mon Aug 19 11:06:48 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id 37553DA0 for ; Mon, 19 Aug 2013 11:06:48 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:1900:2254:206c::16:87]) (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 23FD7251D for ; Mon, 19 Aug 2013 11:06:48 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.7/8.14.7) with ESMTP id r7JB6mIB006101 for ; Mon, 19 Aug 2013 11:06:48 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.7/8.14.7/Submit) id r7JB6lOx006099 for freebsd-numerics@FreeBSD.org; Mon, 19 Aug 2013 11:06:47 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 19 Aug 2013 11:06:47 GMT Message-Id: <201308191106.r7JB6lOx006099@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-numerics@FreeBSD.org Subject: Current problem reports assigned to 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: Mon, 19 Aug 2013 11:06:48 -0000 Note: to view an individual PR, use: http://www.freebsd.org/cgi/query-pr.cgi?pr=(number). The following is a listing of current problems submitted by FreeBSD users. These represent problem reports covering all versions including experimental development code and obsolete releases. S Tracker Resp. Description -------------------------------------------------------------------------------- o stand/175811 numerics libstdc++ needs complex support in order use C99 o bin/170206 numerics [msun] [patch] complex arcsinh, log, etc. o stand/82654 numerics C99 long double math functions are missing 3 problems total. From owner-freebsd-numerics@FreeBSD.ORG Thu Aug 22 19:33:39 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 592BB52A for ; Thu, 22 Aug 2013 19:33:39 +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 3ACB8289B for ; Thu, 22 Aug 2013 19:33:39 +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 r7MJXTRi082870 for ; Thu, 22 Aug 2013 12:33:29 -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 r7MJXTm9082869 for freebsd-numerics@freebsd.org; Thu, 22 Aug 2013 12:33:29 -0700 (PDT) (envelope-from sgk) Date: Thu, 22 Aug 2013 12:33:29 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: tweaks to erff() threshold values Message-ID: <20130822193329.GA81453@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 19:33:39 -0000 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); 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 From owner-freebsd-numerics@FreeBSD.ORG Fri Aug 23 03:26:35 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 DFC76F2D for ; Fri, 23 Aug 2013 03:26:35 +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 A67E0294B for ; Fri, 23 Aug 2013 03:26:35 +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 r7N3QWeA043340 for ; Thu, 22 Aug 2013 20:26:32 -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 r7N3QWsZ043339 for freebsd-numerics@freebsd.org; Thu, 22 Aug 2013 20:26:32 -0700 (PDT) (envelope-from sgk) Date: Thu, 22 Aug 2013 20:26:32 -0700 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: tweaks to erff() threshold values Message-ID: <20130823032632.GA43322@troutmask.apl.washington.edu> References: <20130822193329.GA81453@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130822193329.GA81453@troutmask.apl.washington.edu> 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: Fri, 23 Aug 2013 03:26:35 -0000 On Thu, Aug 22, 2013 at 12:33:29PM -0700, Steve Kargl wrote: > 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(). > The thresholds in erfcf() are also not well-chosen. Here's a new diff fixing both erff() and erfcf(). -- 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); @@ -160,7 +160,7 @@ } if(ix < 0x3f580000) { /* |x|<0.84375 */ - if(ix < 0x23800000) /* |x|<2**-56 */ + if(ix < 0x33800000) /* |x|<2**-24*/ return one-x; z = x*x; r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); @@ -184,7 +184,7 @@ z = erx+P/Q; return one+z; } } - if (ix < 0x41e00000) { /* |x|<28 */ + if (ix < 0x41131d83) { /* |x|<9.194705 */ x = fabsf(x); s = one/(x*x); if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ @@ -193,7 +193,7 @@ S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( sa5+s*(sa6+s*(sa7+s*sa8))))))); } else { /* |x| >= 1/.35 ~ 2.857143 */ - if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */ + if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */ R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( rb5+s*rb6))))); S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( From owner-freebsd-numerics@FreeBSD.ORG Fri Aug 23 11:12:48 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id C9D518D1 for ; Fri, 23 Aug 2013 11:12:48 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail109.syd.optusnet.com.au (mail109.syd.optusnet.com.au [211.29.132.80]) by mx1.freebsd.org (Postfix) with ESMTP id 44E092EBC for ; Fri, 23 Aug 2013 11:12:48 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail109.syd.optusnet.com.au (Postfix) with ESMTPS id CD9D9D62E77; Fri, 23 Aug 2013 21:12:34 +1000 (EST) Date: Fri, 23 Aug 2013 21:12:33 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl Subject: Re: (2nd time) tweaks to erff() threshold values In-Reply-To: <20130822213315.GA6708@troutmask.apl.washington.edu> Message-ID: <20130823202257.Q1593@besplex.bde.org> References: <20130822213315.GA6708@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=DstvpgP+ c=1 sm=1 tr=0 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 a=PO7r1zJSAAAA:8 a=F6uwjYSw4mAA:10 a=kj9zAlcOel0A:10 a=JzwRw_2MAAAA:8 a=bRjWgPqadKYA:10 a=0q7wHxVP66BF9WyAtZ0A:9 a=9vR6l-q_0Ft0s3_N:21 a=6bBav436BQXeBPxM:21 a=CjuIK1q_8ugA:10 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 11:12:48 -0000 On Thu, 22 Aug 2013, Steve Kargl wrote: > (Having problems with freebsd-numerics list. Hopefully, not > a duplicate.) It was 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(). > ... > 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. > + 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). 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)). A recent discussion in comp.std.c made me more aware that underflow is mishandled for most functions in libm. They need to use something like the above to handle it correctly. For example, sin(x) returns x with inexact for denormal x, but it should return x with inexact AND UNDERFLOW for denormal x, since x is only an approximation to sin(x) -- the exact result is slightly smaller in magnitude in x, so is "tiny" when x is, and underflow means tiny plus inexact. The above still fails to raise underflow if 8*x, efx8*x, and their sum are all exact. 8*x is exact if x has 3 lower bits 0, which is quite often, and the others are rarely exact, except on i386 without clang everything is always exact since it is in extra precision. In nonstandard rounding modes, sin(x) shouldn't be rounded to x, but fdlibm and FreeBSD don't try to support such rounding modes for most functions including sin(). If they were supported, then they would increase the complications for underflow. IEEE x54 allows 2 alternatives for underflow. One is when the result is tiny after rounding, and inexact, and the other is when the result is tiny before rounding, and inexact. This causes complications even in round-to-nearest mode. Consider sin(x) on the largest denormal. The result of x is tiny both before and after rounding, and inexact, so underflow should be rasied. Consider sin(x) on the smallest normal. The result of x is normal after rounding, so is not tiny, but before rounding it is smaller than the smallest normal so should probably be considered as "tiny", so underflow should be raised iff it depends on tinyness before rounding. The discussion on comp.std.c was mainly about how extra precision complicates things further. When an expression like (8*x+efx8*x)/8 is evaluated in extra precision, nothing is tiny, so underflow is not raised. Inexact is probably not raised either for this particular expression. Underflow should be raised iff the extra precision is discarded, but the information for raising it has been partly lost. It will only be raised on conversion if the result is "tiny" then (before or after rounding) and is inexact THEN. It is inexact then iff any of the low bits that are discarded is nonzero. But sometimes the extra-precision result makes all these bits zero and discards nonzero bits even further out. The delayed conversion cannot do the right thing since it doesn't know about the bits already discarded. C standards get this wrong by having strict requirements for raising underflow. It is very difficult to raise it exactly when specified. C standards don't have strict requirements for raising inexact because they know about this problem for inexact. But underflow is even harded to decide than inexact. Implementations that support trapping of overflows have additional complications. IEEE x54 specifies that if traps for underflow are enabled, then underflow traps occur iff the result is tiny, while the underflow exception is raised iff the result is tiny and inexact. Thus, underflow traps, if enabled, occur more often than underflow exceptions. The trapped case is simpler but different (if the trap handler handler just returns). C standards don't support trapping of FP exceptions, and are even further from supporting returning from trap handlers. > return x + efx*x; > } > z = x*x; > ... 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; @ } The whole erf implementation is probably not the best way for float precision, but is good for understanding higher precisions. Bruce 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 From owner-freebsd-numerics@FreeBSD.ORG Fri Aug 23 17:15:27 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 13FDA758 for ; Fri, 23 Aug 2013 17:15:27 +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 CE955252C for ; Fri, 23 Aug 2013 17:15:26 +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 r7NHFQqn047750; Fri, 23 Aug 2013 10:15:26 -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 r7NHFQrY047749; Fri, 23 Aug 2013 10:15:26 -0700 (PDT) (envelope-from sgk) Date: Fri, 23 Aug 2013 10:15:26 -0700 From: Steve Kargl To: Bruce Evans Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130823171526.GA47736@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 17:15:27 -0000 On Fri, Aug 23, 2013 at 09:12:33PM +1000, Bruce Evans wrote: > > @ 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 */ To keep the diff between s_erf.c and s_erff.c small, do you want my to update s_erf.c with the last line above? -- Steve 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; }