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

llvm-svn: 103886

1 files changed, 386 insertions, 270 deletions

diff --git a/libcxx/include/random b/libcxx/include/random
index 4045574b73dc..63a44aae443b 100644
--- a/libcxx/include/random
+++ b/libcxx/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)
-            : __t_(__t), __p_(__p) {}
+        explicit param_type(result_type __t = 1, double __p = 0.5);
 
         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:
@@ -3192,16 +3196,55 @@ public:
 };
 
 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());
-    _IntType __r = 0;
-    _IntType __t = __p.t();
-    for (_IntType __i = 0; __i < __t; ++__i)
-        __r += __bd(__g);
-    return __r;
+    if (__pr.__t_ == 0 || __pr.__p_ == 0)
+        return 0;
+    if (__pr.__p_ == 1)
+        return __pr.__t_;
+    uniform_real_distribution<double> __gen;
+    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,34 +3277,29 @@ operator>>(basic_istream<_CharT, _Traits>& __is,
     return __is;
 }
 
-// poisson_distribution
+// exponential_distribution
 
-template<class _IntType = int>
-class poisson_distribution
+template<class _RealType = double>
+class exponential_distribution
 {
 public:
     // types
-    typedef _IntType result_type;
+    typedef _RealType result_type;
 
     class param_type
     {
-        double __mean_;
-        double __sq_;
-        double __alxm_;
-        double __g_;
+        result_type __lambda_;
     public:
-        typedef poisson_distribution distribution_type;
+        typedef exponential_distribution distribution_type;
 
-        explicit param_type(double __mean = 1.0);
+        explicit param_type(result_type __lambda = 1) : __lambda_(__lambda) {}
 
-        double mean() const {return __mean_;}
+        result_type lambda() const {return __lambda_;}
 
         friend bool operator==(const param_type& __x, const param_type& __y)
-            {return __x.__mean_ == __y.__mean_;}
+            {return __x.__lambda_ == __y.__lambda_;}
         friend bool operator!=(const param_type& __x, const param_type& __y)
             {return !(__x == __y);}
-
-        friend class poisson_distribution;
     };
 
 private:
@@ -3269,8 +3307,9 @@ private:
 
 public:
     // constructors and reset functions
-    explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {}
-    explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
+    explicit exponential_distribution(result_type __lambda = 1)
+        : __p_(param_type(__lambda)) {}
+    explicit exponential_distribution(const param_type& __p) : __p_(__p) {}
     void reset() {}
 
     // generating functions
@@ -3279,127 +3318,249 @@ public:
     template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
 
     // property functions
-    double mean() const {return __p_.mean();}
+    result_type lambda() const {return __p_.lambda();}
 
     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();}
+    result_type max() const
+        {return -std::log(1-std::nextafter(result_type(1), result_type(-1))) /
+                 __p_.lambda();}
 
-    friend bool operator==(const poisson_distribution& __x,
-                           const poisson_distribution& __y)
+    friend bool operator==(const exponential_distribution& __x,
+                           const exponential_distribution& __y)
         {return __x.__p_ == __y.__p_;}
-    friend bool operator!=(const poisson_distribution& __x,
-                           const poisson_distribution& __y)
+    friend bool operator!=(const exponential_distribution& __x,
+                           const exponential_distribution& __y)
         {return !(__x == __y);}
 };
 
-template<class _IntType>
-poisson_distribution<_IntType>::param_type::param_type(double __mean)
-    : __mean_(__mean)
+template <class _RealType>
+template<class _URNG>
+_RealType
+exponential_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
 {
-    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);
-    }
+    return -_STD::log
+                  (
+                      result_type(1) -
+                      _STD::generate_canonical<result_type,
+                                       numeric_limits<result_type>::digits>(__g)
+                  )
+                  / __p.lambda();
 }
 
-template <class _IntType>
+template <class _CharT, class _Traits, class _RealType>
+basic_ostream<_CharT, _Traits>&
+operator<<(basic_ostream<_CharT, _Traits>& __os,
+           const exponential_distribution<_RealType>& __x)
+{
+    __save_flags<_CharT, _Traits> _(__os);
+    __os.flags(ios_base::dec | ios_base::left);
+    return __os << __x.lambda();
+}
+
+template <class _CharT, class _Traits, class _RealType>
+basic_istream<_CharT, _Traits>&
+operator>>(basic_istream<_CharT, _Traits>& __is,
+           exponential_distribution<_RealType>& __x)
+{
+    typedef exponential_distribution<_RealType> _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 __lambda;
+    __is >> __lambda;
+    if (!__is.fail())
+        __x.param(param_type(__lambda));
+    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>
-_IntType
-poisson_distribution<_IntType>::operator()(_URNG& __g, const param_type& __p)
+_RealType
+normal_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
 {
-    result_type __x;
-    uniform_real_distribution<double> __gen;
-    if (__p.__mean_ < 12.0)
+    result_type _U;
+    if (_V_hot_)
     {
-        __x = result_type(~0);
-        double __t = 1;
-        do
-        {
-            ++__x;
-            __t *= __gen(__g);
-        } while (__t > __p.__g_);
+        _V_hot_ = false;
+        _U = _V_;
     }
     else
     {
-        double __t;
-        const double __pi = 3.14159265358979323846264338328;
+        uniform_real_distribution<result_type> _Uni(-1, 1);
+        result_type __u;
+        result_type __v;
+        result_type __s;
         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);
+            __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 __x;
+    return _U * __p.stddev() + __p.mean();
 }
 
-template <class _CharT, class _Traits, class _IntType>
+template <class _CharT, class _Traits, class _RT>
 basic_ostream<_CharT, _Traits>&
 operator<<(basic_ostream<_CharT, _Traits>& __os,
-           const poisson_distribution<_IntType>& __x)
+           const normal_distribution<_RT>& __x)
 {
     __save_flags<_CharT, _Traits> _(__os);
     __os.flags(ios_base::dec | ios_base::left);
-    return __os << __x.mean();
+    _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 _IntType>
+template <class _CharT, class _Traits, class _RT>
 basic_istream<_CharT, _Traits>&
 operator>>(basic_istream<_CharT, _Traits>& __is,
-           poisson_distribution<_IntType>& __x)
+           normal_distribution<_RT>& __x)
 {
-    typedef poisson_distribution<_IntType> _Eng;
+    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);
-    double __mean;
-    __is >> __mean;
+    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));
+    {
+        __x.param(param_type(__mean, __stddev));
+        __x._V_hot_ = _V_hot;
+        __x._V_ = _V;
+    }
     return __is;
 }
 
