This is the mail archive of the 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/03/2018 02:49 PM, Michele Pezzutti wrote:


On 01/02/2018 05:43 PM, Michele Pezzutti wrote:

On 01/02/2018 02:28 AM, Ed Smith-Rowland wrote:
I like the patch.

I have a similar one in the tr29124 branch.

Anyway, I got held up and I think it's good to have new folks looking into this.

It looks good except that you need to uglify k.

I looked at the GSL implementation, based on same reference, and their loop is cleaner. What about porting that implementation here? Possible?

My implementation is also using one more term for P than for Q, which is discouraged in GSL, according to their comments.

Ed, do you have any comment about this point?
Regarding the 2nd line, after my observations, usually term stops contributing to P, Q before k > nu/2, so actually, an offset by one is most likely without consequences.
GSL implementation is nevertheless more elegant.

I've seen their implementation.  It's tight.  Feel free to port it.
I assume this is it:
gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result);

You do want to give Q or b one more term so as to make that last factor as small as possible.  They're right.

In principal, these ideas should be ported to I,K but I think (and IIRC GSL agrees) these under,over-flow before they have much effect.
I think I put the full series in there.  They could use a similar clean-up.

But let's get this in there first.


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