Date: Thu, 14 Jan 2010 17:53:17 -0800 (PST) From: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Subject: Re: standards/142803: j0 Bessel function inaccurate near zeros of the function Message-ID: <201001150153.o0F1rHu7048195@troutmask.apl.washington.edu> In-Reply-To: <20100115110443.Y63344@delplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
Bruce Evans wrote: > On Thu, 14 Jan 2010, Steven G. Kargl wrote: > > Bruce Evans wrote: > > >> 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: > > > >>> x my j0f(x) libm j0f(x) MPFR j0 my err libm err > > 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 > > Wonder why this jumps. Below x=10, I use the ascending series. Above x=0, I use an asymptotic approximation (AA). x = 11.79... is the first zero I hit with the AA. > > 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 > > ... > > > > 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. > > The extra precision is almost certainly necessary. Whether double > precision is nearly enough is unclear, but the error near 11.7 suggests > that it is nearly enough except there. The large error might be caused > by that zero alone (among small zeros) being very close to a representable > value. The AA is j0(x) = (P(x) * cos(x+pi/4) + Q(x) * sin(x+pi/4)) where P(x) and Q(x) are infinite, divergent sums. I use the first 11 terms in P(x) and Q(x) to achieve the 'my err' = 75.9. Additional terms cause an increase in 'my err' as does fewer terms. This is probably the limit of double precision. I haven't investigated the intervals around the values I listed. So, there may be larger errors that are yet to be found. BTW, the MPFR website has a document that describes all their algorithms. They claim that the AA can be used for |x| > p*log(2)/2 where p is the precision of x; however, in the mpfr code the criterion is |x| > p/2. http://www.mpfr.org/algo.html > I forgot the mention that the error table in my previous mail is on amd64 > for comparing float precision functions with double precision ones, assuming > that the latter are correct, which they aren't, but they are hopefully > correct enough for this comparision. The errors on i386 are much larger, > due to i386 still using i387 hardware trigonometric which are extremely > inaccurate near zeros, starting at the first zero. Here are both tables: My values are also computed on amd64. I seldomly use i386 for numerical work. A quick change to my test code to use the double precision j0() suggests that it has sufficient precision for the comparison you've done. -- Steve http://troutmask.apl.washington.edu/~kargl/
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201001150153.o0F1rHu7048195>