160e95dc3d
The deserialization functions for random number distributions fail to check the stream state before using the extracted values. In some cases this leads to using indeterminate values to resize a vector, and then filling that vector with indeterminate values. No values that affect control flow should be used without checking that a good value was read from the stream. Additionally, where reasonable to do so, defer modifying any state in the distribution until all values have been successfully read, to avoid modifying some of the distribution's parameters and leaving others unchanged. PR libstdc++/93205 * include/bits/random.h (operator>>): Check stream operation succeeds. * include/bits/random.tcc (operator<<): Remove redundant __ostream_type typedefs. (operator>>): Remove redundant __istream_type typedefs. Check stream operations succeed. (__extract_params): New function to fill a vector from a stream. * testsuite/26_numerics/random/pr60037-neg.cc: Adjust dg-error line. From-SVN: r280061
3317 lines
101 KiB
C++
3317 lines
101 KiB
C++
// random number generation (out of line) -*- C++ -*-
|
|
|
|
// Copyright (C) 2009-2020 Free Software Foundation, Inc.
|
|
//
|
|
// This file is part of the GNU ISO C++ Library. This library is free
|
|
// software; you can redistribute it and/or modify it under the
|
|
// terms of the GNU General Public License as published by the
|
|
// Free Software Foundation; either version 3, or (at your option)
|
|
// any later version.
|
|
|
|
// This library is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
// GNU General Public License for more details.
|
|
|
|
// Under Section 7 of GPL version 3, you are granted additional
|
|
// permissions described in the GCC Runtime Library Exception, version
|
|
// 3.1, as published by the Free Software Foundation.
|
|
|
|
// You should have received a copy of the GNU General Public License and
|
|
// a copy of the GCC Runtime Library Exception along with this program;
|
|
// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
// <http://www.gnu.org/licenses/>.
|
|
|
|
/** @file bits/random.tcc
|
|
* This is an internal header file, included by other library headers.
|
|
* Do not attempt to use it directly. @headername{random}
|
|
*/
|
|
|
|
#ifndef _RANDOM_TCC
|
|
#define _RANDOM_TCC 1
|
|
|
|
#include <numeric> // std::accumulate and std::partial_sum
|
|
|
|
namespace std _GLIBCXX_VISIBILITY(default)
|
|
{
|
|
_GLIBCXX_BEGIN_NAMESPACE_VERSION
|
|
|
|
/*
|
|
* (Further) implementation-space details.
|
|
*/
|
|
namespace __detail
|
|
{
|
|
// General case for x = (ax + c) mod m -- use Schrage's algorithm
|
|
// to avoid integer overflow.
|
|
//
|
|
// Preconditions: a > 0, m > 0.
|
|
//
|
|
// Note: only works correctly for __m % __a < __m / __a.
|
|
template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
|
|
_Tp
|
|
_Mod<_Tp, __m, __a, __c, false, true>::
|
|
__calc(_Tp __x)
|
|
{
|
|
if (__a == 1)
|
|
__x %= __m;
|
|
else
|
|
{
|
|
static const _Tp __q = __m / __a;
|
|
static const _Tp __r = __m % __a;
|
|
|
|
_Tp __t1 = __a * (__x % __q);
|
|
_Tp __t2 = __r * (__x / __q);
|
|
if (__t1 >= __t2)
|
|
__x = __t1 - __t2;
|
|
else
|
|
__x = __m - __t2 + __t1;
|
|
}
|
|
|
|
if (__c != 0)
|
|
{
|
|
const _Tp __d = __m - __x;
|
|
if (__d > __c)
|
|
__x += __c;
|
|
else
|
|
__x = __c - __d;
|
|
}
|
|
return __x;
|
|
}
|
|
|
|
template<typename _InputIterator, typename _OutputIterator,
|
|
typename _Tp>
|
|
_OutputIterator
|
|
__normalize(_InputIterator __first, _InputIterator __last,
|
|
_OutputIterator __result, const _Tp& __factor)
|
|
{
|
|
for (; __first != __last; ++__first, ++__result)
|
|
*__result = *__first / __factor;
|
|
return __result;
|
|
}
|
|
|
|
} // namespace __detail
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
constexpr _UIntType
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::multiplier;
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
constexpr _UIntType
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::increment;
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
constexpr _UIntType
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::modulus;
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
constexpr _UIntType
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::default_seed;
|
|
|
|
/**
|
|
* Seeds the LCR with integral value @p __s, adjusted so that the
|
|
* ring identity is never a member of the convergence set.
|
|
*/
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
void
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::
|
|
seed(result_type __s)
|
|
{
|
|
if ((__detail::__mod<_UIntType, __m>(__c) == 0)
|
|
&& (__detail::__mod<_UIntType, __m>(__s) == 0))
|
|
_M_x = 1;
|
|
else
|
|
_M_x = __detail::__mod<_UIntType, __m>(__s);
|
|
}
|
|
|
|
/**
|
|
* Seeds the LCR engine with a value generated by @p __q.
|
|
*/
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m>
|
|
template<typename _Sseq>
|
|
auto
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>::
|
|
seed(_Sseq& __q)
|
|
-> _If_seed_seq<_Sseq>
|
|
{
|
|
const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
|
|
: std::__lg(__m);
|
|
const _UIntType __k = (__k0 + 31) / 32;
|
|
uint_least32_t __arr[__k + 3];
|
|
__q.generate(__arr + 0, __arr + __k + 3);
|
|
_UIntType __factor = 1u;
|
|
_UIntType __sum = 0u;
|
|
for (size_t __j = 0; __j < __k; ++__j)
|
|
{
|
|
__sum += __arr[__j + 3] * __factor;
|
|
__factor *= __detail::_Shift<_UIntType, 32>::__value;
|
|
}
|
|
seed(__sum);
|
|
}
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const linear_congruential_engine<_UIntType,
|
|
__a, __c, __m>& __lcr)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
|
|
__os.fill(__os.widen(' '));
|
|
|
|
__os << __lcr._M_x;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _UIntType, _UIntType __a, _UIntType __c, _UIntType __m,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
linear_congruential_engine<_UIntType, __a, __c, __m>& __lcr)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec);
|
|
|
|
__is >> __lcr._M_x;
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::word_size;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::state_size;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::shift_size;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::mask_bits;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::xor_mask;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_u;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_d;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_s;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_b;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_t;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_c;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr size_t
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::tempering_l;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
initialization_multiplier;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
constexpr _UIntType
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::default_seed;
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
void
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
seed(result_type __sd)
|
|
{
|
|
_M_x[0] = __detail::__mod<_UIntType,
|
|
__detail::_Shift<_UIntType, __w>::__value>(__sd);
|
|
|
|
for (size_t __i = 1; __i < state_size; ++__i)
|
|
{
|
|
_UIntType __x = _M_x[__i - 1];
|
|
__x ^= __x >> (__w - 2);
|
|
__x *= __f;
|
|
__x += __detail::__mod<_UIntType, __n>(__i);
|
|
_M_x[__i] = __detail::__mod<_UIntType,
|
|
__detail::_Shift<_UIntType, __w>::__value>(__x);
|
|
}
|
|
_M_p = state_size;
|
|
}
|
|
|
|
template<typename _UIntType,
|
|
size_t __w, size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
template<typename _Sseq>
|
|
auto
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
seed(_Sseq& __q)
|
|
-> _If_seed_seq<_Sseq>
|
|
{
|
|
const _UIntType __upper_mask = (~_UIntType()) << __r;
|
|
const size_t __k = (__w + 31) / 32;
|
|
uint_least32_t __arr[__n * __k];
|
|
__q.generate(__arr + 0, __arr + __n * __k);
|
|
|
|
bool __zero = true;
|
|
for (size_t __i = 0; __i < state_size; ++__i)
|
|
{
|
|
_UIntType __factor = 1u;
|
|
_UIntType __sum = 0u;
|
|
for (size_t __j = 0; __j < __k; ++__j)
|
|
{
|
|
__sum += __arr[__k * __i + __j] * __factor;
|
|
__factor *= __detail::_Shift<_UIntType, 32>::__value;
|
|
}
|
|
_M_x[__i] = __detail::__mod<_UIntType,
|
|
__detail::_Shift<_UIntType, __w>::__value>(__sum);
|
|
|
|
if (__zero)
|
|
{
|
|
if (__i == 0)
|
|
{
|
|
if ((_M_x[0] & __upper_mask) != 0u)
|
|
__zero = false;
|
|
}
|
|
else if (_M_x[__i] != 0u)
|
|
__zero = false;
|
|
}
|
|
}
|
|
if (__zero)
|
|
_M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
|
|
_M_p = state_size;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w,
|
|
size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
void
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
_M_gen_rand(void)
|
|
{
|
|
const _UIntType __upper_mask = (~_UIntType()) << __r;
|
|
const _UIntType __lower_mask = ~__upper_mask;
|
|
|
|
for (size_t __k = 0; __k < (__n - __m); ++__k)
|
|
{
|
|
_UIntType __y = ((_M_x[__k] & __upper_mask)
|
|
| (_M_x[__k + 1] & __lower_mask));
|
|
_M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
|
|
^ ((__y & 0x01) ? __a : 0));
|
|
}
|
|
|
|
for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
|
|
{
|
|
_UIntType __y = ((_M_x[__k] & __upper_mask)
|
|
| (_M_x[__k + 1] & __lower_mask));
|
|
_M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
|
|
^ ((__y & 0x01) ? __a : 0));
|
|
}
|
|
|
|
_UIntType __y = ((_M_x[__n - 1] & __upper_mask)
|
|
| (_M_x[0] & __lower_mask));
|
|
_M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
|
|
^ ((__y & 0x01) ? __a : 0));
|
|
_M_p = 0;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w,
|
|
size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
void
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
discard(unsigned long long __z)
|
|
{
|
|
while (__z > state_size - _M_p)
|
|
{
|
|
__z -= state_size - _M_p;
|
|
_M_gen_rand();
|
|
}
|
|
_M_p += __z;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w,
|
|
size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f>
|
|
typename
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::result_type
|
|
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
|
|
__s, __b, __t, __c, __l, __f>::
|
|
operator()()
|
|
{
|
|
// Reload the vector - cost is O(n) amortized over n calls.
|
|
if (_M_p >= state_size)
|
|
_M_gen_rand();
|
|
|
|
// Calculate o(x(i)).
|
|
result_type __z = _M_x[_M_p++];
|
|
__z ^= (__z >> __u) & __d;
|
|
__z ^= (__z << __s) & __b;
|
|
__z ^= (__z << __t) & __c;
|
|
__z ^= (__z >> __l);
|
|
|
|
return __z;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w,
|
|
size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f, typename _CharT, typename _Traits>
|
|
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)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
|
|
__os.fill(__space);
|
|
|
|
for (size_t __i = 0; __i < __n; ++__i)
|
|
__os << __x._M_x[__i] << __space;
|
|
__os << __x._M_p;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w,
|
|
size_t __n, size_t __m, size_t __r,
|
|
_UIntType __a, size_t __u, _UIntType __d, size_t __s,
|
|
_UIntType __b, size_t __t, _UIntType __c, size_t __l,
|
|
_UIntType __f, typename _CharT, typename _Traits>
|
|
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)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
for (size_t __i = 0; __i < __n; ++__i)
|
|
__is >> __x._M_x[__i];
|
|
__is >> __x._M_p;
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
constexpr size_t
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::word_size;
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
constexpr size_t
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::short_lag;
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
constexpr size_t
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::long_lag;
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
constexpr _UIntType
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::default_seed;
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
void
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::
|
|
seed(result_type __value)
|
|
{
|
|
std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
|
|
__lcg(__value == 0u ? default_seed : __value);
|
|
|
|
const size_t __n = (__w + 31) / 32;
|
|
|
|
for (size_t __i = 0; __i < long_lag; ++__i)
|
|
{
|
|
_UIntType __sum = 0u;
|
|
_UIntType __factor = 1u;
|
|
for (size_t __j = 0; __j < __n; ++__j)
|
|
{
|
|
__sum += __detail::__mod<uint_least32_t,
|
|
__detail::_Shift<uint_least32_t, 32>::__value>
|
|
(__lcg()) * __factor;
|
|
__factor *= __detail::_Shift<_UIntType, 32>::__value;
|
|
}
|
|
_M_x[__i] = __detail::__mod<_UIntType,
|
|
__detail::_Shift<_UIntType, __w>::__value>(__sum);
|
|
}
|
|
_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
|
|
_M_p = 0;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
template<typename _Sseq>
|
|
auto
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::
|
|
seed(_Sseq& __q)
|
|
-> _If_seed_seq<_Sseq>
|
|
{
|
|
const size_t __k = (__w + 31) / 32;
|
|
uint_least32_t __arr[__r * __k];
|
|
__q.generate(__arr + 0, __arr + __r * __k);
|
|
|
|
for (size_t __i = 0; __i < long_lag; ++__i)
|
|
{
|
|
_UIntType __sum = 0u;
|
|
_UIntType __factor = 1u;
|
|
for (size_t __j = 0; __j < __k; ++__j)
|
|
{
|
|
__sum += __arr[__k * __i + __j] * __factor;
|
|
__factor *= __detail::_Shift<_UIntType, 32>::__value;
|
|
}
|
|
_M_x[__i] = __detail::__mod<_UIntType,
|
|
__detail::_Shift<_UIntType, __w>::__value>(__sum);
|
|
}
|
|
_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
|
|
_M_p = 0;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
|
|
typename subtract_with_carry_engine<_UIntType, __w, __s, __r>::
|
|
result_type
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>::
|
|
operator()()
|
|
{
|
|
// Derive short lag index from current index.
|
|
long __ps = _M_p - short_lag;
|
|
if (__ps < 0)
|
|
__ps += long_lag;
|
|
|
|
// Calculate new x(i) without overflow or division.
|
|
// NB: Thanks to the requirements for _UIntType, _M_x[_M_p] + _M_carry
|
|
// cannot overflow.
|
|
_UIntType __xi;
|
|
if (_M_x[__ps] >= _M_x[_M_p] + _M_carry)
|
|
{
|
|
__xi = _M_x[__ps] - _M_x[_M_p] - _M_carry;
|
|
_M_carry = 0;
|
|
}
|
|
else
|
|
{
|
|
__xi = (__detail::_Shift<_UIntType, __w>::__value
|
|
- _M_x[_M_p] - _M_carry + _M_x[__ps]);
|
|
_M_carry = 1;
|
|
}
|
|
_M_x[_M_p] = __xi;
|
|
|
|
// Adjust current index to loop around in ring buffer.
|
|
if (++_M_p >= long_lag)
|
|
_M_p = 0;
|
|
|
|
return __xi;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const subtract_with_carry_engine<_UIntType,
|
|
__w, __s, __r>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
|
|
__os.fill(__space);
|
|
|
|
for (size_t __i = 0; __i < __r; ++__i)
|
|
__os << __x._M_x[__i] << __space;
|
|
__os << __x._M_carry << __space << __x._M_p;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _UIntType, size_t __w, size_t __s, size_t __r,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
subtract_with_carry_engine<_UIntType, __w, __s, __r>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
for (size_t __i = 0; __i < __r; ++__i)
|
|
__is >> __x._M_x[__i];
|
|
__is >> __x._M_carry;
|
|
__is >> __x._M_p;
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RandomNumberEngine, size_t __p, size_t __r>
|
|
constexpr size_t
|
|
discard_block_engine<_RandomNumberEngine, __p, __r>::block_size;
|
|
|
|
template<typename _RandomNumberEngine, size_t __p, size_t __r>
|
|
constexpr size_t
|
|
discard_block_engine<_RandomNumberEngine, __p, __r>::used_block;
|
|
|
|
template<typename _RandomNumberEngine, size_t __p, size_t __r>
|
|
typename discard_block_engine<_RandomNumberEngine,
|
|
__p, __r>::result_type
|
|
discard_block_engine<_RandomNumberEngine, __p, __r>::
|
|
operator()()
|
|
{
|
|
if (_M_n >= used_block)
|
|
{
|
|
_M_b.discard(block_size - _M_n);
|
|
_M_n = 0;
|
|
}
|
|
++_M_n;
|
|
return _M_b();
|
|
}
|
|
|
|
template<typename _RandomNumberEngine, size_t __p, size_t __r,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const discard_block_engine<_RandomNumberEngine,
|
|
__p, __r>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
|
|
__os.fill(__space);
|
|
|
|
__os << __x.base() << __space << __x._M_n;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RandomNumberEngine, size_t __p, size_t __r,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
discard_block_engine<_RandomNumberEngine, __p, __r>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
__is >> __x._M_b >> __x._M_n;
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
|
|
typename independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
|
|
result_type
|
|
independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
|
|
operator()()
|
|
{
|
|
typedef typename _RandomNumberEngine::result_type _Eresult_type;
|
|
const _Eresult_type __r
|
|
= (_M_b.max() - _M_b.min() < std::numeric_limits<_Eresult_type>::max()
|
|
? _M_b.max() - _M_b.min() + 1 : 0);
|
|
const unsigned __edig = std::numeric_limits<_Eresult_type>::digits;
|
|
const unsigned __m = __r ? std::__lg(__r) : __edig;
|
|
|
|
typedef typename std::common_type<_Eresult_type, result_type>::type
|
|
__ctype;
|
|
const unsigned __cdig = std::numeric_limits<__ctype>::digits;
|
|
|
|
unsigned __n, __n0;
|
|
__ctype __s0, __s1, __y0, __y1;
|
|
|
|
for (size_t __i = 0; __i < 2; ++__i)
|
|
{
|
|
__n = (__w + __m - 1) / __m + __i;
|
|
__n0 = __n - __w % __n;
|
|
const unsigned __w0 = __w / __n; // __w0 <= __m
|
|
|
|
__s0 = 0;
|
|
__s1 = 0;
|
|
if (__w0 < __cdig)
|
|
{
|
|
__s0 = __ctype(1) << __w0;
|
|
__s1 = __s0 << 1;
|
|
}
|
|
|
|
__y0 = 0;
|
|
__y1 = 0;
|
|
if (__r)
|
|
{
|
|
__y0 = __s0 * (__r / __s0);
|
|
if (__s1)
|
|
__y1 = __s1 * (__r / __s1);
|
|
|
|
if (__r - __y0 <= __y0 / __n)
|
|
break;
|
|
}
|
|
else
|
|
break;
|
|
}
|
|
|
|
result_type __sum = 0;
|
|
for (size_t __k = 0; __k < __n0; ++__k)
|
|
{
|
|
__ctype __u;
|
|
do
|
|
__u = _M_b() - _M_b.min();
|
|
while (__y0 && __u >= __y0);
|
|
__sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
|
|
}
|
|
for (size_t __k = __n0; __k < __n; ++__k)
|
|
{
|
|
__ctype __u;
|
|
do
|
|
__u = _M_b() - _M_b.min();
|
|
while (__y1 && __u >= __y1);
|
|
__sum = __s1 * __sum + (__s1 ? __u % __s1 : __u);
|
|
}
|
|
return __sum;
|
|
}
|
|
|
|
|
|
template<typename _RandomNumberEngine, size_t __k>
|
|
constexpr size_t
|
|
shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
|
|
|
|
template<typename _RandomNumberEngine, size_t __k>
|
|
typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
|
|
shuffle_order_engine<_RandomNumberEngine, __k>::
|
|
operator()()
|
|
{
|
|
size_t __j = __k * ((_M_y - _M_b.min())
|
|
/ (_M_b.max() - _M_b.min() + 1.0L));
|
|
_M_y = _M_v[__j];
|
|
_M_v[__j] = _M_b();
|
|
|
|
return _M_y;
|
|
}
|
|
|
|
template<typename _RandomNumberEngine, size_t __k,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const shuffle_order_engine<_RandomNumberEngine, __k>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
|
|
__os.fill(__space);
|
|
|
|
__os << __x.base();
|
|
for (size_t __i = 0; __i < __k; ++__i)
|
|
__os << __space << __x._M_v[__i];
|
|
__os << __space << __x._M_y;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RandomNumberEngine, size_t __k,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
shuffle_order_engine<_RandomNumberEngine, __k>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
__is >> __x._M_b;
|
|
for (size_t __i = 0; __i < __k; ++__i)
|
|
__is >> __x._M_v[__i];
|
|
__is >> __x._M_y;
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const uniform_int_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
|
|
__os << __x.a() << __space << __x.b();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
uniform_int_distribution<_IntType>& __x)
|
|
{
|
|
using param_type
|
|
= typename uniform_int_distribution<_IntType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_IntType __a, __b;
|
|
if (__is >> __a >> __b)
|
|
__x.param(param_type(__a, __b));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
uniform_real_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
auto __range = __p.b() - __p.a();
|
|
while (__f != __t)
|
|
*__f++ = __aurng() * __range + __p.a();
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const uniform_real_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.a() << __space << __x.b();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
uniform_real_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename uniform_real_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::skipws);
|
|
|
|
_RealType __a, __b;
|
|
if (__is >> __a >> __b)
|
|
__x.param(param_type(__a, __b));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::bernoulli_distribution::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
auto __limit = __p.p() * (__aurng.max() - __aurng.min());
|
|
|
|
while (__f != __t)
|
|
*__f++ = (__aurng() - __aurng.min()) < __limit;
|
|
}
|
|
|
|
template<typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const bernoulli_distribution& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__os.widen(' '));
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
__os << __x.p();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename geometric_distribution<_IntType>::result_type
|
|
geometric_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
// About the epsilon thing see this thread:
|
|
// http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
|
|
const double __naf =
|
|
(1 - std::numeric_limits<double>::epsilon()) / 2;
|
|
// The largest _RealType convertible to _IntType.
|
|
const double __thr =
|
|
std::numeric_limits<_IntType>::max() + __naf;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
double __cand;
|
|
do
|
|
__cand = std::floor(std::log(1.0 - __aurng()) / __param._M_log_1_p);
|
|
while (__cand >= __thr);
|
|
|
|
return result_type(__cand + __naf);
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
geometric_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
// About the epsilon thing see this thread:
|
|
// http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00971.html
|
|
const double __naf =
|
|
(1 - std::numeric_limits<double>::epsilon()) / 2;
|
|
// The largest _RealType convertible to _IntType.
|
|
const double __thr =
|
|
std::numeric_limits<_IntType>::max() + __naf;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
while (__f != __t)
|
|
{
|
|
double __cand;
|
|
do
|
|
__cand = std::floor(std::log(1.0 - __aurng())
|
|
/ __param._M_log_1_p);
|
|
while (__cand >= __thr);
|
|
|
|
*__f++ = __cand + __naf;
|
|
}
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const geometric_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__os.widen(' '));
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
__os << __x.p();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
geometric_distribution<_IntType>& __x)
|
|
{
|
|
using param_type = typename geometric_distribution<_IntType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::skipws);
|
|
|
|
double __p;
|
|
if (__is >> __p)
|
|
__x.param(param_type(__p));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
// This is Leger's algorithm, also in Devroye, Ch. X, Example 1.5.
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename negative_binomial_distribution<_IntType>::result_type
|
|
negative_binomial_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
const double __y = _M_gd(__urng);
|
|
|
|
// XXX Is the constructor too slow?
|
|
std::poisson_distribution<result_type> __poisson(__y);
|
|
return __poisson(__urng);
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename negative_binomial_distribution<_IntType>::result_type
|
|
negative_binomial_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
typedef typename std::gamma_distribution<double>::param_type
|
|
param_type;
|
|
|
|
const double __y =
|
|
_M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
|
|
|
|
std::poisson_distribution<result_type> __poisson(__y);
|
|
return __poisson(__urng);
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
negative_binomial_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
{
|
|
const double __y = _M_gd(__urng);
|
|
|
|
// XXX Is the constructor too slow?
|
|
std::poisson_distribution<result_type> __poisson(__y);
|
|
*__f++ = __poisson(__urng);
|
|
}
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
negative_binomial_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
typename std::gamma_distribution<result_type>::param_type
|
|
__p2(__p.k(), (1.0 - __p.p()) / __p.p());
|
|
|
|
while (__f != __t)
|
|
{
|
|
const double __y = _M_gd(__urng, __p2);
|
|
|
|
std::poisson_distribution<result_type> __poisson(__y);
|
|
*__f++ = __poisson(__urng);
|
|
}
|
|
}
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const negative_binomial_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__os.widen(' '));
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
__os << __x.k() << __space << __x.p()
|
|
<< __space << __x._M_gd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
negative_binomial_distribution<_IntType>& __x)
|
|
{
|
|
using param_type
|
|
= typename negative_binomial_distribution<_IntType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::skipws);
|
|
|
|
_IntType __k;
|
|
double __p;
|
|
if (__is >> __k >> __p >> __x._M_gd)
|
|
__x.param(param_type(__k, __p));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _IntType>
|
|
void
|
|
poisson_distribution<_IntType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
#if _GLIBCXX_USE_C99_MATH_TR1
|
|
if (_M_mean >= 12)
|
|
{
|
|
const double __m = std::floor(_M_mean);
|
|
_M_lm_thr = std::log(_M_mean);
|
|
_M_lfm = std::lgamma(__m + 1);
|
|
_M_sm = std::sqrt(__m);
|
|
|
|
const double __pi_4 = 0.7853981633974483096156608458198757L;
|
|
const double __dx = std::sqrt(2 * __m * std::log(32 * __m
|
|
/ __pi_4));
|
|
_M_d = std::round(std::max<double>(6.0, std::min(__m, __dx)));
|
|
const double __cx = 2 * __m + _M_d;
|
|
_M_scx = std::sqrt(__cx / 2);
|
|
_M_1cx = 1 / __cx;
|
|
|
|
_M_c2b = std::sqrt(__pi_4 * __cx) * std::exp(_M_1cx);
|
|
_M_cb = 2 * __cx * std::exp(-_M_d * _M_1cx * (1 + _M_d / 2))
|
|
/ _M_d;
|
|
}
|
|
else
|
|
#endif
|
|
_M_lm_thr = std::exp(-_M_mean);
|
|
}
|
|
|
|
/**
|
|
* A rejection algorithm when mean >= 12 and a simple method based
|
|
* upon the multiplication of uniform random variates otherwise.
|
|
* NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
|
|
* is defined.
|
|
*
|
|
* Reference:
|
|
* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
|
|
* New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
|
|
*/
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename poisson_distribution<_IntType>::result_type
|
|
poisson_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
#if _GLIBCXX_USE_C99_MATH_TR1
|
|
if (__param.mean() >= 12)
|
|
{
|
|
double __x;
|
|
|
|
// See comments above...
|
|
const double __naf =
|
|
(1 - std::numeric_limits<double>::epsilon()) / 2;
|
|
const double __thr =
|
|
std::numeric_limits<_IntType>::max() + __naf;
|
|
|
|
const double __m = std::floor(__param.mean());
|
|
// sqrt(pi / 2)
|
|
const double __spi_2 = 1.2533141373155002512078826424055226L;
|
|
const double __c1 = __param._M_sm * __spi_2;
|
|
const double __c2 = __param._M_c2b + __c1;
|
|
const double __c3 = __c2 + 1;
|
|
const double __c4 = __c3 + 1;
|
|
// 1 / 78
|
|
const double __178 = 0.0128205128205128205128205128205128L;
|
|
// e^(1 / 78)
|
|
const double __e178 = 1.0129030479320018583185514777512983L;
|
|
const double __c5 = __c4 + __e178;
|
|
const double __c = __param._M_cb + __c5;
|
|
const double __2cx = 2 * (2 * __m + __param._M_d);
|
|
|
|
bool __reject = true;
|
|
do
|
|
{
|
|
const double __u = __c * __aurng();
|
|
const double __e = -std::log(1.0 - __aurng());
|
|
|
|
double __w = 0.0;
|
|
|
|
if (__u <= __c1)
|
|
{
|
|
const double __n = _M_nd(__urng);
|
|
const double __y = -std::abs(__n) * __param._M_sm - 1;
|
|
__x = std::floor(__y);
|
|
__w = -__n * __n / 2;
|
|
if (__x < -__m)
|
|
continue;
|
|
}
|
|
else if (__u <= __c2)
|
|
{
|
|
const double __n = _M_nd(__urng);
|
|
const double __y = 1 + std::abs(__n) * __param._M_scx;
|
|
__x = std::ceil(__y);
|
|
__w = __y * (2 - __y) * __param._M_1cx;
|
|
if (__x > __param._M_d)
|
|
continue;
|
|
}
|
|
else if (__u <= __c3)
|
|
// NB: This case not in the book, nor in the Errata,
|
|
// but should be ok...
|
|
__x = -1;
|
|
else if (__u <= __c4)
|
|
__x = 0;
|
|
else if (__u <= __c5)
|
|
{
|
|
__x = 1;
|
|
// Only in the Errata, see libstdc++/83237.
|
|
__w = __178;
|
|
}
|
|
else
|
|
{
|
|
const double __v = -std::log(1.0 - __aurng());
|
|
const double __y = __param._M_d
|
|
+ __v * __2cx / __param._M_d;
|
|
__x = std::ceil(__y);
|
|
__w = -__param._M_d * __param._M_1cx * (1 + __y / 2);
|
|
}
|
|
|
|
__reject = (__w - __e - __x * __param._M_lm_thr
|
|
> __param._M_lfm - std::lgamma(__x + __m + 1));
|
|
|
|
__reject |= __x + __m >= __thr;
|
|
|
|
} while (__reject);
|
|
|
|
return result_type(__x + __m + __naf);
|
|
}
|
|
else
|
|
#endif
|
|
{
|
|
_IntType __x = 0;
|
|
double __prod = 1.0;
|
|
|
|
do
|
|
{
|
|
__prod *= __aurng();
|
|
__x += 1;
|
|
}
|
|
while (__prod > __param._M_lm_thr);
|
|
|
|
return __x - 1;
|
|
}
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
poisson_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
// We could duplicate everything from operator()...
|
|
while (__f != __t)
|
|
*__f++ = this->operator()(__urng, __param);
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const poisson_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
__os << __x.mean() << __space << __x._M_nd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
poisson_distribution<_IntType>& __x)
|
|
{
|
|
using param_type = typename poisson_distribution<_IntType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::skipws);
|
|
|
|
double __mean;
|
|
if (__is >> __mean >> __x._M_nd)
|
|
__x.param(param_type(__mean));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _IntType>
|
|
void
|
|
binomial_distribution<_IntType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
const double __p12 = _M_p <= 0.5 ? _M_p : 1.0 - _M_p;
|
|
|
|
_M_easy = true;
|
|
|
|
#if _GLIBCXX_USE_C99_MATH_TR1
|
|
if (_M_t * __p12 >= 8)
|
|
{
|
|
_M_easy = false;
|
|
const double __np = std::floor(_M_t * __p12);
|
|
const double __pa = __np / _M_t;
|
|
const double __1p = 1 - __pa;
|
|
|
|
const double __pi_4 = 0.7853981633974483096156608458198757L;
|
|
const double __d1x =
|
|
std::sqrt(__np * __1p * std::log(32 * __np
|
|
/ (81 * __pi_4 * __1p)));
|
|
_M_d1 = std::round(std::max<double>(1.0, __d1x));
|
|
const double __d2x =
|
|
std::sqrt(__np * __1p * std::log(32 * _M_t * __1p
|
|
/ (__pi_4 * __pa)));
|
|
_M_d2 = std::round(std::max<double>(1.0, __d2x));
|
|
|
|
// sqrt(pi / 2)
|
|
const double __spi_2 = 1.2533141373155002512078826424055226L;
|
|
_M_s1 = std::sqrt(__np * __1p) * (1 + _M_d1 / (4 * __np));
|
|
_M_s2 = std::sqrt(__np * __1p) * (1 + _M_d2 / (4 * _M_t * __1p));
|
|
_M_c = 2 * _M_d1 / __np;
|
|
_M_a1 = std::exp(_M_c) * _M_s1 * __spi_2;
|
|
const double __a12 = _M_a1 + _M_s2 * __spi_2;
|
|
const double __s1s = _M_s1 * _M_s1;
|
|
_M_a123 = __a12 + (std::exp(_M_d1 / (_M_t * __1p))
|
|
* 2 * __s1s / _M_d1
|
|
* std::exp(-_M_d1 * _M_d1 / (2 * __s1s)));
|
|
const double __s2s = _M_s2 * _M_s2;
|
|
_M_s = (_M_a123 + 2 * __s2s / _M_d2
|
|
* std::exp(-_M_d2 * _M_d2 / (2 * __s2s)));
|
|
_M_lf = (std::lgamma(__np + 1)
|
|
+ std::lgamma(_M_t - __np + 1));
|
|
_M_lp1p = std::log(__pa / __1p);
|
|
|
|
_M_q = -std::log(1 - (__p12 - __pa) / __1p);
|
|
}
|
|
else
|
|
#endif
|
|
_M_q = -std::log(1 - __p12);
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename binomial_distribution<_IntType>::result_type
|
|
binomial_distribution<_IntType>::
|
|
_M_waiting(_UniformRandomNumberGenerator& __urng,
|
|
_IntType __t, double __q)
|
|
{
|
|
_IntType __x = 0;
|
|
double __sum = 0.0;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
do
|
|
{
|
|
if (__t == __x)
|
|
return __x;
|
|
const double __e = -std::log(1.0 - __aurng());
|
|
__sum += __e / (__t - __x);
|
|
__x += 1;
|
|
}
|
|
while (__sum <= __q);
|
|
|
|
return __x - 1;
|
|
}
|
|
|
|
/**
|
|
* A rejection algorithm when t * p >= 8 and a simple waiting time
|
|
* method - the second in the referenced book - otherwise.
|
|
* NB: The former is available only if _GLIBCXX_USE_C99_MATH_TR1
|
|
* is defined.
|
|
*
|
|
* Reference:
|
|
* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
|
|
* New York, 1986, Ch. X, Sect. 4 (+ Errata!).
|
|
*/
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename binomial_distribution<_IntType>::result_type
|
|
binomial_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
result_type __ret;
|
|
const _IntType __t = __param.t();
|
|
const double __p = __param.p();
|
|
const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
#if _GLIBCXX_USE_C99_MATH_TR1
|
|
if (!__param._M_easy)
|
|
{
|
|
double __x;
|
|
|
|
// See comments above...
|
|
const double __naf =
|
|
(1 - std::numeric_limits<double>::epsilon()) / 2;
|
|
const double __thr =
|
|
std::numeric_limits<_IntType>::max() + __naf;
|
|
|
|
const double __np = std::floor(__t * __p12);
|
|
|
|
// sqrt(pi / 2)
|
|
const double __spi_2 = 1.2533141373155002512078826424055226L;
|
|
const double __a1 = __param._M_a1;
|
|
const double __a12 = __a1 + __param._M_s2 * __spi_2;
|
|
const double __a123 = __param._M_a123;
|
|
const double __s1s = __param._M_s1 * __param._M_s1;
|
|
const double __s2s = __param._M_s2 * __param._M_s2;
|
|
|
|
bool __reject;
|
|
do
|
|
{
|
|
const double __u = __param._M_s * __aurng();
|
|
|
|
double __v;
|
|
|
|
if (__u <= __a1)
|
|
{
|
|
const double __n = _M_nd(__urng);
|
|
const double __y = __param._M_s1 * std::abs(__n);
|
|
__reject = __y >= __param._M_d1;
|
|
if (!__reject)
|
|
{
|
|
const double __e = -std::log(1.0 - __aurng());
|
|
__x = std::floor(__y);
|
|
__v = -__e - __n * __n / 2 + __param._M_c;
|
|
}
|
|
}
|
|
else if (__u <= __a12)
|
|
{
|
|
const double __n = _M_nd(__urng);
|
|
const double __y = __param._M_s2 * std::abs(__n);
|
|
__reject = __y >= __param._M_d2;
|
|
if (!__reject)
|
|
{
|
|
const double __e = -std::log(1.0 - __aurng());
|
|
__x = std::floor(-__y);
|
|
__v = -__e - __n * __n / 2;
|
|
}
|
|
}
|
|
else if (__u <= __a123)
|
|
{
|
|
const double __e1 = -std::log(1.0 - __aurng());
|
|
const double __e2 = -std::log(1.0 - __aurng());
|
|
|
|
const double __y = __param._M_d1
|
|
+ 2 * __s1s * __e1 / __param._M_d1;
|
|
__x = std::floor(__y);
|
|
__v = (-__e2 + __param._M_d1 * (1 / (__t - __np)
|
|
-__y / (2 * __s1s)));
|
|
__reject = false;
|
|
}
|
|
else
|
|
{
|
|
const double __e1 = -std::log(1.0 - __aurng());
|
|
const double __e2 = -std::log(1.0 - __aurng());
|
|
|
|
const double __y = __param._M_d2
|
|
+ 2 * __s2s * __e1 / __param._M_d2;
|
|
__x = std::floor(-__y);
|
|
__v = -__e2 - __param._M_d2 * __y / (2 * __s2s);
|
|
__reject = false;
|
|
}
|
|
|
|
__reject = __reject || __x < -__np || __x > __t - __np;
|
|
if (!__reject)
|
|
{
|
|
const double __lfx =
|
|
std::lgamma(__np + __x + 1)
|
|
+ std::lgamma(__t - (__np + __x) + 1);
|
|
__reject = __v > __param._M_lf - __lfx
|
|
+ __x * __param._M_lp1p;
|
|
}
|
|
|
|
__reject |= __x + __np >= __thr;
|
|
}
|
|
while (__reject);
|
|
|
|
__x += __np + __naf;
|
|
|
|
const _IntType __z = _M_waiting(__urng, __t - _IntType(__x),
|
|
__param._M_q);
|
|
__ret = _IntType(__x) + __z;
|
|
}
|
|
else
|
|
#endif
|
|
__ret = _M_waiting(__urng, __t, __param._M_q);
|
|
|
|
if (__p12 != __p)
|
|
__ret = __t - __ret;
|
|
return __ret;
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
binomial_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
// We could duplicate everything from operator()...
|
|
while (__f != __t)
|
|
*__f++ = this->operator()(__urng, __param);
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const binomial_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
__os << __x.t() << __space << __x.p()
|
|
<< __space << __x._M_nd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _IntType,
|
|
typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
binomial_distribution<_IntType>& __x)
|
|
{
|
|
using param_type = typename binomial_distribution<_IntType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_IntType __t;
|
|
double __p;
|
|
if (__is >> __t >> __p >> __x._M_nd)
|
|
__x.param(param_type(__t, __p));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::exponential_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
while (__f != __t)
|
|
*__f++ = -std::log(result_type(1) - __aurng()) / __p.lambda();
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const exponential_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__os.widen(' '));
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.lambda();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
exponential_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename exponential_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __lambda;
|
|
if (__is >> __lambda)
|
|
__x.param(param_type(__lambda));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
/**
|
|
* Polar method due to Marsaglia.
|
|
*
|
|
* Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
|
|
* New York, 1986, Ch. V, Sect. 4.4.
|
|
*/
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename normal_distribution<_RealType>::result_type
|
|
normal_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
result_type __ret;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
|
|
if (_M_saved_available)
|
|
{
|
|
_M_saved_available = false;
|
|
__ret = _M_saved;
|
|
}
|
|
else
|
|
{
|
|
result_type __x, __y, __r2;
|
|
do
|
|
{
|
|
__x = result_type(2.0) * __aurng() - 1.0;
|
|
__y = result_type(2.0) * __aurng() - 1.0;
|
|
__r2 = __x * __x + __y * __y;
|
|
}
|
|
while (__r2 > 1.0 || __r2 == 0.0);
|
|
|
|
const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
|
|
_M_saved = __x * __mult;
|
|
_M_saved_available = true;
|
|
__ret = __y * __mult;
|
|
}
|
|
|
|
__ret = __ret * __param.stddev() + __param.mean();
|
|
return __ret;
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
normal_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
|
|
if (__f == __t)
|
|
return;
|
|
|
|
if (_M_saved_available)
|
|
{
|
|
_M_saved_available = false;
|
|
*__f++ = _M_saved * __param.stddev() + __param.mean();
|
|
|
|
if (__f == __t)
|
|
return;
|
|
}
|
|
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
|
|
while (__f + 1 < __t)
|
|
{
|
|
result_type __x, __y, __r2;
|
|
do
|
|
{
|
|
__x = result_type(2.0) * __aurng() - 1.0;
|
|
__y = result_type(2.0) * __aurng() - 1.0;
|
|
__r2 = __x * __x + __y * __y;
|
|
}
|
|
while (__r2 > 1.0 || __r2 == 0.0);
|
|
|
|
const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
|
|
*__f++ = __y * __mult * __param.stddev() + __param.mean();
|
|
*__f++ = __x * __mult * __param.stddev() + __param.mean();
|
|
}
|
|
|
|
if (__f != __t)
|
|
{
|
|
result_type __x, __y, __r2;
|
|
do
|
|
{
|
|
__x = result_type(2.0) * __aurng() - 1.0;
|
|
__y = result_type(2.0) * __aurng() - 1.0;
|
|
__r2 = __x * __x + __y * __y;
|
|
}
|
|
while (__r2 > 1.0 || __r2 == 0.0);
|
|
|
|
const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
|
|
_M_saved = __x * __mult;
|
|
_M_saved_available = true;
|
|
*__f = __y * __mult * __param.stddev() + __param.mean();
|
|
}
|
|
}
|
|
|
|
template<typename _RealType>
|
|
bool
|
|
operator==(const std::normal_distribution<_RealType>& __d1,
|
|
const std::normal_distribution<_RealType>& __d2)
|
|
{
|
|
if (__d1._M_param == __d2._M_param
|
|
&& __d1._M_saved_available == __d2._M_saved_available)
|
|
{
|
|
if (__d1._M_saved_available
|
|
&& __d1._M_saved == __d2._M_saved)
|
|
return true;
|
|
else if(!__d1._M_saved_available)
|
|
return true;
|
|
else
|
|
return false;
|
|
}
|
|
else
|
|
return false;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const normal_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.mean() << __space << __x.stddev()
|
|
<< __space << __x._M_saved_available;
|
|
if (__x._M_saved_available)
|
|
__os << __space << __x._M_saved;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
normal_distribution<_RealType>& __x)
|
|
{
|
|
using param_type = typename normal_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
double __mean, __stddev;
|
|
bool __saved_avail;
|
|
if (__is >> __mean >> __stddev >> __saved_avail)
|
|
{
|
|
if (__saved_avail && (__is >> __x._M_saved))
|
|
{
|
|
__x._M_saved_available = __saved_avail;
|
|
__x.param(param_type(__mean, __stddev));
|
|
}
|
|
}
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
lognormal_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
*__f++ = std::exp(__p.s() * _M_nd(__urng) + __p.m());
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const lognormal_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.m() << __space << __x.s()
|
|
<< __space << __x._M_nd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
lognormal_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename lognormal_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __m, __s;
|
|
if (__is >> __m >> __s >> __x._M_nd)
|
|
__x.param(param_type(__m, __s));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::chi_squared_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
*__f++ = 2 * _M_gd(__urng);
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::chi_squared_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const typename
|
|
std::gamma_distribution<result_type>::param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
*__f++ = 2 * _M_gd(__urng, __p);
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const chi_squared_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.n() << __space << __x._M_gd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
chi_squared_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename chi_squared_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __n;
|
|
if (__is >> __n >> __x._M_gd)
|
|
__x.param(param_type(__n));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename cauchy_distribution<_RealType>::result_type
|
|
cauchy_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
_RealType __u;
|
|
do
|
|
__u = __aurng();
|
|
while (__u == 0.5);
|
|
|
|
const _RealType __pi = 3.1415926535897932384626433832795029L;
|
|
return __p.a() + __p.b() * std::tan(__pi * __u);
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
cauchy_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
const _RealType __pi = 3.1415926535897932384626433832795029L;
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
while (__f != __t)
|
|
{
|
|
_RealType __u;
|
|
do
|
|
__u = __aurng();
|
|
while (__u == 0.5);
|
|
|
|
*__f++ = __p.a() + __p.b() * std::tan(__pi * __u);
|
|
}
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const cauchy_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.a() << __space << __x.b();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
cauchy_distribution<_RealType>& __x)
|
|
{
|
|
using param_type = typename cauchy_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __a, __b;
|
|
if (__is >> __a >> __b)
|
|
__x.param(param_type(__a, __b));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::fisher_f_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
*__f++ = ((_M_gd_x(__urng) * n()) / (_M_gd_y(__urng) * m()));
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::fisher_f_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
typedef typename std::gamma_distribution<result_type>::param_type
|
|
param_type;
|
|
param_type __p1(__p.m() / 2);
|
|
param_type __p2(__p.n() / 2);
|
|
while (__f != __t)
|
|
*__f++ = ((_M_gd_x(__urng, __p1) * n())
|
|
/ (_M_gd_y(__urng, __p2) * m()));
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const fisher_f_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.m() << __space << __x.n()
|
|
<< __space << __x._M_gd_x << __space << __x._M_gd_y;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
fisher_f_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename fisher_f_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __m, __n;
|
|
if (__is >> __m >> __n >> __x._M_gd_x >> __x._M_gd_y)
|
|
__x.param(param_type(__m, __n));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::student_t_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
while (__f != __t)
|
|
*__f++ = _M_nd(__urng) * std::sqrt(n() / _M_gd(__urng));
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
std::student_t_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
typename std::gamma_distribution<result_type>::param_type
|
|
__p2(__p.n() / 2, 2);
|
|
while (__f != __t)
|
|
*__f++ = _M_nd(__urng) * std::sqrt(__p.n() / _M_gd(__urng, __p2));
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const student_t_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
student_t_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename student_t_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __n;
|
|
if (__is >> __n >> __x._M_nd >> __x._M_gd)
|
|
__x.param(param_type(__n));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
void
|
|
gamma_distribution<_RealType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
_M_malpha = _M_alpha < 1.0 ? _M_alpha + _RealType(1.0) : _M_alpha;
|
|
|
|
const _RealType __a1 = _M_malpha - _RealType(1.0) / _RealType(3.0);
|
|
_M_a2 = _RealType(1.0) / std::sqrt(_RealType(9.0) * __a1);
|
|
}
|
|
|
|
/**
|
|
* Marsaglia, G. and Tsang, W. W.
|
|
* "A Simple Method for Generating Gamma Variables"
|
|
* ACM Transactions on Mathematical Software, 26, 3, 363-372, 2000.
|
|
*/
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename gamma_distribution<_RealType>::result_type
|
|
gamma_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
|
|
result_type __u, __v, __n;
|
|
const result_type __a1 = (__param._M_malpha
|
|
- _RealType(1.0) / _RealType(3.0));
|
|
|
|
do
|
|
{
|
|
do
|
|
{
|
|
__n = _M_nd(__urng);
|
|
__v = result_type(1.0) + __param._M_a2 * __n;
|
|
}
|
|
while (__v <= 0.0);
|
|
|
|
__v = __v * __v * __v;
|
|
__u = __aurng();
|
|
}
|
|
while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
|
|
&& (std::log(__u) > (0.5 * __n * __n + __a1
|
|
* (1.0 - __v + std::log(__v)))));
|
|
|
|
if (__param.alpha() == __param._M_malpha)
|
|
return __a1 * __v * __param.beta();
|
|
else
|
|
{
|
|
do
|
|
__u = __aurng();
|
|
while (__u == 0.0);
|
|
|
|
return (std::pow(__u, result_type(1.0) / __param.alpha())
|
|
* __a1 * __v * __param.beta());
|
|
}
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
gamma_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
|
|
result_type __u, __v, __n;
|
|
const result_type __a1 = (__param._M_malpha
|
|
- _RealType(1.0) / _RealType(3.0));
|
|
|
|
if (__param.alpha() == __param._M_malpha)
|
|
while (__f != __t)
|
|
{
|
|
do
|
|
{
|
|
do
|
|
{
|
|
__n = _M_nd(__urng);
|
|
__v = result_type(1.0) + __param._M_a2 * __n;
|
|
}
|
|
while (__v <= 0.0);
|
|
|
|
__v = __v * __v * __v;
|
|
__u = __aurng();
|
|
}
|
|
while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
|
|
&& (std::log(__u) > (0.5 * __n * __n + __a1
|
|
* (1.0 - __v + std::log(__v)))));
|
|
|
|
*__f++ = __a1 * __v * __param.beta();
|
|
}
|
|
else
|
|
while (__f != __t)
|
|
{
|
|
do
|
|
{
|
|
do
|
|
{
|
|
__n = _M_nd(__urng);
|
|
__v = result_type(1.0) + __param._M_a2 * __n;
|
|
}
|
|
while (__v <= 0.0);
|
|
|
|
__v = __v * __v * __v;
|
|
__u = __aurng();
|
|
}
|
|
while (__u > result_type(1.0) - 0.0331 * __n * __n * __n * __n
|
|
&& (std::log(__u) > (0.5 * __n * __n + __a1
|
|
* (1.0 - __v + std::log(__v)))));
|
|
|
|
do
|
|
__u = __aurng();
|
|
while (__u == 0.0);
|
|
|
|
*__f++ = (std::pow(__u, result_type(1.0) / __param.alpha())
|
|
* __a1 * __v * __param.beta());
|
|
}
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const gamma_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.alpha() << __space << __x.beta()
|
|
<< __space << __x._M_nd;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
gamma_distribution<_RealType>& __x)
|
|
{
|
|
using param_type = typename gamma_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __alpha_val, __beta_val;
|
|
if (__is >> __alpha_val >> __beta_val >> __x._M_nd)
|
|
__x.param(param_type(__alpha_val, __beta_val));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename weibull_distribution<_RealType>::result_type
|
|
weibull_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
return __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
|
|
result_type(1) / __p.a());
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
weibull_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
auto __inv_a = result_type(1) / __p.a();
|
|
|
|
while (__f != __t)
|
|
*__f++ = __p.b() * std::pow(-std::log(result_type(1) - __aurng()),
|
|
__inv_a);
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const weibull_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.a() << __space << __x.b();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
weibull_distribution<_RealType>& __x)
|
|
{
|
|
using param_type = typename weibull_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __a, __b;
|
|
if (__is >> __a >> __b)
|
|
__x.param(param_type(__a, __b));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename extreme_value_distribution<_RealType>::result_type
|
|
extreme_value_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
return __p.a() - __p.b() * std::log(-std::log(result_type(1)
|
|
- __aurng()));
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
extreme_value_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __p)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
|
|
__aurng(__urng);
|
|
|
|
while (__f != __t)
|
|
*__f++ = __p.a() - __p.b() * std::log(-std::log(result_type(1)
|
|
- __aurng()));
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const extreme_value_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
__os << __x.a() << __space << __x.b();
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
extreme_value_distribution<_RealType>& __x)
|
|
{
|
|
using param_type
|
|
= typename extreme_value_distribution<_RealType>::param_type;
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
_RealType __a, __b;
|
|
if (__is >> __a >> __b)
|
|
__x.param(param_type(__a, __b));
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _IntType>
|
|
void
|
|
discrete_distribution<_IntType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
if (_M_prob.size() < 2)
|
|
{
|
|
_M_prob.clear();
|
|
return;
|
|
}
|
|
|
|
const double __sum = std::accumulate(_M_prob.begin(),
|
|
_M_prob.end(), 0.0);
|
|
// Now normalize the probabilites.
|
|
__detail::__normalize(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
|
|
__sum);
|
|
// Accumulate partial sums.
|
|
_M_cp.reserve(_M_prob.size());
|
|
std::partial_sum(_M_prob.begin(), _M_prob.end(),
|
|
std::back_inserter(_M_cp));
|
|
// Make sure the last cumulative probability is one.
|
|
_M_cp[_M_cp.size() - 1] = 1.0;
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _Func>
|
|
discrete_distribution<_IntType>::param_type::
|
|
param_type(size_t __nw, double __xmin, double __xmax, _Func __fw)
|
|
: _M_prob(), _M_cp()
|
|
{
|
|
const size_t __n = __nw == 0 ? 1 : __nw;
|
|
const double __delta = (__xmax - __xmin) / __n;
|
|
|
|
_M_prob.reserve(__n);
|
|
for (size_t __k = 0; __k < __nw; ++__k)
|
|
_M_prob.push_back(__fw(__xmin + __k * __delta + 0.5 * __delta));
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename discrete_distribution<_IntType>::result_type
|
|
discrete_distribution<_IntType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
if (__param._M_cp.empty())
|
|
return result_type(0);
|
|
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
const double __p = __aurng();
|
|
auto __pos = std::lower_bound(__param._M_cp.begin(),
|
|
__param._M_cp.end(), __p);
|
|
|
|
return __pos - __param._M_cp.begin();
|
|
}
|
|
|
|
template<typename _IntType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
discrete_distribution<_IntType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
|
|
if (__param._M_cp.empty())
|
|
{
|
|
while (__f != __t)
|
|
*__f++ = result_type(0);
|
|
return;
|
|
}
|
|
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
while (__f != __t)
|
|
{
|
|
const double __p = __aurng();
|
|
auto __pos = std::lower_bound(__param._M_cp.begin(),
|
|
__param._M_cp.end(), __p);
|
|
|
|
*__f++ = __pos - __param._M_cp.begin();
|
|
}
|
|
}
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const discrete_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<double>::max_digits10);
|
|
|
|
std::vector<double> __prob = __x.probabilities();
|
|
__os << __prob.size();
|
|
for (auto __dit = __prob.begin(); __dit != __prob.end(); ++__dit)
|
|
__os << __space << *__dit;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
namespace __detail
|
|
{
|
|
template<typename _ValT, typename _CharT, typename _Traits>
|
|
basic_istream<_CharT, _Traits>&
|
|
__extract_params(basic_istream<_CharT, _Traits>& __is,
|
|
vector<_ValT>& __vals, size_t __n)
|
|
{
|
|
__vals.reserve(__n);
|
|
while (__n--)
|
|
{
|
|
_ValT __val;
|
|
if (__is >> __val)
|
|
__vals.push_back(__val);
|
|
else
|
|
break;
|
|
}
|
|
return __is;
|
|
}
|
|
} // namespace __detail
|
|
|
|
template<typename _IntType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
discrete_distribution<_IntType>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
size_t __n;
|
|
if (__is >> __n)
|
|
{
|
|
std::vector<double> __prob_vec;
|
|
if (__detail::__extract_params(__is, __prob_vec, __n))
|
|
__x.param({__prob_vec.begin(), __prob_vec.end()});
|
|
}
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
void
|
|
piecewise_constant_distribution<_RealType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
if (_M_int.size() < 2
|
|
|| (_M_int.size() == 2
|
|
&& _M_int[0] == _RealType(0)
|
|
&& _M_int[1] == _RealType(1)))
|
|
{
|
|
_M_int.clear();
|
|
_M_den.clear();
|
|
return;
|
|
}
|
|
|
|
const double __sum = std::accumulate(_M_den.begin(),
|
|
_M_den.end(), 0.0);
|
|
|
|
__detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
|
|
__sum);
|
|
|
|
_M_cp.reserve(_M_den.size());
|
|
std::partial_sum(_M_den.begin(), _M_den.end(),
|
|
std::back_inserter(_M_cp));
|
|
|
|
// Make sure the last cumulative probability is one.
|
|
_M_cp[_M_cp.size() - 1] = 1.0;
|
|
|
|
for (size_t __k = 0; __k < _M_den.size(); ++__k)
|
|
_M_den[__k] /= _M_int[__k + 1] - _M_int[__k];
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _InputIteratorB, typename _InputIteratorW>
|
|
piecewise_constant_distribution<_RealType>::param_type::
|
|
param_type(_InputIteratorB __bbegin,
|
|
_InputIteratorB __bend,
|
|
_InputIteratorW __wbegin)
|
|
: _M_int(), _M_den(), _M_cp()
|
|
{
|
|
if (__bbegin != __bend)
|
|
{
|
|
for (;;)
|
|
{
|
|
_M_int.push_back(*__bbegin);
|
|
++__bbegin;
|
|
if (__bbegin == __bend)
|
|
break;
|
|
|
|
_M_den.push_back(*__wbegin);
|
|
++__wbegin;
|
|
}
|
|
}
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _Func>
|
|
piecewise_constant_distribution<_RealType>::param_type::
|
|
param_type(initializer_list<_RealType> __bl, _Func __fw)
|
|
: _M_int(), _M_den(), _M_cp()
|
|
{
|
|
_M_int.reserve(__bl.size());
|
|
for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
|
|
_M_int.push_back(*__biter);
|
|
|
|
_M_den.reserve(_M_int.size() - 1);
|
|
for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
|
|
_M_den.push_back(__fw(0.5 * (_M_int[__k + 1] + _M_int[__k])));
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _Func>
|
|
piecewise_constant_distribution<_RealType>::param_type::
|
|
param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
|
|
: _M_int(), _M_den(), _M_cp()
|
|
{
|
|
const size_t __n = __nw == 0 ? 1 : __nw;
|
|
const _RealType __delta = (__xmax - __xmin) / __n;
|
|
|
|
_M_int.reserve(__n + 1);
|
|
for (size_t __k = 0; __k <= __nw; ++__k)
|
|
_M_int.push_back(__xmin + __k * __delta);
|
|
|
|
_M_den.reserve(__n);
|
|
for (size_t __k = 0; __k < __nw; ++__k)
|
|
_M_den.push_back(__fw(_M_int[__k] + 0.5 * __delta));
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename piecewise_constant_distribution<_RealType>::result_type
|
|
piecewise_constant_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
const double __p = __aurng();
|
|
if (__param._M_cp.empty())
|
|
return __p;
|
|
|
|
auto __pos = std::lower_bound(__param._M_cp.begin(),
|
|
__param._M_cp.end(), __p);
|
|
const size_t __i = __pos - __param._M_cp.begin();
|
|
|
|
const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
|
|
|
|
return __param._M_int[__i] + (__p - __pref) / __param._M_den[__i];
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
piecewise_constant_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
if (__param._M_cp.empty())
|
|
{
|
|
while (__f != __t)
|
|
*__f++ = __aurng();
|
|
return;
|
|
}
|
|
|
|
while (__f != __t)
|
|
{
|
|
const double __p = __aurng();
|
|
|
|
auto __pos = std::lower_bound(__param._M_cp.begin(),
|
|
__param._M_cp.end(), __p);
|
|
const size_t __i = __pos - __param._M_cp.begin();
|
|
|
|
const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
|
|
|
|
*__f++ = (__param._M_int[__i]
|
|
+ (__p - __pref) / __param._M_den[__i]);
|
|
}
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const piecewise_constant_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
std::vector<_RealType> __int = __x.intervals();
|
|
__os << __int.size() - 1;
|
|
|
|
for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
|
|
__os << __space << *__xit;
|
|
|
|
std::vector<double> __den = __x.densities();
|
|
for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
|
|
__os << __space << *__dit;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
piecewise_constant_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
size_t __n;
|
|
if (__is >> __n)
|
|
{
|
|
std::vector<_RealType> __int_vec;
|
|
if (__detail::__extract_params(__is, __int_vec, __n + 1))
|
|
{
|
|
std::vector<double> __den_vec;
|
|
if (__detail::__extract_params(__is, __den_vec, __n))
|
|
{
|
|
__x.param({ __int_vec.begin(), __int_vec.end(),
|
|
__den_vec.begin() });
|
|
}
|
|
}
|
|
}
|
|
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _RealType>
|
|
void
|
|
piecewise_linear_distribution<_RealType>::param_type::
|
|
_M_initialize()
|
|
{
|
|
if (_M_int.size() < 2
|
|
|| (_M_int.size() == 2
|
|
&& _M_int[0] == _RealType(0)
|
|
&& _M_int[1] == _RealType(1)
|
|
&& _M_den[0] == _M_den[1]))
|
|
{
|
|
_M_int.clear();
|
|
_M_den.clear();
|
|
return;
|
|
}
|
|
|
|
double __sum = 0.0;
|
|
_M_cp.reserve(_M_int.size() - 1);
|
|
_M_m.reserve(_M_int.size() - 1);
|
|
for (size_t __k = 0; __k < _M_int.size() - 1; ++__k)
|
|
{
|
|
const _RealType __delta = _M_int[__k + 1] - _M_int[__k];
|
|
__sum += 0.5 * (_M_den[__k + 1] + _M_den[__k]) * __delta;
|
|
_M_cp.push_back(__sum);
|
|
_M_m.push_back((_M_den[__k + 1] - _M_den[__k]) / __delta);
|
|
}
|
|
|
|
// Now normalize the densities...
|
|
__detail::__normalize(_M_den.begin(), _M_den.end(), _M_den.begin(),
|
|
__sum);
|
|
// ... and partial sums...
|
|
__detail::__normalize(_M_cp.begin(), _M_cp.end(), _M_cp.begin(), __sum);
|
|
// ... and slopes.
|
|
__detail::__normalize(_M_m.begin(), _M_m.end(), _M_m.begin(), __sum);
|
|
|
|
// Make sure the last cumulative probablility is one.
|
|
_M_cp[_M_cp.size() - 1] = 1.0;
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _InputIteratorB, typename _InputIteratorW>
|
|
piecewise_linear_distribution<_RealType>::param_type::
|
|
param_type(_InputIteratorB __bbegin,
|
|
_InputIteratorB __bend,
|
|
_InputIteratorW __wbegin)
|
|
: _M_int(), _M_den(), _M_cp(), _M_m()
|
|
{
|
|
for (; __bbegin != __bend; ++__bbegin, ++__wbegin)
|
|
{
|
|
_M_int.push_back(*__bbegin);
|
|
_M_den.push_back(*__wbegin);
|
|
}
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _Func>
|
|
piecewise_linear_distribution<_RealType>::param_type::
|
|
param_type(initializer_list<_RealType> __bl, _Func __fw)
|
|
: _M_int(), _M_den(), _M_cp(), _M_m()
|
|
{
|
|
_M_int.reserve(__bl.size());
|
|
_M_den.reserve(__bl.size());
|
|
for (auto __biter = __bl.begin(); __biter != __bl.end(); ++__biter)
|
|
{
|
|
_M_int.push_back(*__biter);
|
|
_M_den.push_back(__fw(*__biter));
|
|
}
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _Func>
|
|
piecewise_linear_distribution<_RealType>::param_type::
|
|
param_type(size_t __nw, _RealType __xmin, _RealType __xmax, _Func __fw)
|
|
: _M_int(), _M_den(), _M_cp(), _M_m()
|
|
{
|
|
const size_t __n = __nw == 0 ? 1 : __nw;
|
|
const _RealType __delta = (__xmax - __xmin) / __n;
|
|
|
|
_M_int.reserve(__n + 1);
|
|
_M_den.reserve(__n + 1);
|
|
for (size_t __k = 0; __k <= __nw; ++__k)
|
|
{
|
|
_M_int.push_back(__xmin + __k * __delta);
|
|
_M_den.push_back(__fw(_M_int[__k] + __delta));
|
|
}
|
|
|
|
_M_initialize();
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _UniformRandomNumberGenerator>
|
|
typename piecewise_linear_distribution<_RealType>::result_type
|
|
piecewise_linear_distribution<_RealType>::
|
|
operator()(_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
|
|
__aurng(__urng);
|
|
|
|
const double __p = __aurng();
|
|
if (__param._M_cp.empty())
|
|
return __p;
|
|
|
|
auto __pos = std::lower_bound(__param._M_cp.begin(),
|
|
__param._M_cp.end(), __p);
|
|
const size_t __i = __pos - __param._M_cp.begin();
|
|
|
|
const double __pref = __i > 0 ? __param._M_cp[__i - 1] : 0.0;
|
|
|
|
const double __a = 0.5 * __param._M_m[__i];
|
|
const double __b = __param._M_den[__i];
|
|
const double __cm = __p - __pref;
|
|
|
|
_RealType __x = __param._M_int[__i];
|
|
if (__a == 0)
|
|
__x += __cm / __b;
|
|
else
|
|
{
|
|
const double __d = __b * __b + 4.0 * __a * __cm;
|
|
__x += 0.5 * (std::sqrt(__d) - __b) / __a;
|
|
}
|
|
|
|
return __x;
|
|
}
|
|
|
|
template<typename _RealType>
|
|
template<typename _ForwardIterator,
|
|
typename _UniformRandomNumberGenerator>
|
|
void
|
|
piecewise_linear_distribution<_RealType>::
|
|
__generate_impl(_ForwardIterator __f, _ForwardIterator __t,
|
|
_UniformRandomNumberGenerator& __urng,
|
|
const param_type& __param)
|
|
{
|
|
__glibcxx_function_requires(_ForwardIteratorConcept<_ForwardIterator>)
|
|
// We could duplicate everything from operator()...
|
|
while (__f != __t)
|
|
*__f++ = this->operator()(__urng, __param);
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_ostream<_CharT, _Traits>&
|
|
operator<<(std::basic_ostream<_CharT, _Traits>& __os,
|
|
const piecewise_linear_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_ostream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __os.flags();
|
|
const _CharT __fill = __os.fill();
|
|
const std::streamsize __precision = __os.precision();
|
|
const _CharT __space = __os.widen(' ');
|
|
__os.flags(__ios_base::scientific | __ios_base::left);
|
|
__os.fill(__space);
|
|
__os.precision(std::numeric_limits<_RealType>::max_digits10);
|
|
|
|
std::vector<_RealType> __int = __x.intervals();
|
|
__os << __int.size() - 1;
|
|
|
|
for (auto __xit = __int.begin(); __xit != __int.end(); ++__xit)
|
|
__os << __space << *__xit;
|
|
|
|
std::vector<double> __den = __x.densities();
|
|
for (auto __dit = __den.begin(); __dit != __den.end(); ++__dit)
|
|
__os << __space << *__dit;
|
|
|
|
__os.flags(__flags);
|
|
__os.fill(__fill);
|
|
__os.precision(__precision);
|
|
return __os;
|
|
}
|
|
|
|
template<typename _RealType, typename _CharT, typename _Traits>
|
|
std::basic_istream<_CharT, _Traits>&
|
|
operator>>(std::basic_istream<_CharT, _Traits>& __is,
|
|
piecewise_linear_distribution<_RealType>& __x)
|
|
{
|
|
using __ios_base = typename basic_istream<_CharT, _Traits>::ios_base;
|
|
|
|
const typename __ios_base::fmtflags __flags = __is.flags();
|
|
__is.flags(__ios_base::dec | __ios_base::skipws);
|
|
|
|
size_t __n;
|
|
if (__is >> __n)
|
|
{
|
|
vector<_RealType> __int_vec;
|
|
if (__detail::__extract_params(__is, __int_vec, __n + 1))
|
|
{
|
|
vector<double> __den_vec;
|
|
if (__detail::__extract_params(__is, __den_vec, __n + 1))
|
|
{
|
|
__x.param({ __int_vec.begin(), __int_vec.end(),
|
|
__den_vec.begin() });
|
|
}
|
|
}
|
|
}
|
|
__is.flags(__flags);
|
|
return __is;
|
|
}
|
|
|
|
|
|
template<typename _IntType>
|
|
seed_seq::seed_seq(std::initializer_list<_IntType> __il)
|
|
{
|
|
for (auto __iter = __il.begin(); __iter != __il.end(); ++__iter)
|
|
_M_v.push_back(__detail::__mod<result_type,
|
|
__detail::_Shift<result_type, 32>::__value>(*__iter));
|
|
}
|
|
|
|
template<typename _InputIterator>
|
|
seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
|
|
{
|
|
for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
|
|
_M_v.push_back(__detail::__mod<result_type,
|
|
__detail::_Shift<result_type, 32>::__value>(*__iter));
|
|
}
|
|
|
|
template<typename _RandomAccessIterator>
|
|
void
|
|
seed_seq::generate(_RandomAccessIterator __begin,
|
|
_RandomAccessIterator __end)
|
|
{
|
|
typedef typename iterator_traits<_RandomAccessIterator>::value_type
|
|
_Type;
|
|
|
|
if (__begin == __end)
|
|
return;
|
|
|
|
std::fill(__begin, __end, _Type(0x8b8b8b8bu));
|
|
|
|
const size_t __n = __end - __begin;
|
|
const size_t __s = _M_v.size();
|
|
const size_t __t = (__n >= 623) ? 11
|
|
: (__n >= 68) ? 7
|
|
: (__n >= 39) ? 5
|
|
: (__n >= 7) ? 3
|
|
: (__n - 1) / 2;
|
|
const size_t __p = (__n - __t) / 2;
|
|
const size_t __q = __p + __t;
|
|
const size_t __m = std::max(size_t(__s + 1), __n);
|
|
|
|
for (size_t __k = 0; __k < __m; ++__k)
|
|
{
|
|
_Type __arg = (__begin[__k % __n]
|
|
^ __begin[(__k + __p) % __n]
|
|
^ __begin[(__k - 1) % __n]);
|
|
_Type __r1 = __arg ^ (__arg >> 27);
|
|
__r1 = __detail::__mod<_Type,
|
|
__detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
|
|
_Type __r2 = __r1;
|
|
if (__k == 0)
|
|
__r2 += __s;
|
|
else if (__k <= __s)
|
|
__r2 += __k % __n + _M_v[__k - 1];
|
|
else
|
|
__r2 += __k % __n;
|
|
__r2 = __detail::__mod<_Type,
|
|
__detail::_Shift<_Type, 32>::__value>(__r2);
|
|
__begin[(__k + __p) % __n] += __r1;
|
|
__begin[(__k + __q) % __n] += __r2;
|
|
__begin[__k % __n] = __r2;
|
|
}
|
|
|
|
for (size_t __k = __m; __k < __m + __n; ++__k)
|
|
{
|
|
_Type __arg = (__begin[__k % __n]
|
|
+ __begin[(__k + __p) % __n]
|
|
+ __begin[(__k - 1) % __n]);
|
|
_Type __r3 = __arg ^ (__arg >> 27);
|
|
__r3 = __detail::__mod<_Type,
|
|
__detail::_Shift<_Type, 32>::__value>(1566083941u * __r3);
|
|
_Type __r4 = __r3 - __k % __n;
|
|
__r4 = __detail::__mod<_Type,
|
|
__detail::_Shift<_Type, 32>::__value>(__r4);
|
|
__begin[(__k + __p) % __n] ^= __r3;
|
|
__begin[(__k + __q) % __n] ^= __r4;
|
|
__begin[__k % __n] = __r4;
|
|
}
|
|
}
|
|
|
|
template<typename _RealType, size_t __bits,
|
|
typename _UniformRandomNumberGenerator>
|
|
_RealType
|
|
generate_canonical(_UniformRandomNumberGenerator& __urng)
|
|
{
|
|
static_assert(std::is_floating_point<_RealType>::value,
|
|
"template argument must be a floating point type");
|
|
|
|
const size_t __b
|
|
= std::min(static_cast<size_t>(std::numeric_limits<_RealType>::digits),
|
|
__bits);
|
|
const long double __r = static_cast<long double>(__urng.max())
|
|
- static_cast<long double>(__urng.min()) + 1.0L;
|
|
const size_t __log2r = std::log(__r) / std::log(2.0L);
|
|
const size_t __m = std::max<size_t>(1UL,
|
|
(__b + __log2r - 1UL) / __log2r);
|
|
_RealType __ret;
|
|
_RealType __sum = _RealType(0);
|
|
_RealType __tmp = _RealType(1);
|
|
for (size_t __k = __m; __k != 0; --__k)
|
|
{
|
|
__sum += _RealType(__urng() - __urng.min()) * __tmp;
|
|
__tmp *= __r;
|
|
}
|
|
__ret = __sum / __tmp;
|
|
if (__builtin_expect(__ret >= _RealType(1), 0))
|
|
{
|
|
#if _GLIBCXX_USE_C99_MATH_TR1
|
|
__ret = std::nextafter(_RealType(1), _RealType(0));
|
|
#else
|
|
__ret = _RealType(1)
|
|
- std::numeric_limits<_RealType>::epsilon() / _RealType(2);
|
|
#endif
|
|
}
|
|
return __ret;
|
|
}
|
|
|
|
_GLIBCXX_END_NAMESPACE_VERSION
|
|
} // namespace
|
|
|
|
#endif
|