// random number generation (out of line) -*- C++ -*-
-// Copyright (C) 2009 Free Software Foundation, Inc.
+// Copyright (C) 2009, 2010, 2011, 2012 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
/** @file bits/random.tcc
* This is an internal header file, included by other library headers.
- * You should not attempt to use it directly.
+ * Do not attempt to use it directly. @headername{random}
*/
-#include <numeric>
-#include <algorithm>
+#ifndef _RANDOM_TCC
+#define _RANDOM_TCC 1
-namespace std
+#include <numeric> // std::accumulate and std::partial_sum
+
+namespace std _GLIBCXX_VISIBILITY(default)
{
/*
* (Further) implementation-space details.
*/
namespace __detail
{
+ _GLIBCXX_BEGIN_NAMESPACE_VERSION
+
// General case for x = (ax + c) mod m -- use Schrage's algorithm to
// avoid integer overflow.
//
//
// Preconditions: a > 0, m > 0.
//
- template<typename _Tp, _Tp __a, _Tp __c, _Tp __m, bool>
+ // XXX FIXME: as-is, only works correctly for __m % __a < __m / __a.
+ //
+ template<typename _Tp, _Tp __m, _Tp __a, _Tp __c, bool>
struct _Mod
{
static _Tp
// Special case for m == 0 -- use unsigned integer overflow as modulo
// operator.
- template<typename _Tp, _Tp __a, _Tp __c, _Tp __m>
- struct _Mod<_Tp, __a, __c, __m, true>
+ template<typename _Tp, _Tp __m, _Tp __a, _Tp __c>
+ struct _Mod<_Tp, __m, __a, __c, true>
{
static _Tp
__calc(_Tp __x)
{ return __a * __x + __c; }
};
+
+ template<typename _InputIterator, typename _OutputIterator,
+ typename _UnaryOperation>
+ _OutputIterator
+ __transform(_InputIterator __first, _InputIterator __last,
+ _OutputIterator __result, _UnaryOperation __unary_op)
+ {
+ for (; __first != __last; ++__first, ++__result)
+ *__result = __unary_op(*__first);
+ return __result;
+ }
+
+ _GLIBCXX_END_NAMESPACE_VERSION
} // namespace __detail
+_GLIBCXX_BEGIN_NAMESPACE_VERSION
+
+ 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.
linear_congruential_engine<_UIntType, __a, __c, __m>::
seed(result_type __s)
{
- if ((__detail::__mod<_UIntType, 1U, 0U, __m>(__c) == 0U)
- && (__detail::__mod<_UIntType, 1U, 0U, __m>(__s) == 0U))
- _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(1U);
+ if ((__detail::__mod<_UIntType, __m>(__c) == 0)
+ && (__detail::__mod<_UIntType, __m>(__s) == 0))
+ _M_x = 1;
else
- _M_x = __detail::__mod<_UIntType, 1U, 0U, __m>(__s);
+ _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>
- void
- linear_congruential_engine<_UIntType, __a, __c, __m>::
- seed(seed_seq& __q)
- {
- const _UIntType __k0 = __m == 0 ? std::numeric_limits<_UIntType>::digits
- : std::__lg(__m);
- const _UIntType __k = (__k0 + 31) / 32;
- _UIntType __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 _Sseq>
+ typename std::enable_if<std::is_class<_Sseq>::value>::type
+ linear_congruential_engine<_UIntType, __a, __c, __m>::
+ seed(_Sseq& __q)
+ {
+ 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>
_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, 1, 0,
+ _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 += __i;
- _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
+ __x += __detail::__mod<_UIntType, __n>(__i);
+ _M_x[__i] = __detail::__mod<_UIntType,
__detail::_Shift<_UIntType, __w>::__value>(__x);
}
_M_p = state_size;
_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,
+ template<typename _Sseq>
+ typename std::enable_if<std::is_class<_Sseq>::value>::type
+ mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
__s, __b, __t, __c, __l, __f>::
- seed(seed_seq& __q)
- {
- const _UIntType __upper_mask = (~_UIntType()) << __r;
- const size_t __k = (__w + 31) / 32;
- _UIntType __arr[__k * __n];
- __q.generate(__arr + 0, __arr + __k * __n);
-
- 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[__i * __k + __j] * __factor;
- __factor *= __detail::_Shift<_UIntType, 32>::__value;
- }
- _M_x[__i] = __detail::__mod<_UIntType, 1U, 0U,
- __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;
- }
- }
+ seed(_Sseq& __q)
+ {
+ 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 - 1U>::__value;
- }
+ _M_x[0] = __detail::_Shift<_UIntType, __w - 1>::__value;
+ }
template<typename _UIntType, size_t __w,
size_t __n, size_t __m, size_t __r,
__os.flags(__ios_base::dec | __ios_base::fixed | __ios_base::left);
__os.fill(__space);
- for (size_t __i = 0; __i < __n - 1; ++__i)
+ for (size_t __i = 0; __i < __n; ++__i)
__os << __x._M_x[__i] << __space;
- __os << __x._M_x[__n - 1];
+ __os << __x._M_p;
__os.flags(__flags);
__os.fill(__fill);
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)
{
- if (__value == 0)
- __value = default_seed;
+ std::linear_congruential_engine<result_type, 40014u, 0u, 2147483563u>
+ __lcg(__value == 0u ? default_seed : __value);
- std::linear_congruential_engine<result_type, 40014U, 0U, 2147483563U>
- __lcg(__value);
-
- // I hope this is right. The "10000" tests work for the ranluxen.
- const size_t __n = (word_size + 31) / 32;
+ const size_t __n = (__w + 31) / 32;
for (size_t __i = 0; __i < long_lag; ++__i)
{
- _UIntType __sum = 0U;
- _UIntType __factor = 1U;
+ _UIntType __sum = 0u;
+ _UIntType __factor = 1u;
for (size_t __j = 0; __j < __n; ++__j)
{
- __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
+ __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, 1, 0,
+ _M_x[__i] = __detail::__mod<_UIntType,
__detail::_Shift<_UIntType, __w>::__value>(__sum);
}
_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
}
template<typename _UIntType, size_t __w, size_t __s, size_t __r>
- void
- subtract_with_carry_engine<_UIntType, __w, __s, __r>::
- seed(seed_seq& __q)
- {
- const size_t __n = (word_size + 31) / 32;
- _UIntType __arr[long_lag + __n];
- __q.generate(__arr + 0, __arr + long_lag + __n);
+ template<typename _Sseq>
+ typename std::enable_if<std::is_class<_Sseq>::value>::type
+ subtract_with_carry_engine<_UIntType, __w, __s, __r>::
+ seed(_Sseq& __q)
+ {
+ 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 < __n; ++__j)
- {
- __sum += __detail::__mod<__detail::_UInt32Type, 1, 0, 0>
- (__arr[__i * __n + __j]) * __factor;
- __factor *= __detail::_Shift<_UIntType, 32>::__value;
- }
- _M_x[__i] = __detail::__mod<_UIntType, 1, 0,
- __detail::_Shift<_UIntType, __w>::__value>(__sum);
- }
- _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
- _M_p = 0;
- }
+ 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>::
for (size_t __i = 0; __i < __r; ++__i)
__os << __x._M_x[__i] << __space;
- __os << __x._M_carry;
+ __os << __x._M_carry << __space << __x._M_p;
__os.flags(__flags);
__os.fill(__fill);
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>::
independent_bits_engine<_RandomNumberEngine, __w, _UIntType>::
operator()()
{
- const long double __r = static_cast<long double>(_M_b.max())
- - static_cast<long double>(_M_b.min()) + 1.0L;
- const result_type __m = std::log(__r) / std::log(2.0L);
- result_type __n, __n0, __y0, __y1, __s0, __s1;
+ 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 result_type __w0 = __w / __n;
- const result_type __w1 = __w0 + 1;
- __s0 = result_type(1) << __w0;
- __s1 = result_type(1) << __w1;
- __y0 = __s0 * (__r / __s0);
- __y1 = __s1 * (__r / __s1);
- if (__r - __y0 <= __y0 / __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)
{
- result_type __u;
+ __ctype __u;
do
__u = _M_b() - _M_b.min();
- while (__u >= __y0);
- __sum = __s0 * __sum + __u % __s0;
+ while (__y0 && __u >= __y0);
+ __sum = __s0 * __sum + (__s0 ? __u % __s0 : __u);
}
for (size_t __k = __n0; __k < __n; ++__k)
{
- result_type __u;
+ __ctype __u;
do
__u = _M_b() - _M_b.min();
- while (__u >= __y1);
- __sum = __s1 * __sum + __u % __s1;
+ 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()()
template<typename _UniformRandomNumberGenerator>
typename uniform_int_distribution<_IntType>::result_type
uniform_int_distribution<_IntType>::
- _M_call(_UniformRandomNumberGenerator& __urng,
- result_type __min, result_type __max, true_type)
+ operator()(_UniformRandomNumberGenerator& __urng,
+ const param_type& __param)
{
- // XXX Must be fixed to work well for *arbitrary* __urng.max(),
- // __urng.min(), __max, __min. Currently works fine only in the
- // most common case __urng.max() - __urng.min() >= __max - __min,
- // with __urng.max() > __urng.min() >= 0.
- typedef typename __gnu_cxx::__add_unsigned<typename
- _UniformRandomNumberGenerator::result_type>::__type __urntype;
- typedef typename __gnu_cxx::__add_unsigned<result_type>::__type
- __utype;
- typedef typename __gnu_cxx::__conditional_type<(sizeof(__urntype)
- > sizeof(__utype)),
- __urntype, __utype>::__type __uctype;
+ typedef typename _UniformRandomNumberGenerator::result_type
+ _Gresult_type;
+ typedef typename std::make_unsigned<result_type>::type __utype;
+ typedef typename std::common_type<_Gresult_type, __utype>::type
+ __uctype;
- result_type __ret;
+ const __uctype __urngmin = __urng.min();
+ const __uctype __urngmax = __urng.max();
+ const __uctype __urngrange = __urngmax - __urngmin;
+ const __uctype __urange
+ = __uctype(__param.b()) - __uctype(__param.a());
- const __urntype __urnmin = __urng.min();
- const __urntype __urnmax = __urng.max();
- const __urntype __urnrange = __urnmax - __urnmin;
- const __uctype __urange = __max - __min;
- const __uctype __udenom = (__urnrange <= __urange
- ? 1 : __urnrange / (__urange + 1));
- do
- __ret = (__urntype(__urng()) - __urnmin) / __udenom;
- while (__ret > __max - __min);
+ __uctype __ret;
+
+ if (__urngrange > __urange)
+ {
+ // downscaling
+ const __uctype __uerange = __urange + 1; // __urange can be zero
+ const __uctype __scaling = __urngrange / __uerange;
+ const __uctype __past = __uerange * __scaling;
+ do
+ __ret = __uctype(__urng()) - __urngmin;
+ while (__ret >= __past);
+ __ret /= __scaling;
+ }
+ else if (__urngrange < __urange)
+ {
+ // upscaling
+ /*
+ Note that every value in [0, urange]
+ can be written uniquely as
+
+ (urngrange + 1) * high + low
+
+ where
+
+ high in [0, urange / (urngrange + 1)]
+
+ and
+
+ low in [0, urngrange].
+ */
+ __uctype __tmp; // wraparound control
+ do
+ {
+ const __uctype __uerngrange = __urngrange + 1;
+ __tmp = (__uerngrange * operator()
+ (__urng, param_type(0, __urange / __uerngrange)));
+ __ret = __tmp + (__uctype(__urng()) - __urngmin);
+ }
+ while (__ret > __urange || __ret < __tmp);
+ }
+ else
+ __ret = __uctype(__urng()) - __urngmin;
- return __ret + __min;
+ return __ret + __param.a();
}
template<typename _IntType, typename _CharT, typename _Traits>
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.a() << __space << __x.b();
const std::streamsize __precision = __os.precision();
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__os.widen(' '));
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
__os << __x.p();
// The largest _RealType convertible to _IntType.
const double __thr =
std::numeric_limits<_IntType>::max() + __naf;
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+ __detail::_Adaptor<_UniformRandomNumberGenerator, double>
__aurng(__urng);
double __cand;
do
- __cand = std::ceil(std::log(__aurng()) / __param._M_log_p);
+ __cand = std::floor(std::log(__aurng()) / __param._M_log_1_p);
while (__cand >= __thr);
return result_type(__cand + __naf);
const std::streamsize __precision = __os.precision();
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__os.widen(' '));
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
__os << __x.p();
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
param_type;
const double __y =
- _M_gd(__urng, param_type(__p.k(), __p.p() / (1.0 - __p.p())));
+ _M_gd(__urng, param_type(__p.k(), (1.0 - __p.p()) / __p.p()));
std::poisson_distribution<result_type> __poisson(__y);
return __poisson(__urng);
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__os.widen(' '));
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
__os << __x.k() << __space << __x.p()
<< __space << __x._M_gd;
* is defined.
*
* Reference:
- * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+ * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
* New York, 1986, Ch. X, Sects. 3.3 & 3.4 (+ Errata!).
*/
template<typename _IntType>
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
__os << __x.mean() << __space << __x._M_nd;
* is defined.
*
* Reference:
- * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+ * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
* New York, 1986, Ch. X, Sect. 4 (+ Errata!).
*/
template<typename _IntType>
{
result_type __ret;
const _IntType __t = __param.t();
- const _IntType __p = __param.p();
+ const double __p = __param.p();
const double __p12 = __p <= 0.5 ? __p : 1.0 - __p;
__detail::_Adaptor<_UniformRandomNumberGenerator, double>
__aurng(__urng);
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
__os << __x.t() << __space << __x.p()
<< __space << __x._M_nd;
const std::streamsize __precision = __os.precision();
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__os.widen(' '));
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.lambda();
/**
* Polar method due to Marsaglia.
*
- * Devroye, L. "Non-Uniform Random Variates Generation." Springer-Verlag,
+ * Devroye, L. Non-Uniform Random Variates Generation. Springer-Verlag,
* New York, 1986, Ch. V, Sect. 4.4.
*/
template<typename _RealType>
return __ret;
}
+ 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 _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.mean() << __space << __x.stddev()
<< __space << __x._M_saved_available;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.m() << __space << __x.s()
<< __space << __x._M_nd;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.n() << __space << __x._M_gd;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.a() << __space << __x.b();
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.m() << __space << __x.n()
<< __space << __x._M_gd_x << __space << __x._M_gd_y;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.n() << __space << __x._M_nd << __space << __x._M_gd;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.alpha() << __space << __x.beta()
<< __space << __x._M_nd;
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.a() << __space << __x.b();
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
__os << __x.a() << __space << __x.b();
if (_M_prob.size() < 2)
{
_M_prob.clear();
- _M_prob.push_back(1.0);
return;
}
const double __sum = std::accumulate(_M_prob.begin(),
_M_prob.end(), 0.0);
// Now normalize the probabilites.
- std::transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
- std::bind2nd(std::divides<double>(), __sum));
+ __detail::__transform(_M_prob.begin(), _M_prob.end(), _M_prob.begin(),
+ std::bind2nd(std::divides<double>(), __sum));
// Accumulate partial sums.
_M_cp.reserve(_M_prob.size());
std::partial_sum(_M_prob.begin(), _M_prob.end(),
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+ if (__param._M_cp.empty())
+ return result_type(0);
+
+ __detail::_Adaptor<_UniformRandomNumberGenerator, double>
__aurng(__urng);
const double __p = __aurng();
const _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<double>::digits10 + 1);
+ __os.precision(std::numeric_limits<double>::max_digits10);
std::vector<double> __prob = __x.probabilities();
__os << __prob.size();
piecewise_constant_distribution<_RealType>::param_type::
_M_initialize()
{
- if (_M_int.size() < 2)
+ if (_M_int.size() < 2
+ || (_M_int.size() == 2
+ && _M_int[0] == _RealType(0)
+ && _M_int[1] == _RealType(1)))
{
_M_int.clear();
- _M_int.reserve(2);
- _M_int.push_back(_RealType(0));
- _M_int.push_back(_RealType(1));
-
_M_den.clear();
- _M_den.push_back(1.0);
-
return;
}
const double __sum = std::accumulate(_M_den.begin(),
_M_den.end(), 0.0);
- std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
- std::bind2nd(std::divides<double>(), __sum));
+ __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
+ std::bind2nd(std::divides<double>(), __sum));
_M_cp.reserve(_M_den.size());
std::partial_sum(_M_den.begin(), _M_den.end(),
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+ __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 _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
std::vector<_RealType> __int = __x.intervals();
__os << __int.size() - 1;
piecewise_linear_distribution<_RealType>::param_type::
_M_initialize()
{
- if (_M_int.size() < 2)
+ 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_int.reserve(2);
- _M_int.push_back(_RealType(0));
- _M_int.push_back(_RealType(1));
-
_M_den.clear();
- _M_den.reserve(2);
- _M_den.push_back(1.0);
- _M_den.push_back(1.0);
-
return;
}
}
// Now normalize the densities...
- std::transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
- std::bind2nd(std::divides<double>(), __sum));
+ __detail::__transform(_M_den.begin(), _M_den.end(), _M_den.begin(),
+ std::bind2nd(std::divides<double>(), __sum));
// ... and partial sums...
- std::transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
- std::bind2nd(std::divides<double>(), __sum));
+ __detail::__transform(_M_cp.begin(), _M_cp.end(), _M_cp.begin(),
+ std::bind2nd(std::divides<double>(), __sum));
// ... and slopes.
- std::transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
- std::bind2nd(std::divides<double>(), __sum));
+ __detail::__transform(_M_m.begin(), _M_m.end(), _M_m.begin(),
+ std::bind2nd(std::divides<double>(), __sum));
// Make sure the last cumulative probablility is one.
_M_cp[_M_cp.size() - 1] = 1.0;
}
operator()(_UniformRandomNumberGenerator& __urng,
const param_type& __param)
{
- __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
+ __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 _CharT __space = __os.widen(' ');
__os.flags(__ios_base::scientific | __ios_base::left);
__os.fill(__space);
- __os.precision(std::numeric_limits<_RealType>::digits10 + 1);
+ __os.precision(std::numeric_limits<_RealType>::max_digits10);
std::vector<_RealType> __int = __x.intervals();
__os << __int.size() - 1;
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, 1, 0,
+ _M_v.push_back(__detail::__mod<result_type,
__detail::_Shift<result_type, 32>::__value>(*__iter));
}
seed_seq::seed_seq(_InputIterator __begin, _InputIterator __end)
{
for (_InputIterator __iter = __begin; __iter != __end; ++__iter)
- _M_v.push_back(__detail::__mod<result_type, 1, 0,
+ _M_v.push_back(__detail::__mod<result_type,
__detail::_Shift<result_type, 32>::__value>(*__iter));
}
if (__begin == __end)
return;
- std::fill(__begin, __end, _Type(0x8b8b8b8bU));
+ std::fill(__begin, __end, _Type(0x8b8b8b8bu));
const size_t __n = __end - __begin;
const size_t __s = _M_v.size();
: (__n - 1) / 2;
const size_t __p = (__n - __t) / 2;
const size_t __q = __p + __t;
- const size_t __m = std::max(__s + 1, __n);
+ 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, 1664525U, 0U,
- __detail::_Shift<_Type, 32>::__value>(__r1);
+ _Type __r1 = __arg ^ (__arg >> 27);
+ __r1 = __detail::__mod<_Type,
+ __detail::_Shift<_Type, 32>::__value>(1664525u * __r1);
_Type __r2 = __r1;
if (__k == 0)
__r2 += __s;
__r2 += __k % __n + _M_v[__k - 1];
else
__r2 += __k % __n;
- __r2 = __detail::__mod<_Type, 1U, 0U,
- __detail::_Shift<_Type, 32>::__value>(__r2);
+ __r2 = __detail::__mod<_Type,
+ __detail::_Shift<_Type, 32>::__value>(__r2);
__begin[(__k + __p) % __n] += __r1;
__begin[(__k + __q) % __n] += __r2;
__begin[__k % __n] = __r2;
_Type __arg = (__begin[__k % __n]
+ __begin[(__k + __p) % __n]
+ __begin[(__k - 1) % __n]);
- _Type __r3 = __arg ^ (__arg << 27);
- __r3 = __detail::__mod<_Type, 1566083941U, 0U,
- __detail::_Shift<_Type, 32>::__value>(__r3);
+ _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, 1U, 0U,
- __detail::_Shift<_Type, 32>::__value>(__r4);
- __begin[(__k + __p) % __n] ^= __r4;
- __begin[(__k + __q) % __n] ^= __r3;
+ __r4 = __detail::__mod<_Type,
+ __detail::_Shift<_Type, 32>::__value>(__r4);
+ __begin[(__k + __p) % __n] ^= __r3;
+ __begin[(__k + __q) % __n] ^= __r4;
__begin[__k % __n] = __r4;
}
}
}
return __sum / __tmp;
}
-}
+
+_GLIBCXX_END_NAMESPACE_VERSION
+} // namespace
+
+#endif
OSDN Git Service
