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