From owner-freebsd-numerics@freebsd.org Thu Sep 6 22:30:35 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 3F94CFE7A14 for ; Thu, 6 Sep 2018 22:30:35 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (Client CN "troutmask", Issuer "troutmask" (not verified)) by mx1.freebsd.org (Postfix) with ESMTPS id D58CC8F52E for ; Thu, 6 Sep 2018 22:30:34 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.15.2/8.15.2) with ESMTPS id w86MUQCV043015 (version=TLSv1.2 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Thu, 6 Sep 2018 15:30:26 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.15.2/8.15.2/Submit) id w86MUQu8043014; Thu, 6 Sep 2018 15:30:26 -0700 (PDT) (envelope-from sgk) Date: Thu, 6 Sep 2018 15:30:26 -0700 From: Steve Kargl To: Bruce Evans 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> Reply-To: sgk@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> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20180906202658.K978@besplex.bde.org> User-Agent: Mutt/1.10.1 (2018-07-13) 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: Thu, 06 Sep 2018 22:30:35 -0000 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