Date: Wed, 5 Sep 2018 22:06:29 +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: <20180905201540.D1142@besplex.bde.org> In-Reply-To: <20180903235724.GA95333@troutmask.apl.washington.edu> References: <20180903235724.GA95333@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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: 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 + .... To work, pzero and qzero must not actually be these nominal functions. A non-tricky implementation would use : z = pzero(x) / qzero(x); where pzero and qzero depend on the subinterval like now and have all the other terms folded into them. We already do this for |x| < 2 (except we spell the rational function r/s and don't fold all of the terms directly into r and s). This might be the best implementation, and I would also try using smaller intervals with polynomial approximations. 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. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180905201540.D1142>