re PR libstdc++/56216 (TR1 bessel functions bomb at x == 0!)

PR libstdc++/56216

From-SVN: r195886

libstdc++-v3/ChangeLog

+2013-02-08  Edward Smith-Rowland  <3dw4rd@verizon.net>
+
+	PR libstdc++/56216
+	* include/tr1/special_function_util.h: Remove spurious const
+	from numeric arguments.
+	* include/tr1/riemann_zeta.tcc: Ditto.
+	* include/tr1/exp_integral.tcc: Ditto.
+	* include/tr1/bessel_function.tcc: Ditto.
+	* include/tr1/hypergeometric.tcc: Ditto.
+	* include/tr1/modified_bessel_func.tcc: Ditto.
+	* include/tr1/poly_laguerre.tcc: Ditto.
+	* include/tr1/gamma.tcc: Ditto.
+	* include/tr1/legendre_function.tcc: Ditto.
+	* include/tr1/poly_hermite.tcc: Ditto.
+	* include/tr1/ell_integral.tcc: Ditto.
+	* include/tr1/bessel_function.tcc (__cyl_bessel_ij_series):
+	If argument is zero return function value.
+	* testsuite/tr1/5_numerical_facilities/special_functions/08_cyl_bessel_i/pr56216.cc:
+	New file.
+
 2013-02-07  Paolo Carlini  <paolo.carlini@oracle.com>
 
 	* testsuite/27_io/basic_ios/pr56193.cc: Tweak.

libstdc++-v3/include/tr1/bessel_function.tcc