-// exponential_distribution
+// poisson_distribution
 
-template<class _RealType = double>
-class exponential_distribution
+template<class _IntType = int>
+class poisson_distribution
 {
 public:
     // types
-    typedef _RealType result_type;
+    typedef _IntType result_type;
 
     class param_type
     {
-        result_type __lambda_;
+        double __mean_;
+        double __s_;
+        double __d_;
+        double __l_;
+        double __omega_;
+        double __c0_;
+        double __c1_;
+        double __c2_;
+        double __c3_;
+        double __c_;
+
     public:
-        typedef exponential_distribution distribution_type;
+        typedef poisson_distribution distribution_type;
 
-        explicit param_type(result_type __lambda = 1) : __lambda_(__lambda) {}
+        explicit param_type(double __mean = 1.0);
 
-        result_type lambda() const {return __lambda_;}
+        double mean() const {return __mean_;}
 
         friend bool operator==(const param_type& __x, const param_type& __y)
-            {return __x.__lambda_ == __y.__lambda_;}
+            {return __x.__mean_ == __y.__mean_;}
         friend bool operator!=(const param_type& __x, const param_type& __y)
             {return !(__x == __y);}
+
+        friend class poisson_distribution;
     };
 
 private:
@@ -3407,9 +3568,8 @@ private:
 
 public:
     // constructors and reset functions
-    explicit exponential_distribution(result_type __lambda = 1)
-        : __p_(param_type(__lambda)) {}
-    explicit exponential_distribution(const param_type& __p) : __p_(__p) {}
+    explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {}
+    explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
     void reset() {}
 
     // generating functions
@@ -3418,62 +3578,166 @@ public:
     template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
 
     // property functions
-    result_type lambda() const {return __p_.lambda();}
+    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 -std::log(1-std::nextafter(result_type(1), result_type(-1))) /
-                 __p_.lambda();}
+    result_type max() const {return numeric_limits<result_type>::max();}
 
-    friend bool operator==(const exponential_distribution& __x,
-                           const exponential_distribution& __y)
+    friend bool operator==(const poisson_distribution& __x,
+                           const poisson_distribution& __y)
         {return __x.__p_ == __y.__p_;}
-    friend bool operator!=(const exponential_distribution& __x,
-                           const exponential_distribution& __y)
+    friend bool operator!=(const poisson_distribution& __x,
+                           const poisson_distribution& __y)
         {return !(__x == __y);}
 };
 
-template <class _RealType>
+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>
-_RealType
-exponential_distribution<_RealType>::operator()(_URNG& __g, const param_type& __p)
+_IntType
+poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr)
 {
-    return -_STD::log
-                  (
-                      result_type(1) -
-                      _STD::generate_canonical<result_type,
-                                       numeric_limits<result_type>::digits>(__g)
-                  )
-                  / __p.lambda();
+    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 _RealType>
+template <class _CharT, class _Traits, class _IntType>
 basic_ostream<_CharT, _Traits>&
 operator<<(basic_ostream<_CharT, _Traits>& __os,
-           const exponential_distribution<_RealType>& __x)
+           const poisson_distribution<_IntType>& __x)
 {
     __save_flags<_CharT, _Traits> _(__os);
     __os.flags(ios_base::dec | ios_base::left);
-    return __os << __x.lambda();
+    return __os << __x.mean();
 }
 
-template <class _CharT, class _Traits, class _RealType>
+template <class _CharT, class _Traits, class _IntType>
 basic_istream<_CharT, _Traits>&
 operator>>(basic_istream<_CharT, _Traits>& __is,
-           exponential_distribution<_RealType>& __x)
+           poisson_distribution<_IntType>& __x)
 {
-    typedef exponential_distribution<_RealType> _Eng;
-    typedef typename _Eng::result_type result_type;
+    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);
-    result_type __lambda;
-    __is >> __lambda;
+    double __mean;
+    __is >> __mean;
     if (!__is.fail())
-        __x.param(param_type(__lambda));
+        __x.param(param_type(__mean));
     return __is;
 }
 
@@ -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