Revisited [rand.dist.bern.bin] and [rand.dist.pois.poisson] with better algorithms

git-svn-id: https://llvm.org/svn/llvm-project/libcxx/trunk@103886 91177308-0d34-0410-b5e6-96231b3b80d8

include/random

@@ -3143,11 +3143,13 @@ public:
     {
         result_type __t_;
         double __p_;
+        double __pr_;
+        double __odds_ratio_;
+        result_type __r0_;
     public:
         typedef binomial_distribution distribution_type;
 
-        explicit param_type(result_type __t = 1, double __p = 0.5)
+        explicit param_type(result_type __t = 1, double __p = 0.5);
-            : __t_(__t), __p_(__p) {}
 
         result_type t() const {return __t_;}
         double p() const {return __p_;}
@@ -3156,6 +3158,8 @@ public:
             {return __x.__t_ == __y.__t_ && __x.__p_ == __y.__p_;}
         friend bool operator!=(const param_type& __x, const param_type& __y)
             {return !(__x == __y);}
+
+        friend class binomial_distribution;
     };
 
 private:
@@ -3191,17 +3195,56 @@ public:
         {return !(__x == __y);}
 };
 
+template<class _IntType>
+binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p)
+    : __t_(__t), __p_(__p)
+{
+    if (0 < __p_ && __p_ < 1)
+    {
+        __r0_ = static_cast<result_type>((__t_ + 1) * __p_);
+        __pr_ = _STD::exp(_STD::lgamma(__t_ + 1.) - _STD::lgamma(__r0_ + 1.) -
+                          _STD::lgamma(__t_ - __r0_ + 1.) + __r0_ * _STD::log(__p_) +
+                          (__t_ - __r0_) * _STD::log(1 - __p_));
+        __odds_ratio_ = __p_ / (1 - __p_);
+    }
+}
+
 template<class _IntType>
 template<class _URNG>
 _IntType
-binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __p)
+binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr)
 {
-    bernoulli_distribution __bd(__p.p());
+    if (__pr.__t_ == 0 || __pr.__p_ == 0)
-    _IntType __r = 0;
+        return 0;
-    _IntType __t = __p.t();
+    if (__pr.__p_ == 1)
-    for (_IntType __i = 0; __i < __t; ++__i)
+        return __pr.__t_;
-        __r += __bd(__g);
+    uniform_real_distribution<double> __gen;
-    return __r;
+    double __u = __gen(__g) - __pr.__pr_;
+    if (__u < 0)
+        return __pr.__r0_;
+    double __pu = __pr.__pr_;
+    double __pd = __pu;
+    result_type __ru = __pr.__r0_;
+    result_type __rd = __ru;
+    while (true)
+    {
+        if (__rd >= 1)
+        {
+            __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1));
+            __u -= __pd;
+            if (__u < 0)
+                return __rd - 1;
+        }
+        --__rd;
+        ++__ru;
+        if (__ru <= __pr.__t_)
+        {
+            __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru;
+            __u -= __pu;
+            if (__u < 0)
+                return __ru;
+        }
+    }
 }
 
 template <class _CharT, class _Traits, class _IntType>
@@ -3234,149 +3277,6 @@ operator>>(basic_istream<_CharT, _Traits>& __is,
     return __is;
 }
 
