From owner-freebsd-numerics@freebsd.org Fri Sep 7 14:07:32 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 206CDFFBDE9 for ; Fri, 7 Sep 2018 14:07:32 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail110.syd.optusnet.com.au (mail110.syd.optusnet.com.au [211.29.132.97]) by mx1.freebsd.org (Postfix) with ESMTP id 66C428AFF1 for ; Fri, 7 Sep 2018 14:07:30 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from [192.168.0.102] (c110-21-101-228.carlnfd1.nsw.optusnet.com.au [110.21.101.228]) by mail110.syd.optusnet.com.au (Postfix) with ESMTPS id 2B74910685F; Sat, 8 Sep 2018 00:07:27 +1000 (AEST) Date: Sat, 8 Sep 2018 00:07:26 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl cc: freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) In-Reply-To: <20180906223026.GA43005@troutmask.apl.washington.edu> Message-ID: <20180907230825.A933@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> <20180906223026.GA43005@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.2 cv=I9sVfJog c=1 sm=1 tr=0 a=PalzARQSbocsUSjMRkwAPg==:117 a=PalzARQSbocsUSjMRkwAPg==:17 a=kj9zAlcOel0A:10 a=x7bEGLp0ZPQA:10 a=xz5PcUCpAAAA:8 a=NdYQpMVuwmpbWQRTCfkA:9 a=CjuIK1q_8ugA:10 a=DwsaIXcYyF8A:10 a=pM-ymjXk_DYA:10 a=IWbigzpecG3ytlOIVP4C:22 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: Fri, 07 Sep 2018 14:07:32 -0000 On Thu, 6 Sep 2018, Steve Kargl wrote: > On Fri, Sep 07, 2018 at 02:02:20AM +1000, Bruce Evans wrote: >> ... >> I now think that no single fudge factor works for both P0 and Q0. P0 >> and Q0 have very different characteristics. Merely choosing the >> ... >> 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. In fact, this must more or less work, and the difficulty is to make it work more than less. Start with any reasonable formula like (u*cc-v*ss). (This is not specific to j0, except it gives an example of a nontrivial formula.) cc and ss are held fixed. The best choice of u is not known, but we if we can get fairly close to it (perhaps even off by a factor of 2), then we can throw away any previous approximation for v and solve for v (in high precision), then approximate that. There must be no singularities for this to work in high precision. We expect to avoid the singularities by some combination of starting with u close enough, keeping the interval small, and special handling of zeros. Here the formula for v is something divided by ss, so we want to avoid zeros of ss. ss is sin(x) - cos(x) which is close to sqrt(2) of the entire interval [2,2.857], so there is no problem with pole singularities. In the next interval, ss has a zero in the middle of the interval, so we should switch the roles of u and v so as to divide by cc instead (it is not so close to -sqrt(2)). Some magic is needed to make the resulting v easy to approximate. Hopefully not much more than a good u and a small interval. Then there is the known zero in the middle of the interval. v should be chosen to at least have the correct signs on the FP values on each side of the zero. I noticed that this method cannot work for jn, since u and v depend on n so there is no way to generate before runtime. jn actually uses a different method and I suspect that it is much more inaccurate. My tests don't cover jn since its n arg is an integer and I only have coverage for 2 real args or 1 complex arg. > I'm beginning to like my use of the product formula. It > allows one to remove the zero from an interval. I did some I haven't seen that. Is it for the final evaluation or for determining the approximation? I don't know of any better way determining approximations near zeros except to multiply by a linear factor to remove the zero. It is becoming clearer why a separate interval is needed for each zero. Polynomials can only give good approximations for one zero at a time except in special cases. > 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. The one in Net/2 (1992) is the same as now, all the way down to style bugs like spelling qzero as pzero in 1 place in the comment about qzero. Bruce