This is the mail archive of the libstdc++@gcc.gnu.org mailing list for the libstdc++ project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]