-// poisson_distribution
-
-template<class _IntType = int>
-class poisson_distribution
-{
-public:
-    // types
-    typedef _IntType result_type;
-
-    class param_type
-    {
-        double __mean_;
-        double __sq_;
-        double __alxm_;
-        double __g_;
-    public:
-        typedef poisson_distribution distribution_type;
-
-        explicit param_type(double __mean = 1.0);
-
-        double mean() const {return __mean_;}
-
-        friend bool operator==(const param_type& __x, const param_type& __y)
-            {return __x.__mean_ == __y.__mean_;}
-        friend bool operator!=(const param_type& __x, const param_type& __y)
-            {return !(__x == __y);}
-
-        friend class poisson_distribution;
-    };
-
-private:
-    param_type __p_;
-
-public:
-    // constructors and reset functions
-    explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {}
-    explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
-    void reset() {}
-
-    // generating functions
-    template<class _URNG> result_type operator()(_URNG& __g)
-        {return (*this)(__g, __p_);}
-    template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
-
-    // property functions
-    double mean() const {return __p_.mean();}
-
-    param_type param() const {return __p_;}
-    void param(const param_type& __p) {__p_ = __p;}
-
-    result_type min() const {return 0;}
-    result_type max() const {return numeric_limits<result_type>::max();}
-
-    friend bool operator==(const poisson_distribution& __x,
-                           const poisson_distribution& __y)
-        {return __x.__p_ == __y.__p_;}
-    friend bool operator!=(const poisson_distribution& __x,
-                           const poisson_distribution& __y)
-        {return !(__x == __y);}
-};
-
-template<class _IntType>
-poisson_distribution<_IntType>::param_type::param_type(double __mean)
-    : __mean_(__mean)
-{
-    if (__mean_ < 12.0)
-    {
-        __g_ = _STD::exp(-__mean_);
-        __alxm_ = 0;
-        __sq_ = 0;
-    }
-    else
-    {
-        __sq_ = _STD::sqrt(2.0 * __mean_);
-        __alxm_ = _STD::log(__mean_);
-        __g_ = __mean_ * __alxm_ - _STD::lgamma(__mean_ + 1);
-    }
-}
-
-template <class _IntType>
-template<class _URNG>
-_IntType
-poisson_distribution<_IntType>::operator()(_URNG& __g, const param_type& __p)
-{
-    result_type __x;
-    uniform_real_distribution<double> __gen;
-    if (__p.__mean_ < 12.0)
-    {
-        __x = result_type(~0);
-        double __t = 1;
-        do
-        {
-            ++__x;
-            __t *= __gen(__g);
-        } while (__t > __p.__g_);
-    }
-    else
-    {
-        double __t;
-        const double __pi = 3.14159265358979323846264338328;
-        do
-        {
-            double _X;
-            double __y;
-            do
-            {
-                __y = _STD::tan(__pi * __gen(__g));
-                _X = __p.__sq_ * __y + __p.__mean_;
-            } while (_X < 0);
-            __x = static_cast<result_type>(_X);
-            __t = 0.9 * (1 + __y * __y) * _STD::exp(__x * __p.__alxm_ -
-                                            _STD::lgamma(__x + 1.0) - __p.__g_);
-        } while (__gen(__g) > __t);
-    }
-    return __x;
-}
-
-template <class _CharT, class _Traits, class _IntType>
-basic_ostream<_CharT, _Traits>&
-operator<<(basic_ostream<_CharT, _Traits>& __os,
-           const poisson_distribution<_IntType>& __x)
-{
-    __save_flags<_CharT, _Traits> _(__os);
-    __os.flags(ios_base::dec | ios_base::left);
-    return __os << __x.mean();
-}
-
-template <class _CharT, class _Traits, class _IntType>
-basic_istream<_CharT, _Traits>&
-operator>>(basic_istream<_CharT, _Traits>& __is,
-           poisson_distribution<_IntType>& __x)
-{
-    typedef poisson_distribution<_IntType> _Eng;
-    typedef typename _Eng::param_type param_type;
-    __save_flags<_CharT, _Traits> _(__is);
-    __is.flags(ios_base::dec | ios_base::skipws);
-    double __mean;
-    __is >> __mean;
-    if (!__is.fail())
-        __x.param(param_type(__mean));
-    return __is;
-}
-
 // exponential_distribution
 
 template<class _RealType = double>
@@ -3477,6 +3377,370 @@ operator>>(basic_istream<_CharT, _Traits>& __is,
     return __is;
 }
 
