Date: Thu, 6 Sep 2018 15:30:26 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Message-ID: <20180906223026.GA43005@troutmask.apl.washington.edu> In-Reply-To: <20180906202658.K978@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>
next in thread | previous in thread | raw e-mail | index | archive | help
On Fri, Sep 07, 2018 at 02:02:20AM +1000, Bruce Evans wrote: > On Wed, 5 Sep 2018, Steve Kargl wrote: > > > On Thu, Sep 06, 2018 at 04:09:05AM +1000, Bruce Evans wrote: > >> On Wed, 5 Sep 2018, Steve Kargl wrote: > >>> 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. > > > > And as I noted, there is no documentation stating the approximations > > pzero(x) and qzero(x) aren't approximations for the asymptotic series > > P0(x) and Q0(x). If you are correct, then pzero(x) and qzero(x) are > > approximations to fudge*P0(x) and fudge*Q0(x). What fudge is and how > > it is determined is not documented. > > The documentation (comments in the source code) is indeed deficient. It > is too verbose about routine general methods but says little about non- > routine things. > > I now think that no single fudge factor works for both P0 and Q0. P0 > and Q0 have very different characteristics. Merely choosing the > truncation point for the asymptotic expansion makes them very different. > One thing that can go wrong is that if you truncate to more than about > 5 terms, a zero point for one or both of P0 and Q0 ends up in the > subinterval being handled. The fudge factor would need to be nearly > infinite to compensate for the zero, and since P0 and Q0 won't have > the zero at the same place, the compensation for one would destroy the > other. > > The infinities show up in attempts to calculate the fudge factor near the > zero. > > When P0 is truncated after the s^14 term, then on [2, 2.857] P0 (1/x) > is strictly increasing from -5.61 to 0.95, with the zero in the middle > giving singularities, but pzero is strictly increasing from 1 - 0.013 to > 1 - 0.007. However, when P0 is truncated after the s^6 term, it is > similar to pzero (strictly increasing from 1 - 0.019 to 1 - 0.009). > pzero is apparently based on the latter P0. > > To find different fudge factors, I would try starting with a common one > and then iteratively adjust to different ones. > > j0 has a zero in the middle of at least the first subinterval, and the > relative error should be minimized for this. I think the choice of > subintervals is related to this. With only one zero per subinterval, > the relative error near it can be minimized without messing up the > relative error near other zeros. The adjustment for this is absolutely > tiny, so it is easier to arrange that it doesn't mess up the absolute > error too much. > I'm beginning to like my use of the product formula. It allows one to remove the zero from an interval. I did some code spelunking and found #ifndef lint static char sccsid[] = "@(#)j0.c 4.2 (Berkeley) 8/21/85"; #endif not lint Found at http://www.retro11.de/ouxr/43bsd/usr/src/usr.lib/libm/j0.c.html The algorithm is completely different than what we have in libm. I won't have time to look at this until the weekend. -- Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180906223026.GA43005>