random (mersenne_twister<>::seed()): Fix per tr1/5.1.4.2, p8.

2006-06-05  Paolo Carlini  <pcarlini@suse.de>

	* include/tr1/random (mersenne_twister<>::seed()): Fix per
	tr1/5.1.4.2, p8.
	* include/tr1/random.tcc (mod_w): Add.
	(mersenne_twister<>::seed(unsigned long)): Fix per tr1/5.1.4.2, p9.
	(mersenne_twister<>::seed(Gen&, false_type)): Adjust to use mod_w.
	* testsuite/tr1/5_numerical_facilies/random/mt19937.cc: Fix
	expected result per tr1/5.1.5, p2.
	* testsuite/tr1/5_numerical_facilies/random/mersenne_twister/
	cons/default.cc: Adjust.

	* include/tr1/random (exponential_distribution<>::operator()()): Fix.

From-SVN: r114412
This commit is contained in:
Paolo Carlini 2006-06-05 21:23:59 +00:00 committed by Paolo Carlini
parent 596475f036
commit 36ac3ed6b7
5 changed files with 58 additions and 54 deletions

View file

@ -1,3 +1,17 @@
2006-06-05 Paolo Carlini <pcarlini@suse.de>
* include/tr1/random (mersenne_twister<>::seed()): Fix per
tr1/5.1.4.2, p8.
* include/tr1/random.tcc (mod_w): Add.
(mersenne_twister<>::seed(unsigned long)): Fix per tr1/5.1.4.2, p9.
(mersenne_twister<>::seed(Gen&, false_type)): Adjust to use mod_w.
* testsuite/tr1/5_numerical_facilies/random/mt19937.cc: Fix
expected result per tr1/5.1.5, p2.
* testsuite/tr1/5_numerical_facilies/random/mersenne_twister/
cons/default.cc: Adjust.
* include/tr1/random (exponential_distribution<>::operator()()): Fix.
2006-06-05 Paolo Carlini <pcarlini@suse.de> 2006-06-05 Paolo Carlini <pcarlini@suse.de>
* include/tr1/random.tcc (Max::value()): Cast 1 to Tp(1) and * include/tr1/random.tcc (Max::value()): Cast 1 to Tp(1) and

View file

