From owner-freebsd-numerics@FreeBSD.ORG  Tue Dec  2 21:43:31 2014
Return-Path: <owner-freebsd-numerics@FreeBSD.ORG>
Delivered-To: freebsd-numerics@freebsd.org
Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115])
 (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits))
 (No client certificate requested)
 by hub.freebsd.org (Postfix) with ESMTPS id CA7D26EC
 for <freebsd-numerics@freebsd.org>; Tue,  2 Dec 2014 21:43:31 +0000 (UTC)
Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu
 [128.95.76.21])
 (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits))
 (Client CN "troutmask", Issuer "troutmask" (not verified))
 by mx1.freebsd.org (Postfix) with ESMTPS id AC03E86D
 for <freebsd-numerics@freebsd.org>; Tue,  2 Dec 2014 21:43:31 +0000 (UTC)
Received: from troutmask.apl.washington.edu (localhost [127.0.0.1])
 by troutmask.apl.washington.edu (8.14.9/8.14.9) with ESMTP id sB2LhPA3094929
 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO)
 for <freebsd-numerics@freebsd.org>; Tue, 2 Dec 2014 13:43:25 -0800 (PST)
 (envelope-from sgk@troutmask.apl.washington.edu)
Received: (from sgk@localhost)
 by troutmask.apl.washington.edu (8.14.9/8.14.9/Submit) id sB2LhPCh094928
 for freebsd-numerics@freebsd.org; Tue, 2 Dec 2014 13:43:25 -0800 (PST)
 (envelope-from sgk)
Date: Tue, 2 Dec 2014 13:43:25 -0800
From: Steve Kargl <sgk@troutmask.apl.washington.edu>
To: freebsd-numerics@freebsd.org
Subject: bug in j0f()
Message-ID: <20141202214325.GA94909@troutmask.apl.washington.edu>
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
User-Agent: Mutt/1.5.23 (2014-03-12)
X-BeenThere: freebsd-numerics@freebsd.org
X-Mailman-Version: 2.1.18-1
Precedence: list
List-Id: "Discussions of high quality implementation of libm functions."
 <freebsd-numerics.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/options/freebsd-numerics>, 
 <mailto:freebsd-numerics-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-numerics/>
List-Post: <mailto:freebsd-numerics@freebsd.org>
List-Help: <mailto:freebsd-numerics-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
 <mailto:freebsd-numerics-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Tue, 02 Dec 2014 21:43:31 -0000

Anyone object to the following patch?

Index: e_j0f.c
===================================================================
--- e_j0f.c	(revision 275211)
+++ e_j0f.c	(working copy)
@@ -62,7 +62,7 @@
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
+		if(ix>0x4b800000) z = (invsqrtpi*cc)/sqrtf(x);
 		else {
 		    u = pzerof(x); v = qzerof(x);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);


Clearly, the value 0x8000000 was a mistake in the original
translation from double to float.  In looking over other
codebases, one finds

OpenBSD
		if(ix>0x80000000U) z = (invsqrtpi*cc)/sqrtf(x);

NetBSD

#ifdef DEAD_CODE
		if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
		else
#endif
		{

glibc
                if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);

OpenBSD is clearly wrong.  NetBSD suppresses a code optimization.
glibc simply copied the test from e_j0.c to e_j0f.c.  The glibc
choice is somewhat too small (although I haven't done any extensive
testing, yet).

prn(0x48000000); /* 0x1p17  = 1.31072000e+05 */
prn(0x4b800000); /* 0x1p24  = 1.67772160e+07 */

-- 
Steve