From owner-freebsd-numerics@FreeBSD.ORG Tue Dec 2 21:43:31 2014 Return-Path: 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 ; 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 ; 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 ; 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 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-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 From owner-freebsd-numerics@FreeBSD.ORG Wed Dec 3 00:09:43 2014 Return-Path: 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 2630BEF4 for ; Wed, 3 Dec 2014 00:09:43 +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 05EE1A4A for ; Wed, 3 Dec 2014 00:09:43 +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 sB309ftS098477 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Tue, 2 Dec 2014 16:09:41 -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 sB309fKR098476 for freebsd-numerics@freebsd.org; Tue, 2 Dec 2014 16:09:41 -0800 (PST) (envelope-from sgk) Date: Tue, 2 Dec 2014 16:09:41 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: bug in j0f() Message-ID: <20141203000941.GA98467@troutmask.apl.washington.edu> References: <20141202214325.GA94909@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20141202214325.GA94909@troutmask.apl.washington.edu> 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 03 Dec 2014 00:09:43 -0000 On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote: > 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); Exhaustive testing in the range 0x1p38 to 0x1p100 indicated at the constant should be 0x54000000. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Wed Dec 3 00:29:09 2014 Return-Path: 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 5F6523A4 for ; Wed, 3 Dec 2014 00:29:09 +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 3E546C1D for ; Wed, 3 Dec 2014 00:29:09 +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 sB30T8LC098601 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Tue, 2 Dec 2014 16:29:08 -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 sB30T8g3098600 for freebsd-numerics@freebsd.org; Tue, 2 Dec 2014 16:29:08 -0800 (PST) (envelope-from sgk) Date: Tue, 2 Dec 2014 16:29:08 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: bug in j0f() Message-ID: <20141203002908.GA98589@troutmask.apl.washington.edu> References: <20141202214325.GA94909@troutmask.apl.washington.edu> <20141203000941.GA98467@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20141203000941.GA98467@troutmask.apl.washington.edu> 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 03 Dec 2014 00:29:09 -0000 On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote: > On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote: > > 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); > > Exhaustive testing in the range 0x1p38 to 0x1p100 > indicated at the constant should be 0x54000000. > Note, a similar wrong condition exists within y0f(). I have not tested y0f(), but propose making a similar change in y0f() as well. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Wed Dec 3 02:32:08 2014 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) (using TLSv1.2 with cipher AECDH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTPS id C6FEB926 for ; Wed, 3 Dec 2014 02:32:08 +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 A47DAAA1 for ; Wed, 3 Dec 2014 02:32:08 +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 sB32W7BQ099069 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO) for ; Tue, 2 Dec 2014 18:32:07 -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 sB32W7o5099068 for freebsd-numerics@freebsd.org; Tue, 2 Dec 2014 18:32:07 -0800 (PST) (envelope-from sgk) Date: Tue, 2 Dec 2014 18:32:07 -0800 From: Steve Kargl To: freebsd-numerics@freebsd.org Subject: Re: bug in j0f() Message-ID: <20141203023207.GA99054@troutmask.apl.washington.edu> References: <20141202214325.GA94909@troutmask.apl.washington.edu> <20141203000941.GA98467@troutmask.apl.washington.edu> <20141203002908.GA98589@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20141203002908.GA98589@troutmask.apl.washington.edu> 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 03 Dec 2014 02:32:08 -0000 On Tue, Dec 02, 2014 at 04:29:08PM -0800, Steve Kargl wrote: > On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote: > > On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote: > > > 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); > > > > Exhaustive testing in the range 0x1p38 to 0x1p100 > > indicated at the constant should be 0x54000000. > > > > Note, a similar wrong condition exists within y0f(). I have > not tested y0f(), but propose making a similar change in y0f() > as well. > While I'm monologuing, I might as well point out that the rational approximations in j0f (and y0f and most likely j1f and y1f) are over-specified. I suspect that the polynomials in the rational approximation can be reduced by one or two terms. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Wed Dec 3 13:08:30 2014 Return-Path: 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 1DBEB712 for ; Wed, 3 Dec 2014 13:08:30 +0000 (UTC) Received: from mail104.syd.optusnet.com.au (mail104.syd.optusnet.com.au [211.29.132.246]) by mx1.freebsd.org (Postfix) with ESMTP id D3C402DE for ; Wed, 3 Dec 2014 13:08:29 +0000 (UTC) Received: from c122-106-147-133.carlnfd1.nsw.optusnet.com.au (c122-106-147-133.carlnfd1.nsw.optusnet.com.au [122.106.147.133]) by mail104.syd.optusnet.com.au (Postfix) with ESMTPS id 2429E42555E; Thu, 4 Dec 2014 00:08:19 +1100 (AEDT) Date: Thu, 4 Dec 2014 00:08:19 +1100 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl Subject: Re: bug in j0f() In-Reply-To: <20141203023207.GA99054@troutmask.apl.washington.edu> Message-ID: <20141203233540.Q44095@besplex.bde.org> References: <20141202214325.GA94909@troutmask.apl.washington.edu> <20141203000941.GA98467@troutmask.apl.washington.edu> <20141203002908.GA98589@troutmask.apl.washington.edu> <20141203023207.GA99054@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.1 cv=dMCfxopb c=1 sm=1 tr=0 a=7NqvjVvQucbO2RlWB8PEog==:117 a=PO7r1zJSAAAA:8 a=kj9zAlcOel0A:10 a=JzwRw_2MAAAA:8 a=cVtEa_Lb80mRl8RVK4oA:9 a=CjuIK1q_8ugA:10 Cc: freebsd-numerics@freebsd.org 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 03 Dec 2014 13:08:30 -0000 On Tue, 2 Dec 2014, Steve Kargl wrote: > On Tue, Dec 02, 2014 at 04:29:08PM -0800, Steve Kargl wrote: >> On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote: >>> On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote: >>>> Anyone object to the following patch? OK (with 0x54000000). >>>> 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); >>> >>> Exhaustive testing in the range 0x1p38 to 0x1p100 >>> indicated at the constant should be 0x54000000. My tests agree. Tested on amd64 and i386. >> Note, a similar wrong condition exists within y0f(). I have >> not tested y0f(), but propose making a similar change in y0f() >> as well. Not so exhaustive testing gave 0x54800000 on amd64. > While I'm monologuing, I might as well point out that the > rational approximations in j0f (and y0f and most likely > j1f and y1f) are over-specified. I suspect that the > polynomials in the rational approximation can be reduced > by one or two terms. Also, the cutoffs of 2**-13 and 2**-27 are the same for both precisions, thus likely to be wrong for float precision. Bruce From owner-freebsd-numerics@FreeBSD.ORG Wed Dec 3 15:49:46 2014 Return-Path: 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 636DA5DD for ; Wed, 3 Dec 2014 15:49:46 +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 40D96C31 for ; Wed, 3 Dec 2014 15:49:46 +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 sB3FnidJ001986 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Wed, 3 Dec 2014 07:49:44 -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 sB3FniT4001985; Wed, 3 Dec 2014 07:49:44 -0800 (PST) (envelope-from sgk) Date: Wed, 3 Dec 2014 07:49:44 -0800 From: Steve Kargl To: Bruce Evans Subject: Re: bug in j0f() Message-ID: <20141203154944.GA1976@troutmask.apl.washington.edu> References: <20141202214325.GA94909@troutmask.apl.washington.edu> <20141203000941.GA98467@troutmask.apl.washington.edu> <20141203002908.GA98589@troutmask.apl.washington.edu> <20141203023207.GA99054@troutmask.apl.washington.edu> <20141203233540.Q44095@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20141203233540.Q44095@besplex.bde.org> User-Agent: Mutt/1.5.23 (2014-03-12) Cc: freebsd-numerics@freebsd.org 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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 03 Dec 2014 15:49:46 -0000 On Thu, Dec 04, 2014 at 12:08:19AM +1100, Bruce Evans wrote: > On Tue, 2 Dec 2014, Steve Kargl wrote: > > > On Tue, Dec 02, 2014 at 04:29:08PM -0800, Steve Kargl wrote: > >> On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote: > >>> On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote: > >>>> Anyone object to the following patch? > > OK (with 0x54000000). > > >>>> 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); > >>> > >>> Exhaustive testing in the range 0x1p38 to 0x1p100 > >>> indicated at the constant should be 0x54000000. > > My tests agree. Tested on amd64 and i386. > > >> Note, a similar wrong condition exists within y0f(). I have > >> not tested y0f(), but propose making a similar change in y0f() > >> as well. > > Not so exhaustive testing gave 0x54800000 on amd64. > Thanks for confirming the values and suggestion the above for y0f(). > > While I'm monologuing, I might as well point out that the > > rational approximations in j0f (and y0f and most likely > > j1f and y1f) are over-specified. I suspect that the > > polynomials in the rational approximation can be reduced > > by one or two terms. > > Also, the cutoffs of 2**-13 and 2**-27 are the same for both precisions, > thus likely to be wrong for float precision. > I haven't gotten that far into analyzing the code. I'll probably get there in a few months. :-) -- Steve