From owner-freebsd-numerics@freebsd.org Wed Sep 5 18:09:17 2018 Return-Path: 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 EA99CFF8AB8 for ; Wed, 5 Sep 2018 18:09:16 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail106.syd.optusnet.com.au (mail106.syd.optusnet.com.au [211.29.132.42]) by mx1.freebsd.org (Postfix) with ESMTP id 695EF8CA2E for ; Wed, 5 Sep 2018 18:09:16 +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 mail106.syd.optusnet.com.au (Postfix) with ESMTPS id 226823CC636; Thu, 6 Sep 2018 04:09:05 +1000 (AEST) Date: Thu, 6 Sep 2018 04:09:05 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl cc: freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) In-Reply-To: <20180905152104.GA26453@troutmask.apl.washington.edu> Message-ID: <20180906034525.A2959@besplex.bde.org> References: <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org> <20180905152104.GA26453@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.2 cv=DZtnkrlW c=1 sm=1 tr=0 a=PalzARQSbocsUSjMRkwAPg==:117 a=PalzARQSbocsUSjMRkwAPg==:17 a=kj9zAlcOel0A:10 a=x7bEGLp0ZPQA:10 a=sdg16lVg9fCeTce6EDsA: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." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Sep 2018 18:09:17 -0000 On Wed, 5 Sep 2018, Steve Kargl wrote: > On Wed, Sep 05, 2018 at 10:06:29PM +1000, 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 >> single formula given by the asymptotic expansion for all the subintervals: > > I've scoured the literature and web for methods of computing > Bessel functions. These functions are important to my real > work. I have not found any paper, webpage, documentation, etc. > that describes what "the related functions" are. They are just the functions in the asymptotic expansion with errors corrected as I discussed. >> ... >> However, we use the asympotic formula for all cases above x = 2. It is >> off by at least 10% unless pzero and qzero are adjusted. But 10% is not >> much, and adjusting only pzero and qzero apparently works to compensate >> for this. To use the same formula for all cases, the complicated factors >> cc, ss and sqrt(x) must not be folded into pzero and qzero; instead these >> factors supply some of the errors which must be compensated for. >> >> To find the adjustments, I would try applying the same scale factor to >> the nominal pzero and qzero. E.g., for j0, the error factor is >> invsqrtpi*(u*cc-v*ss)/sqrt(x) / j0(x). Divide the nominal pzero and qzero >> by this. Apply normal minimax methods to the modfided pzero and qzero. > > Obviously, I had not thought about scaling P0(x) and Q0(x) with a > common factor to force agreement between j0(x) and its approximation. I think I understand why we use the asymptotic formula as a base. If we choose a minimax function for each interval, this would need to have a very high order. Much the same order as the rational function part of the asymptotic expansion. But the latter has a clever grouping of terms using cos and sin and in the internals of P0 and Q0. This makes it faster to calculate and also probably has better accuracy. Bruce