+// normal_distribution
+
+template<class _RealType = double>
+class normal_distribution
+{
+public:
+    // types
+    typedef _RealType result_type;
+
+    class param_type
+    {
+        result_type __mean_;
+        result_type __stddev_;
+    public:
+        typedef normal_distribution distribution_type;
+
+        explicit param_type(result_type __mean = 0, result_type __stddev = 1)
+            : __mean_(__mean), __stddev_(__stddev) {}
+
+        result_type mean() const {return __mean_;}
+        result_type stddev() const {return __stddev_;}
+
+        friend bool operator==(const param_type& __x, const param_type& __y)
+            {return __x.__mean_ == __y.__mean_ && __x.__stddev_ == __y.__stddev_;}
+        friend bool operator!=(const param_type& __x, const param_type& __y)
+            {return !(__x == __y);}
+    };
+
+private:
+    param_type __p_;
+    result_type _V_;
+    bool _V_hot_;
+
+public:
+    // constructors and reset functions
+    explicit normal_distribution(result_type __mean = 0, result_type __stddev = 1)
+        : __p_(param_type(__mean, __stddev)), _V_hot_(false) {}
+    explicit normal_distribution(const param_type& __p)
+        : __p_(__p), _V_hot_(false) {}
+    void reset() {_V_hot_ = false;}
+
+    // generating functions
+    template<class _URNG> result_type operator()(_URNG& __g)
+        {return (*this)(__g, __p_);}
+    template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
+
+    // property functions
+    result_type mean() const {return __p_.mean();}
+    result_type stddev() const {return __p_.stddev();}
+
+    param_type param() const {return __p_;}
+    void param(const param_type& __p) {__p_ = __p;}
+
+    result_type min() const {return -numeric_limits<result_type>::infinity();}
+    result_type max() const {return numeric_limits<result_type>::infinity();}
+
+    friend bool operator==(const normal_distribution& __x,
+                           const normal_distribution& __y)
+        {return __x.__p_ == __y.__p_ && __x._V_hot_ == __y._V_hot_ &&
+                (!__x._V_hot_ || __x._V_ == __y._V_);}
+    friend bool operator!=(const normal_distribution& __x,
+                           const normal_distribution& __y)
+        {return !(__x == __y);}
+
+    template <class _CharT, class _Traits, class _RT>
+    friend
+    basic_ostream<_CharT, _Traits>&
+    operator<<(basic_ostream<_CharT, _Traits>& __os,
+               const normal_distribution<_RT>& __x);
+    
+    template <class _CharT, class _Traits, class _RT>
+    friend
+    basic_istream<_CharT, _Traits>&
+    operator>>(basic_istream<_CharT, _Traits>& __is,
+               normal_distribution<_RT>& __x);
+};
+
+template <class _RealType>
+template<class _URNG>
+_RealType
+normal_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
+{
+    result_type _U;
+    if (_V_hot_)
+    {
+        _V_hot_ = false;
+        _U = _V_;
+    }
+    else
+    {
+        uniform_real_distribution<result_type> _Uni(-1, 1);
+        result_type __u;
+        result_type __v;
+        result_type __s;
+        do
+        {
+            __u = _Uni(__g);
+            __v = _Uni(__g);
+            __s = __u * __u + __v * __v;
+        } while (__s > 1 || __s == 0);
+        result_type _F = _STD::sqrt(-2 * _STD::log(__s) / __s);
+        _V_ = __v * _F;
+        _V_hot_ = true;
+        _U = __u * _F;
+    }
+    return _U * __p.stddev() + __p.mean();
+}
+
+template <class _CharT, class _Traits, class _RT>
+basic_ostream<_CharT, _Traits>&
+operator<<(basic_ostream<_CharT, _Traits>& __os,
+           const normal_distribution<_RT>& __x)
+{
+    __save_flags<_CharT, _Traits> _(__os);
+    __os.flags(ios_base::dec | ios_base::left);
+    _CharT __sp = __os.widen(' ');
+    __os.fill(__sp);
+    __os << __x.mean() << __sp << __x.stddev() << __sp << __x._V_hot_;
+    if (__x._V_hot_)
+        __os << __sp << __x._V_;
+    return __os;
+}
+
+template <class _CharT, class _Traits, class _RT>
+basic_istream<_CharT, _Traits>&
+operator>>(basic_istream<_CharT, _Traits>& __is,
+           normal_distribution<_RT>& __x)
+{
+    typedef normal_distribution<_RT> _Eng;
+    typedef typename _Eng::result_type result_type;
+    typedef typename _Eng::param_type param_type;
+    __save_flags<_CharT, _Traits> _(__is);
+    __is.flags(ios_base::dec | ios_base::skipws);
+    result_type __mean;
+    result_type __stddev;
+    result_type _V = 0;
+    bool _V_hot = false;
+    __is >> __mean >> __stddev >> _V_hot;
+    if (_V_hot)
+        __is >> _V;
+    if (!__is.fail())
+    {
+        __x.param(param_type(__mean, __stddev));
+        __x._V_hot_ = _V_hot;
+        __x._V_ = _V;
+    }
+    return __is;
+}
+
+// poisson_distribution
+
+template<class _IntType = int>
+class poisson_distribution
+{
+public:
+    // types
+    typedef _IntType result_type;
+
+    class param_type
+    {
+        double __mean_;
+        double __s_;
+        double __d_;
+        double __l_;
+        double __omega_;
+        double __c0_;
+        double __c1_;
+        double __c2_;
+        double __c3_;
+        double __c_;
+
+    public:
+        typedef poisson_distribution distribution_type;
+
+        explicit param_type(double __mean = 1.0);
+
+        double mean() const {return __mean_;}
+
+        friend bool operator==(const param_type& __x, const param_type& __y)
+            {return __x.__mean_ == __y.__mean_;}
+        friend bool operator!=(const param_type& __x, const param_type& __y)
+            {return !(__x == __y);}
+
+        friend class poisson_distribution;
+    };
+
+private:
+    param_type __p_;
+
+public:
+    // constructors and reset functions
+    explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {}
+    explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
+    void reset() {}
+
+    // generating functions
+    template<class _URNG> result_type operator()(_URNG& __g)
+        {return (*this)(__g, __p_);}
+    template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
+
+    // property functions
+    double mean() const {return __p_.mean();}
+
+    param_type param() const {return __p_;}
+    void param(const param_type& __p) {__p_ = __p;}
+
+    result_type min() const {return 0;}
+    result_type max() const {return numeric_limits<result_type>::max();}
+
+    friend bool operator==(const poisson_distribution& __x,
+                           const poisson_distribution& __y)
+        {return __x.__p_ == __y.__p_;}
+    friend bool operator!=(const poisson_distribution& __x,
+                           const poisson_distribution& __y)
+        {return !(__x == __y);}
+};
+
+template<class _IntType>
+poisson_distribution<_IntType>::param_type::param_type(double __mean)
+    : __mean_(__mean)
+{
+    if (__mean_ < 10)
+    {
+        __s_ = 0;
+        __d_ = 0;
+        __l_ = _STD::exp(-__mean_);
+        __omega_ = 0;
+        __c3_ = 0;
+        __c2_ = 0;
+        __c1_ = 0;
+        __c0_ = 0;
+        __c_ = 0;
+    }
+    else
+    {
+        __s_ = _STD::sqrt(__mean_);
+        __d_ = 6 * __mean_ * __mean_;
+        __l_ = static_cast<result_type>(__mean_ - 1.1484);
+        __omega_ = .3989423 / __s_;
+        double __b1_ = .4166667E-1 / __mean_;
+        double __b2_ = .3 * __b1_ * __b1_;
+        __c3_ = .1428571 * __b1_ * __b2_;
+        __c2_ = __b2_ - 15. * __c3_;
+        __c1_ = __b1_ - 6. * __b2_ + 45. * __c3_;
+        __c0_ = 1. - __b1_ + 3. * __b2_ - 15. * __c3_;
+        __c_ = .1069 / __mean_;
+    }
+}
+
+template <class _IntType>
+template<class _URNG>
+_IntType
+poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr)
+{
+    result_type __x;
+    uniform_real_distribution<double> __urd;
+    if (__pr.__mean_ <= 10)
+    {
+         __x = 0;
+        for (double __p = __urd(__urng); __p > __pr.__l_; ++__x)
+            __p *= __urd(__urng);
+    }
+    else
+    {
+        double __difmuk;
+        double __g = __pr.__mean_ + __pr.__s_ * normal_distribution<double>()(__urng);
+        double __u;
+        if (__g > 0)
+        {
+            __x = static_cast<result_type>(__g);
+            if (__x >= __pr.__l_)
+                return __x;
+            __difmuk = __pr.__mean_ - __x;
+            __u = __urd(__urng);
+            if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk)
+                return __x;
+        }
+        exponential_distribution<double> __edist;
+        for (bool __using_exp_dist = false; true; __using_exp_dist = true)
+        {
+            double __e;
+            if (__using_exp_dist || __g < 0)
+            {
+                double __t;
+                do
+                {
+                    __e = __edist(__urng);
+                    __u = __urd(__urng);
+                    __u += __u - 1;
+                    __t = 1.8 + (__u < 0 ? -__e : __e);
+                } while (__t <= -.6744);
+                __x = __pr.__mean_ + __pr.__s_ * __t;
+                __difmuk = __pr.__mean_ - __x;
+                __using_exp_dist = true;
+            }
+            double __px;
+            double __py;
+            if (__x < 10)
+            {
+                const result_type __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040,
+                                             40320, 362880};
+                __px = -__pr.__mean_;
+                __py = _STD::pow(__pr.__mean_, (double)__x) / __fac[__x];
+            }
+            else
+            {
+                double __del = .8333333E-1 / __x;
+                __del -= 4.8 * __del * __del * __del;
+                double __v = __difmuk / __x;
+                if (_STD::abs(__v) > 0.25)
+                    __px = __x * _STD::log(1 + __v) - __difmuk - __del;
+                else
+                    __px = __x * __v * __v * (((((((.1250060 * __v + -.1384794) *
+                           __v + .1421878) * __v + -.1661269) * __v + .2000118) *
+                           __v + -.2500068) * __v + .3333333) * __v + -.5) - __del;
+                __py = .3989423 / _STD::sqrt(__x);
+            }
+            double __r = (0.5 - __difmuk) / __pr.__s_;
+            double __r2 = __r * __r;
+            double __fx = -0.5 * __r2;
+            double __fy = __pr.__omega_ * (((__pr.__c3_ * __r2 + __pr.__c2_) *
+                                        __r2 + __pr.__c1_) * __r2 + __pr.__c0_);
+            if (__using_exp_dist)
+            {
+                if (__pr.__c_ * _STD::abs(__u) <= __py * _STD::exp(__px + __e) -
+                                                   __fy * _STD::exp(__fx + __e))
+                    break;
+            }
+            else
+            {
+                if (__fy - __u * __fy <= __py * _STD::exp(__px - __fx))
+                    break;
+            }
+        }
+    }
+    return __x;
+}
+
+template <class _CharT, class _Traits, class _IntType>
+basic_ostream<_CharT, _Traits>&
+operator<<(basic_ostream<_CharT, _Traits>& __os,
+           const poisson_distribution<_IntType>& __x)
+{
+    __save_flags<_CharT, _Traits> _(__os);
+    __os.flags(ios_base::dec | ios_base::left);
+    return __os << __x.mean();
+}
+
+template <class _CharT, class _Traits, class _IntType>
+basic_istream<_CharT, _Traits>&
+operator>>(basic_istream<_CharT, _Traits>& __is,
+           poisson_distribution<_IntType>& __x)
+{
+    typedef poisson_distribution<_IntType> _Eng;
+    typedef typename _Eng::param_type param_type;
+    __save_flags<_CharT, _Traits> _(__is);
+    __is.flags(ios_base::dec | ios_base::skipws);
+    double __mean;
+    __is >> __mean;
+    if (!__is.fail())
+        __x.param(param_type(__mean));
+    return __is;
+}
+
 // gamma_distribution
 
 template<class _RealType = double>
