From owner-freebsd-numerics@freebsd.org  Wed Sep  5 12:19:11 2018
Return-Path: <owner-freebsd-numerics@freebsd.org>
Delivered-To: freebsd-numerics@mailman.ysv.freebsd.org
Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2610:1c1:1:606c::19:1])
 by mailman.ysv.freebsd.org (Postfix) with ESMTP id 422BBFEE745
 for <freebsd-numerics@mailman.ysv.freebsd.org>;
 Wed,  5 Sep 2018 12:19:11 +0000 (UTC)
 (envelope-from brde@optusnet.com.au)
Received: from mail104.syd.optusnet.com.au (mail104.syd.optusnet.com.au
 [211.29.132.246])
 by mx1.freebsd.org (Postfix) with ESMTP id B76887E8C2
 for <freebsd-numerics@freebsd.org>; Wed,  5 Sep 2018 12:19:10 +0000 (UTC)
 (envelope-from brde@optusnet.com.au)
Received: from [192.168.0.102] (c110-21-101-228.carlnfd1.nsw.optusnet.com.au
 [110.21.101.228])
 by mail104.syd.optusnet.com.au (Postfix) with ESMTPS id B410A42AC83;
 Wed,  5 Sep 2018 22:19:07 +1000 (AEST)
Date: Wed, 5 Sep 2018 22:19:06 +1000 (EST)
From: Bruce Evans <brde@optusnet.com.au>
X-X-Sender: bde@besplex.bde.org
cc: Steve Kargl <sgk@troutmask.apl.washington.edu>, 
 freebsd-numerics@freebsd.org
Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2)
In-Reply-To: <20180905201540.D1142@besplex.bde.org>
Message-ID: <20180905221116.X1704@besplex.bde.org>
References: <20180903235724.GA95333@troutmask.apl.washington.edu>
 <20180905201540.D1142@besplex.bde.org>
MIME-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
X-Optus-CM-Score: 0
X-Optus-CM-Analysis: v=2.2 cv=I9sVfJog c=1 sm=1 tr=0
 a=PalzARQSbocsUSjMRkwAPg==:117 a=PalzARQSbocsUSjMRkwAPg==:17
 a=9cW_t1CCXrUA:10 a=kj9zAlcOel0A:10 a=x7bEGLp0ZPQA:10
 a=0GgSnPAGSJmx6d4KNYIA:9 a=CjuIK1q_8ugA:10
X-BeenThere: freebsd-numerics@freebsd.org
X-Mailman-Version: 2.1.27
Precedence: list
List-Id: "Discussions of high quality implementation of libm functions."
 <freebsd-numerics.freebsd.org>
List-Unsubscribe: <https://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: <https://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
 <mailto:freebsd-numerics-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Wed, 05 Sep 2018 12:19:11 -0000

On Wed, 5 Sep 2018, Bruce Evans wrote:

> On Mon, 3 Sep 2018, Steve Kargl wrote:
>
>> Anyone know where the approximations for j0 (and y0) come from?
>
> I think they are ordinary minimax rational approximations for related
> functions.  As you noticed, the asymptotic expansion doesn't work below
> about x = 8 (it is off by about 10% for j0(2).  But we want to use the

... j0(2)).

> single formula given by the asymptotic expansion for all the subintervals:
>
> XX 	/*
> XX 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
> XX 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
> XX 	 */
> XX 		if(ix>0x48000000) z = (invsqrtpi*cc)/sqrt(x);
> XX 		else {
> XX 		    u = pzero(x); v = qzero(x);
> XX 		    z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
> XX 		}
> XX 		return z;
> XX 	}
>
> where pzero(s) is nominally 1 - 9/128 s**2 + ... and qzero(s) is nominally
> -1/8 *s + ....

These polynomials are actually only part of the numerators of pzero and
qzero (s = 1/x already gives a non-polynomial, and even more divisions are
used in pzero = 1 + R/S...).

> To work, pzero and qzero must not actually be these nominal functions.
> ...

Bruce