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


Formatting fixed.

diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc b/libstdc++-v3/include/tr1/bessel_function.tcc
index 7ac733d..5f8fc9f 100644
--- a/libstdc++-v3/include/tr1/bessel_function.tcc
+++ b/libstdc++-v3/include/tr1/bessel_function.tcc
@@ -27,6 +27,10 @@
  *  Do not attempt to use it directly. @headername{tr1/cmath}
  */

+/* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4 specfunc/bessel_j.c
+ * Copyright (C) 1996-2003 Gerard Jungman
+ */
+
 //
 // ISO C++ 14882 TR1: 5.2  Special functions
 //
@@ -358,16 +362,42 @@ namespace tr1
     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 __mu = _Tp(4) * __nu * __nu;
+      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();

On 01/06/2018 10:23 AM, Paolo Carlini wrote:
Hi,

On 05/01/2018 23:46, Michele Pezzutti wrote:
+ __term *= (__k == 0) ? _Tp(1) : -(__mu - (2 * __k - 1) * (2 * __k - 1))
+                        / (__k * __8x);
In such cases, per the Coding Standards, you want an outer level of parentheses wrapping the whole right side expression. See toward the end of https://www.gnu.org/prep/standards/html_node/Formatting.html#Formatting. I see that many other places will need fixing, at some point - this rather straightforward rule is often overlooked leading to brittle formatting of many expressions, probably because it's really obvious in practice only together with Emacs, I'm not sure.

Also - this kind of stylistic nitpicking is partially personal taste - the parentheses around (__k == 0) seem definitely redundant to me, and I don't think you would find many examples of that in our code.

About the Assignement, please be patient. For example, used to be the case that when RMS was traveling couldn't process Assignments, for example. It can be slow for many reasons, it's 100% normal.

Paolo..




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