Date: Sat, 13 Nov 2010 12:56:48 +0000 From: Alexander Best <arundel@freebsd.org> To: Ulrich Spoerlein <uqs@FreeBSD.org> Cc: svn-src-head@freebsd.org, svn-src-all@freebsd.org, src-committers@freebsd.org Subject: Re: svn commit: r215237 - head/lib/msun/src Message-ID: <20101113125648.GA25183@freebsd.org> In-Reply-To: <201011131054.oADAsA7I045096@svn.freebsd.org> References: <201011131054.oADAsA7I045096@svn.freebsd.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sat Nov 13 10, Ulrich Spoerlein wrote: > Author: uqs > Date: Sat Nov 13 10:54:10 2010 > New Revision: 215237 > URL: http://svn.freebsd.org/changeset/base/215237 > > Log: > Fix bug in jn(3) and jnf(3) that led to -inf results thank you very much for fixing this long outstanding issue. you might want to have a look at [1], where two more issues have been reported. cheers. alex [1] http://mailman.oakapple.net/pipermail/numeric-interest/2010-September/thread.html > > Explanation by Steve: > jn[f](n,x) for certain ranges of x uses downward recursion to compute > the value of the function. The recursion sequence that is generated is > proportional to the actual desired value, so a normalization step is > taken. This normalization is j0[f](x) divided by the zeroth sequence > member. As Bruce notes, near the zeros of j0[f](x) the computed value > can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only > the leading decimal digit is correct. The solution to the issue is > fairly straight forward. The zeros of j0(x) and j1(x) never coincide, > so as j0(x) approaches a zero, the normalization constant switches to > j1[f](x) divided by the 2nd sequence member. The expectation is that > j1[f](x) is a more accurately computed value. > > PR: bin/144306 > Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu> > Reviewed by: bde > MFC after: 7 days > > Modified: > head/lib/msun/src/e_jn.c > head/lib/msun/src/e_jnf.c > > Modified: head/lib/msun/src/e_jn.c > ============================================================================== > --- head/lib/msun/src/e_jn.c Sat Nov 13 10:38:06 2010 (r215236) > +++ head/lib/msun/src/e_jn.c Sat Nov 13 10:54:10 2010 (r215237) > @@ -200,7 +200,12 @@ __ieee754_jn(int n, double x) > } > } > } > - b = (t*__ieee754_j0(x)/b); > + z = __ieee754_j0(x); > + w = __ieee754_j1(x); > + if (fabs(z) >= fabs(w)) > + b = (t*z/b); > + else > + b = (t*w/a); > } > } > if(sgn==1) return -b; else return b; > > Modified: head/lib/msun/src/e_jnf.c > ============================================================================== > --- head/lib/msun/src/e_jnf.c Sat Nov 13 10:38:06 2010 (r215236) > +++ head/lib/msun/src/e_jnf.c Sat Nov 13 10:54:10 2010 (r215237) > @@ -152,7 +152,12 @@ __ieee754_jnf(int n, float x) > } > } > } > - b = (t*__ieee754_j0f(x)/b); > + z = __ieee754_j0f(x); > + w = __ieee754_j1f(x); > + if (fabsf(z) >= fabsf(w)) > + b = (t*z/b); > + else > + b = (t*w/a); > } > } > if(sgn==1) return -b; else return b; -- a13x
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20101113125648.GA25183>