This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
Re: Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders
On 01/02/2018 04:41 PM, Michele Pezzutti wrote:
On 01/02/2018 05:59 PM, Ed Smith-Rowland wrote:
OK,
on *third* look summing up to k = nu/2 at minimum will a achieve the
result of not blowing up the asymptotic series:
nu^2 - (2k-1)^2. And it will do that without a check.
This stopping criterion should work even near x=nu which would be the
most difficult case. The sum could go further for larger x but let's
just go with your termination criterion for now. Later, with some
experimentation, we could sum up to nu/2 at a minimum *then* snoop
forward until the terms start drifting up. Or we could just solve
k_max for this case as a function of x. Also, we may never need these
extras.
For what I could see, 'term' stops to contributing to P, Q few steps
before k = nu /2.
Right. I can see these terms rising, sometimes quite high, then falling.
In many experiments, the sums converge before the terms blow up.
In my old impl, I gave up too easily.
Perhaps, in general, one could wait until the terms started decreasing,
*then* if subsequent terms start growing again exit loop.
But it seems I never get to that point before sum convergence.
On the good side it seems that if x >= 20 nu then you can go quite high
in nu.
nu = 1000
x = 20000
Jnu = 0.00540683536187055
Nnu = -0.00162387432151849
Jnua = 0.00540683536187055
Nnua = -0.00162387432151849
Jnua - Jnu = 0
Nnua - Nnu = 0
Jpnu = 0.00162170521572447
Npnu = 0.0054001182451475
Jpnua = 0.00162170521572447
Npnua = 0.0054001182451475
Jpnua - Jpnu = 0
Npnua - Npnu = 0
pi x Wronski / 2 = 1.00000022382762