2009-06-08 Paolo Carlini <paolo.carlini@oracle.com>
authorpaolo <paolo@138bc75d-0d04-0410-961f-82ee72b054a4>
Mon, 8 Jun 2009 14:38:48 +0000 (14:38 +0000)
committerpaolo <paolo@138bc75d-0d04-0410-961f-82ee72b054a4>
Mon, 8 Jun 2009 14:38:48 +0000 (14:38 +0000)
* include/bits/random.tcc (gamma_distribution<>::operator()
(_UniformRandomNumberGenerator&, const param_type&): Redo, using
the Marsaglia/Tsang algorithm.
(gamma_distribution<>::param_type::_M_initialize): Adjust.
(operator<<(basic_ostream<>&, gamma_distribution<>),
operator>>(basic_ostream<>&, gamma_distribution<>): Likewise.

* include/bits/random.tcc(student_t_distribution<>::_M_gaussian):
Remove, just use normal_distribution.
(operator<<(basic_ostream<>&, student_t_distribution<>),
operator>>(basic_ostream<>&, student_t_distribution<>): Adjust.
(linear_congruential_engine<>::operator()()): Move inline.
(lognormal_distribution<>::operator()(_UniformRandomNumberGenerator&,
const param_type&)): Move inline, just use normal_distribution.
(operator<<(basic_ostream<>&, lognormal_distribution<>),
operator>>(basic_ostream<>&, lognormal_distribution<>): Adjust.
(weibull_distribution<>::operator()(_UniformRandomNumberGenerator&,
const param_type&)): Move here, out of line.
(piecewise_constant_distribution<>::param_type::param_type()): Move
inline.
* include/bits/random.h: Adjust, minor tweaks.

git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@148276 138bc75d-0d04-0410-961f-82ee72b054a4

libstdc++-v3/ChangeLog
libstdc++-v3/include/bits/random.h
libstdc++-v3/include/bits/random.tcc

index 9708d40..9d3bf47 100644 (file)
@@ -1,3 +1,27 @@
+2009-06-08  Paolo Carlini  <paolo.carlini@oracle.com>
+
+       * include/bits/random.tcc (gamma_distribution<>::operator()
+       (_UniformRandomNumberGenerator&, const param_type&): Redo, using
+       the Marsaglia/Tsang algorithm.
+       (gamma_distribution<>::param_type::_M_initialize): Adjust.
+       (operator<<(basic_ostream<>&, gamma_distribution<>),
+       operator>>(basic_ostream<>&, gamma_distribution<>): Likewise.
+
+       * include/bits/random.tcc(student_t_distribution<>::_M_gaussian):
+       Remove, just use normal_distribution.
+       (operator<<(basic_ostream<>&, student_t_distribution<>),
+       operator>>(basic_ostream<>&, student_t_distribution<>): Adjust.
+       (linear_congruential_engine<>::operator()()): Move inline.
+       (lognormal_distribution<>::operator()(_UniformRandomNumberGenerator&,
+       const param_type&)): Move inline, just use normal_distribution.
+       (operator<<(basic_ostream<>&, lognormal_distribution<>),
+       operator>>(basic_ostream<>&, lognormal_distribution<>): Adjust.
+       (weibull_distribution<>::operator()(_UniformRandomNumberGenerator&,
+       const param_type&)): Move here, out of line.
+       (piecewise_constant_distribution<>::param_type::param_type()): Move
+       inline.
+       * include/bits/random.h: Adjust, minor tweaks.
+
 2009-06-05  Benjamin Kosnik  <bkoz@redhat.com>
 
        * testsuite/29_atomics/atomic_address/cons/aggregate.cc: Remove xfail.
index b5ccf8a..8a21ae5 100644 (file)
@@ -268,7 +268,11 @@ namespace std
        * @brief Gets the next random number in the sequence.
        */
       result_type
-      operator()();
+      operator()()
+      {
+       _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
+       return _M_x;
+      }
 
       /**
        * @brief Compares two linear congruential random number generator
@@ -1588,12 +1592,7 @@ namespace std
       template<typename _UniformRandomNumberGenerator>
        result_type
        operator()(_UniformRandomNumberGenerator& __urng)
-       {
-         typedef typename _UniformRandomNumberGenerator::result_type
-           _UResult_type;
-         return _M_call(__urng, this->a(), this->b(),
-                        typename is_integral<_UResult_type>::type());
-       }
+        { return this->operator()(__urng, this->param()); }
 
       /**
        * Gets a uniform random number in the range @f$[0, n)@f$.
@@ -1765,11 +1764,7 @@ namespace std
       template<typename _UniformRandomNumberGenerator>
        result_type
        operator()(_UniformRandomNumberGenerator& __urng)
-       {
-         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
-           __aurng(__urng);
-         return (__aurng() * (this->b() - this->a())) + this->a();
-       }
+        { return this->operator()(__urng, this->param()); }
 
       template<typename _UniformRandomNumberGenerator>
        result_type
@@ -2014,12 +2009,12 @@ namespace std
       explicit
       lognormal_distribution(_RealType __m = _RealType(0),
                             _RealType __s = _RealType(1))
-      : _M_param(__m, __s)
+      : _M_param(__m, __s), _M_nd()
       { }
 
       explicit
       lognormal_distribution(const param_type& __p)
-      : _M_param(__p)
+      : _M_param(__p), _M_nd()
       { }
 
       /**
@@ -2027,7 +2022,7 @@ namespace std
        */
       void
       reset()
-      { }
+      { _M_nd.reset(); }
 
       /**
        *
@@ -2077,10 +2072,13 @@ namespace std
       template<typename _UniformRandomNumberGenerator>
        result_type
        operator()(_UniformRandomNumberGenerator& __urng,
-                  const param_type& __p);
+                  const param_type& __p)
+        { return std::exp(__p.s() * _M_nd(__urng) + __p.m()); }
 
     private:
       param_type _M_param;
+
+      normal_distribution<result_type> _M_nd;
     };
 
   /**
@@ -2555,12 +2553,12 @@ namespace std
 
       explicit
       student_t_distribution(_RealType __n = _RealType(1))
-      : _M_param(__n)
+      : _M_param(__n), _M_nd()
       { }
 
       explicit
       student_t_distribution(const param_type& __p)
-      : _M_param(__p)
+      : _M_param(__p), _M_nd()
       { }
 
       /**
@@ -2568,7 +2566,7 @@ namespace std
        */
       void
       reset()
-      { }
+      { _M_nd.reset(); }
 
       /**
        *
@@ -2617,12 +2615,9 @@ namespace std
                   const param_type& __p);
 
     private:
-      template<typename _UniformRandomNumberGenerator>
-       result_type
-       _M_gaussian(_UniformRandomNumberGenerator& __urng,
-                   const result_type __sigma);
-
       param_type _M_param;
+
+      normal_distribution<result_type> _M_nd;
     };
 
   /**
@@ -2761,14 +2756,7 @@ namespace std
     template<typename _UniformRandomNumberGenerator>
       result_type
       operator()(_UniformRandomNumberGenerator& __urng)
-      {
-       __detail::_Adaptor<_UniformRandomNumberGenerator, double>
-         __aurng(__urng);
-       if ((__aurng() - __aurng.min())
-            < this->p() * (__aurng.max() - __aurng.min()))
-         return true;
-       return false;
-      }
+      { return this->operator()(__urng, this->param()); }
 
     template<typename _UniformRandomNumberGenerator>
       result_type
@@ -3539,11 +3527,7 @@ namespace std
       template<typename _UniformRandomNumberGenerator>
        result_type
        operator()(_UniformRandomNumberGenerator& __urng)
-       {
-         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
-           __aurng(__urng);
-         return -std::log(__aurng()) / this->lambda();
-       }
+        { return this->operator()(__urng, this->param()); }
 
       template<typename _UniformRandomNumberGenerator>
        result_type
@@ -3633,8 +3617,7 @@ namespace std
        _RealType _M_alpha;
        _RealType _M_beta;
 
-       // Hosts either lambda of GB or d of modified Vaduva's.
-       _RealType _M_l_d;
+       _RealType _M_malpha, _M_a2;
       };
 
     public:
@@ -3645,21 +3628,20 @@ namespace std
       explicit
       gamma_distribution(_RealType __alpha_val = _RealType(1),
                         _RealType __beta_val = _RealType(1))
-      : _M_param(__alpha_val, __beta_val)
+      : _M_param(__alpha_val, __beta_val), _M_nd()
       { }
 
       explicit
       gamma_distribution(const param_type& __p)
-      : _M_param(__p)
+      : _M_param(__p), _M_nd()
       { }
 
       /**
        * @brief Resets the distribution state.
-       *
-       * Does nothing for the gamma distribution.
        */
       void
-      reset() { }
+      reset()
+      { _M_nd.reset(); }
 
       /**
        * @brief Returns the @f$ \alpha @f$ of the distribution.
@@ -3716,6 +3698,8 @@ namespace std
 
     private:
       param_type _M_param;
+
+      normal_distribution<result_type> _M_nd;
     };
 
   /**
@@ -3854,13 +3838,7 @@ namespace std
       template<typename _UniformRandomNumberGenerator>
        result_type
        operator()(_UniformRandomNumberGenerator& __urng,
-                  const param_type& __p)
-       {
-         __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
-           __aurng(__urng);
-         return __p.b() * std::pow(-std::log(__aurng()),
-                                   result_type(1) / __p.a());
-       }
+                  const param_type& __p);
 
     private:
       param_type _M_param;
@@ -4222,7 +4200,9 @@ namespace std
        typedef piecewise_constant_distribution<_RealType> distribution_type;
        friend class piecewise_constant_distribution<_RealType>;
 
-       param_type();
+       param_type()
+       : _M_int(), _M_den(), _M_cp()
+       { _M_initialize(); }
 
        template<typename _InputIteratorB, typename _InputIteratorW>
          param_type(_InputIteratorB __bfirst,
index b110f99..b933f6d 100644 (file)
@@ -128,19 +128,6 @@ namespace std
       seed(__sum);
     }
 
-  /**
-   * Gets the next generated value in sequence.
-   */
-  template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
-    typename linear_congruential_engine<_UIntType, __a, __c, __m>::
-            result_type
-    linear_congruential_engine<_UIntType, __a, __c, __m>::
-    operator()()
-    {
-      _M_x = __detail::__mod<_UIntType, __a, __c, __m>(_M_x);
-      return _M_x;
-    }
-
   template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
           typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
@@ -556,7 +543,7 @@ namespace std
     {
       const long double __r = static_cast<long double>(_M_b.max())
                            - static_cast<long double>(_M_b.min()) + 1.0L;
-      const result_type __m = std::log10(__r) / std::log10(2.0L);
+      const result_type __m = std::log(__r) / std::log(2.0L);
       result_type __n, __n0, __y0, __y1, __s0, __s1;
       for (size_t __i = 0; __i < 2; ++__i)
        {
@@ -874,17 +861,12 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __p)
       {
-       typename gamma_distribution<>::param_type
-         __gamma_param(__p.k(), 1.0);
-       gamma_distribution<> __gamma(__gamma_param);
+       gamma_distribution<> __gamma(__p.k(), 1.0);
        double __x = __gamma(__urng);
 
-       typename poisson_distribution<result_type>::param_type
-         __poisson_param(__x * __p.p() / (1.0 - __p.p()));
-       poisson_distribution<result_type> __poisson(__poisson_param);
-       result_type __m = __poisson(__urng);
-
-       return __m;
+       poisson_distribution<result_type> __poisson(__x * __p.p()
+                                                   / (1.0 - __p.p()));
+       return __poisson(__urng);
       }
 
   template<typename _IntType, typename _CharT, typename _Traits>
@@ -1510,33 +1492,6 @@ namespace std
     }
 
 
-  template<typename _RealType>
-    template<typename _UniformRandomNumberGenerator>
-      typename lognormal_distribution<_RealType>::result_type
-      lognormal_distribution<_RealType>::
-      operator()(_UniformRandomNumberGenerator& __urng,
-                const param_type& __p)
-      {
-       _RealType __u, __v, __r2, __normal;
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
-         __aurng(__urng);
-
-       do
-         {
-           // Choose x,y in uniform square (-1,-1) to (+1,+1).
-           __u = 2 * __aurng() - 1;
-           __v = 2 * __aurng() - 1;
-
-           // See if it is in the unit circle.
-           __r2 = __u * __u + __v * __v;
-         }
-       while (__r2 > 1 || __r2 == 0);
-
-       __normal = __u * std::sqrt(-2 * std::log(__r2) / __r2);
-
-       return std::exp(__p.s() * __normal + __p.m());
-      }
-
   template<typename _RealType, typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
@@ -1553,7 +1508,8 @@ namespace std
       __os.fill(__space);
       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
 
-      __os << __x.m() << __space << __x.s();
+      __os << __x.m() << __space << __x.s()
+          << __space << __x._M_nd;
 
       __os.flags(__flags);
       __os.fill(__fill);
@@ -1573,7 +1529,7 @@ namespace std
       __is.flags(__ios_base::dec | __ios_base::skipws);
 
       _RealType __m, __s;
-      __is >> __m >> __s;
+      __is >> __m >> __s >> __x._M_nd;
       __x.param(typename lognormal_distribution<_RealType>::
                param_type(__m, __s));
 
@@ -1589,9 +1545,7 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __p)
       {
-       typename gamma_distribution<_RealType>::param_type
-         __gamma_param(__p.n() / 2, 1.0);
-       gamma_distribution<_RealType> __gamma(__gamma_param);
+       gamma_distribution<_RealType> __gamma(__p.n() / 2, 1.0);
        return 2 * __gamma(__urng);
       }
 
@@ -1765,35 +1719,6 @@ namespace std
     }
 
 
-  //
-  //  This could be operator() for a Gaussian distribution.
-  //
-  template<typename _RealType>
-    template<typename _UniformRandomNumberGenerator>
-      typename student_t_distribution<_RealType>::result_type
-      student_t_distribution<_RealType>::
-      _M_gaussian(_UniformRandomNumberGenerator& __urng,
-                 const result_type __sigma)
-      {
-       _RealType __x, __y, __r2;
-       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
-         __aurng(__urng);
-
-       do
-         {
-           // Choose x,y in uniform square (-1,-1) to (+1,+1).
-           __x = 2 * __aurng() - 1;
-           __y = 2 * __aurng() - 1;
-
-           // See if it is in the unit circle.
-           __r2 = __x * __x + __y * __y;
-         }
-       while (__r2 > 1 || __r2 == 0);
-
-       // Box-Muller transform.
-       return __sigma * __y * std::sqrt(-2 * std::log(__r2) / __r2);
-      }
-
   template<typename _RealType>
     template<typename _UniformRandomNumberGenerator>
       typename student_t_distribution<_RealType>::result_type
@@ -1803,10 +1728,8 @@ namespace std
       {
        if (__param.n() <= 2.0)
          {
-           _RealType __y1 = _M_gaussian(__urng, 1.0);
-           typename chi_squared_distribution<_RealType>::param_type
-             __chisq_param(__param.n());
-           chi_squared_distribution<_RealType> __chisq(__chisq_param);
+           _RealType __y1 = _M_nd(__urng);
+           chi_squared_distribution<_RealType> __chisq(__param.n());
            _RealType __y2 = __chisq(__urng);
 
            return __y1 / std::sqrt(__y2 / __param.n());
@@ -1814,13 +1737,12 @@ namespace std
        else
          {
            _RealType __y1, __y2, __z;
+           exponential_distribution<_RealType>
+             __exponential(1.0 / (__param.n() / 2.0 - 1.0));
+
            do
              {
-               __y1 = _M_gaussian(__urng, 1.0);
-               typename exponential_distribution<_RealType>::param_type
-                 __exp_param(1.0 / (__param.n() / 2.0 - 1.0));
-               exponential_distribution<_RealType>
-                 __exponential(__exp_param);
+               __y1 = _M_nd(__urng);
                __y2 = __exponential(__urng);
 
                __z = __y1 * __y1 / (__param.n() - 2.0);
@@ -1850,7 +1772,7 @@ namespace std
       __os.fill(__space);
       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
 
-      __os << __x.n();
+      __os << __x.n() << __space << __x._M_nd;
 
       __os.flags(__flags);
       __os.fill(__fill);
@@ -1870,7 +1792,7 @@ namespace std
       __is.flags(__ios_base::dec | __ios_base::skipws);
 
       _RealType __n;
-      __is >> __n;
+      __is >> __n >> __x._M_nd;
       __x.param(typename student_t_distribution<_RealType>::param_type(__n));
 
       __is.flags(__flags);
@@ -1883,28 +1805,16 @@ namespace std
     gamma_distribution<_RealType>::param_type::
     _M_initialize()
     {
-      if (_M_alpha >= 1)
-       _M_l_d = std::sqrt(2 * _M_alpha - 1);
-      else
-       _M_l_d = (std::pow(_M_alpha, _M_alpha / (1 - _M_alpha))
-                 * (1 - _M_alpha));
+      _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
+
+      const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
+      _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
     }
 
   /**
-   * Cheng's rejection algorithm GB for alpha >= 1 and a modification
-   * of Vaduva's rejection from Weibull algorithm due to Devroye for
-   * alpha < 1.
-   *
-   * References:
-   * Cheng, R. C. "The Generation of Gamma Random Variables with Non-integral
-   * Shape Parameter." Applied Statistics, 26, 71-75, 1977.
-   *
-   * Vaduva, I. "Computer Generation of Gamma Gandom Variables by Rejection
-   * and Composition Procedures." Math. Operationsforschung and Statistik,
-   * Series in Statistics, 8, 545-576, 1977.
-   *
-   * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
-   * New York, 1986, Ch. IX, Sect. 3.4 (+ Errata!).
+   * Marsaglia, G. and Tsang, W. W.
+   * "A Simple Method for Generating Gamma Variables"
+   * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
    */
   template<typename _RealType>
     template<typename _UniformRandomNumberGenerator>
@@ -1913,58 +1823,40 @@ namespace std
       operator()(_UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
       {
-       result_type __x;
        __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
          __aurng(__urng);
 
-       bool __reject;
-       const _RealType __alpha_val = __param.alpha();
-       const _RealType __beta_val = __param.beta();
-       if (__alpha_val >= 1)
-         {
-           // alpha - log(4)
-           const result_type __b = __alpha_val
-             - result_type(1.3862943611198906188344642429163531L);
-           const result_type __c = __alpha_val + __param._M_l_d;
-           const result_type __1l = 1 / __param._M_l_d;
-
-           // 1 + log(9 / 2)
-           const result_type __k = 2.5040773967762740733732583523868748L;
+       result_type __u, __v, __n;
+       const result_type __a1 = (__param._M_malpha
+                                 - _RealType(1.0) / _RealType(3.0));
 
+       do
+         {
            do
              {
-               const result_type __u = __aurng() / __beta_val;
-               const result_type __v = __aurng() / __beta_val;
-
-               const result_type __y = __1l * std::log(__v / (1 - __v));
-               __x = __alpha_val * std::exp(__y);
-
-               const result_type __z = __u * __v * __v;
-               const result_type __r = __b + __c * __y - __x;
-
-               __reject = __r < result_type(4.5) * __z - __k;
-               if (__reject)
-                 __reject = __r < std::log(__z);
+               __n = _M_nd(__urng);
+               __v = result_type(1.0) + __param._M_a2 * __n; 
              }
-           while (__reject);
+           while (__v <= 0.0);
+
+           __v = __v * __v * __v;
+           __u = __aurng();
          }
+       while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n
+              && (std::log(__u) > (0.5 * __n * __n + __a1
+                                   * (1.0 - __v + std::log(__v)))));
+
+       if (__param.alpha() == __param._M_malpha)
+         return __a1 * __v * __param.beta();
        else
          {
-           const result_type __c = 1 / __alpha_val;
-
            do
-             {
-               const result_type __z = -std::log(__aurng() / __beta_val);
-               const result_type __e = -std::log(__aurng() / __beta_val);
-
-               __x = std::pow(__z, __c);
-
-               __reject = __z + __e < __param._M_l_d + __x;
-             }
-           while (__reject);
+             __u = __aurng();
+           while (__u == 0.0);
+           
+           return (std::pow(__u, result_type(1.0) / __param.alpha())
+                   * __a1 * __v * __param.beta());
          }
-
-       return __beta_val * __x;
       }
 
   template<typename _RealType, typename _CharT, typename _Traits>
@@ -1983,7 +1875,8 @@ namespace std
       __os.fill(__space);
       __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
 
-      __os << __x.alpha() << __space << __x.beta();
+      __os << __x.alpha() << __space << __x.beta()
+          << __space << __x._M_nd;
 
       __os.flags(__flags);
       __os.fill(__fill);
@@ -2003,7 +1896,7 @@ namespace std
       __is.flags(__ios_base::dec | __ios_base::skipws);
 
       _RealType __alpha_val, __beta_val;
-      __is >> __alpha_val >> __beta_val;
+      __is >> __alpha_val >> __beta_val >> __x._M_nd;
       __x.param(typename gamma_distribution<_RealType>::
                param_type(__alpha_val, __beta_val));
 
@@ -2012,6 +1905,19 @@ namespace std
     }
 
 
+  template<typename _RealType>
+    template<typename _UniformRandomNumberGenerator>
+      typename weibull_distribution<_RealType>::result_type
+      weibull_distribution<_RealType>::
+      operator()(_UniformRandomNumberGenerator& __urng,
+                const param_type& __p)
+      {
+       __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+         __aurng(__urng);
+       return __p.b() * std::pow(-std::log(__aurng()),
+                                 result_type(1) / __p.a());
+      }
+
   template<typename _RealType, typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
     operator<<(std::basic_ostream<_CharT, _Traits>& __os,
@@ -2266,12 +2172,6 @@ namespace std
     }
 
   template<typename _RealType>
-    piecewise_constant_distribution<_RealType>::param_type::
-    param_type()
-    : _M_int(), _M_den(), _M_cp()
-    { _M_initialize(); }
-
-  template<typename _RealType>
     template<typename _InputIteratorB, typename _InputIteratorW>
       piecewise_constant_distribution<_RealType>::param_type::
       param_type(_InputIteratorB __bbegin,
@@ -2315,8 +2215,7 @@ namespace std
   template<typename _RealType>
     template<typename _Func>
       piecewise_constant_distribution<_RealType>::param_type::
-      param_type(size_t __nw, _RealType __xmin, _RealType __xmax,
-                _Func __fw)
+      param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
       : _M_int(), _M_den(), _M_cp()
       {
        for (size_t __i = 0; __i <= __nw; ++__i)
@@ -2713,7 +2612,7 @@ namespace std
                    __bits);
       const long double __r = static_cast<long double>(__urng.max())
                            - static_cast<long double>(__urng.min()) + 1.0L;
-      const size_t __log2r = std::log10(__r) / std::log10(2.0L);
+      const size_t __log2r = std::log(__r) / std::log(2.0L);
       size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r);
       _RealType __sum = _RealType(0);
       _RealType __tmp = _RealType(1);