@@ -3629,154 +3893,6 @@ operator>>(basic_istream<_CharT, _Traits>& __is,
         __x.param(param_type(__alpha, __beta));
     return __is;
 }
-// normal_distribution
-
-template<class _RealType = double>
-class normal_distribution
-{
-public:
-    // types
-    typedef _RealType result_type;
-
-    class param_type
-    {
-        result_type __mean_;
-        result_type __stddev_;
-    public:
-        typedef normal_distribution distribution_type;
-
-        explicit param_type(result_type __mean = 0, result_type __stddev = 1)
-            : __mean_(__mean), __stddev_(__stddev) {}
-
-        result_type mean() const {return __mean_;}
-        result_type stddev() const {return __stddev_;}
-
-        friend bool operator==(const param_type& __x, const param_type& __y)
-            {return __x.__mean_ == __y.__mean_ && __x.__stddev_ == __y.__stddev_;}
-        friend bool operator!=(const param_type& __x, const param_type& __y)
-            {return !(__x == __y);}
-    };
-
-private:
-    param_type __p_;
-    result_type _V_;
-    bool _V_hot_;
-
-public:
-    // constructors and reset functions
-    explicit normal_distribution(result_type __mean = 0, result_type __stddev = 1)
-        : __p_(param_type(__mean, __stddev)), _V_hot_(false) {}
-    explicit normal_distribution(const param_type& __p)
-        : __p_(__p), _V_hot_(false) {}
-    void reset() {_V_hot_ = false;}
-
-    // generating functions
-    template<class _URNG> result_type operator()(_URNG& __g)
-        {return (*this)(__g, __p_);}
-    template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
-
-    // property functions
-    result_type mean() const {return __p_.mean();}
-    result_type stddev() const {return __p_.stddev();}
-
-    param_type param() const {return __p_;}
-    void param(const param_type& __p) {__p_ = __p;}
-
-    result_type min() const {return -numeric_limits<result_type>::infinity();}
-    result_type max() const {return numeric_limits<result_type>::infinity();}
-
-    friend bool operator==(const normal_distribution& __x,
-                           const normal_distribution& __y)
-        {return __x.__p_ == __y.__p_ && __x._V_hot_ == __y._V_hot_ &&
-                (!__x._V_hot_ || __x._V_ == __y._V_);}
-    friend bool operator!=(const normal_distribution& __x,
-                           const normal_distribution& __y)
-        {return !(__x == __y);}
-
-    template <class _CharT, class _Traits, class _RT>
-    friend
-    basic_ostream<_CharT, _Traits>&
-    operator<<(basic_ostream<_CharT, _Traits>& __os,
-               const normal_distribution<_RT>& __x);
-    
-    template <class _CharT, class _Traits, class _RT>
-    friend
-    basic_istream<_CharT, _Traits>&
-    operator>>(basic_istream<_CharT, _Traits>& __is,
-               normal_distribution<_RT>& __x);
-};
-
-template <class _RealType>
-template<class _URNG>
-_RealType
-normal_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
-{
-    result_type _U;
-    if (_V_hot_)
-    {
-        _V_hot_ = false;
-        _U = _V_;
-    }
-    else
-    {
-        uniform_real_distribution<result_type> _Uni(-1, 1);
-        result_type __u;
-        result_type __v;
-        result_type __s;
-        do
-        {
-            __u = _Uni(__g);
-            __v = _Uni(__g);
-            __s = __u * __u + __v * __v;
-        } while (__s > 1 || __s == 0);
-        result_type _F = _STD::sqrt(-2 * _STD::log(__s) / __s);
-        _V_ = __v * _F;
-        _V_hot_ = true;
-        _U = __u * _F;
-    }
-    return _U * __p.stddev() + __p.mean();
-}
-
-template <class _CharT, class _Traits, class _RT>
-basic_ostream<_CharT, _Traits>&
-operator<<(basic_ostream<_CharT, _Traits>& __os,
-           const normal_distribution<_RT>& __x)
-{
-    __save_flags<_CharT, _Traits> _(__os);
-    __os.flags(ios_base::dec | ios_base::left);
-    _CharT __sp = __os.widen(' ');
-    __os.fill(__sp);
-    __os << __x.mean() << __sp << __x.stddev() << __sp << __x._V_hot_;
-    if (__x._V_hot_)
-        __os << __sp << __x._V_;
-    return __os;
-}
-
-template <class _CharT, class _Traits, class _RT>
-basic_istream<_CharT, _Traits>&
-operator>>(basic_istream<_CharT, _Traits>& __is,
-           normal_distribution<_RT>& __x)
-{
-    typedef normal_distribution<_RT> _Eng;
-    typedef typename _Eng::result_type result_type;
-    typedef typename _Eng::param_type param_type;
-    __save_flags<_CharT, _Traits> _(__is);
-    __is.flags(ios_base::dec | ios_base::skipws);
-    result_type __mean;
-    result_type __stddev;
-    result_type _V = 0;
-    bool _V_hot = false;
-    __is >> __mean >> __stddev >> _V_hot;
-    if (_V_hot)
-        __is >> _V;
-    if (!__is.fail())
-    {
-        __x.param(param_type(__mean, __stddev));
-        __x._V_hot_ = _V_hot;
-        __x._V_ = _V;
-    }
-    return __is;
-}
 
 _LIBCPP_END_NAMESPACE_STD
 

