From owner-freebsd-numerics@freebsd.org Sat Sep 8 13:52:21 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 75F2AFF474C for ; Sat, 8 Sep 2018 13:52:21 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail105.syd.optusnet.com.au (mail105.syd.optusnet.com.au [211.29.132.249]) by mx1.freebsd.org (Postfix) with ESMTP id E87CA94F77 for ; Sat, 8 Sep 2018 13:52:20 +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 mail105.syd.optusnet.com.au (Postfix) with ESMTPS id A6A8A104B9F8; Sat, 8 Sep 2018 23:52:10 +1000 (AEST) Date: Sat, 8 Sep 2018 23:52:10 +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: <20180908044931.GA52882@troutmask.apl.washington.edu> Message-ID: <20180908232603.V883@besplex.bde.org> References: <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org> <20180905152104.GA26453@troutmask.apl.washington.edu> <20180906034525.A2959@besplex.bde.org> <20180905190223.GA27865@troutmask.apl.washington.edu> <20180906202658.K978@besplex.bde.org> <20180906223026.GA43005@troutmask.apl.washington.edu> <20180907230825.A933@besplex.bde.org> <20180907175406.GA49351@troutmask.apl.washington.edu> <20180908044931.GA52882@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=FNpr/6gs c=1 sm=1 tr=0 a=PalzARQSbocsUSjMRkwAPg==:117 a=PalzARQSbocsUSjMRkwAPg==:17 a=kj9zAlcOel0A:10 a=x7bEGLp0ZPQA:10 a=HajKcG0aLoVG3s6pNfgA: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: Sat, 08 Sep 2018 13:52:21 -0000 On Fri, 7 Sep 2018, Steve Kargl wrote: > On Fri, Sep 07, 2018 at 10:54:06AM -0700, Steve Kargl wrote: >> >> I'll probably try a different approach tomorrow. j0(x), >> cc, and ss are smoothly varying functions of x. u and v >> appear to be slowly varying (in small intervals), so we can >> have >> >> j0(x) = sqrt(1/pi) * (cc(x) *u-ss(x) *v) / sqrt(x). >> j0(x+e) = sqrt(1/pi) * (cc(x+e)*u-ss(x+e)*v) / sqrt(x+e). >> >> where e is some small perturbution. Two equation and two >> unknowns is easily solved. Simply need to do this across >> the interval to generate u(x) and v(x). > > I may have have a better approach! > > j0(x) = sqrt(1/pi) * (cc(x) * u - ss(x) * v) / sqrt(x) > y0(x) = sqrt(1/pi) * (cc(x) * u + ss(x) * v) / sqrt(x) > > j0(x) + y0(x) = 2 * sqrt(1/pi) * cc(x) * u > j0(x) - y0(x) = - 2 * sqrt(1/pi) * ss(x) * v Divided by sqrt(x) on the RHS. > So, we have > > u = (j0 + y0) / (2 * sqrt(1/pi) * cc) > v = (y0 - j0) / (s * sqrt(1/pi) * ss) Multiplied by sqrt(x) on the RHS. Also, change s back to 2. > Thus, we can find u = r/s or u = 1 + r/s and v = r/s. Good idea. I checked that the u and v calculated by this work right in pari (they give a precision almost as high as the precision used in pari when they are used to implement j0 using libm's formula in pari) (57 decimal digits with 100-digit precision in pari, even though I used the limit formula A&S 9.1.2 for y0 because I couldn't find y0 in pari (use j_nu and y_nu with nu = 1e-50 to get the 57 decimal digits of accuracy). (I just noticed that libm only supports integral nu, presumably because its recursion method only works for integers. Current versions of pari document not very good accuracy for small args and/or small nu, and I'm using an 11 year old version, but nu = 1e-50 works perfectly.) However, the u and v calculated by this have no resemblance to libm's pzero and qzero (with the latter checked by pari to give a precision of about 16 decimal digits when used to implement j0). pzero / u - 1 is between -0.5 and -3.1. qzero / v - 1 is between -1.2 and -1.0. So libm doesn't seem to be doing anything better than fixing a vague approximation to the correct u and then solving for p. Better use the simple product method. Bruce