From owner-freebsd-numerics@FreeBSD.ORG Mon Aug 26 14:27:43 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 7E191B23 for ; Mon, 26 Aug 2013 14:27:43 +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 4311A27D3 for ; Mon, 26 Aug 2013 14:27:43 +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 B4BAAD64236; Tue, 27 Aug 2013 00:27:31 +1000 (EST) Date: Tue, 27 Aug 2013 00:27:30 +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: <20130823165744.GA47369@troutmask.apl.washington.edu> Message-ID: <20130826235215.X2165@besplex.bde.org> 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; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=YYGEuWhf 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=EmHTJ-zS--LAX7py56IA:9 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: Mon, 26 Aug 2013 14:27:43 -0000 Catching up with old mail... On Fri, 23 Aug 2013, Steve Kargl wrote: > On Fri, Aug 23, 2013 at 09:12:33PM +1000, Bruce Evans wrote: > ... >> 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()? Yes. I didn't try them because I neglected to test erfc*() before so I didn't have anything to compare with. I thought that erfc*() was so similar to erf*() that it didn't need testing. But testing shows that it did need testing since it is quite broken. I get errors like the following, and the changes in your later mail don't affect this much: (oops, the change of the threshold from 28 to 9.194705 in this mail did affect this and gave very large errors of up to 2048 starting at the new threshold. The extra errors seem to be the result of underflow to 0 (even on i386 with some extra precision). With that backed out): Very large errors start about here: % x = 0x407b9000 3.93066406 % erfc(x) = 0x3e5d2abbd6aee678 0x32e955df 2.71638191e-08 % erfcf(x) = 0x3e5d2abdc0000000 0x32e955ee 2.71638463e-08 % err = +0x1e9511988 15.29115 This was the largest one found: % x = 0x40ffd000 7.99414062 % erfc(x) = 0x39ef467e0a9d6031 0xf7a33f0 1.23359547e-29 % erfcf(x) = 0x39ef4685e0000000 0xf7a342f 1.23360018e-29 % err = +0x7d5629fcf 62.66829 This is due to a large cancelation for subtracting 1. The fdlibm algorithm is no good. erfc() is much better in old BSD libm. It is almost the same as fdlibm except it avoids this problem: fdlibm erfc(): % z = x; % SET_LOW_WORD(z,0); % r = __ieee754_exp(-z*z-0.5625)* % __ieee754_exp((z-x)*(z+x)+R/S); % if(hx>0) return r/x; else return two-r/x; Old BSD erfc() (from FreeBSD-1.1.5): % /* return exp(-x^2 - lsqrtPI_hi + R + y)/x; */ % s = ((R + y) - lsqrtPI_hi) + z; % y = (((z-s) - lsqrtPI_hi) + R) + y; % r = exp__D(s, y)/x; % if (x>0) % return r; % else % return two-r; The extra-precision function exp__D() gives some chance of a not very large error. I think it only gives about 10 more bits of accuracy. I don't quite see how it can work for very large x when the cancelation is almost complete. Replacing the expf()'s by exp()'s in the erfc() reduced the error similarly, but it is still above 1 ulp. The old BSD erf*() also has more comments about what the error actually is, but these don't clearly distinguish erfc() from erf(). The old BSD erf*()'s file name is not spelled with an s_. > ... Bruce