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/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




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