@ -493,7 +493,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
void void
seed() seed()
{ seed(0UL); } { seed(5489UL); }
void void
seed(unsigned long value); seed(unsigned long value);
@ -1560,7 +1560,7 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
template<class _UniformRandomNumberGenerator> template<class _UniformRandomNumberGenerator>
result_type result_type
operator()(_UniformRandomNumberGenerator& __urng) operator()(_UniformRandomNumberGenerator& __urng)
{ return std::log(__urng) / _M_lambda; } { return -std::log(__urng()) / _M_lambda; }
/** /**
* Inserts a %exponential_distribution random number distribution * Inserts a %exponential_distribution random number distribution

View file

@ -113,6 +113,11 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
{ return x; } { return x; }
}; };
template<typename _Tp, _Tp w>
inline _Tp
mod_w(_Tp x)
{ return Mod_w<_Tp, w, w == std::numeric_limits<_Tp>::digits>::calc(x); }
// Selector to return the maximum value possible that will fit in // Selector to return the maximum value possible that will fit in
// @p w bits of @p _Tp. // @p w bits of @p _Tp.
template<typename _Tp, _Tp w, bool> template<typename _Tp, _Tp w, bool>
@ -229,27 +234,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>:: mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
seed(unsigned long value) seed(unsigned long value)
{ {
if (value == 0) _M_x[0] = _Private::mod_w<_UInt, w>(value);
value = 4357;
#if 0
// @todo handle case numeric_limits<_UInt>::digits > 32
if (std::numeric_limits<_UInt>::digits > 32)
{
std::tr1::linear_congruential<unsigned long, 69069, 0, 2**32> lcg(value);
seed(lcg);
}
else
{
std::tr1::linear_congruential<unsigned long, 69069, 0, 0> lcg(value);
seed(lcg);
}
#else
std::tr1::linear_congruential<unsigned long, 69069, 0, 0> lcg(value);
seed(lcg);
#endif
}
for (int i = 1; i < n; ++i)
{
_UInt x = _M_x[i - 1];
x ^= x >> (w - 2);
x *= 1812433253ul;
x += i;
_M_x[i] = _Private::mod_w<_UInt, w>(x);
}
_M_p = n;
}
template<class _UInt, int w, int n, int m, int r, template<class _UInt, int w, int n, int m, int r,
_UInt a, int u, int s, _UInt a, int u, int s,
@ -259,13 +255,9 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>:: mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
seed(Gen& gen, false_type) seed(Gen& gen, false_type)
{ {
using _Private::Mod_w; for (int i = 0; i < n; ++i)
using std::numeric_limits; _M_x[i] = _Private::mod_w<_UInt, w>(gen());
_M_p = n;
for (int i = 0; i < state_size; ++i)
_M_x[i] = Mod_w<_UInt, w,
w == numeric_limits<_UInt>::digits>::calc(gen());
_M_p = state_size + 1;
} }
template<class _UInt, int w, int n, int m, int r, template<class _UInt, int w, int n, int m, int r,
@ -281,7 +273,6 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
return Max_w<_UInt, w, w == numeric_limits<_UInt>::digits>::value(); return Max_w<_UInt, w, w == numeric_limits<_UInt>::digits>::value();
} }
template<class _UInt, int w, int n, int m, int r, template<class _UInt, int w, int n, int m, int r,
_UInt a, int u, int s, _UInt a, int u, int s,
_UInt b, int t, _UInt c, int l> _UInt b, int t, _UInt c, int l>
@ -290,8 +281,8 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>:: mersenne_twister<_UInt, w, n, m, r, a, u, s, b, t, c, l>::
operator()() operator()()
{ {
// reload the vector - cost is O(n) amortized over n calls. // Reload the vector - cost is O(n) amortized over n calls.
if (_M_p >= state_size) if (_M_p >= n)
{ {
const _UInt upper_mask = (~_UInt()) << r; const _UInt upper_mask = (~_UInt()) << r;
const _UInt lower_mask = ~upper_mask; const _UInt lower_mask = ~upper_mask;
@ -311,14 +302,14 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
_M_p = 0; _M_p = 0;
} }
// Calculate x(i) // Calculate o(x(i)).
result_type y = _M_x[_M_p++]; result_type z = _M_x[_M_p++];
y ^= (y >> u); z ^= (z >> u);
y ^= (y << s) & b; z ^= (z << s) & b;
y ^= (y << t) & c; z ^= (z << t) & c;
y ^= (y >> l); z ^= (z >> l);
return y; return z;
} }
@ -329,15 +320,14 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
{ {
std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563> std::tr1::linear_congruential<unsigned long, 40014, 0, 2147483563>
lcg(__value); lcg(__value);
for (int i = 0; i < long_lag; ++i) for (int i = 0; i < long_lag; ++i)
_M_x[i] = _Private::mod<_IntType, 1, 0, modulus>(lcg()); _M_x[i] = _Private::mod<_IntType, 1, 0, modulus>(lcg());
_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
_M_p = 0; _M_p = 0;
} }
// //
// This implementation differs from the tr1 spec because the tr1 spec refused // This implementation differs from the tr1 spec because the tr1 spec refused
// to make any sense to me: the exponent of the factor in the spec goes from // to make any sense to me: the exponent of the factor in the spec goes from
@ -367,18 +357,18 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
_M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0; _M_carry = (_M_x[long_lag - 1] == 0) ? 1 : 0;
_M_p = 0; _M_p = 0;
} }
template<typename _IntType, _IntType __m, int __s, int __r> template<typename _IntType, _IntType __m, int __s, int __r>
typename subtract_with_carry<_IntType, __m, __s, __r>::result_type typename subtract_with_carry<_IntType, __m, __s, __r>::result_type
subtract_with_carry<_IntType, __m, __s, __r>:: subtract_with_carry<_IntType, __m, __s, __r>::
operator()() operator()()
{ {
// derive short lag index from current index // Derive short lag index from current index.
int ps = _M_p - short_lag; int ps = _M_p - short_lag;
if (ps < 0) ps += long_lag; if (ps < 0)
ps += long_lag;
// calculate new x(i) without overflow or division
// Calculate new x(i) without overflow or division.
_IntType xi; _IntType xi;
if (_M_x[ps] >= _M_x[_M_p] + _M_carry) if (_M_x[ps] >= _M_x[_M_p] + _M_carry)
{ {
@ -391,15 +381,15 @@ _GLIBCXX_BEGIN_NAMESPACE(tr1)
_M_carry = 1; _M_carry = 1;
} }
_M_x[_M_p++] = xi; _M_x[_M_p++] = xi;
// adjust current index to loop around in ring buffer // Adjust current index to loop around in ring buffer.
if (_M_p >= long_lag) if (_M_p >= long_lag)
_M_p = 0; _M_p = 0;
return xi; return xi;
} }
template<class _E, int __p, int __r> template<class _E, int __p, int __r>
typename discard_block<_E, __p, __r>::result_type typename discard_block<_E, __p, __r>::result_type
discard_block<_E, __p, __r>:: discard_block<_E, __p, __r>::

View file

@ -38,7 +38,7 @@ test01()
VERIFY( x.min() == 0 ); VERIFY( x.min() == 0 );
VERIFY( x.max() == 4294967295ul ); VERIFY( x.max() == 4294967295ul );
VERIFY( x() == 4290933890ul ); VERIFY( x() == 3499211612ul );
} }
int main() int main()

View file

@ -33,7 +33,7 @@ test01()
for (int i = 0; i < 9999; ++i) for (int i = 0; i < 9999; ++i)
a(); a();
VERIFY( a() == 3346425566ul ); VERIFY( a() == 4123659995ul );
} }
int main() int main()