On 01/04/2018 06:16 PM, Ed Smith-Rowland wrote:
On 01/03/2018 02:49 PM, Michele Pezzutti wrote:
Hi.
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:
int
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.
Here it is.
diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc
b/libstdc++-v3/include/tr1/bessel_function.tcc
index 7ac733d..7842350 100644
--- a/libstdc++-v3/include/tr1/bessel_function.tcc
+++ b/libstdc++-v3/include/tr1/bessel_function.tcc
@@ -353,21 +353,47 @@ namespace tr1
* @param __x The argument of the Bessel functions.
* @param __Jnu The output Bessel function of the first kind.
* @param __Nnu The output Neumann function (Bessel function
of the second kind).
+ *
+ * Adapted for libstdc++ from GNU GSL version 2.4 specfunc/bessel_j.c
+ * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 Gerard
Jungman
*/
template <typename _Tp>
void
__cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu)
{
const _Tp __mu = _Tp(4) * __nu * __nu;
- const _Tp __mum1 = __mu - _Tp(1);
- const _Tp __mum9 = __mu - _Tp(9);
- const _Tp __mum25 = __mu - _Tp(25);
- const _Tp __mum49 = __mu - _Tp(49);
- const _Tp __xx = _Tp(64) * __x * __x;
- const _Tp __P = _Tp(1) - __mum1 * __mum9 / (_Tp(2) * __xx)
- * (_Tp(1) - __mum25 * __mum49 / (_Tp(12) * __xx));
- const _Tp __Q = __mum1 / (_Tp(8) * __x)
- * (_Tp(1) - __mum9 * __mum25 / (_Tp(6) * __xx));
+ const _Tp __8x = _Tp(8) * __x;
+
+ _Tp __P = _Tp(0);
+ _Tp __Q = _Tp(0);
+
+ _Tp k = _Tp(0);
+ _Tp __term = _Tp(1);
+
+ int __epsP = 0;
+ int __epsQ = 0;
+
+ _Tp __eps = std::numeric_limits<_Tp>::epsilon();
+
+ do
+ {
+ __term *= (k == 0) ? _Tp(1) : -(__mu - (2 * k - 1) * (2 * k -
1)) / (k * __8x);
+ __epsP = std::abs(__term) < std::abs(__eps * __P);
+ __P += __term;
+
+ k++;
+
+ __term *= (__mu - (2 * k - 1) * (2 * k - 1)) / (k * __8x);
+ __epsQ = std::abs(__term) < std::abs(__eps * __Q);
+ __Q += __term;
+
+ if (__epsP && __epsQ && k > __nu / 2.)
+ break;
+
+ k++;
+ }
+ while (k < 1000);
+
const _Tp __chi = __x - (__nu + _Tp(0.5L))
* __numeric_constants<_Tp>::__pi_2();
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.
Ed