test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp

@@ -14,22 +14,58 @@
 // template<class _URNG> result_type operator()(_URNG& g);
 
 #include <random>
+#include <numeric>
+#include <vector>
 #include <cassert>
 
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
+
 int main()
 {
     {
         typedef std::bernoulli_distribution D;
-        typedef std::minstd_rand0 G;
+        typedef std::minstd_rand G;
         G g;
         D d(.75);
-        int count = 0;
+        const int N = 100000;
-        for (int i = 0; i < 10000; ++i)
+        std::vector<D::result_type> u;
-        {
+        for (int i = 0; i < N; ++i)
-            bool u = d(g);
+            u.push_back(d(g));
-            if (u)
+        double mean = std::accumulate(u.begin(), u.end(),
-                ++count;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(count > 7400);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.p();
+        double x_var = d.p()*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+    }
+    {
+        typedef std::bernoulli_distribution D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(.25);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.p();
+        double x_var = d.p()*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
 }

test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp

@@ -14,24 +14,62 @@
 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
 
 #include <random>
+#include <numeric>
+#include <vector>
 #include <cassert>
 
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
+
 int main()
 {
     {
         typedef std::bernoulli_distribution D;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
+        typedef std::minstd_rand G;
         G g;
         D d(.75);
         P p(.25);
-        int count = 0;
+        const int N = 100000;
-        for (int i = 0; i < 10000; ++i)
+        std::vector<D::result_type> u;
-        {
+        for (int i = 0; i < N; ++i)
-            bool u = d(g, p);
+            u.push_back(d(g, p));
-            if (u)
+        double mean = std::accumulate(u.begin(), u.end(),
-                ++count;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(count < 2600);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = p.p();
+        double x_var = p.p()*(1-p.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+    }
+    {
+        typedef std::bernoulli_distribution D;
+        typedef D::param_type P;
+        typedef std::minstd_rand G;
+        G g;
+        D d(.25);
+        P p(.75);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g, p));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = p.p();
+        double x_var = p.p()*(1-p.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
 }

test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp

@@ -15,53 +15,218 @@
 // template<class _URNG> result_type operator()(_URNG& g);
 
 #include <random>
+#include <numeric>
+#include <vector>
 #include <cassert>
 
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
+
 int main()
 {
     {
         typedef std::binomial_distribution<> D;
-        typedef D::param_type P;
+        typedef std::minstd_rand G;
-        typedef std::minstd_rand0 G;
         G g;
-        D d(16, .25);
+        D d(5, .75);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g));
-            D::result_type u = d(g);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 4);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef D::param_type P;
+        typedef std::minstd_rand G;
-        typedef std::minstd_rand0 G;
         G g;
-        D d(16, .5);
+        D d(30, .03125);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g));
-            D::result_type u = d(g);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 8);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
-        typedef D::param_type P;
+        typedef std::minstd_rand G;
-        typedef std::minstd_rand0 G;
         G g;
-        D d(16, .75);
+        D d(40, .25);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g));
-            D::result_type u = d(g);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 12);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(40, 0);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(mean == x_mean);
+        assert(var == x_var);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(40, 1);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(mean == x_mean);
+        assert(var == x_var);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(400, 0.5);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(1, 0.5);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(0, 0.005);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(mean == x_mean);
+        assert(var == x_var);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(0, 0);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(mean == x_mean);
+        assert(var == x_var);
+    }
+    {
+        typedef std::binomial_distribution<> D;
+        typedef std::minstd_rand G;
+        G g;
+        D d(0, 1);
+        const int N = 100000;
+        std::vector<D::result_type> u;
+        for (int i = 0; i < N; ++i)
+            u.push_back(d(g));
+        double mean = std::accumulate(u.begin(), u.end(),
+                                              double(0)) / u.size();
+        double var = 0;
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = d.t() * d.p();
+        double x_var = x_mean*(1-d.p());
+        assert(mean == x_mean);
+        assert(var == x_var);
     }
 }

