From owner-svn-src-all@FreeBSD.ORG Sat Nov 13 10:54:10 2010 Return-Path: Delivered-To: svn-src-all@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 8455B106564A; Sat, 13 Nov 2010 10:54:10 +0000 (UTC) (envelope-from uqs@FreeBSD.org) Received: from svn.freebsd.org (svn.freebsd.org [IPv6:2001:4f8:fff6::2c]) by mx1.freebsd.org (Postfix) with ESMTP id 57D358FC08; Sat, 13 Nov 2010 10:54:10 +0000 (UTC) Received: from svn.freebsd.org (localhost [127.0.0.1]) by svn.freebsd.org (8.14.3/8.14.3) with ESMTP id oADAsAHs045099; Sat, 13 Nov 2010 10:54:10 GMT (envelope-from uqs@svn.freebsd.org) Received: (from uqs@localhost) by svn.freebsd.org (8.14.3/8.14.3/Submit) id oADAsA7I045096; Sat, 13 Nov 2010 10:54:10 GMT (envelope-from uqs@svn.freebsd.org) Message-Id: <201011131054.oADAsA7I045096@svn.freebsd.org> From: Ulrich Spoerlein Date: Sat, 13 Nov 2010 10:54:10 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org X-SVN-Group: head MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cc: Subject: svn commit: r215237 - head/lib/msun/src X-BeenThere: svn-src-all@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "SVN commit messages for the entire src tree \(except for " user" and " projects" \)" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 13 Nov 2010 10:54:10 -0000 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 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 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;