Date: Thu, 22 Aug 2013 20:26:32 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: Re: tweaks to erff() threshold values Message-ID: <20130823032632.GA43322@troutmask.apl.washington.edu> In-Reply-To: <20130822193329.GA81453@troutmask.apl.washington.edu> References: <20130822193329.GA81453@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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*(
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130823032632.GA43322>