Date: Tue, 27 Aug 2013 00:27:30 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-numerics@FreeBSD.org Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130826235215.X2165@besplex.bde.org> In-Reply-To: <20130823165744.GA47369@troutmask.apl.washington.edu> References: <20130822213315.GA6708@troutmask.apl.washington.edu> <20130823202257.Q1593@besplex.bde.org> <20130823165744.GA47369@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130826235215.X2165>