test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp

@@ -15,56 +15,84 @@
 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
 
 #include <random>
+#include <numeric>
+#include <vector>
 #include <cassert>
 
+template <class T>
+inline
+T
+sqr(T x)
+{
+    return x * x;
+}
+
 int main()
 {
     {
         typedef std::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
+        typedef std::minstd_rand G;
         G g;
         D d(16, .75);
-        P p(16, .25);
+        P p(5, .75);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g, p));
-            D::result_type u = d(g, p);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 4);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = p.t() * p.p();
+        double x_var = x_mean*(1-p.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
+        typedef std::minstd_rand G;
         G g;
         D d(16, .75);
-        P p(16, .5);
+        P p(30, .03125);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g, p));
-            D::result_type u = d(g, p);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 8);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = p.t() * p.p();
+        double x_var = x_mean*(1-p.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
     {
         typedef std::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand0 G;
+        typedef std::minstd_rand G;
         G g;
         D d(16, .75);
-        P p(16, .75);
+        P p(40, .25);
-        int count = 0;
+        const int N = 100000;
-        int r = 0;
+        std::vector<D::result_type> u;
-        for (int i = 0; i < 100; ++i)
+        for (int i = 0; i < N; ++i)
-        {
+            u.push_back(d(g, p));
-            D::result_type u = d(g, p);
+        double mean = std::accumulate(u.begin(), u.end(),
-            r += u;
+                                              double(0)) / u.size();
-        }
+        double var = 0;
-        assert(int(r/100. + .5) == 12);
+        for (int i = 0; i < u.size(); ++i)
+            var += sqr(u[i] - mean);
+        var /= u.size();
+        double x_mean = p.t() * p.p();
+        double x_var = x_mean*(1-p.p());
+        assert(std::abs(mean - x_mean) / x_mean < 0.01);
+        assert(std::abs(var - x_var) / x_var < 0.01);
     }
 }

test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval.pass.cpp

@@ -72,7 +72,7 @@ int main()
         typedef std::minstd_rand G;
         G g;
         D d(20);
-        const int N = 10000;
+        const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
             u.push_back(d(g));

test/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.poisson/eval_param.pass.cpp

@@ -78,7 +78,7 @@ int main()
         G g;
         D d(2);
         P p(20);
-        const int N = 10000;
+        const int N = 100000;
         std::vector<double> u;
         for (int i = 0; i < N; ++i)
             u.push_back(d(g, p));