Fix Bug 83566 - cyl_bessel_j returns wrong result for x>1000 for high orders
Michele Pezzutti
mpezz@tiscali.it
Fri Jan 5 22:46:00 GMT 2018
On 01/05/2018 05:50 PM, Jonathan Wakely wrote:
> On 05/01/18 09:54 -0500, Ed Smith-Rowland wrote:
>> This looks good to me. The only nit is that you have to uglify the k
>> loop variable.
>
> Right, it needs to be __k.
>
> I'm also not sure we can put the copyright notice there, rather than
> at the top of the file. We can collapse the years to 1996-2003 so I
> suggest adding to the top of the file something like:
>
> /* __cyl_bessel_jn_asymp adapted from GNU GSL version 2.4
> specfunc/bessel_j.c
> * Copyright (C) 1996-2003 Gerard Jungman
> */
>
diff --git a/libstdc++-v3/include/tr1/bessel_function.tcc
b/libstdc++-v3/include/tr1/bessel_function.tcc
index 7ac733d..6f5ff34 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,40 @@ 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();
>> (I hope we can un-uglify our variables when modules come!)
>>
>> Do you have copyright assignment in place and all that?
>
> That appears to be underway.
No news since I sent questionnaire though.
>
> Thanks for working on this.
>
Welcome.
More information about the Libstdc++
mailing list