1e4b17023SJohn Marino // random number generation (out of line) -*- C++ -*- 2e4b17023SJohn Marino 35ce9237cSJohn Marino // Copyright (C) 2009-2012 Free Software Foundation, Inc. 4e4b17023SJohn Marino // 5e4b17023SJohn Marino // This file is part of the GNU ISO C++ Library. This library is free 6e4b17023SJohn Marino // software; you can redistribute it and/or modify it under the 7e4b17023SJohn Marino // terms of the GNU General Public License as published by the 8e4b17023SJohn Marino // Free Software Foundation; either version 3, or (at your option) 9e4b17023SJohn Marino // any later version. 10e4b17023SJohn Marino 11e4b17023SJohn Marino // This library is distributed in the hope that it will be useful, 12e4b17023SJohn Marino // but WITHOUT ANY WARRANTY; without even the implied warranty of 13e4b17023SJohn Marino // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14e4b17023SJohn Marino // GNU General Public License for more details. 15e4b17023SJohn Marino 16e4b17023SJohn Marino // Under Section 7 of GPL version 3, you are granted additional 17e4b17023SJohn Marino // permissions described in the GCC Runtime Library Exception, version 18e4b17023SJohn Marino // 3.1, as published by the Free Software Foundation. 19e4b17023SJohn Marino 20e4b17023SJohn Marino // You should have received a copy of the GNU General Public License and 21e4b17023SJohn Marino // a copy of the GCC Runtime Library Exception along with this program; 22e4b17023SJohn Marino // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 23e4b17023SJohn Marino // <http://www.gnu.org/licenses/>. 24e4b17023SJohn Marino 25e4b17023SJohn Marino /** @file bits/random.tcc 26e4b17023SJohn Marino * This is an internal header file, included by other library headers. 27e4b17023SJohn Marino * Do not attempt to use it directly. @headername{random} 28e4b17023SJohn Marino */ 29e4b17023SJohn Marino 30e4b17023SJohn Marino #ifndef _RANDOM_TCC 31e4b17023SJohn Marino #define _RANDOM_TCC 1 32e4b17023SJohn Marino 33e4b17023SJohn Marino #include <numeric> // std::accumulate and std::partial_sum 34e4b17023SJohn Marino 35e4b17023SJohn Marino namespace std _GLIBCXX_VISIBILITY(default) 36e4b17023SJohn Marino { 37e4b17023SJohn Marino /* 38e4b17023SJohn Marino * (Further) implementation-space details. 39e4b17023SJohn Marino */ 40e4b17023SJohn Marino namespace __detail 41e4b17023SJohn Marino { 42e4b17023SJohn Marino _GLIBCXX_BEGIN_NAMESPACE_VERSION 43e4b17023SJohn Marino 44e4b17023SJohn Marino // General case for x = (ax + c) mod m -- use Schrage's algorithm to 45e4b17023SJohn Marino // avoid integer overflow. 46e4b17023SJohn Marino // 47e4b17023SJohn Marino // Because a and c are compile-time integral constants the compiler 48e4b17023SJohn Marino // kindly elides any unreachable paths. 49e4b17023SJohn Marino // 50e4b17023SJohn Marino // Preconditions: a > 0, m > 0. 51e4b17023SJohn Marino // 52e4b17023SJohn Marino // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a. 53e4b17023SJohn Marino // 54e4b17023SJohn Marino template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool> 55e4b17023SJohn Marino struct _Mod 56e4b17023SJohn Marino { 57e4b17023SJohn Marino static _Tp __calcstd::__detail::_Mod58e4b17023SJohn Marino __calc(_Tp __x) 59e4b17023SJohn Marino { 60e4b17023SJohn Marino if (__a == 1) 61e4b17023SJohn Marino __x %= __m; 62e4b17023SJohn Marino else 63e4b17023SJohn Marino { 64e4b17023SJohn Marino static const _Tp __q = __m / __a; 65e4b17023SJohn Marino static const _Tp __r = __m % __a; 66e4b17023SJohn Marino 67e4b17023SJohn Marino _Tp __t1 = __a * (__x % __q); 68e4b17023SJohn Marino _Tp __t2 = __r * (__x / __q); 69e4b17023SJohn Marino if (__t1 >= __t2) 70e4b17023SJohn Marino __x = __t1 - __t2; 71e4b17023SJohn Marino else 72e4b17023SJohn Marino __x = __m - __t2 + __t1; 73e4b17023SJohn Marino } 74e4b17023SJohn Marino 75e4b17023SJohn Marino if (__c != 0) 76e4b17023SJohn Marino { 77e4b17023SJohn Marino const _Tp __d = __m - __x; 78e4b17023SJohn Marino if (__d > __c) 79e4b17023SJohn Marino __x += __c; 80e4b17023SJohn Marino else 81e4b17023SJohn Marino __x = __c - __d; 82e4b17023SJohn Marino } 83e4b17023SJohn Marino return __x; 84e4b17023SJohn Marino } 85e4b17023SJohn Marino }; 86e4b17023SJohn Marino 87e4b17023SJohn Marino // Special case for m == 0 -- use unsigned integer overflow as modulo 88e4b17023SJohn Marino // operator. 89e4b17023SJohn Marino template<typename _Tp, _Tp __m, _Tp __a, _Tp __c> 90e4b17023SJohn Marino struct _Mod<_Tp, __m, __a, __c, true> 91e4b17023SJohn Marino { 92e4b17023SJohn Marino static _Tp __calcstd::__detail::_Mod93e4b17023SJohn Marino __calc(_Tp __x) 94e4b17023SJohn Marino { return __a * __x + __c; } 95e4b17023SJohn Marino }; 96e4b17023SJohn Marino 97e4b17023SJohn Marino template<typename _InputIterator, typename _OutputIterator, 98e4b17023SJohn Marino typename _UnaryOperation> 99e4b17023SJohn Marino _OutputIterator 100e4b17023SJohn Marino __transform(_InputIterator __first, _InputIterator __last, 101e4b17023SJohn Marino _OutputIterator __result, _UnaryOperation __unary_op) 102e4b17023SJohn Marino { 103e4b17023SJohn Marino for (; __first != __last; ++__first, ++__result) 104e4b17023SJohn Marino *__result = __unary_op(*__first); 105e4b17023SJohn Marino return __result; 106e4b17023SJohn Marino } 107e4b17023SJohn Marino 108e4b17023SJohn Marino _GLIBCXX_END_NAMESPACE_VERSION 109e4b17023SJohn Marino } // namespace __detail 110e4b17023SJohn Marino 111e4b17023SJohn Marino _GLIBCXX_BEGIN_NAMESPACE_VERSION 112e4b17023SJohn Marino 113e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 114e4b17023SJohn Marino constexpr _UIntType 115e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier; 116e4b17023SJohn Marino 117e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 118e4b17023SJohn Marino constexpr _UIntType 119e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>::increment; 120e4b17023SJohn Marino 121e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 122e4b17023SJohn Marino constexpr _UIntType 123e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>::modulus; 124e4b17023SJohn Marino 125e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 126e4b17023SJohn Marino constexpr _UIntType 127e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed; 128e4b17023SJohn Marino 129e4b17023SJohn Marino /** 130e4b17023SJohn Marino * Seeds the LCR with integral value @p __s, adjusted so that the 131e4b17023SJohn Marino * ring identity is never a member of the convergence set. 132e4b17023SJohn Marino */ 133e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 134e4b17023SJohn Marino void 135e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>:: seed(result_type __s)136e4b17023SJohn Marino seed(result_type __s) 137e4b17023SJohn Marino { 138e4b17023SJohn Marino if ((__detail::__mod<_UIntType, __m>(__c) == 0) 139e4b17023SJohn Marino && (__detail::__mod<_UIntType, __m>(__s) == 0)) 140e4b17023SJohn Marino _M_x = 1; 141e4b17023SJohn Marino else 142e4b17023SJohn Marino _M_x = __detail::__mod<_UIntType, __m>(__s); 143e4b17023SJohn Marino } 144e4b17023SJohn Marino 145e4b17023SJohn Marino /** 146e4b17023SJohn Marino * Seeds the LCR engine with a value generated by @p __q. 147e4b17023SJohn Marino */ 148e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m> 149e4b17023SJohn Marino template<typename _Sseq> 150e4b17023SJohn Marino typename std::enable_if<std::is_class<_Sseq>::value>::type 151e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>:: seed(_Sseq & __q)152e4b17023SJohn Marino seed(_Sseq& __q) 153e4b17023SJohn Marino { 154e4b17023SJohn Marino const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits 155e4b17023SJohn Marino : std::__lg(__m); 156e4b17023SJohn Marino const _UIntType __k = (__k0 + 31) / 32; 157e4b17023SJohn Marino uint_least32_t __arr[__k + 3]; 158e4b17023SJohn Marino __q.generate(__arr + 0, __arr + __k + 3); 159e4b17023SJohn Marino _UIntType __factor = 1u; 160e4b17023SJohn Marino _UIntType __sum = 0u; 161e4b17023SJohn Marino for (size_t __j = 0; __j < __k; ++__j) 162e4b17023SJohn Marino { 163e4b17023SJohn Marino __sum += __arr[__j + 3] * __factor; 164e4b17023SJohn Marino __factor *= __detail::_Shift<_UIntType, 32>::__value; 165e4b17023SJohn Marino } 166e4b17023SJohn Marino seed(__sum); 167e4b17023SJohn Marino } 168e4b17023SJohn Marino 169e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 170e4b17023SJohn Marino typename _CharT, typename _Traits> 171e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const linear_congruential_engine<_UIntType,__a,__c,__m> & __lcr)172e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 173e4b17023SJohn Marino const linear_congruential_engine<_UIntType, 174e4b17023SJohn Marino __a, __c, __m>& __lcr) 175e4b17023SJohn Marino { 176e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 177e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 178e4b17023SJohn Marino 179e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 180e4b17023SJohn Marino const _CharT __fill = __os.fill(); 181e4b17023SJohn Marino __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 182e4b17023SJohn Marino __os.fill(__os.widen(' ')); 183e4b17023SJohn Marino 184e4b17023SJohn Marino __os << __lcr._M_x; 185e4b17023SJohn Marino 186e4b17023SJohn Marino __os.flags(__flags); 187e4b17023SJohn Marino __os.fill(__fill); 188e4b17023SJohn Marino return __os; 189e4b17023SJohn Marino } 190e4b17023SJohn Marino 191e4b17023SJohn Marino template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m, 192e4b17023SJohn Marino typename _CharT, typename _Traits> 193e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,linear_congruential_engine<_UIntType,__a,__c,__m> & __lcr)194e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 195e4b17023SJohn Marino linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr) 196e4b17023SJohn Marino { 197e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 198e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 199e4b17023SJohn Marino 200e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 201e4b17023SJohn Marino __is.flags(__ios_base::dec); 202e4b17023SJohn Marino 203e4b17023SJohn Marino __is >> __lcr._M_x; 204e4b17023SJohn Marino 205e4b17023SJohn Marino __is.flags(__flags); 206e4b17023SJohn Marino return __is; 207e4b17023SJohn Marino } 208e4b17023SJohn Marino 209e4b17023SJohn Marino 210e4b17023SJohn Marino template<typename _UIntType, 211e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 212e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 213e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 214e4b17023SJohn Marino _UIntType __f> 215e4b17023SJohn Marino constexpr size_t 216e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 217e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::word_size; 218e4b17023SJohn Marino 219e4b17023SJohn Marino template<typename _UIntType, 220e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 221e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 222e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 223e4b17023SJohn Marino _UIntType __f> 224e4b17023SJohn Marino constexpr size_t 225e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 226e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::state_size; 227e4b17023SJohn Marino 228e4b17023SJohn Marino template<typename _UIntType, 229e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 230e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 231e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 232e4b17023SJohn Marino _UIntType __f> 233e4b17023SJohn Marino constexpr size_t 234e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 235e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::shift_size; 236e4b17023SJohn Marino 237e4b17023SJohn Marino template<typename _UIntType, 238e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 239e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 240e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 241e4b17023SJohn Marino _UIntType __f> 242e4b17023SJohn Marino constexpr size_t 243e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 244e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::mask_bits; 245e4b17023SJohn Marino 246e4b17023SJohn Marino template<typename _UIntType, 247e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 248e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 249e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 250e4b17023SJohn Marino _UIntType __f> 251e4b17023SJohn Marino constexpr _UIntType 252e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 253e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::xor_mask; 254e4b17023SJohn Marino 255e4b17023SJohn Marino template<typename _UIntType, 256e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 257e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 258e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 259e4b17023SJohn Marino _UIntType __f> 260e4b17023SJohn Marino constexpr size_t 261e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 262e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_u; 263e4b17023SJohn Marino 264e4b17023SJohn Marino template<typename _UIntType, 265e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 266e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 267e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 268e4b17023SJohn Marino _UIntType __f> 269e4b17023SJohn Marino constexpr _UIntType 270e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 271e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_d; 272e4b17023SJohn Marino 273e4b17023SJohn Marino template<typename _UIntType, 274e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 275e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 276e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 277e4b17023SJohn Marino _UIntType __f> 278e4b17023SJohn Marino constexpr size_t 279e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 280e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_s; 281e4b17023SJohn Marino 282e4b17023SJohn Marino template<typename _UIntType, 283e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 284e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 285e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 286e4b17023SJohn Marino _UIntType __f> 287e4b17023SJohn Marino constexpr _UIntType 288e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 289e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_b; 290e4b17023SJohn Marino 291e4b17023SJohn Marino template<typename _UIntType, 292e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 293e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 294e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 295e4b17023SJohn Marino _UIntType __f> 296e4b17023SJohn Marino constexpr size_t 297e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 298e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_t; 299e4b17023SJohn Marino 300e4b17023SJohn Marino template<typename _UIntType, 301e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 302e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 303e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 304e4b17023SJohn Marino _UIntType __f> 305e4b17023SJohn Marino constexpr _UIntType 306e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 307e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_c; 308e4b17023SJohn Marino 309e4b17023SJohn Marino template<typename _UIntType, 310e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 311e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 312e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 313e4b17023SJohn Marino _UIntType __f> 314e4b17023SJohn Marino constexpr size_t 315e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 316e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::tempering_l; 317e4b17023SJohn Marino 318e4b17023SJohn Marino template<typename _UIntType, 319e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 320e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 321e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 322e4b17023SJohn Marino _UIntType __f> 323e4b17023SJohn Marino constexpr _UIntType 324e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 325e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>:: 326e4b17023SJohn Marino initialization_multiplier; 327e4b17023SJohn Marino 328e4b17023SJohn Marino template<typename _UIntType, 329e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 330e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 331e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 332e4b17023SJohn Marino _UIntType __f> 333e4b17023SJohn Marino constexpr _UIntType 334e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 335e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::default_seed; 336e4b17023SJohn Marino 337e4b17023SJohn Marino template<typename _UIntType, 338e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 339e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 340e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 341e4b17023SJohn Marino _UIntType __f> 342e4b17023SJohn Marino void 343e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 344e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>:: seed(result_type __sd)345e4b17023SJohn Marino seed(result_type __sd) 346e4b17023SJohn Marino { 347e4b17023SJohn Marino _M_x[0] = __detail::__mod<_UIntType, 348e4b17023SJohn Marino __detail::_Shift<_UIntType, __w>::__value>(__sd); 349e4b17023SJohn Marino 350e4b17023SJohn Marino for (size_t __i = 1; __i < state_size; ++__i) 351e4b17023SJohn Marino { 352e4b17023SJohn Marino _UIntType __x = _M_x[__i - 1]; 353e4b17023SJohn Marino __x ^= __x >> (__w - 2); 354e4b17023SJohn Marino __x *= __f; 355e4b17023SJohn Marino __x += __detail::__mod<_UIntType, __n>(__i); 356e4b17023SJohn Marino _M_x[__i] = __detail::__mod<_UIntType, 357e4b17023SJohn Marino __detail::_Shift<_UIntType, __w>::__value>(__x); 358e4b17023SJohn Marino } 359e4b17023SJohn Marino _M_p = state_size; 360e4b17023SJohn Marino } 361e4b17023SJohn Marino 362e4b17023SJohn Marino template<typename _UIntType, 363e4b17023SJohn Marino size_t __w, size_t __n, size_t __m, size_t __r, 364e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 365e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 366e4b17023SJohn Marino _UIntType __f> 367e4b17023SJohn Marino template<typename _Sseq> 368e4b17023SJohn Marino typename std::enable_if<std::is_class<_Sseq>::value>::type 369e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 370e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>:: seed(_Sseq & __q)371e4b17023SJohn Marino seed(_Sseq& __q) 372e4b17023SJohn Marino { 373e4b17023SJohn Marino const _UIntType __upper_mask = (~_UIntType()) << __r; 374e4b17023SJohn Marino const size_t __k = (__w + 31) / 32; 375e4b17023SJohn Marino uint_least32_t __arr[__n * __k]; 376e4b17023SJohn Marino __q.generate(__arr + 0, __arr + __n * __k); 377e4b17023SJohn Marino 378e4b17023SJohn Marino bool __zero = true; 379e4b17023SJohn Marino for (size_t __i = 0; __i < state_size; ++__i) 380e4b17023SJohn Marino { 381e4b17023SJohn Marino _UIntType __factor = 1u; 382e4b17023SJohn Marino _UIntType __sum = 0u; 383e4b17023SJohn Marino for (size_t __j = 0; __j < __k; ++__j) 384e4b17023SJohn Marino { 385e4b17023SJohn Marino __sum += __arr[__k * __i + __j] * __factor; 386e4b17023SJohn Marino __factor *= __detail::_Shift<_UIntType, 32>::__value; 387e4b17023SJohn Marino } 388e4b17023SJohn Marino _M_x[__i] = __detail::__mod<_UIntType, 389e4b17023SJohn Marino __detail::_Shift<_UIntType, __w>::__value>(__sum); 390e4b17023SJohn Marino 391e4b17023SJohn Marino if (__zero) 392e4b17023SJohn Marino { 393e4b17023SJohn Marino if (__i == 0) 394e4b17023SJohn Marino { 395e4b17023SJohn Marino if ((_M_x[0] & __upper_mask) != 0u) 396e4b17023SJohn Marino __zero = false; 397e4b17023SJohn Marino } 398e4b17023SJohn Marino else if (_M_x[__i] != 0u) 399e4b17023SJohn Marino __zero = false; 400e4b17023SJohn Marino } 401e4b17023SJohn Marino } 402e4b17023SJohn Marino if (__zero) 403e4b17023SJohn Marino _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value; 4045ce9237cSJohn Marino _M_p = state_size; 405e4b17023SJohn Marino } 406e4b17023SJohn Marino 407e4b17023SJohn Marino template<typename _UIntType, size_t __w, 408e4b17023SJohn Marino size_t __n, size_t __m, size_t __r, 409e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 410e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 411e4b17023SJohn Marino _UIntType __f> 412e4b17023SJohn Marino typename 413e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 414e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>::result_type 415e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d, 416e4b17023SJohn Marino __s, __b, __t, __c, __l, __f>:: operator ()()417e4b17023SJohn Marino operator()() 418e4b17023SJohn Marino { 419e4b17023SJohn Marino // Reload the vector - cost is O(n) amortized over n calls. 420e4b17023SJohn Marino if (_M_p >= state_size) 421e4b17023SJohn Marino { 422e4b17023SJohn Marino const _UIntType __upper_mask = (~_UIntType()) << __r; 423e4b17023SJohn Marino const _UIntType __lower_mask = ~__upper_mask; 424e4b17023SJohn Marino 425e4b17023SJohn Marino for (size_t __k = 0; __k < (__n - __m); ++__k) 426e4b17023SJohn Marino { 427e4b17023SJohn Marino _UIntType __y = ((_M_x[__k] & __upper_mask) 428e4b17023SJohn Marino | (_M_x[__k + 1] & __lower_mask)); 429e4b17023SJohn Marino _M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1) 430e4b17023SJohn Marino ^ ((__y & 0x01) ? __a : 0)); 431e4b17023SJohn Marino } 432e4b17023SJohn Marino 433e4b17023SJohn Marino for (size_t __k = (__n - __m); __k < (__n - 1); ++__k) 434e4b17023SJohn Marino { 435e4b17023SJohn Marino _UIntType __y = ((_M_x[__k] & __upper_mask) 436e4b17023SJohn Marino | (_M_x[__k + 1] & __lower_mask)); 437e4b17023SJohn Marino _M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1) 438e4b17023SJohn Marino ^ ((__y & 0x01) ? __a : 0)); 439e4b17023SJohn Marino } 440e4b17023SJohn Marino 441e4b17023SJohn Marino _UIntType __y = ((_M_x[__n - 1] & __upper_mask) 442e4b17023SJohn Marino | (_M_x[0] & __lower_mask)); 443e4b17023SJohn Marino _M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1) 444e4b17023SJohn Marino ^ ((__y & 0x01) ? __a : 0)); 445e4b17023SJohn Marino _M_p = 0; 446e4b17023SJohn Marino } 447e4b17023SJohn Marino 448e4b17023SJohn Marino // Calculate o(x(i)). 449e4b17023SJohn Marino result_type __z = _M_x[_M_p++]; 450e4b17023SJohn Marino __z ^= (__z >> __u) & __d; 451e4b17023SJohn Marino __z ^= (__z << __s) & __b; 452e4b17023SJohn Marino __z ^= (__z << __t) & __c; 453e4b17023SJohn Marino __z ^= (__z >> __l); 454e4b17023SJohn Marino 455e4b17023SJohn Marino return __z; 456e4b17023SJohn Marino } 457e4b17023SJohn Marino 458e4b17023SJohn Marino template<typename _UIntType, size_t __w, 459e4b17023SJohn Marino size_t __n, size_t __m, size_t __r, 460e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 461e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 462e4b17023SJohn Marino _UIntType __f, typename _CharT, typename _Traits> 463e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const mersenne_twister_engine<_UIntType,__w,__n,__m,__r,__a,__u,__d,__s,__b,__t,__c,__l,__f> & __x)464e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 465e4b17023SJohn Marino const mersenne_twister_engine<_UIntType, __w, __n, __m, 466e4b17023SJohn Marino __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 467e4b17023SJohn Marino { 468e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 469e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 470e4b17023SJohn Marino 471e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 472e4b17023SJohn Marino const _CharT __fill = __os.fill(); 473e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 474e4b17023SJohn Marino __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 475e4b17023SJohn Marino __os.fill(__space); 476e4b17023SJohn Marino 477e4b17023SJohn Marino for (size_t __i = 0; __i < __n; ++__i) 478e4b17023SJohn Marino __os << __x._M_x[__i] << __space; 479e4b17023SJohn Marino __os << __x._M_p; 480e4b17023SJohn Marino 481e4b17023SJohn Marino __os.flags(__flags); 482e4b17023SJohn Marino __os.fill(__fill); 483e4b17023SJohn Marino return __os; 484e4b17023SJohn Marino } 485e4b17023SJohn Marino 486e4b17023SJohn Marino template<typename _UIntType, size_t __w, 487e4b17023SJohn Marino size_t __n, size_t __m, size_t __r, 488e4b17023SJohn Marino _UIntType __a, size_t __u, _UIntType __d, size_t __s, 489e4b17023SJohn Marino _UIntType __b, size_t __t, _UIntType __c, size_t __l, 490e4b17023SJohn Marino _UIntType __f, typename _CharT, typename _Traits> 491e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,mersenne_twister_engine<_UIntType,__w,__n,__m,__r,__a,__u,__d,__s,__b,__t,__c,__l,__f> & __x)492e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 493e4b17023SJohn Marino mersenne_twister_engine<_UIntType, __w, __n, __m, 494e4b17023SJohn Marino __r, __a, __u, __d, __s, __b, __t, __c, __l, __f>& __x) 495e4b17023SJohn Marino { 496e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 497e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 498e4b17023SJohn Marino 499e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 500e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 501e4b17023SJohn Marino 502e4b17023SJohn Marino for (size_t __i = 0; __i < __n; ++__i) 503e4b17023SJohn Marino __is >> __x._M_x[__i]; 504e4b17023SJohn Marino __is >> __x._M_p; 505e4b17023SJohn Marino 506e4b17023SJohn Marino __is.flags(__flags); 507e4b17023SJohn Marino return __is; 508e4b17023SJohn Marino } 509e4b17023SJohn Marino 510e4b17023SJohn Marino 511e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 512e4b17023SJohn Marino constexpr size_t 513e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size; 514e4b17023SJohn Marino 515e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 516e4b17023SJohn Marino constexpr size_t 517e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag; 518e4b17023SJohn Marino 519e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 520e4b17023SJohn Marino constexpr size_t 521e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag; 522e4b17023SJohn Marino 523e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 524e4b17023SJohn Marino constexpr _UIntType 525e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed; 526e4b17023SJohn Marino 527e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 528e4b17023SJohn Marino void 529e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>:: seed(result_type __value)530e4b17023SJohn Marino seed(result_type __value) 531e4b17023SJohn Marino { 532e4b17023SJohn Marino std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u> 533e4b17023SJohn Marino __lcg(__value == 0u ? default_seed : __value); 534e4b17023SJohn Marino 535e4b17023SJohn Marino const size_t __n = (__w + 31) / 32; 536e4b17023SJohn Marino 537e4b17023SJohn Marino for (size_t __i = 0; __i < long_lag; ++__i) 538e4b17023SJohn Marino { 539e4b17023SJohn Marino _UIntType __sum = 0u; 540e4b17023SJohn Marino _UIntType __factor = 1u; 541e4b17023SJohn Marino for (size_t __j = 0; __j < __n; ++__j) 542e4b17023SJohn Marino { 543e4b17023SJohn Marino __sum += __detail::__mod<uint_least32_t, 544e4b17023SJohn Marino __detail::_Shift<uint_least32_t, 32>::__value> 545e4b17023SJohn Marino (__lcg()) * __factor; 546e4b17023SJohn Marino __factor *= __detail::_Shift<_UIntType, 32>::__value; 547e4b17023SJohn Marino } 548e4b17023SJohn Marino _M_x[__i] = __detail::__mod<_UIntType, 549e4b17023SJohn Marino __detail::_Shift<_UIntType, __w>::__value>(__sum); 550e4b17023SJohn Marino } 551e4b17023SJohn Marino _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 552e4b17023SJohn Marino _M_p = 0; 553e4b17023SJohn Marino } 554e4b17023SJohn Marino 555e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 556e4b17023SJohn Marino template<typename _Sseq> 557e4b17023SJohn Marino typename std::enable_if<std::is_class<_Sseq>::value>::type 558e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>:: seed(_Sseq & __q)559e4b17023SJohn Marino seed(_Sseq& __q) 560e4b17023SJohn Marino { 561e4b17023SJohn Marino const size_t __k = (__w + 31) / 32; 562e4b17023SJohn Marino uint_least32_t __arr[__r * __k]; 563e4b17023SJohn Marino __q.generate(__arr + 0, __arr + __r * __k); 564e4b17023SJohn Marino 565e4b17023SJohn Marino for (size_t __i = 0; __i < long_lag; ++__i) 566e4b17023SJohn Marino { 567e4b17023SJohn Marino _UIntType __sum = 0u; 568e4b17023SJohn Marino _UIntType __factor = 1u; 569e4b17023SJohn Marino for (size_t __j = 0; __j < __k; ++__j) 570e4b17023SJohn Marino { 571e4b17023SJohn Marino __sum += __arr[__k * __i + __j] * __factor; 572e4b17023SJohn Marino __factor *= __detail::_Shift<_UIntType, 32>::__value; 573e4b17023SJohn Marino } 574e4b17023SJohn Marino _M_x[__i] = __detail::__mod<_UIntType, 575e4b17023SJohn Marino __detail::_Shift<_UIntType, __w>::__value>(__sum); 576e4b17023SJohn Marino } 577e4b17023SJohn Marino _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; 578e4b17023SJohn Marino _M_p = 0; 579e4b17023SJohn Marino } 580e4b17023SJohn Marino 581e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r> 582e4b17023SJohn Marino typename subtract_with_carry_engine<_UIntType, __w, __s, __r>:: 583e4b17023SJohn Marino result_type 584e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>:: operator ()()585e4b17023SJohn Marino operator()() 586e4b17023SJohn Marino { 587e4b17023SJohn Marino // Derive short lag index from current index. 588e4b17023SJohn Marino long __ps = _M_p - short_lag; 589e4b17023SJohn Marino if (__ps < 0) 590e4b17023SJohn Marino __ps += long_lag; 591e4b17023SJohn Marino 592e4b17023SJohn Marino // Calculate new x(i) without overflow or division. 593e4b17023SJohn Marino // NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry 594e4b17023SJohn Marino // cannot overflow. 595e4b17023SJohn Marino _UIntType __xi; 596e4b17023SJohn Marino if (_M_x[__ps] >= _M_x[_M_p] + _M_carry) 597e4b17023SJohn Marino { 598e4b17023SJohn Marino __xi = _M_x[__ps] - _M_x[_M_p] - _M_carry; 599e4b17023SJohn Marino _M_carry = 0; 600e4b17023SJohn Marino } 601e4b17023SJohn Marino else 602e4b17023SJohn Marino { 603e4b17023SJohn Marino __xi = (__detail::_Shift<_UIntType, __w>::__value 604e4b17023SJohn Marino - _M_x[_M_p] - _M_carry + _M_x[__ps]); 605e4b17023SJohn Marino _M_carry = 1; 606e4b17023SJohn Marino } 607e4b17023SJohn Marino _M_x[_M_p] = __xi; 608e4b17023SJohn Marino 609e4b17023SJohn Marino // Adjust current index to loop around in ring buffer. 610e4b17023SJohn Marino if (++_M_p >= long_lag) 611e4b17023SJohn Marino _M_p = 0; 612e4b17023SJohn Marino 613e4b17023SJohn Marino return __xi; 614e4b17023SJohn Marino } 615e4b17023SJohn Marino 616e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r, 617e4b17023SJohn Marino typename _CharT, typename _Traits> 618e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const subtract_with_carry_engine<_UIntType,__w,__s,__r> & __x)619e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 620e4b17023SJohn Marino const subtract_with_carry_engine<_UIntType, 621e4b17023SJohn Marino __w, __s, __r>& __x) 622e4b17023SJohn Marino { 623e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 624e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 625e4b17023SJohn Marino 626e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 627e4b17023SJohn Marino const _CharT __fill = __os.fill(); 628e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 629e4b17023SJohn Marino __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 630e4b17023SJohn Marino __os.fill(__space); 631e4b17023SJohn Marino 632e4b17023SJohn Marino for (size_t __i = 0; __i < __r; ++__i) 633e4b17023SJohn Marino __os << __x._M_x[__i] << __space; 634e4b17023SJohn Marino __os << __x._M_carry << __space << __x._M_p; 635e4b17023SJohn Marino 636e4b17023SJohn Marino __os.flags(__flags); 637e4b17023SJohn Marino __os.fill(__fill); 638e4b17023SJohn Marino return __os; 639e4b17023SJohn Marino } 640e4b17023SJohn Marino 641e4b17023SJohn Marino template<typename _UIntType, size_t __w, size_t __s, size_t __r, 642e4b17023SJohn Marino typename _CharT, typename _Traits> 643e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,subtract_with_carry_engine<_UIntType,__w,__s,__r> & __x)644e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 645e4b17023SJohn Marino subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x) 646e4b17023SJohn Marino { 647e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __istream_type; 648e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 649e4b17023SJohn Marino 650e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 651e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 652e4b17023SJohn Marino 653e4b17023SJohn Marino for (size_t __i = 0; __i < __r; ++__i) 654e4b17023SJohn Marino __is >> __x._M_x[__i]; 655e4b17023SJohn Marino __is >> __x._M_carry; 656e4b17023SJohn Marino __is >> __x._M_p; 657e4b17023SJohn Marino 658e4b17023SJohn Marino __is.flags(__flags); 659e4b17023SJohn Marino return __is; 660e4b17023SJohn Marino } 661e4b17023SJohn Marino 662e4b17023SJohn Marino 663e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __p, size_t __r> 664e4b17023SJohn Marino constexpr size_t 665e4b17023SJohn Marino discard_block_engine<_RandomNumberEngine, __p, __r>::block_size; 666e4b17023SJohn Marino 667e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __p, size_t __r> 668e4b17023SJohn Marino constexpr size_t 669e4b17023SJohn Marino discard_block_engine<_RandomNumberEngine, __p, __r>::used_block; 670e4b17023SJohn Marino 671e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __p, size_t __r> 672e4b17023SJohn Marino typename discard_block_engine<_RandomNumberEngine, 673e4b17023SJohn Marino __p, __r>::result_type 674e4b17023SJohn Marino discard_block_engine<_RandomNumberEngine, __p, __r>:: operator ()()675e4b17023SJohn Marino operator()() 676e4b17023SJohn Marino { 677e4b17023SJohn Marino if (_M_n >= used_block) 678e4b17023SJohn Marino { 679e4b17023SJohn Marino _M_b.discard(block_size - _M_n); 680e4b17023SJohn Marino _M_n = 0; 681e4b17023SJohn Marino } 682e4b17023SJohn Marino ++_M_n; 683e4b17023SJohn Marino return _M_b(); 684e4b17023SJohn Marino } 685e4b17023SJohn Marino 686e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __p, size_t __r, 687e4b17023SJohn Marino typename _CharT, typename _Traits> 688e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const discard_block_engine<_RandomNumberEngine,__p,__r> & __x)689e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 690e4b17023SJohn Marino const discard_block_engine<_RandomNumberEngine, 691e4b17023SJohn Marino __p, __r>& __x) 692e4b17023SJohn Marino { 693e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 694e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 695e4b17023SJohn Marino 696e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 697e4b17023SJohn Marino const _CharT __fill = __os.fill(); 698e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 699e4b17023SJohn Marino __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 700e4b17023SJohn Marino __os.fill(__space); 701e4b17023SJohn Marino 702e4b17023SJohn Marino __os << __x.base() << __space << __x._M_n; 703e4b17023SJohn Marino 704e4b17023SJohn Marino __os.flags(__flags); 705e4b17023SJohn Marino __os.fill(__fill); 706e4b17023SJohn Marino return __os; 707e4b17023SJohn Marino } 708e4b17023SJohn Marino 709e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __p, size_t __r, 710e4b17023SJohn Marino typename _CharT, typename _Traits> 711e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,discard_block_engine<_RandomNumberEngine,__p,__r> & __x)712e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 713e4b17023SJohn Marino discard_block_engine<_RandomNumberEngine, __p, __r>& __x) 714e4b17023SJohn Marino { 715e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 716e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 717e4b17023SJohn Marino 718e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 719e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 720e4b17023SJohn Marino 721e4b17023SJohn Marino __is >> __x._M_b >> __x._M_n; 722e4b17023SJohn Marino 723e4b17023SJohn Marino __is.flags(__flags); 724e4b17023SJohn Marino return __is; 725e4b17023SJohn Marino } 726e4b17023SJohn Marino 727e4b17023SJohn Marino 728e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __w, typename _UIntType> 729e4b17023SJohn Marino typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: 730e4b17023SJohn Marino result_type 731e4b17023SJohn Marino independent_bits_engine<_RandomNumberEngine, __w, _UIntType>:: operator ()()732e4b17023SJohn Marino operator()() 733e4b17023SJohn Marino { 734e4b17023SJohn Marino typedef typename _RandomNumberEngine::result_type _Eresult_type; 735e4b17023SJohn Marino const _Eresult_type __r 736e4b17023SJohn Marino = (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max() 737e4b17023SJohn Marino ? _M_b.max() - _M_b.min() + 1 : 0); 738e4b17023SJohn Marino const unsigned __edig = std::numeric_limits<_Eresult_type>::digits; 739e4b17023SJohn Marino const unsigned __m = __r ? std::__lg(__r) : __edig; 740e4b17023SJohn Marino 741e4b17023SJohn Marino typedef typename std::common_type<_Eresult_type, result_type>::type 742e4b17023SJohn Marino __ctype; 743e4b17023SJohn Marino const unsigned __cdig = std::numeric_limits<__ctype>::digits; 744e4b17023SJohn Marino 745e4b17023SJohn Marino unsigned __n, __n0; 746e4b17023SJohn Marino __ctype __s0, __s1, __y0, __y1; 747e4b17023SJohn Marino 748e4b17023SJohn Marino for (size_t __i = 0; __i < 2; ++__i) 749e4b17023SJohn Marino { 750e4b17023SJohn Marino __n = (__w + __m - 1) / __m + __i; 751e4b17023SJohn Marino __n0 = __n - __w % __n; 752e4b17023SJohn Marino const unsigned __w0 = __w / __n; // __w0 <= __m 753e4b17023SJohn Marino 754e4b17023SJohn Marino __s0 = 0; 755e4b17023SJohn Marino __s1 = 0; 756e4b17023SJohn Marino if (__w0 < __cdig) 757e4b17023SJohn Marino { 758e4b17023SJohn Marino __s0 = __ctype(1) << __w0; 759e4b17023SJohn Marino __s1 = __s0 << 1; 760e4b17023SJohn Marino } 761e4b17023SJohn Marino 762e4b17023SJohn Marino __y0 = 0; 763e4b17023SJohn Marino __y1 = 0; 764e4b17023SJohn Marino if (__r) 765e4b17023SJohn Marino { 766e4b17023SJohn Marino __y0 = __s0 * (__r / __s0); 767e4b17023SJohn Marino if (__s1) 768e4b17023SJohn Marino __y1 = __s1 * (__r / __s1); 769e4b17023SJohn Marino 770e4b17023SJohn Marino if (__r - __y0 <= __y0 / __n) 771e4b17023SJohn Marino break; 772e4b17023SJohn Marino } 773e4b17023SJohn Marino else 774e4b17023SJohn Marino break; 775e4b17023SJohn Marino } 776e4b17023SJohn Marino 777e4b17023SJohn Marino result_type __sum = 0; 778e4b17023SJohn Marino for (size_t __k = 0; __k < __n0; ++__k) 779e4b17023SJohn Marino { 780e4b17023SJohn Marino __ctype __u; 781e4b17023SJohn Marino do 782e4b17023SJohn Marino __u = _M_b() - _M_b.min(); 783e4b17023SJohn Marino while (__y0 && __u >= __y0); 784e4b17023SJohn Marino __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u); 785e4b17023SJohn Marino } 786e4b17023SJohn Marino for (size_t __k = __n0; __k < __n; ++__k) 787e4b17023SJohn Marino { 788e4b17023SJohn Marino __ctype __u; 789e4b17023SJohn Marino do 790e4b17023SJohn Marino __u = _M_b() - _M_b.min(); 791e4b17023SJohn Marino while (__y1 && __u >= __y1); 792e4b17023SJohn Marino __sum = __s1 * __sum + (__s1 ? __u % __s1 : __u); 793e4b17023SJohn Marino } 794e4b17023SJohn Marino return __sum; 795e4b17023SJohn Marino } 796e4b17023SJohn Marino 797e4b17023SJohn Marino 798e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __k> 799e4b17023SJohn Marino constexpr size_t 800e4b17023SJohn Marino shuffle_order_engine<_RandomNumberEngine, __k>::table_size; 801e4b17023SJohn Marino 802e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __k> 803e4b17023SJohn Marino typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type 804e4b17023SJohn Marino shuffle_order_engine<_RandomNumberEngine, __k>:: operator ()()805e4b17023SJohn Marino operator()() 806e4b17023SJohn Marino { 807e4b17023SJohn Marino size_t __j = __k * ((_M_y - _M_b.min()) 808e4b17023SJohn Marino / (_M_b.max() - _M_b.min() + 1.0L)); 809e4b17023SJohn Marino _M_y = _M_v[__j]; 810e4b17023SJohn Marino _M_v[__j] = _M_b(); 811e4b17023SJohn Marino 812e4b17023SJohn Marino return _M_y; 813e4b17023SJohn Marino } 814e4b17023SJohn Marino 815e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __k, 816e4b17023SJohn Marino typename _CharT, typename _Traits> 817e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const shuffle_order_engine<_RandomNumberEngine,__k> & __x)818e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 819e4b17023SJohn Marino const shuffle_order_engine<_RandomNumberEngine, __k>& __x) 820e4b17023SJohn Marino { 821e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 822e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 823e4b17023SJohn Marino 824e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 825e4b17023SJohn Marino const _CharT __fill = __os.fill(); 826e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 827e4b17023SJohn Marino __os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left); 828e4b17023SJohn Marino __os.fill(__space); 829e4b17023SJohn Marino 830e4b17023SJohn Marino __os << __x.base(); 831e4b17023SJohn Marino for (size_t __i = 0; __i < __k; ++__i) 832e4b17023SJohn Marino __os << __space << __x._M_v[__i]; 833e4b17023SJohn Marino __os << __space << __x._M_y; 834e4b17023SJohn Marino 835e4b17023SJohn Marino __os.flags(__flags); 836e4b17023SJohn Marino __os.fill(__fill); 837e4b17023SJohn Marino return __os; 838e4b17023SJohn Marino } 839e4b17023SJohn Marino 840e4b17023SJohn Marino template<typename _RandomNumberEngine, size_t __k, 841e4b17023SJohn Marino typename _CharT, typename _Traits> 842e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,shuffle_order_engine<_RandomNumberEngine,__k> & __x)843e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 844e4b17023SJohn Marino shuffle_order_engine<_RandomNumberEngine, __k>& __x) 845e4b17023SJohn Marino { 846e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 847e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 848e4b17023SJohn Marino 849e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 850e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 851e4b17023SJohn Marino 852e4b17023SJohn Marino __is >> __x._M_b; 853e4b17023SJohn Marino for (size_t __i = 0; __i < __k; ++__i) 854e4b17023SJohn Marino __is >> __x._M_v[__i]; 855e4b17023SJohn Marino __is >> __x._M_y; 856e4b17023SJohn Marino 857e4b17023SJohn Marino __is.flags(__flags); 858e4b17023SJohn Marino return __is; 859e4b17023SJohn Marino } 860e4b17023SJohn Marino 861e4b17023SJohn Marino 862e4b17023SJohn Marino template<typename _IntType> 863e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 864e4b17023SJohn Marino typename uniform_int_distribution<_IntType>::result_type 865e4b17023SJohn Marino uniform_int_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)866e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 867e4b17023SJohn Marino const param_type& __param) 868e4b17023SJohn Marino { 869e4b17023SJohn Marino typedef typename _UniformRandomNumberGenerator::result_type 870e4b17023SJohn Marino _Gresult_type; 871e4b17023SJohn Marino typedef typename std::make_unsigned<result_type>::type __utype; 872e4b17023SJohn Marino typedef typename std::common_type<_Gresult_type, __utype>::type 873e4b17023SJohn Marino __uctype; 874e4b17023SJohn Marino 875e4b17023SJohn Marino const __uctype __urngmin = __urng.min(); 876e4b17023SJohn Marino const __uctype __urngmax = __urng.max(); 877e4b17023SJohn Marino const __uctype __urngrange = __urngmax - __urngmin; 878e4b17023SJohn Marino const __uctype __urange 879e4b17023SJohn Marino = __uctype(__param.b()) - __uctype(__param.a()); 880e4b17023SJohn Marino 881e4b17023SJohn Marino __uctype __ret; 882e4b17023SJohn Marino 883e4b17023SJohn Marino if (__urngrange > __urange) 884e4b17023SJohn Marino { 885e4b17023SJohn Marino // downscaling 886e4b17023SJohn Marino const __uctype __uerange = __urange + 1; // __urange can be zero 887e4b17023SJohn Marino const __uctype __scaling = __urngrange / __uerange; 888e4b17023SJohn Marino const __uctype __past = __uerange * __scaling; 889e4b17023SJohn Marino do 890e4b17023SJohn Marino __ret = __uctype(__urng()) - __urngmin; 891e4b17023SJohn Marino while (__ret >= __past); 892e4b17023SJohn Marino __ret /= __scaling; 893e4b17023SJohn Marino } 894e4b17023SJohn Marino else if (__urngrange < __urange) 895e4b17023SJohn Marino { 896e4b17023SJohn Marino // upscaling 897e4b17023SJohn Marino /* 898e4b17023SJohn Marino Note that every value in [0, urange] 899e4b17023SJohn Marino can be written uniquely as 900e4b17023SJohn Marino 901e4b17023SJohn Marino (urngrange + 1) * high + low 902e4b17023SJohn Marino 903e4b17023SJohn Marino where 904e4b17023SJohn Marino 905e4b17023SJohn Marino high in [0, urange / (urngrange + 1)] 906e4b17023SJohn Marino 907e4b17023SJohn Marino and 908e4b17023SJohn Marino 909e4b17023SJohn Marino low in [0, urngrange]. 910e4b17023SJohn Marino */ 911e4b17023SJohn Marino __uctype __tmp; // wraparound control 912e4b17023SJohn Marino do 913e4b17023SJohn Marino { 914e4b17023SJohn Marino const __uctype __uerngrange = __urngrange + 1; 915e4b17023SJohn Marino __tmp = (__uerngrange * operator() 916e4b17023SJohn Marino (__urng, param_type(0, __urange / __uerngrange))); 917e4b17023SJohn Marino __ret = __tmp + (__uctype(__urng()) - __urngmin); 918e4b17023SJohn Marino } 919e4b17023SJohn Marino while (__ret > __urange || __ret < __tmp); 920e4b17023SJohn Marino } 921e4b17023SJohn Marino else 922e4b17023SJohn Marino __ret = __uctype(__urng()) - __urngmin; 923e4b17023SJohn Marino 924e4b17023SJohn Marino return __ret + __param.a(); 925e4b17023SJohn Marino } 926e4b17023SJohn Marino 927e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 928e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const uniform_int_distribution<_IntType> & __x)929e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 930e4b17023SJohn Marino const uniform_int_distribution<_IntType>& __x) 931e4b17023SJohn Marino { 932e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 933e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 934e4b17023SJohn Marino 935e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 936e4b17023SJohn Marino const _CharT __fill = __os.fill(); 937e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 938e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 939e4b17023SJohn Marino __os.fill(__space); 940e4b17023SJohn Marino 941e4b17023SJohn Marino __os << __x.a() << __space << __x.b(); 942e4b17023SJohn Marino 943e4b17023SJohn Marino __os.flags(__flags); 944e4b17023SJohn Marino __os.fill(__fill); 945e4b17023SJohn Marino return __os; 946e4b17023SJohn Marino } 947e4b17023SJohn Marino 948e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 949e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,uniform_int_distribution<_IntType> & __x)950e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 951e4b17023SJohn Marino uniform_int_distribution<_IntType>& __x) 952e4b17023SJohn Marino { 953e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 954e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 955e4b17023SJohn Marino 956e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 957e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 958e4b17023SJohn Marino 959e4b17023SJohn Marino _IntType __a, __b; 960e4b17023SJohn Marino __is >> __a >> __b; 961e4b17023SJohn Marino __x.param(typename uniform_int_distribution<_IntType>:: 962e4b17023SJohn Marino param_type(__a, __b)); 963e4b17023SJohn Marino 964e4b17023SJohn Marino __is.flags(__flags); 965e4b17023SJohn Marino return __is; 966e4b17023SJohn Marino } 967e4b17023SJohn Marino 968e4b17023SJohn Marino 969e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 970e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const uniform_real_distribution<_RealType> & __x)971e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 972e4b17023SJohn Marino const uniform_real_distribution<_RealType>& __x) 973e4b17023SJohn Marino { 974e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 975e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 976e4b17023SJohn Marino 977e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 978e4b17023SJohn Marino const _CharT __fill = __os.fill(); 979e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 980e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 981e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 982e4b17023SJohn Marino __os.fill(__space); 983e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 984e4b17023SJohn Marino 985e4b17023SJohn Marino __os << __x.a() << __space << __x.b(); 986e4b17023SJohn Marino 987e4b17023SJohn Marino __os.flags(__flags); 988e4b17023SJohn Marino __os.fill(__fill); 989e4b17023SJohn Marino __os.precision(__precision); 990e4b17023SJohn Marino return __os; 991e4b17023SJohn Marino } 992e4b17023SJohn Marino 993e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 994e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,uniform_real_distribution<_RealType> & __x)995e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 996e4b17023SJohn Marino uniform_real_distribution<_RealType>& __x) 997e4b17023SJohn Marino { 998e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 999e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1000e4b17023SJohn Marino 1001e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1002e4b17023SJohn Marino __is.flags(__ios_base::skipws); 1003e4b17023SJohn Marino 1004e4b17023SJohn Marino _RealType __a, __b; 1005e4b17023SJohn Marino __is >> __a >> __b; 1006e4b17023SJohn Marino __x.param(typename uniform_real_distribution<_RealType>:: 1007e4b17023SJohn Marino param_type(__a, __b)); 1008e4b17023SJohn Marino 1009e4b17023SJohn Marino __is.flags(__flags); 1010e4b17023SJohn Marino return __is; 1011e4b17023SJohn Marino } 1012e4b17023SJohn Marino 1013e4b17023SJohn Marino 1014e4b17023SJohn Marino template<typename _CharT, typename _Traits> 1015e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const bernoulli_distribution & __x)1016e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1017e4b17023SJohn Marino const bernoulli_distribution& __x) 1018e4b17023SJohn Marino { 1019e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1020e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1021e4b17023SJohn Marino 1022e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1023e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1024e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1025e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1026e4b17023SJohn Marino __os.fill(__os.widen(' ')); 1027e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 1028e4b17023SJohn Marino 1029e4b17023SJohn Marino __os << __x.p(); 1030e4b17023SJohn Marino 1031e4b17023SJohn Marino __os.flags(__flags); 1032e4b17023SJohn Marino __os.fill(__fill); 1033e4b17023SJohn Marino __os.precision(__precision); 1034e4b17023SJohn Marino return __os; 1035e4b17023SJohn Marino } 1036e4b17023SJohn Marino 1037e4b17023SJohn Marino 1038e4b17023SJohn Marino template<typename _IntType> 1039e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1040e4b17023SJohn Marino typename geometric_distribution<_IntType>::result_type 1041e4b17023SJohn Marino geometric_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1042e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1043e4b17023SJohn Marino const param_type& __param) 1044e4b17023SJohn Marino { 1045e4b17023SJohn Marino // About the epsilon thing see this thread: 1046e4b17023SJohn Marino // http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html 1047e4b17023SJohn Marino const double __naf = 1048e4b17023SJohn Marino (1 - std::numeric_limits<double>::epsilon()) / 2; 1049e4b17023SJohn Marino // The largest _RealType convertible to _IntType. 1050e4b17023SJohn Marino const double __thr = 1051e4b17023SJohn Marino std::numeric_limits<_IntType>::max() + __naf; 1052e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1053e4b17023SJohn Marino __aurng(__urng); 1054e4b17023SJohn Marino 1055e4b17023SJohn Marino double __cand; 1056e4b17023SJohn Marino do 10575ce9237cSJohn Marino __cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p); 1058e4b17023SJohn Marino while (__cand >= __thr); 1059e4b17023SJohn Marino 1060e4b17023SJohn Marino return result_type(__cand + __naf); 1061e4b17023SJohn Marino } 1062e4b17023SJohn Marino 1063e4b17023SJohn Marino template<typename _IntType, 1064e4b17023SJohn Marino typename _CharT, typename _Traits> 1065e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const geometric_distribution<_IntType> & __x)1066e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1067e4b17023SJohn Marino const geometric_distribution<_IntType>& __x) 1068e4b17023SJohn Marino { 1069e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1070e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1071e4b17023SJohn Marino 1072e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1073e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1074e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1075e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1076e4b17023SJohn Marino __os.fill(__os.widen(' ')); 1077e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 1078e4b17023SJohn Marino 1079e4b17023SJohn Marino __os << __x.p(); 1080e4b17023SJohn Marino 1081e4b17023SJohn Marino __os.flags(__flags); 1082e4b17023SJohn Marino __os.fill(__fill); 1083e4b17023SJohn Marino __os.precision(__precision); 1084e4b17023SJohn Marino return __os; 1085e4b17023SJohn Marino } 1086e4b17023SJohn Marino 1087e4b17023SJohn Marino template<typename _IntType, 1088e4b17023SJohn Marino typename _CharT, typename _Traits> 1089e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,geometric_distribution<_IntType> & __x)1090e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1091e4b17023SJohn Marino geometric_distribution<_IntType>& __x) 1092e4b17023SJohn Marino { 1093e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1094e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1095e4b17023SJohn Marino 1096e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1097e4b17023SJohn Marino __is.flags(__ios_base::skipws); 1098e4b17023SJohn Marino 1099e4b17023SJohn Marino double __p; 1100e4b17023SJohn Marino __is >> __p; 1101e4b17023SJohn Marino __x.param(typename geometric_distribution<_IntType>::param_type(__p)); 1102e4b17023SJohn Marino 1103e4b17023SJohn Marino __is.flags(__flags); 1104e4b17023SJohn Marino return __is; 1105e4b17023SJohn Marino } 1106e4b17023SJohn Marino 1107e4b17023SJohn Marino // This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5. 1108e4b17023SJohn Marino template<typename _IntType> 1109e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1110e4b17023SJohn Marino typename negative_binomial_distribution<_IntType>::result_type 1111e4b17023SJohn Marino negative_binomial_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng)1112e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng) 1113e4b17023SJohn Marino { 1114e4b17023SJohn Marino const double __y = _M_gd(__urng); 1115e4b17023SJohn Marino 1116e4b17023SJohn Marino // XXX Is the constructor too slow? 1117e4b17023SJohn Marino std::poisson_distribution<result_type> __poisson(__y); 1118e4b17023SJohn Marino return __poisson(__urng); 1119e4b17023SJohn Marino } 1120e4b17023SJohn Marino 1121e4b17023SJohn Marino template<typename _IntType> 1122e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1123e4b17023SJohn Marino typename negative_binomial_distribution<_IntType>::result_type 1124e4b17023SJohn Marino negative_binomial_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1125e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1126e4b17023SJohn Marino const param_type& __p) 1127e4b17023SJohn Marino { 1128*95d28233SJohn Marino typedef typename std::gamma_distribution<double>::param_type 1129e4b17023SJohn Marino param_type; 1130e4b17023SJohn Marino 1131e4b17023SJohn Marino const double __y = 1132e4b17023SJohn Marino _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p())); 1133e4b17023SJohn Marino 1134e4b17023SJohn Marino std::poisson_distribution<result_type> __poisson(__y); 1135e4b17023SJohn Marino return __poisson(__urng); 1136e4b17023SJohn Marino } 1137e4b17023SJohn Marino 1138e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 1139e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const negative_binomial_distribution<_IntType> & __x)1140e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1141e4b17023SJohn Marino const negative_binomial_distribution<_IntType>& __x) 1142e4b17023SJohn Marino { 1143e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1144e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1145e4b17023SJohn Marino 1146e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1147e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1148e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1149e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1150e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1151e4b17023SJohn Marino __os.fill(__os.widen(' ')); 1152e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 1153e4b17023SJohn Marino 1154e4b17023SJohn Marino __os << __x.k() << __space << __x.p() 1155e4b17023SJohn Marino << __space << __x._M_gd; 1156e4b17023SJohn Marino 1157e4b17023SJohn Marino __os.flags(__flags); 1158e4b17023SJohn Marino __os.fill(__fill); 1159e4b17023SJohn Marino __os.precision(__precision); 1160e4b17023SJohn Marino return __os; 1161e4b17023SJohn Marino } 1162e4b17023SJohn Marino 1163e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 1164e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,negative_binomial_distribution<_IntType> & __x)1165e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1166e4b17023SJohn Marino negative_binomial_distribution<_IntType>& __x) 1167e4b17023SJohn Marino { 1168e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1169e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1170e4b17023SJohn Marino 1171e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1172e4b17023SJohn Marino __is.flags(__ios_base::skipws); 1173e4b17023SJohn Marino 1174e4b17023SJohn Marino _IntType __k; 1175e4b17023SJohn Marino double __p; 1176e4b17023SJohn Marino __is >> __k >> __p >> __x._M_gd; 1177e4b17023SJohn Marino __x.param(typename negative_binomial_distribution<_IntType>:: 1178e4b17023SJohn Marino param_type(__k, __p)); 1179e4b17023SJohn Marino 1180e4b17023SJohn Marino __is.flags(__flags); 1181e4b17023SJohn Marino return __is; 1182e4b17023SJohn Marino } 1183e4b17023SJohn Marino 1184e4b17023SJohn Marino 1185e4b17023SJohn Marino template<typename _IntType> 1186e4b17023SJohn Marino void 1187e4b17023SJohn Marino poisson_distribution<_IntType>::param_type:: _M_initialize()1188e4b17023SJohn Marino _M_initialize() 1189e4b17023SJohn Marino { 1190e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 1191e4b17023SJohn Marino if (_M_mean >= 12) 1192e4b17023SJohn Marino { 1193e4b17023SJohn Marino const double __m = std::floor(_M_mean); 1194e4b17023SJohn Marino _M_lm_thr = std::log(_M_mean); 1195e4b17023SJohn Marino _M_lfm = std::lgamma(__m + 1); 1196e4b17023SJohn Marino _M_sm = std::sqrt(__m); 1197e4b17023SJohn Marino 1198e4b17023SJohn Marino const double __pi_4 = 0.7853981633974483096156608458198757L; 1199e4b17023SJohn Marino const double __dx = std::sqrt(2 * __m * std::log(32 * __m 1200e4b17023SJohn Marino / __pi_4)); 1201e4b17023SJohn Marino _M_d = std::round(std::max(6.0, std::min(__m, __dx))); 1202e4b17023SJohn Marino const double __cx = 2 * __m + _M_d; 1203e4b17023SJohn Marino _M_scx = std::sqrt(__cx / 2); 1204e4b17023SJohn Marino _M_1cx = 1 / __cx; 1205e4b17023SJohn Marino 1206e4b17023SJohn Marino _M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx); 1207e4b17023SJohn Marino _M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2)) 1208e4b17023SJohn Marino / _M_d; 1209e4b17023SJohn Marino } 1210e4b17023SJohn Marino else 1211e4b17023SJohn Marino #endif 1212e4b17023SJohn Marino _M_lm_thr = std::exp(-_M_mean); 1213e4b17023SJohn Marino } 1214e4b17023SJohn Marino 1215e4b17023SJohn Marino /** 1216e4b17023SJohn Marino * A rejection algorithm when mean >= 12 and a simple method based 1217e4b17023SJohn Marino * upon the multiplication of uniform random variates otherwise. 1218e4b17023SJohn Marino * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1219e4b17023SJohn Marino * is defined. 1220e4b17023SJohn Marino * 1221e4b17023SJohn Marino * Reference: 1222e4b17023SJohn Marino * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1223e4b17023SJohn Marino * New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!). 1224e4b17023SJohn Marino */ 1225e4b17023SJohn Marino template<typename _IntType> 1226e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1227e4b17023SJohn Marino typename poisson_distribution<_IntType>::result_type 1228e4b17023SJohn Marino poisson_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1229e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1230e4b17023SJohn Marino const param_type& __param) 1231e4b17023SJohn Marino { 1232e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1233e4b17023SJohn Marino __aurng(__urng); 1234e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 1235e4b17023SJohn Marino if (__param.mean() >= 12) 1236e4b17023SJohn Marino { 1237e4b17023SJohn Marino double __x; 1238e4b17023SJohn Marino 1239e4b17023SJohn Marino // See comments above... 1240e4b17023SJohn Marino const double __naf = 1241e4b17023SJohn Marino (1 - std::numeric_limits<double>::epsilon()) / 2; 1242e4b17023SJohn Marino const double __thr = 1243e4b17023SJohn Marino std::numeric_limits<_IntType>::max() + __naf; 1244e4b17023SJohn Marino 1245e4b17023SJohn Marino const double __m = std::floor(__param.mean()); 1246e4b17023SJohn Marino // sqrt(pi / 2) 1247e4b17023SJohn Marino const double __spi_2 = 1.2533141373155002512078826424055226L; 1248e4b17023SJohn Marino const double __c1 = __param._M_sm * __spi_2; 1249e4b17023SJohn Marino const double __c2 = __param._M_c2b + __c1; 1250e4b17023SJohn Marino const double __c3 = __c2 + 1; 1251e4b17023SJohn Marino const double __c4 = __c3 + 1; 1252e4b17023SJohn Marino // e^(1 / 78) 1253e4b17023SJohn Marino const double __e178 = 1.0129030479320018583185514777512983L; 1254e4b17023SJohn Marino const double __c5 = __c4 + __e178; 1255e4b17023SJohn Marino const double __c = __param._M_cb + __c5; 1256e4b17023SJohn Marino const double __2cx = 2 * (2 * __m + __param._M_d); 1257e4b17023SJohn Marino 1258e4b17023SJohn Marino bool __reject = true; 1259e4b17023SJohn Marino do 1260e4b17023SJohn Marino { 1261e4b17023SJohn Marino const double __u = __c * __aurng(); 12625ce9237cSJohn Marino const double __e = -std::log(1.0 - __aurng()); 1263e4b17023SJohn Marino 1264e4b17023SJohn Marino double __w = 0.0; 1265e4b17023SJohn Marino 1266e4b17023SJohn Marino if (__u <= __c1) 1267e4b17023SJohn Marino { 1268e4b17023SJohn Marino const double __n = _M_nd(__urng); 1269e4b17023SJohn Marino const double __y = -std::abs(__n) * __param._M_sm - 1; 1270e4b17023SJohn Marino __x = std::floor(__y); 1271e4b17023SJohn Marino __w = -__n * __n / 2; 1272e4b17023SJohn Marino if (__x < -__m) 1273e4b17023SJohn Marino continue; 1274e4b17023SJohn Marino } 1275e4b17023SJohn Marino else if (__u <= __c2) 1276e4b17023SJohn Marino { 1277e4b17023SJohn Marino const double __n = _M_nd(__urng); 1278e4b17023SJohn Marino const double __y = 1 + std::abs(__n) * __param._M_scx; 1279e4b17023SJohn Marino __x = std::ceil(__y); 1280e4b17023SJohn Marino __w = __y * (2 - __y) * __param._M_1cx; 1281e4b17023SJohn Marino if (__x > __param._M_d) 1282e4b17023SJohn Marino continue; 1283e4b17023SJohn Marino } 1284e4b17023SJohn Marino else if (__u <= __c3) 1285e4b17023SJohn Marino // NB: This case not in the book, nor in the Errata, 1286e4b17023SJohn Marino // but should be ok... 1287e4b17023SJohn Marino __x = -1; 1288e4b17023SJohn Marino else if (__u <= __c4) 1289e4b17023SJohn Marino __x = 0; 1290e4b17023SJohn Marino else if (__u <= __c5) 1291e4b17023SJohn Marino __x = 1; 1292e4b17023SJohn Marino else 1293e4b17023SJohn Marino { 12945ce9237cSJohn Marino const double __v = -std::log(1.0 - __aurng()); 1295e4b17023SJohn Marino const double __y = __param._M_d 1296e4b17023SJohn Marino + __v * __2cx / __param._M_d; 1297e4b17023SJohn Marino __x = std::ceil(__y); 1298e4b17023SJohn Marino __w = -__param._M_d * __param._M_1cx * (1 + __y / 2); 1299e4b17023SJohn Marino } 1300e4b17023SJohn Marino 1301e4b17023SJohn Marino __reject = (__w - __e - __x * __param._M_lm_thr 1302e4b17023SJohn Marino > __param._M_lfm - std::lgamma(__x + __m + 1)); 1303e4b17023SJohn Marino 1304e4b17023SJohn Marino __reject |= __x + __m >= __thr; 1305e4b17023SJohn Marino 1306e4b17023SJohn Marino } while (__reject); 1307e4b17023SJohn Marino 1308e4b17023SJohn Marino return result_type(__x + __m + __naf); 1309e4b17023SJohn Marino } 1310e4b17023SJohn Marino else 1311e4b17023SJohn Marino #endif 1312e4b17023SJohn Marino { 1313e4b17023SJohn Marino _IntType __x = 0; 1314e4b17023SJohn Marino double __prod = 1.0; 1315e4b17023SJohn Marino 1316e4b17023SJohn Marino do 1317e4b17023SJohn Marino { 1318e4b17023SJohn Marino __prod *= __aurng(); 1319e4b17023SJohn Marino __x += 1; 1320e4b17023SJohn Marino } 1321e4b17023SJohn Marino while (__prod > __param._M_lm_thr); 1322e4b17023SJohn Marino 1323e4b17023SJohn Marino return __x - 1; 1324e4b17023SJohn Marino } 1325e4b17023SJohn Marino } 1326e4b17023SJohn Marino 1327e4b17023SJohn Marino template<typename _IntType, 1328e4b17023SJohn Marino typename _CharT, typename _Traits> 1329e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const poisson_distribution<_IntType> & __x)1330e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1331e4b17023SJohn Marino const poisson_distribution<_IntType>& __x) 1332e4b17023SJohn Marino { 1333e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1334e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1335e4b17023SJohn Marino 1336e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1337e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1338e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1339e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1340e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1341e4b17023SJohn Marino __os.fill(__space); 1342e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 1343e4b17023SJohn Marino 1344e4b17023SJohn Marino __os << __x.mean() << __space << __x._M_nd; 1345e4b17023SJohn Marino 1346e4b17023SJohn Marino __os.flags(__flags); 1347e4b17023SJohn Marino __os.fill(__fill); 1348e4b17023SJohn Marino __os.precision(__precision); 1349e4b17023SJohn Marino return __os; 1350e4b17023SJohn Marino } 1351e4b17023SJohn Marino 1352e4b17023SJohn Marino template<typename _IntType, 1353e4b17023SJohn Marino typename _CharT, typename _Traits> 1354e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,poisson_distribution<_IntType> & __x)1355e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1356e4b17023SJohn Marino poisson_distribution<_IntType>& __x) 1357e4b17023SJohn Marino { 1358e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1359e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1360e4b17023SJohn Marino 1361e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1362e4b17023SJohn Marino __is.flags(__ios_base::skipws); 1363e4b17023SJohn Marino 1364e4b17023SJohn Marino double __mean; 1365e4b17023SJohn Marino __is >> __mean >> __x._M_nd; 1366e4b17023SJohn Marino __x.param(typename poisson_distribution<_IntType>::param_type(__mean)); 1367e4b17023SJohn Marino 1368e4b17023SJohn Marino __is.flags(__flags); 1369e4b17023SJohn Marino return __is; 1370e4b17023SJohn Marino } 1371e4b17023SJohn Marino 1372e4b17023SJohn Marino 1373e4b17023SJohn Marino template<typename _IntType> 1374e4b17023SJohn Marino void 1375e4b17023SJohn Marino binomial_distribution<_IntType>::param_type:: _M_initialize()1376e4b17023SJohn Marino _M_initialize() 1377e4b17023SJohn Marino { 1378e4b17023SJohn Marino const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p; 1379e4b17023SJohn Marino 1380e4b17023SJohn Marino _M_easy = true; 1381e4b17023SJohn Marino 1382e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 1383e4b17023SJohn Marino if (_M_t * __p12 >= 8) 1384e4b17023SJohn Marino { 1385e4b17023SJohn Marino _M_easy = false; 1386e4b17023SJohn Marino const double __np = std::floor(_M_t * __p12); 1387e4b17023SJohn Marino const double __pa = __np / _M_t; 1388e4b17023SJohn Marino const double __1p = 1 - __pa; 1389e4b17023SJohn Marino 1390e4b17023SJohn Marino const double __pi_4 = 0.7853981633974483096156608458198757L; 1391e4b17023SJohn Marino const double __d1x = 1392e4b17023SJohn Marino std::sqrt(__np * __1p * std::log(32 * __np 1393e4b17023SJohn Marino / (81 * __pi_4 * __1p))); 1394e4b17023SJohn Marino _M_d1 = std::round(std::max(1.0, __d1x)); 1395e4b17023SJohn Marino const double __d2x = 1396e4b17023SJohn Marino std::sqrt(__np * __1p * std::log(32 * _M_t * __1p 1397e4b17023SJohn Marino / (__pi_4 * __pa))); 1398e4b17023SJohn Marino _M_d2 = std::round(std::max(1.0, __d2x)); 1399e4b17023SJohn Marino 1400e4b17023SJohn Marino // sqrt(pi / 2) 1401e4b17023SJohn Marino const double __spi_2 = 1.2533141373155002512078826424055226L; 1402e4b17023SJohn Marino _M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np)); 1403e4b17023SJohn Marino _M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p)); 1404e4b17023SJohn Marino _M_c = 2 * _M_d1 / __np; 1405e4b17023SJohn Marino _M_a1 = std::exp(_M_c) * _M_s1 * __spi_2; 1406e4b17023SJohn Marino const double __a12 = _M_a1 + _M_s2 * __spi_2; 1407e4b17023SJohn Marino const double __s1s = _M_s1 * _M_s1; 1408e4b17023SJohn Marino _M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p)) 1409e4b17023SJohn Marino * 2 * __s1s / _M_d1 1410e4b17023SJohn Marino * std::exp(-_M_d1 * _M_d1 / (2 * __s1s))); 1411e4b17023SJohn Marino const double __s2s = _M_s2 * _M_s2; 1412e4b17023SJohn Marino _M_s = (_M_a123 + 2 * __s2s / _M_d2 1413e4b17023SJohn Marino * std::exp(-_M_d2 * _M_d2 / (2 * __s2s))); 1414e4b17023SJohn Marino _M_lf = (std::lgamma(__np + 1) 1415e4b17023SJohn Marino + std::lgamma(_M_t - __np + 1)); 1416e4b17023SJohn Marino _M_lp1p = std::log(__pa / __1p); 1417e4b17023SJohn Marino 1418e4b17023SJohn Marino _M_q = -std::log(1 - (__p12 - __pa) / __1p); 1419e4b17023SJohn Marino } 1420e4b17023SJohn Marino else 1421e4b17023SJohn Marino #endif 1422e4b17023SJohn Marino _M_q = -std::log(1 - __p12); 1423e4b17023SJohn Marino } 1424e4b17023SJohn Marino 1425e4b17023SJohn Marino template<typename _IntType> 1426e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1427e4b17023SJohn Marino typename binomial_distribution<_IntType>::result_type 1428e4b17023SJohn Marino binomial_distribution<_IntType>:: _M_waiting(_UniformRandomNumberGenerator & __urng,_IntType __t)1429e4b17023SJohn Marino _M_waiting(_UniformRandomNumberGenerator& __urng, _IntType __t) 1430e4b17023SJohn Marino { 1431e4b17023SJohn Marino _IntType __x = 0; 1432e4b17023SJohn Marino double __sum = 0.0; 1433e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1434e4b17023SJohn Marino __aurng(__urng); 1435e4b17023SJohn Marino 1436e4b17023SJohn Marino do 1437e4b17023SJohn Marino { 14385ce9237cSJohn Marino const double __e = -std::log(1.0 - __aurng()); 1439e4b17023SJohn Marino __sum += __e / (__t - __x); 1440e4b17023SJohn Marino __x += 1; 1441e4b17023SJohn Marino } 1442e4b17023SJohn Marino while (__sum <= _M_param._M_q); 1443e4b17023SJohn Marino 1444e4b17023SJohn Marino return __x - 1; 1445e4b17023SJohn Marino } 1446e4b17023SJohn Marino 1447e4b17023SJohn Marino /** 1448e4b17023SJohn Marino * A rejection algorithm when t * p >= 8 and a simple waiting time 1449e4b17023SJohn Marino * method - the second in the referenced book - otherwise. 1450e4b17023SJohn Marino * NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1 1451e4b17023SJohn Marino * is defined. 1452e4b17023SJohn Marino * 1453e4b17023SJohn Marino * Reference: 1454e4b17023SJohn Marino * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1455e4b17023SJohn Marino * New York, 1986, Ch. X, Sect. 4 (+ Errata!). 1456e4b17023SJohn Marino */ 1457e4b17023SJohn Marino template<typename _IntType> 1458e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1459e4b17023SJohn Marino typename binomial_distribution<_IntType>::result_type 1460e4b17023SJohn Marino binomial_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1461e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1462e4b17023SJohn Marino const param_type& __param) 1463e4b17023SJohn Marino { 1464e4b17023SJohn Marino result_type __ret; 1465e4b17023SJohn Marino const _IntType __t = __param.t(); 1466e4b17023SJohn Marino const double __p = __param.p(); 1467e4b17023SJohn Marino const double __p12 = __p <= 0.5 ? __p : 1.0 - __p; 1468e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 1469e4b17023SJohn Marino __aurng(__urng); 1470e4b17023SJohn Marino 1471e4b17023SJohn Marino #if _GLIBCXX_USE_C99_MATH_TR1 1472e4b17023SJohn Marino if (!__param._M_easy) 1473e4b17023SJohn Marino { 1474e4b17023SJohn Marino double __x; 1475e4b17023SJohn Marino 1476e4b17023SJohn Marino // See comments above... 1477e4b17023SJohn Marino const double __naf = 1478e4b17023SJohn Marino (1 - std::numeric_limits<double>::epsilon()) / 2; 1479e4b17023SJohn Marino const double __thr = 1480e4b17023SJohn Marino std::numeric_limits<_IntType>::max() + __naf; 1481e4b17023SJohn Marino 1482e4b17023SJohn Marino const double __np = std::floor(__t * __p12); 1483e4b17023SJohn Marino 1484e4b17023SJohn Marino // sqrt(pi / 2) 1485e4b17023SJohn Marino const double __spi_2 = 1.2533141373155002512078826424055226L; 1486e4b17023SJohn Marino const double __a1 = __param._M_a1; 1487e4b17023SJohn Marino const double __a12 = __a1 + __param._M_s2 * __spi_2; 1488e4b17023SJohn Marino const double __a123 = __param._M_a123; 1489e4b17023SJohn Marino const double __s1s = __param._M_s1 * __param._M_s1; 1490e4b17023SJohn Marino const double __s2s = __param._M_s2 * __param._M_s2; 1491e4b17023SJohn Marino 1492e4b17023SJohn Marino bool __reject; 1493e4b17023SJohn Marino do 1494e4b17023SJohn Marino { 1495e4b17023SJohn Marino const double __u = __param._M_s * __aurng(); 1496e4b17023SJohn Marino 1497e4b17023SJohn Marino double __v; 1498e4b17023SJohn Marino 1499e4b17023SJohn Marino if (__u <= __a1) 1500e4b17023SJohn Marino { 1501e4b17023SJohn Marino const double __n = _M_nd(__urng); 1502e4b17023SJohn Marino const double __y = __param._M_s1 * std::abs(__n); 1503e4b17023SJohn Marino __reject = __y >= __param._M_d1; 1504e4b17023SJohn Marino if (!__reject) 1505e4b17023SJohn Marino { 15065ce9237cSJohn Marino const double __e = -std::log(1.0 - __aurng()); 1507e4b17023SJohn Marino __x = std::floor(__y); 1508e4b17023SJohn Marino __v = -__e - __n * __n / 2 + __param._M_c; 1509e4b17023SJohn Marino } 1510e4b17023SJohn Marino } 1511e4b17023SJohn Marino else if (__u <= __a12) 1512e4b17023SJohn Marino { 1513e4b17023SJohn Marino const double __n = _M_nd(__urng); 1514e4b17023SJohn Marino const double __y = __param._M_s2 * std::abs(__n); 1515e4b17023SJohn Marino __reject = __y >= __param._M_d2; 1516e4b17023SJohn Marino if (!__reject) 1517e4b17023SJohn Marino { 15185ce9237cSJohn Marino const double __e = -std::log(1.0 - __aurng()); 1519e4b17023SJohn Marino __x = std::floor(-__y); 1520e4b17023SJohn Marino __v = -__e - __n * __n / 2; 1521e4b17023SJohn Marino } 1522e4b17023SJohn Marino } 1523e4b17023SJohn Marino else if (__u <= __a123) 1524e4b17023SJohn Marino { 15255ce9237cSJohn Marino const double __e1 = -std::log(1.0 - __aurng()); 15265ce9237cSJohn Marino const double __e2 = -std::log(1.0 - __aurng()); 1527e4b17023SJohn Marino 1528e4b17023SJohn Marino const double __y = __param._M_d1 1529e4b17023SJohn Marino + 2 * __s1s * __e1 / __param._M_d1; 1530e4b17023SJohn Marino __x = std::floor(__y); 1531e4b17023SJohn Marino __v = (-__e2 + __param._M_d1 * (1 / (__t - __np) 1532e4b17023SJohn Marino -__y / (2 * __s1s))); 1533e4b17023SJohn Marino __reject = false; 1534e4b17023SJohn Marino } 1535e4b17023SJohn Marino else 1536e4b17023SJohn Marino { 15375ce9237cSJohn Marino const double __e1 = -std::log(1.0 - __aurng()); 15385ce9237cSJohn Marino const double __e2 = -std::log(1.0 - __aurng()); 1539e4b17023SJohn Marino 1540e4b17023SJohn Marino const double __y = __param._M_d2 1541e4b17023SJohn Marino + 2 * __s2s * __e1 / __param._M_d2; 1542e4b17023SJohn Marino __x = std::floor(-__y); 1543e4b17023SJohn Marino __v = -__e2 - __param._M_d2 * __y / (2 * __s2s); 1544e4b17023SJohn Marino __reject = false; 1545e4b17023SJohn Marino } 1546e4b17023SJohn Marino 1547e4b17023SJohn Marino __reject = __reject || __x < -__np || __x > __t - __np; 1548e4b17023SJohn Marino if (!__reject) 1549e4b17023SJohn Marino { 1550e4b17023SJohn Marino const double __lfx = 1551e4b17023SJohn Marino std::lgamma(__np + __x + 1) 1552e4b17023SJohn Marino + std::lgamma(__t - (__np + __x) + 1); 1553e4b17023SJohn Marino __reject = __v > __param._M_lf - __lfx 1554e4b17023SJohn Marino + __x * __param._M_lp1p; 1555e4b17023SJohn Marino } 1556e4b17023SJohn Marino 1557e4b17023SJohn Marino __reject |= __x + __np >= __thr; 1558e4b17023SJohn Marino } 1559e4b17023SJohn Marino while (__reject); 1560e4b17023SJohn Marino 1561e4b17023SJohn Marino __x += __np + __naf; 1562e4b17023SJohn Marino 1563e4b17023SJohn Marino const _IntType __z = _M_waiting(__urng, __t - _IntType(__x)); 1564e4b17023SJohn Marino __ret = _IntType(__x) + __z; 1565e4b17023SJohn Marino } 1566e4b17023SJohn Marino else 1567e4b17023SJohn Marino #endif 1568e4b17023SJohn Marino __ret = _M_waiting(__urng, __t); 1569e4b17023SJohn Marino 1570e4b17023SJohn Marino if (__p12 != __p) 1571e4b17023SJohn Marino __ret = __t - __ret; 1572e4b17023SJohn Marino return __ret; 1573e4b17023SJohn Marino } 1574e4b17023SJohn Marino 1575e4b17023SJohn Marino template<typename _IntType, 1576e4b17023SJohn Marino typename _CharT, typename _Traits> 1577e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const binomial_distribution<_IntType> & __x)1578e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1579e4b17023SJohn Marino const binomial_distribution<_IntType>& __x) 1580e4b17023SJohn Marino { 1581e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1582e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1583e4b17023SJohn Marino 1584e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1585e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1586e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1587e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1588e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1589e4b17023SJohn Marino __os.fill(__space); 1590e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 1591e4b17023SJohn Marino 1592e4b17023SJohn Marino __os << __x.t() << __space << __x.p() 1593e4b17023SJohn Marino << __space << __x._M_nd; 1594e4b17023SJohn Marino 1595e4b17023SJohn Marino __os.flags(__flags); 1596e4b17023SJohn Marino __os.fill(__fill); 1597e4b17023SJohn Marino __os.precision(__precision); 1598e4b17023SJohn Marino return __os; 1599e4b17023SJohn Marino } 1600e4b17023SJohn Marino 1601e4b17023SJohn Marino template<typename _IntType, 1602e4b17023SJohn Marino typename _CharT, typename _Traits> 1603e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,binomial_distribution<_IntType> & __x)1604e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1605e4b17023SJohn Marino binomial_distribution<_IntType>& __x) 1606e4b17023SJohn Marino { 1607e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1608e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1609e4b17023SJohn Marino 1610e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1611e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1612e4b17023SJohn Marino 1613e4b17023SJohn Marino _IntType __t; 1614e4b17023SJohn Marino double __p; 1615e4b17023SJohn Marino __is >> __t >> __p >> __x._M_nd; 1616e4b17023SJohn Marino __x.param(typename binomial_distribution<_IntType>:: 1617e4b17023SJohn Marino param_type(__t, __p)); 1618e4b17023SJohn Marino 1619e4b17023SJohn Marino __is.flags(__flags); 1620e4b17023SJohn Marino return __is; 1621e4b17023SJohn Marino } 1622e4b17023SJohn Marino 1623e4b17023SJohn Marino 1624e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1625e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const exponential_distribution<_RealType> & __x)1626e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1627e4b17023SJohn Marino const exponential_distribution<_RealType>& __x) 1628e4b17023SJohn Marino { 1629e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1630e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1631e4b17023SJohn Marino 1632e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1633e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1634e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1635e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1636e4b17023SJohn Marino __os.fill(__os.widen(' ')); 1637e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1638e4b17023SJohn Marino 1639e4b17023SJohn Marino __os << __x.lambda(); 1640e4b17023SJohn Marino 1641e4b17023SJohn Marino __os.flags(__flags); 1642e4b17023SJohn Marino __os.fill(__fill); 1643e4b17023SJohn Marino __os.precision(__precision); 1644e4b17023SJohn Marino return __os; 1645e4b17023SJohn Marino } 1646e4b17023SJohn Marino 1647e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1648e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,exponential_distribution<_RealType> & __x)1649e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1650e4b17023SJohn Marino exponential_distribution<_RealType>& __x) 1651e4b17023SJohn Marino { 1652e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1653e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1654e4b17023SJohn Marino 1655e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1656e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1657e4b17023SJohn Marino 1658e4b17023SJohn Marino _RealType __lambda; 1659e4b17023SJohn Marino __is >> __lambda; 1660e4b17023SJohn Marino __x.param(typename exponential_distribution<_RealType>:: 1661e4b17023SJohn Marino param_type(__lambda)); 1662e4b17023SJohn Marino 1663e4b17023SJohn Marino __is.flags(__flags); 1664e4b17023SJohn Marino return __is; 1665e4b17023SJohn Marino } 1666e4b17023SJohn Marino 1667e4b17023SJohn Marino 1668e4b17023SJohn Marino /** 1669e4b17023SJohn Marino * Polar method due to Marsaglia. 1670e4b17023SJohn Marino * 1671e4b17023SJohn Marino * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag, 1672e4b17023SJohn Marino * New York, 1986, Ch. V, Sect. 4.4. 1673e4b17023SJohn Marino */ 1674e4b17023SJohn Marino template<typename _RealType> 1675e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1676e4b17023SJohn Marino typename normal_distribution<_RealType>::result_type 1677e4b17023SJohn Marino normal_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)1678e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1679e4b17023SJohn Marino const param_type& __param) 1680e4b17023SJohn Marino { 1681e4b17023SJohn Marino result_type __ret; 1682e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1683e4b17023SJohn Marino __aurng(__urng); 1684e4b17023SJohn Marino 1685e4b17023SJohn Marino if (_M_saved_available) 1686e4b17023SJohn Marino { 1687e4b17023SJohn Marino _M_saved_available = false; 1688e4b17023SJohn Marino __ret = _M_saved; 1689e4b17023SJohn Marino } 1690e4b17023SJohn Marino else 1691e4b17023SJohn Marino { 1692e4b17023SJohn Marino result_type __x, __y, __r2; 1693e4b17023SJohn Marino do 1694e4b17023SJohn Marino { 1695e4b17023SJohn Marino __x = result_type(2.0) * __aurng() - 1.0; 1696e4b17023SJohn Marino __y = result_type(2.0) * __aurng() - 1.0; 1697e4b17023SJohn Marino __r2 = __x * __x + __y * __y; 1698e4b17023SJohn Marino } 1699e4b17023SJohn Marino while (__r2 > 1.0 || __r2 == 0.0); 1700e4b17023SJohn Marino 1701e4b17023SJohn Marino const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2); 1702e4b17023SJohn Marino _M_saved = __x * __mult; 1703e4b17023SJohn Marino _M_saved_available = true; 1704e4b17023SJohn Marino __ret = __y * __mult; 1705e4b17023SJohn Marino } 1706e4b17023SJohn Marino 1707e4b17023SJohn Marino __ret = __ret * __param.stddev() + __param.mean(); 1708e4b17023SJohn Marino return __ret; 1709e4b17023SJohn Marino } 1710e4b17023SJohn Marino 1711e4b17023SJohn Marino template<typename _RealType> 1712e4b17023SJohn Marino bool operator ==(const std::normal_distribution<_RealType> & __d1,const std::normal_distribution<_RealType> & __d2)1713e4b17023SJohn Marino operator==(const std::normal_distribution<_RealType>& __d1, 1714e4b17023SJohn Marino const std::normal_distribution<_RealType>& __d2) 1715e4b17023SJohn Marino { 1716e4b17023SJohn Marino if (__d1._M_param == __d2._M_param 1717e4b17023SJohn Marino && __d1._M_saved_available == __d2._M_saved_available) 1718e4b17023SJohn Marino { 1719e4b17023SJohn Marino if (__d1._M_saved_available 1720e4b17023SJohn Marino && __d1._M_saved == __d2._M_saved) 1721e4b17023SJohn Marino return true; 1722e4b17023SJohn Marino else if(!__d1._M_saved_available) 1723e4b17023SJohn Marino return true; 1724e4b17023SJohn Marino else 1725e4b17023SJohn Marino return false; 1726e4b17023SJohn Marino } 1727e4b17023SJohn Marino else 1728e4b17023SJohn Marino return false; 1729e4b17023SJohn Marino } 1730e4b17023SJohn Marino 1731e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1732e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const normal_distribution<_RealType> & __x)1733e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1734e4b17023SJohn Marino const normal_distribution<_RealType>& __x) 1735e4b17023SJohn Marino { 1736e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1737e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1738e4b17023SJohn Marino 1739e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1740e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1741e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1742e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1743e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1744e4b17023SJohn Marino __os.fill(__space); 1745e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1746e4b17023SJohn Marino 1747e4b17023SJohn Marino __os << __x.mean() << __space << __x.stddev() 1748e4b17023SJohn Marino << __space << __x._M_saved_available; 1749e4b17023SJohn Marino if (__x._M_saved_available) 1750e4b17023SJohn Marino __os << __space << __x._M_saved; 1751e4b17023SJohn Marino 1752e4b17023SJohn Marino __os.flags(__flags); 1753e4b17023SJohn Marino __os.fill(__fill); 1754e4b17023SJohn Marino __os.precision(__precision); 1755e4b17023SJohn Marino return __os; 1756e4b17023SJohn Marino } 1757e4b17023SJohn Marino 1758e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1759e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,normal_distribution<_RealType> & __x)1760e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1761e4b17023SJohn Marino normal_distribution<_RealType>& __x) 1762e4b17023SJohn Marino { 1763e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1764e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1765e4b17023SJohn Marino 1766e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1767e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1768e4b17023SJohn Marino 1769e4b17023SJohn Marino double __mean, __stddev; 1770e4b17023SJohn Marino __is >> __mean >> __stddev 1771e4b17023SJohn Marino >> __x._M_saved_available; 1772e4b17023SJohn Marino if (__x._M_saved_available) 1773e4b17023SJohn Marino __is >> __x._M_saved; 1774e4b17023SJohn Marino __x.param(typename normal_distribution<_RealType>:: 1775e4b17023SJohn Marino param_type(__mean, __stddev)); 1776e4b17023SJohn Marino 1777e4b17023SJohn Marino __is.flags(__flags); 1778e4b17023SJohn Marino return __is; 1779e4b17023SJohn Marino } 1780e4b17023SJohn Marino 1781e4b17023SJohn Marino 1782e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1783e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const lognormal_distribution<_RealType> & __x)1784e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1785e4b17023SJohn Marino const lognormal_distribution<_RealType>& __x) 1786e4b17023SJohn Marino { 1787e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1788e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1789e4b17023SJohn Marino 1790e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1791e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1792e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1793e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1794e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1795e4b17023SJohn Marino __os.fill(__space); 1796e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1797e4b17023SJohn Marino 1798e4b17023SJohn Marino __os << __x.m() << __space << __x.s() 1799e4b17023SJohn Marino << __space << __x._M_nd; 1800e4b17023SJohn Marino 1801e4b17023SJohn Marino __os.flags(__flags); 1802e4b17023SJohn Marino __os.fill(__fill); 1803e4b17023SJohn Marino __os.precision(__precision); 1804e4b17023SJohn Marino return __os; 1805e4b17023SJohn Marino } 1806e4b17023SJohn Marino 1807e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1808e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,lognormal_distribution<_RealType> & __x)1809e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1810e4b17023SJohn Marino lognormal_distribution<_RealType>& __x) 1811e4b17023SJohn Marino { 1812e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1813e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1814e4b17023SJohn Marino 1815e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1816e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1817e4b17023SJohn Marino 1818e4b17023SJohn Marino _RealType __m, __s; 1819e4b17023SJohn Marino __is >> __m >> __s >> __x._M_nd; 1820e4b17023SJohn Marino __x.param(typename lognormal_distribution<_RealType>:: 1821e4b17023SJohn Marino param_type(__m, __s)); 1822e4b17023SJohn Marino 1823e4b17023SJohn Marino __is.flags(__flags); 1824e4b17023SJohn Marino return __is; 1825e4b17023SJohn Marino } 1826e4b17023SJohn Marino 1827e4b17023SJohn Marino 1828e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1829e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const chi_squared_distribution<_RealType> & __x)1830e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1831e4b17023SJohn Marino const chi_squared_distribution<_RealType>& __x) 1832e4b17023SJohn Marino { 1833e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1834e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1835e4b17023SJohn Marino 1836e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1837e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1838e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1839e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1840e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1841e4b17023SJohn Marino __os.fill(__space); 1842e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1843e4b17023SJohn Marino 1844e4b17023SJohn Marino __os << __x.n() << __space << __x._M_gd; 1845e4b17023SJohn Marino 1846e4b17023SJohn Marino __os.flags(__flags); 1847e4b17023SJohn Marino __os.fill(__fill); 1848e4b17023SJohn Marino __os.precision(__precision); 1849e4b17023SJohn Marino return __os; 1850e4b17023SJohn Marino } 1851e4b17023SJohn Marino 1852e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1853e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,chi_squared_distribution<_RealType> & __x)1854e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1855e4b17023SJohn Marino chi_squared_distribution<_RealType>& __x) 1856e4b17023SJohn Marino { 1857e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1858e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1859e4b17023SJohn Marino 1860e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1861e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1862e4b17023SJohn Marino 1863e4b17023SJohn Marino _RealType __n; 1864e4b17023SJohn Marino __is >> __n >> __x._M_gd; 1865e4b17023SJohn Marino __x.param(typename chi_squared_distribution<_RealType>:: 1866e4b17023SJohn Marino param_type(__n)); 1867e4b17023SJohn Marino 1868e4b17023SJohn Marino __is.flags(__flags); 1869e4b17023SJohn Marino return __is; 1870e4b17023SJohn Marino } 1871e4b17023SJohn Marino 1872e4b17023SJohn Marino 1873e4b17023SJohn Marino template<typename _RealType> 1874e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 1875e4b17023SJohn Marino typename cauchy_distribution<_RealType>::result_type 1876e4b17023SJohn Marino cauchy_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)1877e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 1878e4b17023SJohn Marino const param_type& __p) 1879e4b17023SJohn Marino { 1880e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 1881e4b17023SJohn Marino __aurng(__urng); 1882e4b17023SJohn Marino _RealType __u; 1883e4b17023SJohn Marino do 1884e4b17023SJohn Marino __u = __aurng(); 1885e4b17023SJohn Marino while (__u == 0.5); 1886e4b17023SJohn Marino 1887e4b17023SJohn Marino const _RealType __pi = 3.1415926535897932384626433832795029L; 1888e4b17023SJohn Marino return __p.a() + __p.b() * std::tan(__pi * __u); 1889e4b17023SJohn Marino } 1890e4b17023SJohn Marino 1891e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1892e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const cauchy_distribution<_RealType> & __x)1893e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1894e4b17023SJohn Marino const cauchy_distribution<_RealType>& __x) 1895e4b17023SJohn Marino { 1896e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1897e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1898e4b17023SJohn Marino 1899e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1900e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1901e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1902e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1903e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1904e4b17023SJohn Marino __os.fill(__space); 1905e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1906e4b17023SJohn Marino 1907e4b17023SJohn Marino __os << __x.a() << __space << __x.b(); 1908e4b17023SJohn Marino 1909e4b17023SJohn Marino __os.flags(__flags); 1910e4b17023SJohn Marino __os.fill(__fill); 1911e4b17023SJohn Marino __os.precision(__precision); 1912e4b17023SJohn Marino return __os; 1913e4b17023SJohn Marino } 1914e4b17023SJohn Marino 1915e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1916e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,cauchy_distribution<_RealType> & __x)1917e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1918e4b17023SJohn Marino cauchy_distribution<_RealType>& __x) 1919e4b17023SJohn Marino { 1920e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1921e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1922e4b17023SJohn Marino 1923e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1924e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1925e4b17023SJohn Marino 1926e4b17023SJohn Marino _RealType __a, __b; 1927e4b17023SJohn Marino __is >> __a >> __b; 1928e4b17023SJohn Marino __x.param(typename cauchy_distribution<_RealType>:: 1929e4b17023SJohn Marino param_type(__a, __b)); 1930e4b17023SJohn Marino 1931e4b17023SJohn Marino __is.flags(__flags); 1932e4b17023SJohn Marino return __is; 1933e4b17023SJohn Marino } 1934e4b17023SJohn Marino 1935e4b17023SJohn Marino 1936e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1937e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const fisher_f_distribution<_RealType> & __x)1938e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1939e4b17023SJohn Marino const fisher_f_distribution<_RealType>& __x) 1940e4b17023SJohn Marino { 1941e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1942e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1943e4b17023SJohn Marino 1944e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1945e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1946e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1947e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1948e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1949e4b17023SJohn Marino __os.fill(__space); 1950e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1951e4b17023SJohn Marino 1952e4b17023SJohn Marino __os << __x.m() << __space << __x.n() 1953e4b17023SJohn Marino << __space << __x._M_gd_x << __space << __x._M_gd_y; 1954e4b17023SJohn Marino 1955e4b17023SJohn Marino __os.flags(__flags); 1956e4b17023SJohn Marino __os.fill(__fill); 1957e4b17023SJohn Marino __os.precision(__precision); 1958e4b17023SJohn Marino return __os; 1959e4b17023SJohn Marino } 1960e4b17023SJohn Marino 1961e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1962e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,fisher_f_distribution<_RealType> & __x)1963e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 1964e4b17023SJohn Marino fisher_f_distribution<_RealType>& __x) 1965e4b17023SJohn Marino { 1966e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 1967e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 1968e4b17023SJohn Marino 1969e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 1970e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 1971e4b17023SJohn Marino 1972e4b17023SJohn Marino _RealType __m, __n; 1973e4b17023SJohn Marino __is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y; 1974e4b17023SJohn Marino __x.param(typename fisher_f_distribution<_RealType>:: 1975e4b17023SJohn Marino param_type(__m, __n)); 1976e4b17023SJohn Marino 1977e4b17023SJohn Marino __is.flags(__flags); 1978e4b17023SJohn Marino return __is; 1979e4b17023SJohn Marino } 1980e4b17023SJohn Marino 1981e4b17023SJohn Marino 1982e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 1983e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const student_t_distribution<_RealType> & __x)1984e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 1985e4b17023SJohn Marino const student_t_distribution<_RealType>& __x) 1986e4b17023SJohn Marino { 1987e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 1988e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 1989e4b17023SJohn Marino 1990e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 1991e4b17023SJohn Marino const _CharT __fill = __os.fill(); 1992e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 1993e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 1994e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 1995e4b17023SJohn Marino __os.fill(__space); 1996e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 1997e4b17023SJohn Marino 1998e4b17023SJohn Marino __os << __x.n() << __space << __x._M_nd << __space << __x._M_gd; 1999e4b17023SJohn Marino 2000e4b17023SJohn Marino __os.flags(__flags); 2001e4b17023SJohn Marino __os.fill(__fill); 2002e4b17023SJohn Marino __os.precision(__precision); 2003e4b17023SJohn Marino return __os; 2004e4b17023SJohn Marino } 2005e4b17023SJohn Marino 2006e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2007e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,student_t_distribution<_RealType> & __x)2008e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2009e4b17023SJohn Marino student_t_distribution<_RealType>& __x) 2010e4b17023SJohn Marino { 2011e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2012e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2013e4b17023SJohn Marino 2014e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2015e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2016e4b17023SJohn Marino 2017e4b17023SJohn Marino _RealType __n; 2018e4b17023SJohn Marino __is >> __n >> __x._M_nd >> __x._M_gd; 2019e4b17023SJohn Marino __x.param(typename student_t_distribution<_RealType>::param_type(__n)); 2020e4b17023SJohn Marino 2021e4b17023SJohn Marino __is.flags(__flags); 2022e4b17023SJohn Marino return __is; 2023e4b17023SJohn Marino } 2024e4b17023SJohn Marino 2025e4b17023SJohn Marino 2026e4b17023SJohn Marino template<typename _RealType> 2027e4b17023SJohn Marino void 2028e4b17023SJohn Marino gamma_distribution<_RealType>::param_type:: _M_initialize()2029e4b17023SJohn Marino _M_initialize() 2030e4b17023SJohn Marino { 2031e4b17023SJohn Marino _M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha; 2032e4b17023SJohn Marino 2033e4b17023SJohn Marino const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0); 2034e4b17023SJohn Marino _M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1); 2035e4b17023SJohn Marino } 2036e4b17023SJohn Marino 2037e4b17023SJohn Marino /** 2038e4b17023SJohn Marino * Marsaglia, G. and Tsang, W. W. 2039e4b17023SJohn Marino * "A Simple Method for Generating Gamma Variables" 2040e4b17023SJohn Marino * ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000. 2041e4b17023SJohn Marino */ 2042e4b17023SJohn Marino template<typename _RealType> 2043e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2044e4b17023SJohn Marino typename gamma_distribution<_RealType>::result_type 2045e4b17023SJohn Marino gamma_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2046e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2047e4b17023SJohn Marino const param_type& __param) 2048e4b17023SJohn Marino { 2049e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2050e4b17023SJohn Marino __aurng(__urng); 2051e4b17023SJohn Marino 2052e4b17023SJohn Marino result_type __u, __v, __n; 2053e4b17023SJohn Marino const result_type __a1 = (__param._M_malpha 2054e4b17023SJohn Marino - _RealType(1.0) / _RealType(3.0)); 2055e4b17023SJohn Marino 2056e4b17023SJohn Marino do 2057e4b17023SJohn Marino { 2058e4b17023SJohn Marino do 2059e4b17023SJohn Marino { 2060e4b17023SJohn Marino __n = _M_nd(__urng); 2061e4b17023SJohn Marino __v = result_type(1.0) + __param._M_a2 * __n; 2062e4b17023SJohn Marino } 2063e4b17023SJohn Marino while (__v <= 0.0); 2064e4b17023SJohn Marino 2065e4b17023SJohn Marino __v = __v * __v * __v; 2066e4b17023SJohn Marino __u = __aurng(); 2067e4b17023SJohn Marino } 2068e4b17023SJohn Marino while (__u > result_type(1.0) - 0.331 * __n * __n * __n * __n 2069e4b17023SJohn Marino && (std::log(__u) > (0.5 * __n * __n + __a1 2070e4b17023SJohn Marino * (1.0 - __v + std::log(__v))))); 2071e4b17023SJohn Marino 2072e4b17023SJohn Marino if (__param.alpha() == __param._M_malpha) 2073e4b17023SJohn Marino return __a1 * __v * __param.beta(); 2074e4b17023SJohn Marino else 2075e4b17023SJohn Marino { 2076e4b17023SJohn Marino do 2077e4b17023SJohn Marino __u = __aurng(); 2078e4b17023SJohn Marino while (__u == 0.0); 2079e4b17023SJohn Marino 2080e4b17023SJohn Marino return (std::pow(__u, result_type(1.0) / __param.alpha()) 2081e4b17023SJohn Marino * __a1 * __v * __param.beta()); 2082e4b17023SJohn Marino } 2083e4b17023SJohn Marino } 2084e4b17023SJohn Marino 2085e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2086e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const gamma_distribution<_RealType> & __x)2087e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2088e4b17023SJohn Marino const gamma_distribution<_RealType>& __x) 2089e4b17023SJohn Marino { 2090e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2091e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2092e4b17023SJohn Marino 2093e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2094e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2095e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2096e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2097e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2098e4b17023SJohn Marino __os.fill(__space); 2099e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 2100e4b17023SJohn Marino 2101e4b17023SJohn Marino __os << __x.alpha() << __space << __x.beta() 2102e4b17023SJohn Marino << __space << __x._M_nd; 2103e4b17023SJohn Marino 2104e4b17023SJohn Marino __os.flags(__flags); 2105e4b17023SJohn Marino __os.fill(__fill); 2106e4b17023SJohn Marino __os.precision(__precision); 2107e4b17023SJohn Marino return __os; 2108e4b17023SJohn Marino } 2109e4b17023SJohn Marino 2110e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2111e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,gamma_distribution<_RealType> & __x)2112e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2113e4b17023SJohn Marino gamma_distribution<_RealType>& __x) 2114e4b17023SJohn Marino { 2115e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2116e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2117e4b17023SJohn Marino 2118e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2119e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2120e4b17023SJohn Marino 2121e4b17023SJohn Marino _RealType __alpha_val, __beta_val; 2122e4b17023SJohn Marino __is >> __alpha_val >> __beta_val >> __x._M_nd; 2123e4b17023SJohn Marino __x.param(typename gamma_distribution<_RealType>:: 2124e4b17023SJohn Marino param_type(__alpha_val, __beta_val)); 2125e4b17023SJohn Marino 2126e4b17023SJohn Marino __is.flags(__flags); 2127e4b17023SJohn Marino return __is; 2128e4b17023SJohn Marino } 2129e4b17023SJohn Marino 2130e4b17023SJohn Marino 2131e4b17023SJohn Marino template<typename _RealType> 2132e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2133e4b17023SJohn Marino typename weibull_distribution<_RealType>::result_type 2134e4b17023SJohn Marino weibull_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)2135e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2136e4b17023SJohn Marino const param_type& __p) 2137e4b17023SJohn Marino { 2138e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2139e4b17023SJohn Marino __aurng(__urng); 21405ce9237cSJohn Marino return __p.b() * std::pow(-std::log(result_type(1) - __aurng()), 2141e4b17023SJohn Marino result_type(1) / __p.a()); 2142e4b17023SJohn Marino } 2143e4b17023SJohn Marino 2144e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2145e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const weibull_distribution<_RealType> & __x)2146e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2147e4b17023SJohn Marino const weibull_distribution<_RealType>& __x) 2148e4b17023SJohn Marino { 2149e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2150e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2151e4b17023SJohn Marino 2152e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2153e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2154e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2155e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2156e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2157e4b17023SJohn Marino __os.fill(__space); 2158e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 2159e4b17023SJohn Marino 2160e4b17023SJohn Marino __os << __x.a() << __space << __x.b(); 2161e4b17023SJohn Marino 2162e4b17023SJohn Marino __os.flags(__flags); 2163e4b17023SJohn Marino __os.fill(__fill); 2164e4b17023SJohn Marino __os.precision(__precision); 2165e4b17023SJohn Marino return __os; 2166e4b17023SJohn Marino } 2167e4b17023SJohn Marino 2168e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2169e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,weibull_distribution<_RealType> & __x)2170e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2171e4b17023SJohn Marino weibull_distribution<_RealType>& __x) 2172e4b17023SJohn Marino { 2173e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2174e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2175e4b17023SJohn Marino 2176e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2177e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2178e4b17023SJohn Marino 2179e4b17023SJohn Marino _RealType __a, __b; 2180e4b17023SJohn Marino __is >> __a >> __b; 2181e4b17023SJohn Marino __x.param(typename weibull_distribution<_RealType>:: 2182e4b17023SJohn Marino param_type(__a, __b)); 2183e4b17023SJohn Marino 2184e4b17023SJohn Marino __is.flags(__flags); 2185e4b17023SJohn Marino return __is; 2186e4b17023SJohn Marino } 2187e4b17023SJohn Marino 2188e4b17023SJohn Marino 2189e4b17023SJohn Marino template<typename _RealType> 2190e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2191e4b17023SJohn Marino typename extreme_value_distribution<_RealType>::result_type 2192e4b17023SJohn Marino extreme_value_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __p)2193e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2194e4b17023SJohn Marino const param_type& __p) 2195e4b17023SJohn Marino { 2196e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, result_type> 2197e4b17023SJohn Marino __aurng(__urng); 21985ce9237cSJohn Marino return __p.a() - __p.b() * std::log(-std::log(result_type(1) 21995ce9237cSJohn Marino - __aurng())); 2200e4b17023SJohn Marino } 2201e4b17023SJohn Marino 2202e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2203e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const extreme_value_distribution<_RealType> & __x)2204e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2205e4b17023SJohn Marino const extreme_value_distribution<_RealType>& __x) 2206e4b17023SJohn Marino { 2207e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2208e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2209e4b17023SJohn Marino 2210e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2211e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2212e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2213e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2214e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2215e4b17023SJohn Marino __os.fill(__space); 2216e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 2217e4b17023SJohn Marino 2218e4b17023SJohn Marino __os << __x.a() << __space << __x.b(); 2219e4b17023SJohn Marino 2220e4b17023SJohn Marino __os.flags(__flags); 2221e4b17023SJohn Marino __os.fill(__fill); 2222e4b17023SJohn Marino __os.precision(__precision); 2223e4b17023SJohn Marino return __os; 2224e4b17023SJohn Marino } 2225e4b17023SJohn Marino 2226e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2227e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,extreme_value_distribution<_RealType> & __x)2228e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2229e4b17023SJohn Marino extreme_value_distribution<_RealType>& __x) 2230e4b17023SJohn Marino { 2231e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2232e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2233e4b17023SJohn Marino 2234e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2235e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2236e4b17023SJohn Marino 2237e4b17023SJohn Marino _RealType __a, __b; 2238e4b17023SJohn Marino __is >> __a >> __b; 2239e4b17023SJohn Marino __x.param(typename extreme_value_distribution<_RealType>:: 2240e4b17023SJohn Marino param_type(__a, __b)); 2241e4b17023SJohn Marino 2242e4b17023SJohn Marino __is.flags(__flags); 2243e4b17023SJohn Marino return __is; 2244e4b17023SJohn Marino } 2245e4b17023SJohn Marino 2246e4b17023SJohn Marino 2247e4b17023SJohn Marino template<typename _IntType> 2248e4b17023SJohn Marino void 2249e4b17023SJohn Marino discrete_distribution<_IntType>::param_type:: _M_initialize()2250e4b17023SJohn Marino _M_initialize() 2251e4b17023SJohn Marino { 2252e4b17023SJohn Marino if (_M_prob.size() < 2) 2253e4b17023SJohn Marino { 2254e4b17023SJohn Marino _M_prob.clear(); 2255e4b17023SJohn Marino return; 2256e4b17023SJohn Marino } 2257e4b17023SJohn Marino 2258e4b17023SJohn Marino const double __sum = std::accumulate(_M_prob.begin(), 2259e4b17023SJohn Marino _M_prob.end(), 0.0); 2260e4b17023SJohn Marino // Now normalize the probabilites. 2261e4b17023SJohn Marino __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(), 2262e4b17023SJohn Marino std::bind2nd(std::divides<double>(), __sum)); 2263e4b17023SJohn Marino // Accumulate partial sums. 2264e4b17023SJohn Marino _M_cp.reserve(_M_prob.size()); 2265e4b17023SJohn Marino std::partial_sum(_M_prob.begin(), _M_prob.end(), 2266e4b17023SJohn Marino std::back_inserter(_M_cp)); 2267e4b17023SJohn Marino // Make sure the last cumulative probability is one. 2268e4b17023SJohn Marino _M_cp[_M_cp.size() - 1] = 1.0; 2269e4b17023SJohn Marino } 2270e4b17023SJohn Marino 2271e4b17023SJohn Marino template<typename _IntType> 2272e4b17023SJohn Marino template<typename _Func> 2273e4b17023SJohn Marino discrete_distribution<_IntType>::param_type:: param_type(size_t __nw,double __xmin,double __xmax,_Func __fw)2274e4b17023SJohn Marino param_type(size_t __nw, double __xmin, double __xmax, _Func __fw) 2275e4b17023SJohn Marino : _M_prob(), _M_cp() 2276e4b17023SJohn Marino { 2277e4b17023SJohn Marino const size_t __n = __nw == 0 ? 1 : __nw; 2278e4b17023SJohn Marino const double __delta = (__xmax - __xmin) / __n; 2279e4b17023SJohn Marino 2280e4b17023SJohn Marino _M_prob.reserve(__n); 2281e4b17023SJohn Marino for (size_t __k = 0; __k < __nw; ++__k) 2282e4b17023SJohn Marino _M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta)); 2283e4b17023SJohn Marino 2284e4b17023SJohn Marino _M_initialize(); 2285e4b17023SJohn Marino } 2286e4b17023SJohn Marino 2287e4b17023SJohn Marino template<typename _IntType> 2288e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2289e4b17023SJohn Marino typename discrete_distribution<_IntType>::result_type 2290e4b17023SJohn Marino discrete_distribution<_IntType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2291e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2292e4b17023SJohn Marino const param_type& __param) 2293e4b17023SJohn Marino { 2294e4b17023SJohn Marino if (__param._M_cp.empty()) 2295e4b17023SJohn Marino return result_type(0); 2296e4b17023SJohn Marino 2297e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2298e4b17023SJohn Marino __aurng(__urng); 2299e4b17023SJohn Marino 2300e4b17023SJohn Marino const double __p = __aurng(); 2301e4b17023SJohn Marino auto __pos = std::lower_bound(__param._M_cp.begin(), 2302e4b17023SJohn Marino __param._M_cp.end(), __p); 2303e4b17023SJohn Marino 2304e4b17023SJohn Marino return __pos - __param._M_cp.begin(); 2305e4b17023SJohn Marino } 2306e4b17023SJohn Marino 2307e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 2308e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const discrete_distribution<_IntType> & __x)2309e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2310e4b17023SJohn Marino const discrete_distribution<_IntType>& __x) 2311e4b17023SJohn Marino { 2312e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2313e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2314e4b17023SJohn Marino 2315e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2316e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2317e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2318e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2319e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2320e4b17023SJohn Marino __os.fill(__space); 2321e4b17023SJohn Marino __os.precision(std::numeric_limits<double>::max_digits10); 2322e4b17023SJohn Marino 2323e4b17023SJohn Marino std::vector<double> __prob = __x.probabilities(); 2324e4b17023SJohn Marino __os << __prob.size(); 2325e4b17023SJohn Marino for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit) 2326e4b17023SJohn Marino __os << __space << *__dit; 2327e4b17023SJohn Marino 2328e4b17023SJohn Marino __os.flags(__flags); 2329e4b17023SJohn Marino __os.fill(__fill); 2330e4b17023SJohn Marino __os.precision(__precision); 2331e4b17023SJohn Marino return __os; 2332e4b17023SJohn Marino } 2333e4b17023SJohn Marino 2334e4b17023SJohn Marino template<typename _IntType, typename _CharT, typename _Traits> 2335e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,discrete_distribution<_IntType> & __x)2336e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2337e4b17023SJohn Marino discrete_distribution<_IntType>& __x) 2338e4b17023SJohn Marino { 2339e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2340e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2341e4b17023SJohn Marino 2342e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2343e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2344e4b17023SJohn Marino 2345e4b17023SJohn Marino size_t __n; 2346e4b17023SJohn Marino __is >> __n; 2347e4b17023SJohn Marino 2348e4b17023SJohn Marino std::vector<double> __prob_vec; 2349e4b17023SJohn Marino __prob_vec.reserve(__n); 2350e4b17023SJohn Marino for (; __n != 0; --__n) 2351e4b17023SJohn Marino { 2352e4b17023SJohn Marino double __prob; 2353e4b17023SJohn Marino __is >> __prob; 2354e4b17023SJohn Marino __prob_vec.push_back(__prob); 2355e4b17023SJohn Marino } 2356e4b17023SJohn Marino 2357e4b17023SJohn Marino __x.param(typename discrete_distribution<_IntType>:: 2358e4b17023SJohn Marino param_type(__prob_vec.begin(), __prob_vec.end())); 2359e4b17023SJohn Marino 2360e4b17023SJohn Marino __is.flags(__flags); 2361e4b17023SJohn Marino return __is; 2362e4b17023SJohn Marino } 2363e4b17023SJohn Marino 2364e4b17023SJohn Marino 2365e4b17023SJohn Marino template<typename _RealType> 2366e4b17023SJohn Marino void 2367e4b17023SJohn Marino piecewise_constant_distribution<_RealType>::param_type:: _M_initialize()2368e4b17023SJohn Marino _M_initialize() 2369e4b17023SJohn Marino { 2370e4b17023SJohn Marino if (_M_int.size() < 2 2371e4b17023SJohn Marino || (_M_int.size() == 2 2372e4b17023SJohn Marino && _M_int[0] == _RealType(0) 2373e4b17023SJohn Marino && _M_int[1] == _RealType(1))) 2374e4b17023SJohn Marino { 2375e4b17023SJohn Marino _M_int.clear(); 2376e4b17023SJohn Marino _M_den.clear(); 2377e4b17023SJohn Marino return; 2378e4b17023SJohn Marino } 2379e4b17023SJohn Marino 2380e4b17023SJohn Marino const double __sum = std::accumulate(_M_den.begin(), 2381e4b17023SJohn Marino _M_den.end(), 0.0); 2382e4b17023SJohn Marino 2383e4b17023SJohn Marino __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(), 2384e4b17023SJohn Marino std::bind2nd(std::divides<double>(), __sum)); 2385e4b17023SJohn Marino 2386e4b17023SJohn Marino _M_cp.reserve(_M_den.size()); 2387e4b17023SJohn Marino std::partial_sum(_M_den.begin(), _M_den.end(), 2388e4b17023SJohn Marino std::back_inserter(_M_cp)); 2389e4b17023SJohn Marino 2390e4b17023SJohn Marino // Make sure the last cumulative probability is one. 2391e4b17023SJohn Marino _M_cp[_M_cp.size() - 1] = 1.0; 2392e4b17023SJohn Marino 2393e4b17023SJohn Marino for (size_t __k = 0; __k < _M_den.size(); ++__k) 2394e4b17023SJohn Marino _M_den[__k] /= _M_int[__k + 1] - _M_int[__k]; 2395e4b17023SJohn Marino } 2396e4b17023SJohn Marino 2397e4b17023SJohn Marino template<typename _RealType> 2398e4b17023SJohn Marino template<typename _InputIteratorB, typename _InputIteratorW> 2399e4b17023SJohn Marino piecewise_constant_distribution<_RealType>::param_type:: param_type(_InputIteratorB __bbegin,_InputIteratorB __bend,_InputIteratorW __wbegin)2400e4b17023SJohn Marino param_type(_InputIteratorB __bbegin, 2401e4b17023SJohn Marino _InputIteratorB __bend, 2402e4b17023SJohn Marino _InputIteratorW __wbegin) 2403e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp() 2404e4b17023SJohn Marino { 2405e4b17023SJohn Marino if (__bbegin != __bend) 2406e4b17023SJohn Marino { 2407e4b17023SJohn Marino for (;;) 2408e4b17023SJohn Marino { 2409e4b17023SJohn Marino _M_int.push_back(*__bbegin); 2410e4b17023SJohn Marino ++__bbegin; 2411e4b17023SJohn Marino if (__bbegin == __bend) 2412e4b17023SJohn Marino break; 2413e4b17023SJohn Marino 2414e4b17023SJohn Marino _M_den.push_back(*__wbegin); 2415e4b17023SJohn Marino ++__wbegin; 2416e4b17023SJohn Marino } 2417e4b17023SJohn Marino } 2418e4b17023SJohn Marino 2419e4b17023SJohn Marino _M_initialize(); 2420e4b17023SJohn Marino } 2421e4b17023SJohn Marino 2422e4b17023SJohn Marino template<typename _RealType> 2423e4b17023SJohn Marino template<typename _Func> 2424e4b17023SJohn Marino piecewise_constant_distribution<_RealType>::param_type:: param_type(initializer_list<_RealType> __bl,_Func __fw)2425e4b17023SJohn Marino param_type(initializer_list<_RealType> __bl, _Func __fw) 2426e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp() 2427e4b17023SJohn Marino { 2428e4b17023SJohn Marino _M_int.reserve(__bl.size()); 2429e4b17023SJohn Marino for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 2430e4b17023SJohn Marino _M_int.push_back(*__biter); 2431e4b17023SJohn Marino 2432e4b17023SJohn Marino _M_den.reserve(_M_int.size() - 1); 2433e4b17023SJohn Marino for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 2434e4b17023SJohn Marino _M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k]))); 2435e4b17023SJohn Marino 2436e4b17023SJohn Marino _M_initialize(); 2437e4b17023SJohn Marino } 2438e4b17023SJohn Marino 2439e4b17023SJohn Marino template<typename _RealType> 2440e4b17023SJohn Marino template<typename _Func> 2441e4b17023SJohn Marino piecewise_constant_distribution<_RealType>::param_type:: param_type(size_t __nw,_RealType __xmin,_RealType __xmax,_Func __fw)2442e4b17023SJohn Marino param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 2443e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp() 2444e4b17023SJohn Marino { 2445e4b17023SJohn Marino const size_t __n = __nw == 0 ? 1 : __nw; 2446e4b17023SJohn Marino const _RealType __delta = (__xmax - __xmin) / __n; 2447e4b17023SJohn Marino 2448e4b17023SJohn Marino _M_int.reserve(__n + 1); 2449e4b17023SJohn Marino for (size_t __k = 0; __k <= __nw; ++__k) 2450e4b17023SJohn Marino _M_int.push_back(__xmin + __k * __delta); 2451e4b17023SJohn Marino 2452e4b17023SJohn Marino _M_den.reserve(__n); 2453e4b17023SJohn Marino for (size_t __k = 0; __k < __nw; ++__k) 2454e4b17023SJohn Marino _M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta)); 2455e4b17023SJohn Marino 2456e4b17023SJohn Marino _M_initialize(); 2457e4b17023SJohn Marino } 2458e4b17023SJohn Marino 2459e4b17023SJohn Marino template<typename _RealType> 2460e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2461e4b17023SJohn Marino typename piecewise_constant_distribution<_RealType>::result_type 2462e4b17023SJohn Marino piecewise_constant_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2463e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2464e4b17023SJohn Marino const param_type& __param) 2465e4b17023SJohn Marino { 2466e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2467e4b17023SJohn Marino __aurng(__urng); 2468e4b17023SJohn Marino 2469e4b17023SJohn Marino const double __p = __aurng(); 2470e4b17023SJohn Marino if (__param._M_cp.empty()) 2471e4b17023SJohn Marino return __p; 2472e4b17023SJohn Marino 2473e4b17023SJohn Marino auto __pos = std::lower_bound(__param._M_cp.begin(), 2474e4b17023SJohn Marino __param._M_cp.end(), __p); 2475e4b17023SJohn Marino const size_t __i = __pos - __param._M_cp.begin(); 2476e4b17023SJohn Marino 2477e4b17023SJohn Marino const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2478e4b17023SJohn Marino 2479e4b17023SJohn Marino return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i]; 2480e4b17023SJohn Marino } 2481e4b17023SJohn Marino 2482e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2483e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const piecewise_constant_distribution<_RealType> & __x)2484e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2485e4b17023SJohn Marino const piecewise_constant_distribution<_RealType>& __x) 2486e4b17023SJohn Marino { 2487e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2488e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2489e4b17023SJohn Marino 2490e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2491e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2492e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2493e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2494e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2495e4b17023SJohn Marino __os.fill(__space); 2496e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 2497e4b17023SJohn Marino 2498e4b17023SJohn Marino std::vector<_RealType> __int = __x.intervals(); 2499e4b17023SJohn Marino __os << __int.size() - 1; 2500e4b17023SJohn Marino 2501e4b17023SJohn Marino for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 2502e4b17023SJohn Marino __os << __space << *__xit; 2503e4b17023SJohn Marino 2504e4b17023SJohn Marino std::vector<double> __den = __x.densities(); 2505e4b17023SJohn Marino for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 2506e4b17023SJohn Marino __os << __space << *__dit; 2507e4b17023SJohn Marino 2508e4b17023SJohn Marino __os.flags(__flags); 2509e4b17023SJohn Marino __os.fill(__fill); 2510e4b17023SJohn Marino __os.precision(__precision); 2511e4b17023SJohn Marino return __os; 2512e4b17023SJohn Marino } 2513e4b17023SJohn Marino 2514e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2515e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,piecewise_constant_distribution<_RealType> & __x)2516e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2517e4b17023SJohn Marino piecewise_constant_distribution<_RealType>& __x) 2518e4b17023SJohn Marino { 2519e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2520e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2521e4b17023SJohn Marino 2522e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2523e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2524e4b17023SJohn Marino 2525e4b17023SJohn Marino size_t __n; 2526e4b17023SJohn Marino __is >> __n; 2527e4b17023SJohn Marino 2528e4b17023SJohn Marino std::vector<_RealType> __int_vec; 2529e4b17023SJohn Marino __int_vec.reserve(__n + 1); 2530e4b17023SJohn Marino for (size_t __i = 0; __i <= __n; ++__i) 2531e4b17023SJohn Marino { 2532e4b17023SJohn Marino _RealType __int; 2533e4b17023SJohn Marino __is >> __int; 2534e4b17023SJohn Marino __int_vec.push_back(__int); 2535e4b17023SJohn Marino } 2536e4b17023SJohn Marino 2537e4b17023SJohn Marino std::vector<double> __den_vec; 2538e4b17023SJohn Marino __den_vec.reserve(__n); 2539e4b17023SJohn Marino for (size_t __i = 0; __i < __n; ++__i) 2540e4b17023SJohn Marino { 2541e4b17023SJohn Marino double __den; 2542e4b17023SJohn Marino __is >> __den; 2543e4b17023SJohn Marino __den_vec.push_back(__den); 2544e4b17023SJohn Marino } 2545e4b17023SJohn Marino 2546e4b17023SJohn Marino __x.param(typename piecewise_constant_distribution<_RealType>:: 2547e4b17023SJohn Marino param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 2548e4b17023SJohn Marino 2549e4b17023SJohn Marino __is.flags(__flags); 2550e4b17023SJohn Marino return __is; 2551e4b17023SJohn Marino } 2552e4b17023SJohn Marino 2553e4b17023SJohn Marino 2554e4b17023SJohn Marino template<typename _RealType> 2555e4b17023SJohn Marino void 2556e4b17023SJohn Marino piecewise_linear_distribution<_RealType>::param_type:: _M_initialize()2557e4b17023SJohn Marino _M_initialize() 2558e4b17023SJohn Marino { 2559e4b17023SJohn Marino if (_M_int.size() < 2 2560e4b17023SJohn Marino || (_M_int.size() == 2 2561e4b17023SJohn Marino && _M_int[0] == _RealType(0) 2562e4b17023SJohn Marino && _M_int[1] == _RealType(1) 2563e4b17023SJohn Marino && _M_den[0] == _M_den[1])) 2564e4b17023SJohn Marino { 2565e4b17023SJohn Marino _M_int.clear(); 2566e4b17023SJohn Marino _M_den.clear(); 2567e4b17023SJohn Marino return; 2568e4b17023SJohn Marino } 2569e4b17023SJohn Marino 2570e4b17023SJohn Marino double __sum = 0.0; 2571e4b17023SJohn Marino _M_cp.reserve(_M_int.size() - 1); 2572e4b17023SJohn Marino _M_m.reserve(_M_int.size() - 1); 2573e4b17023SJohn Marino for (size_t __k = 0; __k < _M_int.size() - 1; ++__k) 2574e4b17023SJohn Marino { 2575e4b17023SJohn Marino const _RealType __delta = _M_int[__k + 1] - _M_int[__k]; 2576e4b17023SJohn Marino __sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta; 2577e4b17023SJohn Marino _M_cp.push_back(__sum); 2578e4b17023SJohn Marino _M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta); 2579e4b17023SJohn Marino } 2580e4b17023SJohn Marino 2581e4b17023SJohn Marino // Now normalize the densities... 2582e4b17023SJohn Marino __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(), 2583e4b17023SJohn Marino std::bind2nd(std::divides<double>(), __sum)); 2584e4b17023SJohn Marino // ... and partial sums... 2585e4b17023SJohn Marino __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), 2586e4b17023SJohn Marino std::bind2nd(std::divides<double>(), __sum)); 2587e4b17023SJohn Marino // ... and slopes. 2588e4b17023SJohn Marino __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(), 2589e4b17023SJohn Marino std::bind2nd(std::divides<double>(), __sum)); 2590e4b17023SJohn Marino // Make sure the last cumulative probablility is one. 2591e4b17023SJohn Marino _M_cp[_M_cp.size() - 1] = 1.0; 2592e4b17023SJohn Marino } 2593e4b17023SJohn Marino 2594e4b17023SJohn Marino template<typename _RealType> 2595e4b17023SJohn Marino template<typename _InputIteratorB, typename _InputIteratorW> 2596e4b17023SJohn Marino piecewise_linear_distribution<_RealType>::param_type:: param_type(_InputIteratorB __bbegin,_InputIteratorB __bend,_InputIteratorW __wbegin)2597e4b17023SJohn Marino param_type(_InputIteratorB __bbegin, 2598e4b17023SJohn Marino _InputIteratorB __bend, 2599e4b17023SJohn Marino _InputIteratorW __wbegin) 2600e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp(), _M_m() 2601e4b17023SJohn Marino { 2602e4b17023SJohn Marino for (; __bbegin != __bend; ++__bbegin, ++__wbegin) 2603e4b17023SJohn Marino { 2604e4b17023SJohn Marino _M_int.push_back(*__bbegin); 2605e4b17023SJohn Marino _M_den.push_back(*__wbegin); 2606e4b17023SJohn Marino } 2607e4b17023SJohn Marino 2608e4b17023SJohn Marino _M_initialize(); 2609e4b17023SJohn Marino } 2610e4b17023SJohn Marino 2611e4b17023SJohn Marino template<typename _RealType> 2612e4b17023SJohn Marino template<typename _Func> 2613e4b17023SJohn Marino piecewise_linear_distribution<_RealType>::param_type:: param_type(initializer_list<_RealType> __bl,_Func __fw)2614e4b17023SJohn Marino param_type(initializer_list<_RealType> __bl, _Func __fw) 2615e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp(), _M_m() 2616e4b17023SJohn Marino { 2617e4b17023SJohn Marino _M_int.reserve(__bl.size()); 2618e4b17023SJohn Marino _M_den.reserve(__bl.size()); 2619e4b17023SJohn Marino for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter) 2620e4b17023SJohn Marino { 2621e4b17023SJohn Marino _M_int.push_back(*__biter); 2622e4b17023SJohn Marino _M_den.push_back(__fw(*__biter)); 2623e4b17023SJohn Marino } 2624e4b17023SJohn Marino 2625e4b17023SJohn Marino _M_initialize(); 2626e4b17023SJohn Marino } 2627e4b17023SJohn Marino 2628e4b17023SJohn Marino template<typename _RealType> 2629e4b17023SJohn Marino template<typename _Func> 2630e4b17023SJohn Marino piecewise_linear_distribution<_RealType>::param_type:: param_type(size_t __nw,_RealType __xmin,_RealType __xmax,_Func __fw)2631e4b17023SJohn Marino param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw) 2632e4b17023SJohn Marino : _M_int(), _M_den(), _M_cp(), _M_m() 2633e4b17023SJohn Marino { 2634e4b17023SJohn Marino const size_t __n = __nw == 0 ? 1 : __nw; 2635e4b17023SJohn Marino const _RealType __delta = (__xmax - __xmin) / __n; 2636e4b17023SJohn Marino 2637e4b17023SJohn Marino _M_int.reserve(__n + 1); 2638e4b17023SJohn Marino _M_den.reserve(__n + 1); 2639e4b17023SJohn Marino for (size_t __k = 0; __k <= __nw; ++__k) 2640e4b17023SJohn Marino { 2641e4b17023SJohn Marino _M_int.push_back(__xmin + __k * __delta); 2642e4b17023SJohn Marino _M_den.push_back(__fw(_M_int[__k] + __delta)); 2643e4b17023SJohn Marino } 2644e4b17023SJohn Marino 2645e4b17023SJohn Marino _M_initialize(); 2646e4b17023SJohn Marino } 2647e4b17023SJohn Marino 2648e4b17023SJohn Marino template<typename _RealType> 2649e4b17023SJohn Marino template<typename _UniformRandomNumberGenerator> 2650e4b17023SJohn Marino typename piecewise_linear_distribution<_RealType>::result_type 2651e4b17023SJohn Marino piecewise_linear_distribution<_RealType>:: operator ()(_UniformRandomNumberGenerator & __urng,const param_type & __param)2652e4b17023SJohn Marino operator()(_UniformRandomNumberGenerator& __urng, 2653e4b17023SJohn Marino const param_type& __param) 2654e4b17023SJohn Marino { 2655e4b17023SJohn Marino __detail::_Adaptor<_UniformRandomNumberGenerator, double> 2656e4b17023SJohn Marino __aurng(__urng); 2657e4b17023SJohn Marino 2658e4b17023SJohn Marino const double __p = __aurng(); 2659e4b17023SJohn Marino if (__param._M_cp.empty()) 2660e4b17023SJohn Marino return __p; 2661e4b17023SJohn Marino 2662e4b17023SJohn Marino auto __pos = std::lower_bound(__param._M_cp.begin(), 2663e4b17023SJohn Marino __param._M_cp.end(), __p); 2664e4b17023SJohn Marino const size_t __i = __pos - __param._M_cp.begin(); 2665e4b17023SJohn Marino 2666e4b17023SJohn Marino const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0; 2667e4b17023SJohn Marino 2668e4b17023SJohn Marino const double __a = 0.5 * __param._M_m[__i]; 2669e4b17023SJohn Marino const double __b = __param._M_den[__i]; 2670e4b17023SJohn Marino const double __cm = __p - __pref; 2671e4b17023SJohn Marino 2672e4b17023SJohn Marino _RealType __x = __param._M_int[__i]; 2673e4b17023SJohn Marino if (__a == 0) 2674e4b17023SJohn Marino __x += __cm / __b; 2675e4b17023SJohn Marino else 2676e4b17023SJohn Marino { 2677e4b17023SJohn Marino const double __d = __b * __b + 4.0 * __a * __cm; 2678e4b17023SJohn Marino __x += 0.5 * (std::sqrt(__d) - __b) / __a; 2679e4b17023SJohn Marino } 2680e4b17023SJohn Marino 2681e4b17023SJohn Marino return __x; 2682e4b17023SJohn Marino } 2683e4b17023SJohn Marino 2684e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2685e4b17023SJohn Marino std::basic_ostream<_CharT, _Traits>& operator <<(std::basic_ostream<_CharT,_Traits> & __os,const piecewise_linear_distribution<_RealType> & __x)2686e4b17023SJohn Marino operator<<(std::basic_ostream<_CharT, _Traits>& __os, 2687e4b17023SJohn Marino const piecewise_linear_distribution<_RealType>& __x) 2688e4b17023SJohn Marino { 2689e4b17023SJohn Marino typedef std::basic_ostream<_CharT, _Traits> __ostream_type; 2690e4b17023SJohn Marino typedef typename __ostream_type::ios_base __ios_base; 2691e4b17023SJohn Marino 2692e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __os.flags(); 2693e4b17023SJohn Marino const _CharT __fill = __os.fill(); 2694e4b17023SJohn Marino const std::streamsize __precision = __os.precision(); 2695e4b17023SJohn Marino const _CharT __space = __os.widen(' '); 2696e4b17023SJohn Marino __os.flags(__ios_base::scientific | __ios_base::left); 2697e4b17023SJohn Marino __os.fill(__space); 2698e4b17023SJohn Marino __os.precision(std::numeric_limits<_RealType>::max_digits10); 2699e4b17023SJohn Marino 2700e4b17023SJohn Marino std::vector<_RealType> __int = __x.intervals(); 2701e4b17023SJohn Marino __os << __int.size() - 1; 2702e4b17023SJohn Marino 2703e4b17023SJohn Marino for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit) 2704e4b17023SJohn Marino __os << __space << *__xit; 2705e4b17023SJohn Marino 2706e4b17023SJohn Marino std::vector<double> __den = __x.densities(); 2707e4b17023SJohn Marino for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit) 2708e4b17023SJohn Marino __os << __space << *__dit; 2709e4b17023SJohn Marino 2710e4b17023SJohn Marino __os.flags(__flags); 2711e4b17023SJohn Marino __os.fill(__fill); 2712e4b17023SJohn Marino __os.precision(__precision); 2713e4b17023SJohn Marino return __os; 2714e4b17023SJohn Marino } 2715e4b17023SJohn Marino 2716e4b17023SJohn Marino template<typename _RealType, typename _CharT, typename _Traits> 2717e4b17023SJohn Marino std::basic_istream<_CharT, _Traits>& operator >>(std::basic_istream<_CharT,_Traits> & __is,piecewise_linear_distribution<_RealType> & __x)2718e4b17023SJohn Marino operator>>(std::basic_istream<_CharT, _Traits>& __is, 2719e4b17023SJohn Marino piecewise_linear_distribution<_RealType>& __x) 2720e4b17023SJohn Marino { 2721e4b17023SJohn Marino typedef std::basic_istream<_CharT, _Traits> __istream_type; 2722e4b17023SJohn Marino typedef typename __istream_type::ios_base __ios_base; 2723e4b17023SJohn Marino 2724e4b17023SJohn Marino const typename __ios_base::fmtflags __flags = __is.flags(); 2725e4b17023SJohn Marino __is.flags(__ios_base::dec | __ios_base::skipws); 2726e4b17023SJohn Marino 2727e4b17023SJohn Marino size_t __n; 2728e4b17023SJohn Marino __is >> __n; 2729e4b17023SJohn Marino 2730e4b17023SJohn Marino std::vector<_RealType> __int_vec; 2731e4b17023SJohn Marino __int_vec.reserve(__n + 1); 2732e4b17023SJohn Marino for (size_t __i = 0; __i <= __n; ++__i) 2733e4b17023SJohn Marino { 2734e4b17023SJohn Marino _RealType __int; 2735e4b17023SJohn Marino __is >> __int; 2736e4b17023SJohn Marino __int_vec.push_back(__int); 2737e4b17023SJohn Marino } 2738e4b17023SJohn Marino 2739e4b17023SJohn Marino std::vector<double> __den_vec; 2740e4b17023SJohn Marino __den_vec.reserve(__n + 1); 2741e4b17023SJohn Marino for (size_t __i = 0; __i <= __n; ++__i) 2742e4b17023SJohn Marino { 2743e4b17023SJohn Marino double __den; 2744e4b17023SJohn Marino __is >> __den; 2745e4b17023SJohn Marino __den_vec.push_back(__den); 2746e4b17023SJohn Marino } 2747e4b17023SJohn Marino 2748e4b17023SJohn Marino __x.param(typename piecewise_linear_distribution<_RealType>:: 2749e4b17023SJohn Marino param_type(__int_vec.begin(), __int_vec.end(), __den_vec.begin())); 2750e4b17023SJohn Marino 2751e4b17023SJohn Marino __is.flags(__flags); 2752e4b17023SJohn Marino return __is; 2753e4b17023SJohn Marino } 2754e4b17023SJohn Marino 2755e4b17023SJohn Marino 2756e4b17023SJohn Marino template<typename _IntType> seed_seq(std::initializer_list<_IntType> __il)2757e4b17023SJohn Marino seed_seq::seed_seq(std::initializer_list<_IntType> __il) 2758e4b17023SJohn Marino { 2759e4b17023SJohn Marino for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter) 2760e4b17023SJohn Marino _M_v.push_back(__detail::__mod<result_type, 2761e4b17023SJohn Marino __detail::_Shift<result_type, 32>::__value>(*__iter)); 2762e4b17023SJohn Marino } 2763e4b17023SJohn Marino 2764e4b17023SJohn Marino template<typename _InputIterator> seed_seq(_InputIterator __begin,_InputIterator __end)2765e4b17023SJohn Marino seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end) 2766e4b17023SJohn Marino { 2767e4b17023SJohn Marino for (_InputIterator __iter = __begin; __iter != __end; ++__iter) 2768e4b17023SJohn Marino _M_v.push_back(__detail::__mod<result_type, 2769e4b17023SJohn Marino __detail::_Shift<result_type, 32>::__value>(*__iter)); 2770e4b17023SJohn Marino } 2771e4b17023SJohn Marino 2772e4b17023SJohn Marino template<typename _RandomAccessIterator> 2773e4b17023SJohn Marino void generate(_RandomAccessIterator __begin,_RandomAccessIterator __end)2774e4b17023SJohn Marino seed_seq::generate(_RandomAccessIterator __begin, 2775e4b17023SJohn Marino _RandomAccessIterator __end) 2776e4b17023SJohn Marino { 2777e4b17023SJohn Marino typedef typename iterator_traits<_RandomAccessIterator>::value_type 2778e4b17023SJohn Marino _Type; 2779e4b17023SJohn Marino 2780e4b17023SJohn Marino if (__begin == __end) 2781e4b17023SJohn Marino return; 2782e4b17023SJohn Marino 2783e4b17023SJohn Marino std::fill(__begin, __end, _Type(0x8b8b8b8bu)); 2784e4b17023SJohn Marino 2785e4b17023SJohn Marino const size_t __n = __end - __begin; 2786e4b17023SJohn Marino const size_t __s = _M_v.size(); 2787e4b17023SJohn Marino const size_t __t = (__n >= 623) ? 11 2788e4b17023SJohn Marino : (__n >= 68) ? 7 2789e4b17023SJohn Marino : (__n >= 39) ? 5 2790e4b17023SJohn Marino : (__n >= 7) ? 3 2791e4b17023SJohn Marino : (__n - 1) / 2; 2792e4b17023SJohn Marino const size_t __p = (__n - __t) / 2; 2793e4b17023SJohn Marino const size_t __q = __p + __t; 2794e4b17023SJohn Marino const size_t __m = std::max(size_t(__s + 1), __n); 2795e4b17023SJohn Marino 2796e4b17023SJohn Marino for (size_t __k = 0; __k < __m; ++__k) 2797e4b17023SJohn Marino { 2798e4b17023SJohn Marino _Type __arg = (__begin[__k % __n] 2799e4b17023SJohn Marino ^ __begin[(__k + __p) % __n] 2800e4b17023SJohn Marino ^ __begin[(__k - 1) % __n]); 2801e4b17023SJohn Marino _Type __r1 = __arg ^ (__arg >> 27); 2802e4b17023SJohn Marino __r1 = __detail::__mod<_Type, 2803e4b17023SJohn Marino __detail::_Shift<_Type, 32>::__value>(1664525u * __r1); 2804e4b17023SJohn Marino _Type __r2 = __r1; 2805e4b17023SJohn Marino if (__k == 0) 2806e4b17023SJohn Marino __r2 += __s; 2807e4b17023SJohn Marino else if (__k <= __s) 2808e4b17023SJohn Marino __r2 += __k % __n + _M_v[__k - 1]; 2809e4b17023SJohn Marino else 2810e4b17023SJohn Marino __r2 += __k % __n; 2811e4b17023SJohn Marino __r2 = __detail::__mod<_Type, 2812e4b17023SJohn Marino __detail::_Shift<_Type, 32>::__value>(__r2); 2813e4b17023SJohn Marino __begin[(__k + __p) % __n] += __r1; 2814e4b17023SJohn Marino __begin[(__k + __q) % __n] += __r2; 2815e4b17023SJohn Marino __begin[__k % __n] = __r2; 2816e4b17023SJohn Marino } 2817e4b17023SJohn Marino 2818e4b17023SJohn Marino for (size_t __k = __m; __k < __m + __n; ++__k) 2819e4b17023SJohn Marino { 2820e4b17023SJohn Marino _Type __arg = (__begin[__k % __n] 2821e4b17023SJohn Marino + __begin[(__k + __p) % __n] 2822e4b17023SJohn Marino + __begin[(__k - 1) % __n]); 2823e4b17023SJohn Marino _Type __r3 = __arg ^ (__arg >> 27); 2824e4b17023SJohn Marino __r3 = __detail::__mod<_Type, 2825e4b17023SJohn Marino __detail::_Shift<_Type, 32>::__value>(1566083941u * __r3); 2826e4b17023SJohn Marino _Type __r4 = __r3 - __k % __n; 2827e4b17023SJohn Marino __r4 = __detail::__mod<_Type, 2828e4b17023SJohn Marino __detail::_Shift<_Type, 32>::__value>(__r4); 2829e4b17023SJohn Marino __begin[(__k + __p) % __n] ^= __r3; 2830e4b17023SJohn Marino __begin[(__k + __q) % __n] ^= __r4; 2831e4b17023SJohn Marino __begin[__k % __n] = __r4; 2832e4b17023SJohn Marino } 2833e4b17023SJohn Marino } 2834e4b17023SJohn Marino 2835e4b17023SJohn Marino template<typename _RealType, size_t __bits, 2836e4b17023SJohn Marino typename _UniformRandomNumberGenerator> 2837e4b17023SJohn Marino _RealType 2838e4b17023SJohn Marino generate_canonical(_UniformRandomNumberGenerator& __urng) 2839e4b17023SJohn Marino { 2840e4b17023SJohn Marino const size_t __b 2841e4b17023SJohn Marino = std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits), 2842e4b17023SJohn Marino __bits); 2843e4b17023SJohn Marino const long double __r = static_cast<long double>(__urng.max()) 2844e4b17023SJohn Marino - static_cast<long double>(__urng.min()) + 1.0L; 2845e4b17023SJohn Marino const size_t __log2r = std::log(__r) / std::log(2.0L); 2846e4b17023SJohn Marino size_t __k = std::max<size_t>(1UL, (__b + __log2r - 1UL) / __log2r); 2847e4b17023SJohn Marino _RealType __sum = _RealType(0); 2848e4b17023SJohn Marino _RealType __tmp = _RealType(1); 2849e4b17023SJohn Marino for (; __k != 0; --__k) 2850e4b17023SJohn Marino { 2851e4b17023SJohn Marino __sum += _RealType(__urng() - __urng.min()) * __tmp; 2852e4b17023SJohn Marino __tmp *= __r; 2853e4b17023SJohn Marino } 2854e4b17023SJohn Marino return __sum / __tmp; 2855e4b17023SJohn Marino } 2856e4b17023SJohn Marino 2857e4b17023SJohn Marino _GLIBCXX_END_NAMESPACE_VERSION 2858e4b17023SJohn Marino } // namespace 2859e4b17023SJohn Marino 2860e4b17023SJohn Marino #endif 2861