From owner-freebsd-standards@FreeBSD.ORG Thu Jan 14 22:49:39 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id CA601106566B; Thu, 14 Jan 2010 22:49:39 +0000 (UTC) (envelope-from kargl@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id ABC7C8FC16; Thu, 14 Jan 2010 22:49:39 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.3/8.14.3) with ESMTP id o0EMnk1t047429; Thu, 14 Jan 2010 14:49:46 -0800 (PST) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.14.3/8.14.3/Submit) id o0EMnjLS047428; Thu, 14 Jan 2010 14:49:45 -0800 (PST) (envelope-from kargl) From: "Steven G. Kargl" Message-Id: <201001142249.o0EMnjLS047428@troutmask.apl.washington.edu> In-Reply-To: <20100114200057.C62558@delplex.bde.org> To: Bruce Evans Date: Thu, 14 Jan 2010 14:49:45 -0800 (PST) X-Mailer: ELM [version 2.4ME+ PL124d (25)] MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset="US-ASCII" Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 14 Jan 2010 22:49:39 -0000 Bruce Evans wrote: > On Wed, 13 Jan 2010, Steven G. Kargl wrote: > > >> Description: > > > > The j0 bessel function supplied by libm is fairly inaccurate at > > arguments at and near a zero of the function. Here's the first > > 30 zeros computed by j0f, my implementation of j0f, a 4000-bit > > significand computed via MPFR (the assumed exact value), and the > > relative absolute error. > > This is a very hard and relatively unimportant problem. Yes, it is very hard, but apparently you do not use bessel functions in your everyday life. :) I only discover this issue because I need bessel functions of complex arguments and I found my routines have issues in the vicinity of zeros. So, I decided to look at the libm routines. > > x my j0f(x) libm j0f(x) MPFR j0 my err libm err > > 2.404825 5.6434434E-08 5.9634296E-08 5.6434400E-08 1.64 152824.59 > > 5.520078 2.4476664E-08 2.4153294E-08 2.4476659E-08 0.31 18878.52 > > 8.653728 1.0355323E-07 1.0359805E-07 1.0355306E-07 6.36 1694.47 > > 11.791534 -2.4511966E-09 -3.5193941E-09 -3.5301714E-09 78243.14 781.53 > > Hmm. I forgot to mention that 'my err' and 'libm err' are in units of epsilon (ie, FLT_EPSILON for j0f). > > Note, my j0f(x) currently uses double precision to accumulate intermediate > > values. Below x = 10, I use the ascending series to compute the value. > > Above x = 10, I'm using an asymptotic approximation. I haven't investigated > > whether additional terms in an asymptotic approximation would pull 'my err' > > for x = 11, 14, 18, and 21 closer to the exact value. > (snip) > Anyway, if you can get anywhere near < 10 ulp error near all zeros using > only an asymptotic method, then that would be good. Then the asymptotic > method would also be capable of locating the zeros very accurately. But > I would be very surprised if this worked. I know of nothing similar for > reducing mod Pi for trigonometric functions, which seems a simpler problem. > I would expect it to at best involve thousands of binary digits in the > tables for the asymptotic method, and corresponding thousands of digits > of precision in the calculation (4000 as for mfpr enough for the 2**100th > zero?). The 4000-bit setting for mpfr was a hold over from testing mpfr_j0 against my ascending series implementation of j0 with mpfr primitives. As few as 128-bits is sufficient to achieve the following: 2.404825 5.6434398E-08 5.9634296E-08 5.6434400E-08 0.05 152824.59 5.520078 2.4476657E-08 2.4153294E-08 2.4476659E-08 0.10 18878.52 8.653728 1.0355303E-07 1.0359805E-07 1.0355306E-07 0.86 1694.47 11.791534 -3.5291243E-09 -3.5193941E-09 -3.5301714E-09 75.93 781.53 14.930918 -6.4815082E-09 -6.3911618E-09 -6.4815052E-09 0.23 6722.88 18.071064 5.1532352E-09 5.3149818E-09 5.1532318E-09 0.23 10910.50 21.211637 -1.5023349E-07 -1.5002509E-07 -1.5023348E-07 2.70 56347.01 24.352472 1.2524569E-07 1.2516310E-07 1.2524570E-07 0.28 2834.53 27.493479 5.4331110E-08 5.4263626E-08 5.4331104E-08 0.29 3261.75 30.634607 1.2205545E-07 1.2203689E-07 1.2205546E-07 0.09 645.39 33.775822 -2.0213095E-07 -2.0206903E-07 -2.0213095E-07 0.27 6263.95 36.917099 8.4751576E-08 8.4749573E-08 8.4751581E-08 0.18 82.59 40.058426 -1.7484838E-08 -1.7475532E-08 -1.7484840E-08 0.12 767.56 43.199791 -9.2091398E-08 -9.2135146E-08 -9.2091406E-08 2.47 13530.51 46.341187 2.1663259E-07 2.1664336E-07 2.1663259E-07 0.16 268.90 49.482609 -1.2502527E-07 -1.2504157E-07 -1.2502526E-07 2.69 23512.60 52.624050 1.8706569E-07 1.8707487E-07 1.8706569E-07 0.01 251.43 55.765511 -2.0935557E-08 -2.0932896E-08 -2.0935556E-08 0.10 227.04 58.906982 1.5637660E-07 1.5634730E-07 1.5637661E-07 0.28 892.23 62.048470 3.5779891E-08 3.5787338E-08 3.5779899E-08 0.42 402.61 As I suspected by adding additional terms to the asymptotic approximation and performing all computations with double precision, reduces 'my err' (5th column). The value at x=11.7... is the best I can get. The asymptotic approximations contain divergent series and additional terms do not help. -- Steve http://troutmask.apl.washington.edu/~kargl/