re PR libstdc++/29520 (tr1: discrete_distributions vs large floating point values)

2006-10-29  Paolo Carlini  <pcarlini@suse.de>

	PR libstdc++/29520
	* include/tr1/random (geometric_distribution<>::
	operator()(_UniformRandomNumberGenerator&)): Only declare.
	* include/tr1/random.tcc (geometric_distribution<>::
	operator()(_UniformRandomNumberGenerator&),
	poisson_distribution<>::operator()(_UniformRandomNumberGenerator&),
	binomial_distribution<>::operator()(_UniformRandomNumberGenerator&)):
	Reject candidate floating point values not convertible to the
	result_type.

From-SVN: r118135

diff --git a/libstdc++-v3/ChangeLog b/libstdc++-v3/ChangeLog
index cb13016..a6a0589 100644 (file)
--- a/libstdc++-v3/ChangeLog
+++ b/libstdc++-v3/ChangeLog
@@ -1,3 +1,15 @@
+2006-10-29  Paolo Carlini  <pcarlini@suse.de>
+       
+       PR libstdc++/29520
+       * include/tr1/random (geometric_distribution<>::
+       operator()(_UniformRandomNumberGenerator&)): Only declare.
+       * include/tr1/random.tcc (geometric_distribution<>::
+       operator()(_UniformRandomNumberGenerator&),
+       poisson_distribution<>::operator()(_UniformRandomNumberGenerator&),
+       binomial_distribution<>::operator()(_UniformRandomNumberGenerator&)):
+       Reject candidate floating point values not convertible to the
+       result_type.
+
 2006-10-28  Paolo Carlini  <pcarlini@suse.de>
        
        * include/tr1/array (array<>::_M_at): New.

diff --git a/libstdc++-v3/include/tr1/random b/libstdc++-v3/include/tr1/random
index 9d1aa84..9114d80 100644 (file)
--- a/libstdc++-v3/include/tr1/random
+++ b/libstdc++-v3/include/tr1/random
@@ -1759,8 +1759,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
 
       template<class _UniformRandomNumberGenerator>
         result_type
-        operator()(_UniformRandomNumberGenerator& __urng)
-        { return result_type(std::ceil(std::log(__urng()) / _M_log_p)); }
+        operator()(_UniformRandomNumberGenerator& __urng);
 
       /**
        * Inserts a %geometric_distribution random number distribution

diff --git a/libstdc++-v3/include/tr1/random.tcc b/libstdc++-v3/include/tr1/random.tcc
index 8624f12..67d79ed 100644 (file)
--- a/libstdc++-v3/include/tr1/random.tcc
+++ b/libstdc++-v3/include/tr1/random.tcc
@@ -768,6 +768,28 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
     }
 
 
+  template<typename _IntType, typename _RealType>
+    template<class _UniformRandomNumberGenerator>
+      typename geometric_distribution<_IntType, _RealType>::result_type
+      geometric_distribution<_IntType, _RealType>::
+      operator()(_UniformRandomNumberGenerator& __urng)
+      {
+       // About the epsilon thing see this thread:
+        // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
+       const _RealType __naf =
+         (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
+       // The largest _RealType convertible to _IntType.
+       const _RealType __thr =
+         std::numeric_limits<_IntType>::max() + __naf;
+
+       _RealType __cand;
+       do
+         __cand = std::ceil(std::log(__urng()) / _M_log_p);
+       while (__cand >= __thr);
+
+       return result_type(__cand + __naf);
+      }
+
   template<typename _IntType, typename _RealType,
           typename _CharT, typename _Traits>
     std::basic_ostream<_CharT, _Traits>&
@@ -841,6 +863,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
          {
            _RealType __x;
 
+           // See comments above...
+           const _RealType __naf =
+             (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
+           const _RealType __thr =
+             std::numeric_limits<_IntType>::max() + __naf;
+
            const _RealType __m = std::floor(_M_mean);
            // sqrt(pi / 2)
            const _RealType __spi_2 = 1.2533141373155002512078826424055226L;
@@ -899,9 +927,11 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
                __reject = (__w - __e - __x * _M_lm_thr
                            > _M_lfm - std::tr1::lgamma(__x + __m + 1));
 
+               __reject |= __x + __m >= __thr;
+
              } while (__reject);
 
-           return _IntType(__x + __m + 0.5);
+           return result_type(__x + __m + __naf);
          }
        else
 #endif
@@ -1055,6 +1085,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
          {
            _RealType __x;
 
+           // See comments above...
+           const _RealType __naf =
+             (1 - std::numeric_limits<_RealType>::epsilon()) / 2;
+           const _RealType __thr =
+             std::numeric_limits<_IntType>::max() + __naf;
+
            const _RealType __np = std::floor(_M_t * __p12);
            const _RealType __pa = __np / _M_t;
 
@@ -1127,10 +1163,12 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
                      + std::tr1::lgamma(_M_t - (__np + __x) + 1);
                    __reject = __v > _M_lf - __lfx + __x * _M_lp1p;
                  }
+
+               __reject |= __x + __np >= __thr;
              }
            while (__reject);
 
-           __x += __np + 0.5;
+           __x += __np + __naf;
 
            const _IntType __z = _M_waiting(__urng, _M_t - _IntType(__x)); 
            __ret = _IntType(__x) + __z;