From owner-freebsd-numerics@freebsd.org Wed Sep 5 15:21:12 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 7F88CFF3A35 for ; Wed, 5 Sep 2018 15:21:12 +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 1FB0C85A2D for ; Wed, 5 Sep 2018 15:21:12 +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 w85FL4gX026575 (version=TLSv1.2 cipher=DHE-RSA-AES256-GCM-SHA384 bits=256 verify=NO); Wed, 5 Sep 2018 08:21:04 -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 w85FL4lG026574; Wed, 5 Sep 2018 08:21:04 -0700 (PDT) (envelope-from sgk) Date: Wed, 5 Sep 2018 08:21:04 -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: <20180905152104.GA26453@troutmask.apl.washington.edu> Reply-To: sgk@troutmask.apl.washington.edu References: <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20180905201540.D1142@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: Wed, 05 Sep 2018 15:21:12 -0000 On Wed, Sep 05, 2018 at 10:06:29PM +1000, Bruce Evans wrote: > 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: 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. > 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. If |x| > p/2*log(2) with p the precision of x, P0(x) and Q0(x) can be used directly. It is the 2 <= |x| < p/2*log(2) range that is the issue, but your last paragraph below may be what I have been missing. > > 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. Obviously, I had not thought about scaling P0(x) and Q0(x) with a common factor to force agreement between j0(x) and its approximation. -- Steve