Date: Sat, 8 Sep 2018 23:52:10 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Message-ID: <20180908232603.V883@besplex.bde.org> In-Reply-To: <20180908044931.GA52882@troutmask.apl.washington.edu> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180908232603.V883>