@@ -86,8 +86,8 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __gamma_temme(const _Tp __mu,
-                   _Tp & __gam1, _Tp & __gam2, _Tp & __gampl, _Tp & __gammi)
+    __gamma_temme(_Tp __mu,
+                  _Tp & __gam1, _Tp & __gam2, _Tp & __gampl, _Tp & __gammi)
     {
 #if _GLIBCXX_USE_C99_MATH_TR1
       __gampl = _Tp(1) / std::tr1::tgamma(_Tp(1) + __mu);
@@ -124,7 +124,7 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __bessel_jn(const _Tp __nu, const _Tp __x,
+    __bessel_jn(_Tp __nu, _Tp __x,
                 _Tp & __Jnu, _Tp & __Nnu, _Tp & __Jpnu, _Tp & __Npnu)
     {
       if (__x == _Tp(0))
@@ -349,11 +349,8 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __cyl_bessel_jn_asymp(const _Tp __nu, const _Tp __x,
-                          _Tp & __Jnu, _Tp & __Nnu)
+    __cyl_bessel_jn_asymp(_Tp __nu, _Tp __x, _Tp & __Jnu, _Tp & __Nnu)
     {
-      const _Tp __coef = std::sqrt(_Tp(2)
-                             / (__numeric_constants<_Tp>::__pi() * __x));
       const _Tp __mu   = _Tp(4) * __nu * __nu;
       const _Tp __mum1 = __mu - _Tp(1);
       const _Tp __mum9 = __mu - _Tp(9);
@@ -370,6 +367,8 @@ namespace tr1
       const _Tp __c = std::cos(__chi);
       const _Tp __s = std::sin(__chi);
 
+      const _Tp __coef = std::sqrt(_Tp(2)
+                             / (__numeric_constants<_Tp>::__pi() * __x));
       __Jnu = __coef * (__c * __P - __s * __Q);
       __Nnu = __coef * (__s * __P + __c * __Q);
 
@@ -406,10 +405,18 @@ namespace tr1
      */
     template <typename _Tp>
     _Tp
-    __cyl_bessel_ij_series(const _Tp __nu, const _Tp __x, const _Tp __sgn,
-                           const unsigned int __max_iter)
+    __cyl_bessel_ij_series(_Tp __nu, _Tp __x, _Tp __sgn,
+                           unsigned int __max_iter)
     {
-
+      if (__x == _Tp(0))
+	{
+          if (__nu == _Tp(0))
+            return _Tp(1);
+          else if (__nu == _Tp(1))
+            return _Tp(0);
+          else
+            return _Tp(0);
+	}
       const _Tp __x2 = __x / _Tp(2);
       _Tp __fact = __nu * std::log(__x2);
 #if _GLIBCXX_USE_C99_MATH_TR1
@@ -450,7 +457,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __cyl_bessel_j(const _Tp __nu, const _Tp __x)
+    __cyl_bessel_j(_Tp __nu, _Tp __x)
     {
       if (__nu < _Tp(0) || __x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "
@@ -492,7 +499,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __cyl_neumann_n(const _Tp __nu, const _Tp __x)
+    __cyl_neumann_n(_Tp __nu, _Tp __x)
     {
       if (__nu < _Tp(0) || __x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "
@@ -529,7 +536,7 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __sph_bessel_jn(const unsigned int __n, const _Tp __x,
+    __sph_bessel_jn(unsigned int __n, _Tp __x,
                     _Tp & __j_n, _Tp & __n_n, _Tp & __jp_n, _Tp & __np_n)
     {
       const _Tp __nu = _Tp(__n) + _Tp(0.5L);
@@ -564,7 +571,7 @@ namespace tr1
      */
     template <typename _Tp>
     _Tp
-    __sph_bessel(const unsigned int __n, const _Tp __x)
+    __sph_bessel(unsigned int __n, _Tp __x)
     {
       if (__x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "
@@ -602,7 +609,7 @@ namespace tr1
      */
     template <typename _Tp>
     _Tp
-    __sph_neumann(const unsigned int __n, const _Tp __x)
+    __sph_neumann(unsigned int __n, _Tp __x)
     {
       if (__x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "

libstdc++-v3/include/tr1/ell_integral.tcc

@@ -70,7 +70,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z)
+    __ellint_rf(_Tp __x, _Tp __y, _Tp __z)
     {
       const _Tp __min = std::numeric_limits<_Tp>::min();
       const _Tp __max = std::numeric_limits<_Tp>::max();
@@ -149,7 +149,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __comp_ellint_1_series(const _Tp __k)
+    __comp_ellint_1_series(_Tp __k)
     {
 
       const _Tp __kk = __k * __k;
@@ -187,7 +187,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __comp_ellint_1(const _Tp __k)
+    __comp_ellint_1(_Tp __k)
     {
 
       if (__isnan(__k))
@@ -215,7 +215,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_1(const _Tp __k, const _Tp __phi)
+    __ellint_1(_Tp __k, _Tp __phi)
     {
 
       if (__isnan(__k) || __isnan(__phi))
@@ -262,7 +262,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __comp_ellint_2_series(const _Tp __k)
+    __comp_ellint_2_series(_Tp __k)
     {
 
       const _Tp __kk = __k * __k;
@@ -310,7 +310,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z)
+    __ellint_rd(_Tp __x, _Tp __y, _Tp __z)
     {
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
       const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6));
@@ -398,7 +398,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __comp_ellint_2(const _Tp __k)
+    __comp_ellint_2(_Tp __k)
     {
 
       if (__isnan(__k))
@@ -432,7 +432,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_2(const _Tp __k, const _Tp __phi)
+    __ellint_2(_Tp __k, _Tp __phi)
     {
 
       if (__isnan(__k) || __isnan(__phi))
@@ -491,7 +491,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_rc(const _Tp __x, const _Tp __y)
+    __ellint_rc(_Tp __x, _Tp __y)
     {
       const _Tp __min = std::numeric_limits<_Tp>::min();
       const _Tp __max = std::numeric_limits<_Tp>::max();
@@ -562,7 +562,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p)
+    __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p)
     {
       const _Tp __min = std::numeric_limits<_Tp>::min();
       const _Tp __max = std::numeric_limits<_Tp>::max();
@@ -666,7 +666,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __comp_ellint_3(const _Tp __k, const _Tp __nu)
+    __comp_ellint_3(_Tp __k, _Tp __nu)
     {
 
       if (__isnan(__k) || __isnan(__nu))
@@ -706,7 +706,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi)
+    __ellint_3(_Tp __k, _Tp __nu, _Tp __phi)
     {
 
       if (__isnan(__k) || __isnan(__nu) || __isnan(__phi))

libstdc++-v3/include/tr1/exp_integral.tcc

@@ -58,7 +58,7 @@ namespace tr1
   {
   _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
-    template<typename _Tp> _Tp __expint_E1(const _Tp);
+    template<typename _Tp> _Tp __expint_E1(_Tp);
 
     /**
      *   @brief Return the exponential integral @f$ E_1(x) @f$
@@ -75,7 +75,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_E1_series(const _Tp __x)
+    __expint_E1_series(_Tp __x)
     {
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
       _Tp __term = _Tp(1);
@@ -112,7 +112,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_E1_asymp(const _Tp __x)
+    __expint_E1_asymp(_Tp __x)
     {
       _Tp __term = _Tp(1);
       _Tp __esum = _Tp(1);
@@ -149,7 +149,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_En_series(const unsigned int __n, const _Tp __x)
+    __expint_En_series(unsigned int __n, _Tp __x)
     {
       const unsigned int __max_iter = 100;
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
@@ -195,7 +195,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_En_cont_frac(const unsigned int __n, const _Tp __x)
+    __expint_En_cont_frac(unsigned int __n, _Tp __x)
     {
       const unsigned int __max_iter = 100;
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
@@ -240,7 +240,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_En_recursion(const unsigned int __n, const _Tp __x)
+    __expint_En_recursion(unsigned int __n, _Tp __x)
     {
       _Tp __En;
       _Tp __E1 = __expint_E1(__x);
@@ -284,7 +284,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_Ei_series(const _Tp __x)
+    __expint_Ei_series(_Tp __x)
     {
       _Tp __term = _Tp(1);
       _Tp __sum = _Tp(0);
@@ -315,7 +315,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_Ei_asymp(const _Tp __x)
+    __expint_Ei_asymp(_Tp __x)
     {
       _Tp __term = _Tp(1);
       _Tp __sum = _Tp(1);
@@ -348,7 +348,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_Ei(const _Tp __x)
+    __expint_Ei(_Tp __x)
     {
       if (__x < _Tp(0))
         return -__expint_E1(-__x);
@@ -372,7 +372,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_E1(const _Tp __x)
+    __expint_E1(_Tp __x)
     {
       if (__x < _Tp(0))
         return -__expint_Ei(-__x);
@@ -402,7 +402,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_asymp(const unsigned int __n, const _Tp __x)
+    __expint_asymp(unsigned int __n, _Tp __x)
     {
       _Tp __term = _Tp(1);
       _Tp __sum = _Tp(1);
@@ -436,7 +436,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint_large_n(const unsigned int __n, const _Tp __x)
+    __expint_large_n(unsigned int __n, _Tp __x)
     {
       const _Tp __xpn = __x + __n;
       const _Tp __xpn2 = __xpn * __xpn;
@@ -470,7 +470,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __expint(const unsigned int __n, const _Tp __x)
+    __expint(unsigned int __n, _Tp __x)
     {
       //  Return NaN on NaN input.
       if (__isnan(__x))
@@ -510,7 +510,7 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __expint(const _Tp __x)
+    __expint(_Tp __x)
     {
       if (__isnan(__x))
         return std::numeric_limits<_Tp>::quiet_NaN();

libstdc++-v3/include/tr1/gamma.tcc

@@ -67,7 +67,8 @@ namespace tr1
      *   @return  The Bernoulli number of order n.
      */
     template <typename _Tp>
-    _Tp __bernoulli_series(unsigned int __n)
+    _Tp
+    __bernoulli_series(unsigned int __n)
     {
 
       static const _Tp __num[28] = {
@@ -130,10 +131,8 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __bernoulli(const int __n)
-    {
-      return __bernoulli_series<_Tp>(__n);
-    }
+    __bernoulli(int __n)
+    { return __bernoulli_series<_Tp>(__n); }
 
 
     /**
@@ -146,7 +145,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __log_gamma_bernoulli(const _Tp __x)
+    __log_gamma_bernoulli(_Tp __x)
     {
       _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
                + _Tp(0.5L) * std::log(_Tp(2)
@@ -174,7 +173,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __log_gamma_lanczos(const _Tp __x)
+    __log_gamma_lanczos(_Tp __x)
     {
       const _Tp __xm1 = __x - _Tp(1);
 
@@ -218,7 +217,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __log_gamma(const _Tp __x)
+    __log_gamma(_Tp __x)
     {
       if (__x > _Tp(0.5L))
         return __log_gamma_lanczos(__x);
@@ -245,7 +244,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __log_gamma_sign(const _Tp __x)
+    __log_gamma_sign(_Tp __x)
     {
       if (__x > _Tp(0))
         return _Tp(1);
@@ -276,7 +275,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __log_bincoef(const unsigned int __n, const unsigned int __k)
+    __log_bincoef(unsigned int __n, unsigned int __k)
     {
       //  Max e exponent before overflow.
       static const _Tp __max_bincoeff
@@ -307,7 +306,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __bincoef(const unsigned int __n, const unsigned int __k)
+    __bincoef(unsigned int __n, unsigned int __k)
     {
       //  Max e exponent before overflow.
       static const _Tp __max_bincoeff
@@ -330,10 +329,8 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __gamma(const _Tp __x)
-    {
-      return std::exp(__log_gamma(__x));
-    }
+    __gamma(_Tp __x)
+    { return std::exp(__log_gamma(__x)); }
 
 
     /**
@@ -351,7 +348,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __psi_series(const _Tp __x)
+    __psi_series(_Tp __x)
     {
       _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
       const unsigned int __max_iter = 100000;
@@ -381,7 +378,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __psi_asymp(const _Tp __x)
+    __psi_asymp(_Tp __x)
     {
       _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
       const _Tp __xx = __x * __x;
@@ -412,7 +409,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __psi(const _Tp __x)
+    __psi(_Tp __x)
     {
       const int __n = static_cast<int>(__x + 0.5L);
       const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
@@ -441,7 +438,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __psi(const unsigned int __n, const _Tp __x)
+    __psi(unsigned int __n, _Tp __x)
     {
       if (__x <= _Tp(0))
         std::__throw_domain_error(__N("Argument out of range "

libstdc++-v3/include/tr1/hypergeometric.tcc

@@ -75,7 +75,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __conf_hyperg_series(const _Tp __a, const _Tp __c, const _Tp __x)
+    __conf_hyperg_series(_Tp __a, _Tp __c, _Tp __x)
     {
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
 
@@ -112,7 +112,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __conf_hyperg_luke(const _Tp __a, const _Tp __c, const _Tp __xin)
+    __conf_hyperg_luke(_Tp __a, _Tp __c, _Tp __xin)
     {
       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));
       const int __nmax = 20000;
@@ -218,8 +218,8 @@ namespace tr1
      *   @return  The confluent hypergeometric function.
      */
     template<typename _Tp>
-    inline _Tp
-    __conf_hyperg(const _Tp __a, const _Tp __c, const _Tp __x)
+    _Tp
+    __conf_hyperg(_Tp __a, _Tp __c, _Tp __x)
     {
 #if _GLIBCXX_USE_C99_MATH_TR1
       const _Tp __c_nint = std::tr1::nearbyint(__c);
@@ -263,8 +263,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __hyperg_series(const _Tp __a, const _Tp __b,
-                    const _Tp __c, const _Tp __x)
+    __hyperg_series(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
     {
       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
 
@@ -297,8 +296,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __hyperg_luke(const _Tp __a, const _Tp __b, const _Tp __c,
-                  const _Tp __xin)
+    __hyperg_luke(_Tp __a, _Tp __b, _Tp __c, _Tp __xin)
     {
       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));
       const int __nmax = 20000;
@@ -432,8 +430,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __hyperg_reflect(const _Tp __a, const _Tp __b, const _Tp __c,
-                     const _Tp __x)
+    __hyperg_reflect(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
     {
       const _Tp __d = __c - __a - __b;
       const int __intd  = std::floor(__d + _Tp(0.5L));
@@ -722,8 +719,8 @@ namespace tr1
      *   @return  The confluent hypergeometric function.
      */
     template<typename _Tp>
-    inline _Tp
-    __hyperg(const _Tp __a, const _Tp __b, const _Tp __c, const _Tp __x)
+    _Tp
+    __hyperg(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
     {
 #if _GLIBCXX_USE_C99_MATH_TR1
       const _Tp __a_nint = std::tr1::nearbyint(__a);

libstdc++-v3/include/tr1/legendre_function.tcc

@@ -72,7 +72,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __poly_legendre_p(const unsigned int __l, const _Tp __x)
+    __poly_legendre_p(unsigned int __l, _Tp __x)
     {
 
       if ((__x < _Tp(-1)) || (__x > _Tp(+1)))
@@ -129,8 +129,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __assoc_legendre_p(const unsigned int __l, const unsigned int __m,
-                       const _Tp __x)
+    __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x)
     {
 
       if (__x < _Tp(-1) || __x > _Tp(+1))
@@ -209,8 +208,7 @@ namespace tr1
      */
     template <typename _Tp>
     _Tp
-    __sph_legendre(const unsigned int __l, const unsigned int __m,
-                   const _Tp __theta)
+    __sph_legendre(unsigned int __l, unsigned int __m, _Tp __theta)
     {
       if (__isnan(__theta))
         return std::numeric_limits<_Tp>::quiet_NaN();

libstdc++-v3/include/tr1/modified_bessel_func.tcc

@@ -77,7 +77,7 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __bessel_ik(const _Tp __nu, const _Tp __x,
+    __bessel_ik(_Tp __nu, _Tp __x,
                 _Tp & __Inu, _Tp & __Knu, _Tp & __Ipnu, _Tp & __Kpnu)
     {
       if (__x == _Tp(0))
@@ -132,7 +132,7 @@ namespace tr1
         }
       if (__i > __max_iter)
         std::__throw_runtime_error(__N("Argument x too large "
-                                       "in __bessel_jn; "
+                                       "in __bessel_ik; "
                                        "try asymptotic expansion."));
       _Tp __Inul = __fp_min;
       _Tp __Ipnul = __h * __Inul;
@@ -185,7 +185,7 @@ namespace tr1
             }
           if (__i > __max_iter)
             std::__throw_runtime_error(__N("Bessel k series failed to converge "
-                                           "in __bessel_jn."));
+                                           "in __bessel_ik."));
           __Kmu = __sum;
           __Knu1 = __sum1 * __xi2;
         }
@@ -221,7 +221,7 @@ namespace tr1
             }
           if (__i > __max_iter)
             std::__throw_runtime_error(__N("Steed's method failed "
-                                           "in __bessel_jn."));
+                                           "in __bessel_ik."));
           __h = __a1 * __h;
           __Kmu = std::sqrt(__numeric_constants<_Tp>::__pi() / (_Tp(2) * __x))
                 * std::exp(-__x) / __s;
@@ -261,7 +261,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __cyl_bessel_i(const _Tp __nu, const _Tp __x)
+    __cyl_bessel_i(_Tp __nu, _Tp __x)
     {
       if (__nu < _Tp(0) || __x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "
@@ -297,7 +297,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __cyl_bessel_k(const _Tp __nu, const _Tp __x)
+    __cyl_bessel_k(_Tp __nu, _Tp __x)
     {
       if (__nu < _Tp(0) || __x < _Tp(0))
         std::__throw_domain_error(__N("Bad argument "
@@ -331,7 +331,7 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __sph_bessel_ik(const unsigned int __n, const _Tp __x,
+    __sph_bessel_ik(unsigned int __n, _Tp __x,
                     _Tp & __i_n, _Tp & __k_n, _Tp & __ip_n, _Tp & __kp_n)
     {
       const _Tp __nu = _Tp(__n) + _Tp(0.5L);
@@ -366,8 +366,7 @@ namespace tr1
      */
     template <typename _Tp>
     void
-    __airy(const _Tp __x,
-           _Tp & __Ai, _Tp & __Bi, _Tp & __Aip, _Tp & __Bip)
+    __airy(_Tp __x, _Tp & __Ai, _Tp & __Bi, _Tp & __Aip, _Tp & __Bip)
     {
       const _Tp __absx = std::abs(__x);
       const _Tp __rootx = std::sqrt(__absx);

libstdc++-v3/include/tr1/poly_hermite.tcc

@@ -66,7 +66,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __poly_hermite_recursion(const unsigned int __n, const _Tp __x)
+    __poly_hermite_recursion(unsigned int __n, _Tp __x)
     {
       //  Compute H_0.
       _Tp __H_0 = 1;
@@ -108,7 +108,7 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __poly_hermite(const unsigned int __n, const _Tp __x)
+    __poly_hermite(unsigned int __n, _Tp __x)
     {
       if (__isnan(__x))
         return std::numeric_limits<_Tp>::quiet_NaN();

libstdc++-v3/include/tr1/poly_laguerre.tcc

@@ -67,8 +67,7 @@ namespace tr1
      */
     template<typename _Tpa, typename _Tp>
     _Tp
-    __poly_laguerre_large_n(const unsigned __n, const _Tpa __alpha1,
-                            const _Tp __x)
+    __poly_laguerre_large_n(unsigned __n, _Tpa __alpha1, _Tp __x)
     {
       const _Tp __a = -_Tp(__n);
       const _Tp __b = _Tp(__alpha1) + _Tp(1);
@@ -122,8 +121,7 @@ namespace tr1
      */
     template<typename _Tpa, typename _Tp>
     _Tp
-    __poly_laguerre_hyperg(const unsigned int __n, const _Tpa __alpha1,
-			   const _Tp __x)
+    __poly_laguerre_hyperg(unsigned int __n, _Tpa __alpha1, _Tp __x)
     {
       const _Tp __b = _Tp(__alpha1) + _Tp(1);
       const _Tp __mx = -__x;
@@ -179,8 +177,7 @@ namespace tr1
      */
     template<typename _Tpa, typename _Tp>
     _Tp
-    __poly_laguerre_recursion(const unsigned int __n,
-                              const _Tpa __alpha1, const _Tp __x)
+    __poly_laguerre_recursion(unsigned int __n, _Tpa __alpha1, _Tp __x)
     {
       //   Compute l_0.
       _Tp __l_0 = _Tp(1);
@@ -238,9 +235,8 @@ namespace tr1
      *           degree @f$ \alpha @f$, and argument x.
      */
     template<typename _Tpa, typename _Tp>
-    inline _Tp
-    __poly_laguerre(const unsigned int __n, const _Tpa __alpha1,
-                    const _Tp __x)
+    _Tp
+    __poly_laguerre(unsigned int __n, _Tpa __alpha1, _Tp __x)
     {
       if (__x < _Tp(0))
         std::__throw_domain_error(__N("Negative argument "
@@ -292,11 +288,8 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __assoc_laguerre(const unsigned int __n, const unsigned int __m,
-                     const _Tp __x)
-    {
-      return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x);
-    }
+    __assoc_laguerre(unsigned int __n, unsigned int __m, _Tp __x)
+    { return __poly_laguerre<unsigned int, _Tp>(__n, __m, __x); }
 
 
     /**
@@ -315,10 +308,8 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __laguerre(const unsigned int __n, const _Tp __x)
-    {
-      return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x);
-    }
+    __laguerre(unsigned int __n, _Tp __x)
+    { return __poly_laguerre<unsigned int, _Tp>(__n, 0, __x); }
 
   _GLIBCXX_END_NAMESPACE_VERSION
   } // namespace std::tr1::__detail

libstdc++-v3/include/tr1/riemann_zeta.tcc

@@ -70,7 +70,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __riemann_zeta_sum(const _Tp __s)
+    __riemann_zeta_sum(_Tp __s)
     {
       //  A user shouldn't get to this.
       if (__s < _Tp(1))
@@ -107,7 +107,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __riemann_zeta_alt(const _Tp __s)
+    __riemann_zeta_alt(_Tp __s)
     {
       _Tp __sgn = _Tp(1);
       _Tp __zeta = _Tp(0);
@@ -149,7 +149,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __riemann_zeta_glob(const _Tp __s)
+    __riemann_zeta_glob(_Tp __s)
     {
       _Tp __zeta = _Tp(0);
 
@@ -244,7 +244,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __riemann_zeta_product(const _Tp __s)
+    __riemann_zeta_product(_Tp __s)
     {
       static const _Tp __prime[] = {
         _Tp(2), _Tp(3), _Tp(5), _Tp(7), _Tp(11), _Tp(13), _Tp(17), _Tp(19),
@@ -285,7 +285,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __riemann_zeta(const _Tp __s)
+    __riemann_zeta(_Tp __s)
     {
       if (__isnan(__s))
         return std::numeric_limits<_Tp>::quiet_NaN();
@@ -357,7 +357,7 @@ namespace tr1
      */
     template<typename _Tp>
     _Tp
-    __hurwitz_zeta_glob(const _Tp __a, const _Tp __s)
+    __hurwitz_zeta_glob(_Tp __a, _Tp __s)
     {
       _Tp __zeta = _Tp(0);
 
@@ -422,10 +422,8 @@ namespace tr1
      */
     template<typename _Tp>
     inline _Tp
-    __hurwitz_zeta(const _Tp __a, const _Tp __s)
-    {
-      return __hurwitz_zeta_glob(__a, __s);
-    }
+    __hurwitz_zeta(_Tp __a, _Tp __s)
+    { return __hurwitz_zeta_glob(__a, __s); }
 
   _GLIBCXX_END_NAMESPACE_VERSION
   } // namespace std::tr1::__detail

libstdc++-v3/include/tr1/special_function_util.h

@@ -107,30 +107,22 @@ namespace tr1
     /// out of intrinsics, this will disappear completely in favor of
     /// std::isnan.
     template<typename _Tp>
-    inline bool __isnan(const _Tp __x)
-    {
-      return std::isnan(__x);
-    }
+    inline bool __isnan(_Tp __x)
+    { return std::isnan(__x); }
 
 #else
 
     template<typename _Tp>
     inline bool __isnan(const _Tp __x)
-    {
-      return __builtin_isnan(__x);
-    }
+    { return __builtin_isnan(__x); }
 
     template<>
-    inline bool __isnan<float>(const float __x)
-    {
-      return __builtin_isnanf(__x);
-    }
+    inline bool __isnan<float>(float __x)
+    { return __builtin_isnanf(__x); }
 
     template<>
-    inline bool __isnan<long double>(const long double __x)
-    {
-      return __builtin_isnanl(__x);
-    }
+    inline bool __isnan<long double>(long double __x)
+    { return __builtin_isnanl(__x); }
 
 #endif
 

libstdc++-v3/testsuite/tr1/5_numerical_facilities/special_functions/08_cyl_bessel_i/pr56216.cc

+// 2013-02-08  Edward Smith-Rowland <3dw4rd@verizon.net>
+//
+// Copyright (C) 2013 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library.  This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3.  If not see
+// <http://www.gnu.org/licenses/>.
+
+// PR libstdc++/56216 - Crash of Bessel functions at x==0!
+
+#include <testsuite_hooks.h>
+#include <tr1/cmath>
+
+void
+test01()
+{
+  bool test __attribute__((unused)) = true;
+
+  double j0 = std::tr1::cyl_bessel_j(0.0, 0.0);
+  double i0 = std::tr1::cyl_bessel_i(0.0, 0.0);
+  double j1 = std::tr1::cyl_bessel_j(1.0, 0.0);
+  double i1 = std::tr1::cyl_bessel_i(1.0, 0.0);
+
+  VERIFY(j0 == 1.0);
+  VERIFY(i0 == 1.0);
+  VERIFY(j1 == 0.0);
+  VERIFY(i1 == 0.0);
+}
+
+int
+main()
+{
+  test01();
+
+  return 0;
+}