boost/random
2
0
Fork 0
mirror of https://github.com/boostorg/random.git synced 2026-01-19 04:22:17 +00:00
Code Issues Projects Releases Wiki Activity

This commit was manufactured by cvs2svn to create tag

Browse Source 
'perforce_2_4_merge_1'.

[SVN r13112]
This commit is contained in:
nobody
2002-03-06 14:13:30 +00:00
parent eb910ff396
commit 9e0bca7c70
38 changed files with 0 additions and 7530 deletions
Download Patch File Download Diff File Expand all files Collapse all files

147
histogram.cpp
View File

@@ -1,147 +0,0 @@
/* boost histogram.cpp graphical verification of distribution functions
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* This test program allows to visibly examine the results of the
* distribution functions.
*/
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <cmath>
#include <string>
#include <boost/random.hpp>
void plot_histogram(const std::vector<int>& slots, int samples,
double from, double to)
{
int m = *std::max_element(slots.begin(), slots.end());
const int nRows = 20;
std::cout.setf(std::ios::fixed|std::ios::left);
std::cout.precision(5);
for(int r = 0; r < nRows; r++) {
double y = ((nRows - r) * double(m))/(nRows * samples);
std::cout << std::setw(10) << y << " ";
for(unsigned int col = 0; col < slots.size(); col++) {
char out = ' ';
if(slots[col]/double(samples) >= y)
out = 'x';
std::cout << out;
}
std::cout << std::endl;
}
std::cout << std::setw(12) << " "
<< std::setw(10) << from;
std::cout.setf(std::ios::right, std::ios::adjustfield);
std::cout << std::setw(slots.size()-10) << to << std::endl;
}
// I am not sure whether these two should be in the library as well
// maintain sum of NumberGenerator results
template<class NumberGenerator,
class Sum = typename NumberGenerator::result_type>
class sum_result
{
public:
typedef NumberGenerator base_type;
typedef typename base_type::result_type result_type;
explicit sum_result(const base_type & g) : gen(g), _sum(0) { }
result_type operator()() { result_type r = gen(); _sum += r; return r; }
base_type & base() { return gen; }
Sum sum() const { return _sum; }
void reset() { _sum = 0; }
private:
base_type gen;
Sum _sum;
};
// maintain square sum of NumberGenerator results
template<class NumberGenerator,
class Sum = typename NumberGenerator::result_type>
class squaresum_result
{
public:
typedef NumberGenerator base_type;
typedef typename base_type::result_type result_type;
explicit squaresum_result(const base_type & g) : gen(g), _sum(0) { }
result_type operator()() { result_type r = gen(); _sum += r*r; return r; }
base_type & base() { return gen; }
Sum squaresum() const { return _sum; }
void reset() { _sum = 0; }
private:
base_type gen;
Sum _sum;
};
template<class RNG>
void histogram(RNG base, int samples, double from, double to,
const std::string & name)
{
typedef squaresum_result<sum_result<RNG, double>, double > SRNG;
SRNG gen((sum_result<RNG, double>(base)));
const int nSlots = 60;
std::vector<int> slots(nSlots,0);
for(int i = 0; i < samples; i++) {
double val = gen();
if(val < from || val >= to) // early check avoids overflow
continue;
int slot = int((val-from)/(to-from) * nSlots);
if(slot < 0 || slot > (int)slots.size())
continue;
slots[slot]++;
}
std::cout << name << std::endl;
plot_histogram(slots, samples, from, to);
double mean = gen.base().sum() / samples;
std::cout << "mean: " << mean
<< " sigma: " << std::sqrt(gen.squaresum()/samples-mean*mean)
<< "\n" << std::endl;
}
template<class PRNG>
void histograms()
{
PRNG rng;
using namespace boost;
histogram(uniform_01<PRNG>(rng), 100000, -0.5, 1.5, "uniform_01(0,1)");
histogram(bernoulli_distribution<PRNG>(rng, 0.2), 100000, -0.5, 1.5,
"bernoulli(0.2)");
histogram(triangle_distribution<PRNG>(rng, 1, 2, 8), 100000, 0, 10,
"triangle(1,2,8)");
histogram(geometric_distribution<PRNG>(rng, 5.0/6.0), 100000, 0, 10,
"geometric(5/6)");
histogram(exponential_distribution<PRNG>(rng, 0.3), 100000, 0, 10,
"exponential(0.3)");
histogram(cauchy_distribution<PRNG>(rng), 100000, -5, 5,
"cauchy");
histogram(lognormal_distribution<PRNG>(rng, 3, 2), 100000, 0, 10,
"lognormal");
histogram(normal_distribution<PRNG>(rng), 100000, -3, 3,
"normal");
histogram(normal_distribution<PRNG>(rng, 0.5, 0.5), 100000, -3, 3,
"normal(0.5, 0.5)");
}
int main()
{
histograms<boost::mt19937>();
}

94
include/boost/random/additive_combine.hpp
View File

@@ -1,94 +0,0 @@
/* boost random/additive_combine.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_ADDITIVE_COMBINE_HPP
#define BOOST_RANDOM_ADDITIVE_COMBINE_HPP
#include <iostream>
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/random/linear_congruential.hpp>
namespace boost {
namespace random {
// L'Ecuyer 1988
template<class MLCG1, class MLCG2,
#ifndef BOOST_NO_DEPENDENT_TYPES_IN_TEMPLATE_VALUE_PARAMETERS
typename MLCG1::result_type
#else
int32_t
#endif
val>
class additive_combine
{
public:
typedef MLCG1 first_base;
typedef MLCG2 second_base;
typedef typename MLCG1::result_type result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = 1;
static const result_type max_value = MLCG1::max_value-1;
#else
enum { has_fixed_range = false };
#endif
result_type min() const { return 1; }
result_type max() const { return _mlcg1.max()-1; }
additive_combine() : _mlcg1(), _mlcg2() { }
additive_combine(typename MLCG1::result_type seed1,
typename MLCG2::result_type seed2)
: _mlcg1(seed1), _mlcg2(seed2) { }
result_type operator()() {
result_type z = _mlcg1() - _mlcg2();
if(z < 1)
z += MLCG1::modulus-1;
return z;
}
bool validation(result_type x) const { return val == x; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, const additive_combine& r)
{ os << r._mlcg1 << " " << r._mlcg2; return os; }
friend std::istream& operator>>(std::istream& is, additive_combine& r)
{ is >> r._mlcg1 >> std::ws >> r._mlcg2; return is; }
friend bool operator==(const additive_combine& x, const additive_combine& y)
{ return x._mlcg1 == y._mlcg1 && x._mlcg2 == y._mlcg2; }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const additive_combine& rhs) const
{ return _mlcg1 == rhs._mlcg1 && _mlcg2 == rhs._mlcg2; }
#endif
private:
MLCG1 _mlcg1;
MLCG2 _mlcg2;
};
} // namespace random
typedef random::additive_combine<
random::linear_congruential<int32_t, 40014, 0, 2147483563, 0>,
random::linear_congruential<int32_t, 40692, 0, 2147483399, 0>,
2060321752> ecuyer1988;
} // namespace boost
#endif // BOOST_RANDOM_ADDITIVE_COMBINE_HPP

64
include/boost/random/bernoulli_distribution.hpp
View File

@@ -1,64 +0,0 @@
/* boost random/bernoulli_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP
#define BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP
#include <cassert>
namespace boost {
// Bernoulli distribution: p(true) = p, p(false) = 1-p (boolean)
template<class UniformRandomNumberGenerator>
class bernoulli_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef bool result_type;
bernoulli_distribution(base_type & rng, double p)
: _rng(rng),
_threshold(static_cast<base_result>
(p * (_rng.max() - _rng.min())) + _rng.min())
{
// for p == 0, we can only set _threshold = 0, which is not enough
assert(p > 0);
}
// compiler-generated copy ctor is fine
// assignment is disallowed because there is a reference member
result_type operator()() { return _rng() <= _threshold; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const bernoulli_distribution& x,
const bernoulli_distribution& y)
{ return x._threshold == y._threshold && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const bernoulli_distribution& rhs) const
{ return _threshold == rhs._threshold && _rng == rhs._rng; }
#endif
private:
typedef typename base_type::result_type base_result;
base_type & _rng;
const base_result _threshold;
};
} // namespace boost
#endif // BOOST_RANDOM_BERNOULLI_DISTRIBUTION_HPP

76
include/boost/random/cauchy_distribution.hpp
View File

@@ -1,76 +0,0 @@
/* boost random/cauchy_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP
#define BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP
#include <cmath>
#include <boost/random/uniform_01.hpp>
namespace boost {
#if defined(__GNUC__) && (__GNUC__ < 3)
// Special gcc workaround: gcc 2.95.x ignores using-declarations
// in template classes (confirmed by gcc author Martin v. Loewis)
using std::tan;
#endif
// Cauchy distribution: p(x) = sigma/(pi*(sigma**2 + (x-median)**2))
template<class UniformRandomNumberGenerator, class RealType = double>
class cauchy_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
cauchy_distribution(base_type & rng, result_type median = 0,
result_type sigma = 1)
: _rng(rng), _median(median), _sigma(sigma) { }
// compiler-generated copy constructor is fine
// uniform_01 cannot be assigned, neither can this class
result_type operator()()
{
const double pi = 3.14159265358979323846;
#ifndef BOOST_NO_STDC_NAMESPACE
using std::tan;
#endif
return _median + _sigma * tan(pi*(_rng()-0.5));
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const cauchy_distribution& x,
const cauchy_distribution& y)
{
return x._median == y._median && x._sigma == y._sigma && x._rng == y._rng;
}
#else
// Use a member function
bool operator==(const cauchy_distribution& rhs) const
{
return _median == rhs._median && _sigma == rhs._sigma && _rng == rhs._rng;
}
#endif
private:
uniform_01<base_type, result_type> _rng;
result_type _median, _sigma;
};
} // namespace boost
#endif // BOOST_RANDOM_CAUCHY_DISTRIBUTION_HPP

224
include/boost/random/detail/const_mod.hpp
View File

@@ -1,224 +0,0 @@
/* boost random/detail/const_mod.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_CONST_MOD_HPP
#define BOOST_RANDOM_CONST_MOD_HPP
#include <cassert>
#include <boost/static_assert.hpp>
#include <boost/cstdint.hpp>
#include <boost/integer_traits.hpp>
namespace boost {
namespace random {
/*
* Some random number generators require modular arithmetic. Put
* everything we need here.
* IntType must be an integral type.
*/
namespace detail {
template<bool is_signed>
struct do_add
{ };
template<>
struct do_add<true>
{
template<class IntType>
static IntType add(IntType m, IntType x, IntType c)
{
x += (c-m);
if(x < 0)
x += m;
return x;
}
};
template<>
struct do_add<false>
{
template<class IntType>
static IntType add(IntType, IntType, IntType)
{
// difficult
assert(!"const_mod::add with c too large");
return 0;
}
};
} // namespace detail
template<class IntType, IntType m>
class const_mod
{
public:
static IntType add(IntType x, IntType c)
{
if(c == 0)
return x;
else if(c <= traits::const_max - m) // i.e. m+c < max
return add_small(x, c);
else
return detail::do_add<traits::is_signed>::add(m, x, c);
}
static IntType mult(IntType a, IntType x)
{
if(a == 1)
return x;
else if(m <= traits::const_max/a) // i.e. a*m <= max
return mult_small(a, x);
else if(traits::is_signed && (m%a < m/a))
return mult_schrage(a, x);
else {
// difficult
assert(!"const_mod::mult with a too large");
return 0;
}
}
static IntType mult_add(IntType a, IntType x, IntType c)
{
if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max
return (a*x+c) % m;
else
return add(mult(a, x), c);
}
static IntType invert(IntType x)
{ return x == 0 ? 0 : invert_euclidian(x); }
private:
typedef integer_traits<IntType> traits;
const_mod(); // don't instantiate
static IntType add_small(IntType x, IntType c)
{
x += c;
if(x >= m)
x -= m;
return x;
}
static IntType mult_small(IntType a, IntType x)
{
return a*x % m;
}
static IntType mult_schrage(IntType a, IntType value)
{
const IntType q = m / a;
const IntType r = m % a;
assert(r < q); // check that overflow cannot happen
value = a*(value%q) - r*(value/q);
while(value <= 0)
value += m;
return value;
}
// invert c in the finite field (mod m) (m must be prime)
static IntType invert_euclidian(IntType c)
{
// we are interested in the gcd factor for c, because this is our inverse
BOOST_STATIC_ASSERT(m > 0);
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
#endif
assert(c > 0);
IntType l1 = 0;
IntType l2 = 1;
IntType n = c;
IntType p = m;
for(;;) {
IntType q = p / n;
l1 -= q * l2; // this requires a signed IntType!
p -= q * n;
if(p == 0)
return (l2 < 1 ? l2 + m : l2);
IntType q2 = n / p;
l2 -= q2 * l1;
n -= q2 * p;
if(n == 0)
return (l1 < 1 ? l1 + m : l1);
}
}
};
// The modulus is exactly the word size: rely on machine overflow handling.
// Due to a GCC bug, we cannot partially specialize in the presence of
// template value parameters.
template<>
class const_mod<unsigned int, 0>
{
typedef unsigned int IntType;
public:
static IntType add(IntType x, IntType c) { return x+c; }
static IntType mult(IntType a, IntType x) { return a*x; }
static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
// m is not prime, thus invert is not useful
private: // don't instantiate
const_mod();
};
template<>
class const_mod<unsigned long, 0>
{
typedef unsigned long IntType;
public:
static IntType add(IntType x, IntType c) { return x+c; }
static IntType mult(IntType a, IntType x) { return a*x; }
static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
// m is not prime, thus invert is not useful
private: // don't instantiate
const_mod();
};
// the modulus is some power of 2: rely partly on machine overflow handling
// we only specialize for rand48 at the moment
#ifndef BOOST_NO_INT64_T
template<>
class const_mod<uint64_t, uint64_t(1) << 48>
{
typedef uint64_t IntType;
public:
static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
static IntType mult(IntType a, IntType x) { return mod(a*x); }
static IntType mult_add(IntType a, IntType x, IntType c)
{ return mod(a*x+c); }
static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
// m is not prime, thus invert is not useful
private: // don't instantiate
const_mod();
};
#endif /* !BOOST_NO_INT64_T */
} // namespace random
} // namespace boost
#endif // BOOST_RANDOM_CONST_MOD_HPP

34
include/boost/random/detail/iterator_mixin.hpp
View File

@@ -1,34 +0,0 @@
#ifndef BOOST_ITERATOR_MIXIN_HPP
#define BOOST_ITERATOR_MIXIN_HPP
#include <boost/operators.hpp>
namespace boost {
// must be in boost namespace, otherwise the inline friend trick fails
template<class Generator, class ResultType>
class generator_iterator_mixin_adapter
: incrementable<Generator>, equality_comparable<Generator>
{
public:
typedef std::input_iterator_tag iterator_category;
typedef ResultType value_type;
typedef std::ptrdiff_t difference_type;
typedef const value_type * pointer;
typedef const value_type & reference;
Generator& operator++() { v = cast()(); return cast(); }
const value_type& operator*() const { return v; }
protected:
// instantiate from derived classes only
generator_iterator_mixin_adapter() { }
void iterator_init() { operator++(); }
private:
Generator & cast() { return static_cast<Generator&>(*this); }
value_type v;
};
} // namespace boost
#endif // BOOST_ITERATOR_MIXIN_HPP

79
include/boost/random/detail/signed_unsigned_compare.hpp
View File

@@ -1,79 +0,0 @@
#ifndef BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_COMPARE
#define BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_COMPARE
#include <boost/limits.hpp>
namespace boost {
namespace random {
/*
* Correctly compare two numbers whose types possibly differ in signedness.
* See boost::numeric_cast<> for the general idea.
* Most "if" statements involve only compile-time constants, so the
* optimizing compiler can do its job easily.
*
* With most compilers, the straightforward implementation produces a
* bunch of (legitimate) warnings. Some template magic helps, though.
*/
namespace detail {
template<bool signed1, bool signed2>
struct do_compare
{ };
template<>
struct do_compare<false, false>
{
// cast to the larger type is automatic with built-in types
template<class T1, class T2>
static bool equal(T1 x, T2 y) { return x == y; }
template<class T1, class T2>
static bool lessthan(T1 x, T2 y) { return x < y; }
};
template<>
struct do_compare<true, true> : do_compare<false, false>
{ };
template<>
struct do_compare<true, false>
{
template<class T1, class T2>
static bool equal(T1 x, T2 y) { return x >= 0 && static_cast<T2>(x) == y; }
template<class T1, class T2>
static bool lessthan(T1 x, T2 y) { return x < 0 || static_cast<T2>(x) < y; }
};
template<>
struct do_compare<false, true>
{
template<class T1, class T2>
static bool equal(T1 x, T2 y) { return y >= 0 && x == static_cast<T1>(y); }
template<class T1, class T2>
static bool lessthan(T1 x, T2 y) { return y >= 0 && x < static_cast<T1>(y); }
};
} // namespace detail
template<class T1, class T2>
int equal_signed_unsigned(T1 x, T2 y)
{
typedef std::numeric_limits<T1> x_traits;
typedef std::numeric_limits<T2> y_traits;
return detail::do_compare<x_traits::is_signed, y_traits::is_signed>::equal(x, y);
}
template<class T1, class T2>
int lessthan_signed_unsigned(T1 x, T2 y)
{
typedef std::numeric_limits<T1> x_traits;
typedef std::numeric_limits<T2> y_traits;
return detail::do_compare<x_traits::is_signed, y_traits::is_signed>::lessthan(x, y);
}
} // namespace random
} // namespace boost
#endif // BOOST_RANDOM_DETAIL_SIGNED_UNSIGNED_COMPARE

66
include/boost/random/exponential_distribution.hpp
View File

@@ -1,66 +0,0 @@
/* boost random/exponential_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP
#include <cmath>
#include <cassert>
#include <boost/random/uniform_01.hpp>
namespace boost {
// exponential distribution: p(x) = lambda * exp(-lambda * x)
template<class UniformRandomNumberGenerator, class RealType = double>
class exponential_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
exponential_distribution(base_type& rng, result_type lambda)
: _rng(rng), _lambda(lambda) { assert(lambda > 0); }
// compiler-generated copy ctor is fine
// uniform_01 cannot be assigned, neither can this class
result_type operator()()
{
#ifndef BOOST_NO_STDC_NAMESPACE
using std::log;
#endif
return -1.0 / _lambda * log(1-_rng());
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const exponential_distribution& x,
const exponential_distribution& y)
{ return x._lambda == y._lambda && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const exponential_distribution& rhs) const
{ return _lambda == rhs._lambda && _rng == rhs._rng; }
#endif
private:
uniform_01<base_type, RealType> _rng;
const result_type _lambda;
};
} // namespace boost
#endif // BOOST_RANDOM_EXPONENTIAL_DISTRIBUTION_HPP

81
include/boost/random/geometric_distribution.hpp
View File

@@ -1,81 +0,0 @@
/* boost random/geometric_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP
#define BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP
#include <cmath> // std::log
#include <cassert>
#include <boost/random/uniform_01.hpp>
namespace boost {
#if defined(__GNUC__) && (__GNUC__ < 3)
// Special gcc workaround: gcc 2.95.x ignores using-declarations
// in template classes (confirmed by gcc author Martin v. Loewis)
using std::log;
#endif
// geometric distribution: p(i) = (1-p) * pow(p, i-1) (integer)
template<class UniformRandomNumberGenerator, class IntType = int>
class geometric_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
geometric_distribution(base_type & rng, double p)
: _rng(rng)
{
assert(0.0 < p && p < 1.0);
#ifndef BOOST_NO_STDC_NAMESPACE
using std::log;
#endif
_log_p = log(p);
}
// compiler-generated copy ctor is fine
// uniform_01 cannot be assigned, neither can this class
result_type operator()()
{
#ifndef BOOST_NO_STDC_NAMESPACE
using std::log;
#endif
return IntType (log(1-_rng()) / _log_p) + 1;
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const geometric_distribution& x,
const geometric_distribution& y)
{ return x._log_p == y._log_p && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const geometric_distribution& rhs) const
{ return _log_p == rhs._log_p && _rng == rhs._rng; }
#endif
private:
uniform_01<base_type> _rng;
typename uniform_01<base_type>::result_type _log_p;
};
} // namespace boost
#endif // BOOST_RANDOM_GEOMETRIC_DISTRIBUTION_HPP

108
include/boost/random/inversive_congruential.hpp
View File

@@ -1,108 +0,0 @@
/* boost random/inversive_congruential.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP
#define BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP
#include <iostream>
#include <cassert>
#include <boost/static_assert.hpp>
#include <boost/random/detail/const_mod.hpp>
namespace boost {
namespace random {
// Eichenauer and Lehn 1986
template<class IntType, IntType a, IntType b, IntType p, IntType val>
class inversive_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = (b == 0 ? 1 : 0);
static const result_type max_value = p-1;
#else
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
#endif
BOOST_STATIC_CONSTANT(result_type, multiplier = a);
BOOST_STATIC_CONSTANT(result_type, increment = b);
BOOST_STATIC_CONSTANT(result_type, modulus = p);
result_type min() const { return b == 0 ? 1 : 0; }
result_type max() const { return p-1; }
explicit inversive_congruential(IntType y0 = 1) : value(y0)
{
BOOST_STATIC_ASSERT(b >= 0);
BOOST_STATIC_ASSERT(p > 1);
BOOST_STATIC_ASSERT(a >= 1);
if(b == 0)
assert(y0 > 0);
}
void seed(IntType y0) { value = y0; if(b == 0) assert(y0 > 0); }
IntType operator()()
{
typedef const_mod<IntType, p> do_mod;
value = do_mod::mult_add(a, do_mod::invert(value), b);
return value;
}
bool validation(result_type x) const { return val == x; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, inversive_congruential x)
{ os << x.value; return os; }
friend std::istream& operator>>(std::istream& is, inversive_congruential& x)
{ is >> x.value; return is; }
friend bool operator==(inversive_congruential x, inversive_congruential y)
{ return x.value == y.value; }
#else
// Use a member function; Streamable concept not supported.
bool operator==(inversive_congruential rhs) const
{ return value == rhs.value; }
#endif
private:
IntType value;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const bool inversive_congruential<IntType, a, b, p, val>::has_fixed_range;
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::min_value;
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::max_value;
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::multiplier;
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::increment;
template<class IntType, IntType a, IntType b, IntType p, IntType val>
const typename inversive_congruential<IntType, a, b, p, val>::result_type inversive_congruential<IntType, a, b, p, val>::modulus;
#endif
} // namespace random
typedef random::inversive_congruential<int32_t, 9102, 2147483647-36884165,
2147483647, 0> hellekalek1995;
} // namespace boost
#endif // BOOST_RANDOM_INVERSIVE_CONGRUENTIAL_HPP

216
include/boost/random/lagged_fibonacci.hpp
View File

@@ -1,216 +0,0 @@
/* boost random/lagged_fibonacci.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
#include <iostream>
#include <algorithm> // std::max
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/cstdint.hpp>
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_01.hpp>
namespace boost {
namespace random {
// lagged Fibonacci generator for the range [0..1)
// contributed by Matthias Troyer
// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
template<class T, unsigned int p, unsigned int q>
struct fibonacci_validation
{
BOOST_STATIC_CONSTANT(bool, is_specialized = false);
static T value() { return 0; }
static T tolerance() { return 0; }
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class T, unsigned int p, unsigned int q>
const bool fibonacci_validation<T, p, q>::is_specialized;
#endif
#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
template<> \
struct fibonacci_validation<T, P, Q> \
{ \
BOOST_STATIC_CONSTANT(bool, is_specialized = true); \
static T value() { return V; } \
static T tolerance() \
{ return std::max(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
};
// (The extra static_cast<T> in the std::max call above is actually
// unnecessary except for HP aCC 1.30, which claims that
// numeric_limits<double>::epsilon() doesn't actually return a double.)
BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
#undef BOOST_RANDOM_FIBONACCI_VAL
template<class FloatType, unsigned int p, unsigned int q>
class lagged_fibonacci
{
public:
typedef FloatType result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
result_type min() const { return 0.0; }
result_type max() const { return 1.0; }
lagged_fibonacci() { seed(); }
explicit lagged_fibonacci(uint32_t value) { seed(value); }
template<class Generator>
explicit lagged_fibonacci(Generator & gen) { seed(gen); }
// compiler-generated copy ctor and assignment operator are fine
void seed(uint32_t value = 331u)
{
minstd_rand0 intgen(value);
seed(intgen);
}
// For GCC, moving this function out-of-line prevents inlining, which may
// reduce overall object code size. However, MSVC does not grok
// out-of-line template member functions.
template<class Generator>
void seed(Generator & gen)
{
uniform_01<Generator, FloatType> gen01(gen);
// I could have used std::generate_n, but it takes "gen" by value
for(unsigned int j = 0; j < long_lag; ++j)
x[j] = gen01();
i = long_lag;
}
result_type operator()()
{
if(i >= long_lag)
fill();
return x[i++];
}
bool validation(result_type x) const
{
result_type v = fibonacci_validation<result_type, p, q>::value();
result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
// std::abs is a source of trouble: sometimes, it's not overloaded
// for double, plus the usual namespace std noncompliance -> avoid it
// using std::abs;
// return abs(x - v) < 5 * epsilon
return x > v - epsilon && x < v + epsilon;
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, const lagged_fibonacci& f)
{
os << f.i << " ";
std::streamsize prec =
os.precision(std::numeric_limits<FloatType>::digits10);
for(unsigned int i = 0; i < long_lag; ++i)
os << f.x[i] << " ";
os.precision(prec);
return os;
}
friend std::istream& operator>>(std::istream& is, lagged_fibonacci& f)
{
is >> f.i >> std::ws;
for(unsigned int i = 0; i < long_lag; ++i)
is >> f.x[i] >> std::ws;
return is;
}
friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
{ return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const lagged_fibonacci& rhs) const
{ return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
#endif
private:
void fill();
unsigned int i;
FloatType x[long_lag];
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class FloatType, unsigned int p, unsigned int q>
const bool lagged_fibonacci<FloatType, p, q>::has_fixed_range;
template<class FloatType, unsigned int p, unsigned int q>
const unsigned int lagged_fibonacci<FloatType, p, q>::long_lag;
template<class FloatType, unsigned int p, unsigned int q>
const unsigned int lagged_fibonacci<FloatType, p, q>::short_lag;
#endif
template<class FloatType, unsigned int p, unsigned int q>
void lagged_fibonacci<FloatType, p, q>::fill()
{
// two loops to avoid costly modulo operations
{ // extra scope for MSVC brokenness w.r.t. for scope
for(unsigned int j = 0; j < short_lag; ++j) {
FloatType t = x[j] + x[j+(long_lag-short_lag)];
if(t >= 1.0)
t -= 1.0;
x[j] = t;
}
}
for(unsigned int j = short_lag; j < long_lag; ++j) {
FloatType t = x[j] + x[j-short_lag];
if(t >= 1.0)
t -= 1.0;
x[j] = t;
}
i = 0;
}
} // namespace random
typedef random::lagged_fibonacci<double, 607, 273> lagged_fibonacci607;
typedef random::lagged_fibonacci<double, 1279, 418> lagged_fibonacci1279;
typedef random::lagged_fibonacci<double, 2281, 1252> lagged_fibonacci2281;
typedef random::lagged_fibonacci<double, 3217, 576> lagged_fibonacci3217;
typedef random::lagged_fibonacci<double, 4423, 2098> lagged_fibonacci4423;
typedef random::lagged_fibonacci<double, 9689, 5502> lagged_fibonacci9689;
typedef random::lagged_fibonacci<double, 19937, 9842> lagged_fibonacci19937;
typedef random::lagged_fibonacci<double, 23209, 13470> lagged_fibonacci23209;
typedef random::lagged_fibonacci<double, 44497, 21034> lagged_fibonacci44497;
// It is possible to partially specialize uniform_01<> on lagged_fibonacci<>
// to help the compiler generate efficient code. For GCC, this seems useless,
// because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now.
} // namespace boost
#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP

153
include/boost/random/linear_congruential.hpp
View File

@@ -1,153 +0,0 @@
/* boost random/linear_congruential.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
#define BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP
#include <iostream>
#include <cassert>
#include <boost/config.hpp>
#include <boost/random/detail/const_mod.hpp>
namespace boost {
namespace random {
// compile-time configurable linear congruential generator
template<class IntType, IntType a, IntType c, IntType m, IntType val>
class linear_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = ( c == 0 ? 1 : 0 );
static const result_type max_value = m-1;
#else
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
#endif
BOOST_STATIC_CONSTANT(IntType, multiplier = a);
BOOST_STATIC_CONSTANT(IntType, increment = c);
BOOST_STATIC_CONSTANT(IntType, modulus = m);
result_type min() const { return c == 0 ? 1 : 0; }
result_type max() const { return m-1; }
explicit linear_congruential(IntType x0 = 1)
: _x(x0)
{
assert(c || x0); /* if c == 0 and x(0) == 0 then x(n) = 0 for all n */
// overflow check
// disabled because it gives spurious "divide by zero" gcc warnings
// assert(m == 0 || (a*(m-1)+c) % m == (c < a ? c-a+m : c-a));
}
// compiler-generated copy constructor and assignment operator are fine
void seed(IntType x0) { assert(c || x0); _x = x0; }
IntType operator()()
{
_x = const_mod<IntType, m>::mult_add(a, _x, c);
return _x;
}
bool validation(IntType x) const { return val == x; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os,
const linear_congruential& lcg)
{ os << lcg._x; return os; }
friend std::istream& operator>>(std::istream& is, linear_congruential& lcg)
{ is >> lcg._x; return is; }
friend bool operator==(const linear_congruential& x,
const linear_congruential& y)
{ return x._x == y._x; }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const linear_congruential& rhs) const
{ return _x == rhs._x; }
#endif
private:
IntType _x;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class IntType, IntType a, IntType c, IntType m, IntType val>
const bool linear_congruential<IntType, a, c, m, val>::has_fixed_range;
template<class IntType, IntType a, IntType c, IntType m, IntType val>
const typename linear_congruential<IntType, a, c, m, val>::result_type linear_congruential<IntType, a, c, m, val>::min_value;
template<class IntType, IntType a, IntType c, IntType m, IntType val>
const typename linear_congruential<IntType, a, c, m, val>::result_type linear_congruential<IntType, a, c, m, val>::max_value;
#endif
} // namespace random
// validation values from the publications
typedef random::linear_congruential<int32_t, 16807, 0, 2147483647,
1043618065> minstd_rand0;
typedef random::linear_congruential<int32_t, 48271, 0, 2147483647,
399268537> minstd_rand;
#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
// emulate the lrand48() C library function; requires support for uint64_t
class rand48
{
public:
typedef int32_t result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const int32_t min_value = 0;
static const int32_t max_value = integer_traits<int32_t>::const_max;
#else
enum { has_fixed_range = false };
#endif
int32_t min() const { return 0; }
int32_t max() const { return std::numeric_limits<int32_t>::max(); }
explicit rand48(int32_t x0 = 1) : lcf(cnv(x0)) { }
explicit rand48(uint64_t x0) : lcf(x0) { }
// compiler-generated copy ctor and assignment operator are fine
void seed(int32_t x0) { lcf.seed(cnv(x0)); }
void seed(uint64_t x0) { lcf.seed(x0); }
int32_t operator()() { return lcf() >> 17; }
// by experiment from lrand48()
bool validation(int32_t x) const { return x == 1993516219; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, const rand48& r)
{ os << r.lcf; return os; }
friend std::istream& operator>>(std::istream& is, rand48& r)
{ is >> r.lcf; return is; }
friend bool operator==(const rand48& x, const rand48& y)
{ return x.lcf == y.lcf; }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const rand48& rhs) const
{ return lcf == rhs.lcf; }
#endif
private:
random::linear_congruential<uint64_t,
uint64_t(0xDEECE66DUL) | (uint64_t(0x5) << 32), // xxxxULL is not portable
0xB, uint64_t(1)<<48, /* unknown */ 0> lcf;
static uint64_t cnv(int32_t x)
{ return (static_cast<uint64_t>(x) << 16) | 0x330e; }
};
#endif /* !BOOST_NO_INT64_T && !BOOST_NO_INTEGRAL_INT64_T */
} // namespace boost
#endif // BOOST_RANDOM_LINEAR_CONGRUENTIAL_HPP

83
include/boost/random/lognormal_distribution.hpp
View File

@@ -1,83 +0,0 @@
/* boost random/lognormal_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP
#include <cmath> // std::exp, std::sqrt
#include <cassert>
#include <boost/random/normal_distribution.hpp>
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std {
using ::log;
using ::sqrt;
}
#endif
namespace boost {
#if defined(__GNUC__) && (__GNUC__ < 3)
// Special gcc workaround: gcc 2.95.x ignores using-declarations
// in template classes (confirmed by gcc author Martin v. Loewis)
using std::sqrt;
using std::exp;
#endif
template<class UniformRandomNumberGenerator, class RealType = double>
class lognormal_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
lognormal_distribution(base_type & rng, result_type mean,
result_type sigma)
: _rng(rng, std::log(mean*mean/std::sqrt(sigma*sigma + mean*mean)),
std::sqrt(std::log(sigma*sigma/mean/mean+1)))
{
assert(mean > 0);
}
// compiler-generated copy constructor is fine
// normal_distribution cannot be assigned, neither can this class
result_type operator()()
{
#ifndef BOOST_NO_STDC_NAMESPACE
// allow for Koenig lookup
using std::exp;
#endif
return exp(_rng());
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const lognormal_distribution& x,
const lognormal_distribution& y)
{ return x._rng == y._rng; }
#else
// Use a member function
bool operator==(const lognormal_distribution& rhs) const
{ return _rng == rhs._rng; }
#endif
private:
normal_distribution<base_type, result_type> _rng;
};
} // namespace boost
#endif // BOOST_RANDOM_LOGNORMAL_DISTRIBUTION_HPP

244
include/boost/random/mersenne_twister.hpp
View File

@@ -1,244 +0,0 @@
/* boost random/mersenne_twister.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP
#define BOOST_RANDOM_MERSENNE_TWISTER_HPP
#include <iostream>
#include <algorithm> // std::copy
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/integer_traits.hpp>
#include <boost/cstdint.hpp>
#include <boost/random/linear_congruential.hpp>
namespace boost {
namespace random {
// http://www.math.keio.ac.jp/matumoto/emt.html
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
class mersenne_twister
{
public:
typedef DataType result_type;
BOOST_STATIC_CONSTANT(int, state_size = n);
BOOST_STATIC_CONSTANT(int, shift_size = m);
BOOST_STATIC_CONSTANT(int, mask_bits = r);
BOOST_STATIC_CONSTANT(DataType, parameter_a = a);
BOOST_STATIC_CONSTANT(int, output_u = u);
BOOST_STATIC_CONSTANT(int, output_s = s);
BOOST_STATIC_CONSTANT(DataType, output_b = b);
BOOST_STATIC_CONSTANT(int, output_t = t);
BOOST_STATIC_CONSTANT(DataType, output_c = c);
BOOST_STATIC_CONSTANT(int, output_l = l);
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = integer_traits<result_type>::const_min;
static const result_type max_value = integer_traits<result_type>::const_max;
#else
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
#endif
result_type min() const { return std::numeric_limits<result_type>::min(); }
result_type max() const { return std::numeric_limits<result_type>::max(); }
mersenne_twister() { seed(); }
#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
// Work around overload resolution problem (Gennadiy E. Rozental)
explicit mersenne_twister(const DataType& value)
#else
explicit mersenne_twister(DataType value)
#endif
{ seed(value); }
template<class Generator>
explicit mersenne_twister(Generator & gen) { seed(gen); }
// compiler-generated copy ctor and assignment operator are fine
void seed() { seed(DataType(4357)); }
#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
// Work around overload resolution problem (Gennadiy E. Rozental)
void seed(const DataType& value)
#else
void seed(DataType value)
#endif
{
random::linear_congruential<uint32_t, 69069, 0, 0, /* unknown */ 0>
gen(value);
seed(gen);
}
// For GCC, moving this function out-of-line prevents inlining, which may
// reduce overall object code size. However, MSVC does not grok
// out-of-line definitions of member function templates.
template<class Generator>
void seed(Generator & gen)
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_signed);
#endif
// I could have used std::generate_n, but it takes "gen" by value
for(int j = 0; j < n; j++)
x[j] = gen();
i = n;
}
result_type operator()();
bool validation(result_type v) const { return val == v; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, const mersenne_twister& mt)
{
os << mt.i << " ";
std::copy(mt.x, mt.x+n, std::ostream_iterator<data_type>(os, " "));
return os;
}
friend std::istream& operator>>(std::istream& is, mersenne_twister& mt)
{
is >> mt.i >> std::ws;
for(int i = 0; i < mt.state_size; ++i)
is >> mt.x[i] >> std::ws;
return is;
}
friend bool operator==(const mersenne_twister& x, const mersenne_twister& y)
{ return x.i == y.i && std::equal(x.x, x.x+n, y.x); }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const mersenne_twister& rhs) const
{ return i == rhs.i && std::equal(x, x+n, rhs.x); }
#endif
private:
typedef DataType data_type;
void twist();
int i;
data_type x[n];
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const bool mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const typename mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::result_type mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::min_value;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const typename mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::result_type mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::max_value;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::state_size;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::shift_size;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::mask_bits;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const DataType mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::parameter_a;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_u;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_s;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const DataType mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_b;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_t;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const DataType mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_c;
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
const int mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::output_l;
#endif
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
void mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::twist()
{
const data_type upper_mask = (~0u) << r;
const data_type lower_mask = ~upper_mask;
/*
for(int j = 0; j < n; j++) {
// Step 2
data_type y = (x[j] & upper_mask) | (x[(j+1)%n] & lower_mask);
// Step 3
x[j] = x[(j+m)%n] ^ (y >> 1) ^ (y&1 ? a : 0);
}
*/
// split loop to avoid costly modulo operations
{ // extra scope for MSVC brokenness w.r.t. for scope
for(int j = 0; j < n-m; j++) {
data_type y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j+m] ^ (y >> 1) ^ (y&1 ? a : 0);
}
}
for(int j = n-m; j < n-1; j++) {
data_type y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j+m-n] ^ (y >> 1) ^ (y&1 ? a : 0);
}
// last iteration
data_type y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
x[n-1] = x[n-1+m-n] ^ (y >> 1) ^ (y&1 ? a : 0);
i = 0;
}
template<class DataType, int n, int m, int r, DataType a, int u,
int s, DataType b, int t, DataType c, int l, DataType val>
inline typename mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::result_type
mersenne_twister<DataType,n,m,r,a,u,s,b,t,c,l,val>::operator()()
{
if(i >= n)
twist();
// Step 4
data_type z = x[i];
++i;
z ^= (z >> u);
z ^= ((z << s) & b);
z ^= ((z << t) & c);
z ^= (z >> l);
return z;
}
} // namespace random
typedef random::mersenne_twister<uint32_t,351,175,19,0xccab8ee7,11,
7,0x31b6ab00,15,0xffe50000,17, /* unknown */ 0> mt11213b;
// validation by experiment from mt19937.c
typedef random::mersenne_twister<uint32_t,624,397,31,0x9908b0df,11,
7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
} // namespace boost
#endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP

96
include/boost/random/normal_distribution.hpp
View File

@@ -1,96 +0,0 @@
/* boost random/normal_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
#define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
#include <cmath>
#include <cassert>
#include <boost/random/uniform_01.hpp>
namespace boost {
// deterministic polar method, uses trigonometric functions
template<class UniformRandomNumberGenerator, class RealType = double>
class normal_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
explicit normal_distribution(base_type & rng, const result_type& mean = 0,
const result_type& sigma = 1)
: _rng(rng), _mean(mean), _sigma(sigma), _valid(false)
{
assert(sigma >= 0);
}
// compiler-generated copy constructor is NOT fine, need to purge cache
normal_distribution(const normal_distribution& other)
: _rng(other._rng), _mean(other._mean), _sigma(other._sigma), _valid(false)
{
}
// uniform_01 cannot be assigned, neither can this class
result_type operator()()
{
#ifndef BOOST_NO_STDC_NAMESPACE
// allow for Koenig lookup
using std::sqrt; using std::log; using std::sin; using std::cos;
#endif
if(!_valid) {
_r1 = _rng();
_r2 = _rng();
_cached_rho = sqrt(-2 * log(1.0-_r2));
_valid = true;
} else {
_valid = false;
}
// Can we have a boost::mathconst please?
const double pi = 3.14159265358979323846;
return _cached_rho * (_valid ? cos(2*pi*_r1) : sin(2*pi*_r1)) * _sigma +
_mean;
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const normal_distribution& x,
const normal_distribution& y)
{
return x._mean == y._mean && x._sigma == y._sigma &&
x._valid == y._valid && x._rng == y._rng;
}
#else
// Use a member function
bool operator==(const normal_distribution& rhs) const
{
return _mean == rhs._mean && _sigma == rhs._sigma &&
_valid == rhs._valid && _rng == rhs._rng;
}
#endif
private:
uniform_01<base_type, RealType> _rng;
const result_type _mean, _sigma;
result_type _r1, _r2, _cached_rho;
bool _valid;
};
} // namespace boost
#endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP

62
include/boost/random/random_number_generator.hpp
View File

@@ -1,62 +0,0 @@
/* boost random/random_number_generator.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP
#define BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/detail/iterator_mixin.hpp>
#include <boost/random/detail/signed_unsigned_compare.hpp>
namespace boost {
// a model for RandomNumberGenerator std:25.2.11 [lib.alg.random.shuffle]
template<class UniformRandomNumberGenerator, class IntType = long>
class random_number_generator
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType argument_type;
typedef IntType result_type;
random_number_generator(base_type & rng) : _rng(rng) {
// MSVC requires the typedef workaround
typedef typename base_type::result_type base_result_type;
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(std::numeric_limits<base_result_type>::is_integer);
BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
#endif
}
// compiler-generated copy ctor is fine
// assignment is disallowed because there is a reference member
result_type operator()(argument_type n) {
return uniform_int<base_type>(_rng, 0, n-1)();
}
private:
base_type & _rng;
};
} // namespace boost
#endif // BOOST_RANDOM_RANDOM_NUMBER_GENERATOR_HPP

158
include/boost/random/shuffle_output.hpp
View File

@@ -1,158 +0,0 @@
/* boost random/shuffle_output.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_SHUFFLE_OUTPUT_HPP
#define BOOST_RANDOM_SHUFFLE_OUTPUT_HPP
#include <iostream>
#include <algorithm> // std::copy
#include <cassert>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/cstdint.hpp>
#include <boost/random/linear_congruential.hpp>
namespace boost {
namespace random {
// Carter Bays and S.D. Durham 1979
template<class UniformRandomNumberGenerator, int k,
class IntType = typename UniformRandomNumberGenerator::result_type,
#ifndef BOOST_NO_DEPENDENT_TYPES_IN_TEMPLATE_VALUE_PARAMETERS
typename UniformRandomNumberGenerator::result_type
#else
uint32_t
#endif
val = 0>
class shuffle_output
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef typename base_type::result_type result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
BOOST_STATIC_CONSTANT(int, buffer_size = k);
shuffle_output() : _rng() { init(); }
#if defined(BOOST_MSVC) && _MSC_VER <= 1200
// MSVC does not implicitly generate the copy constructor here
shuffle_output(const shuffle_output & x)
: _rng(x._rng), y(x.y) { std::copy(x.v, x.v+k, v); }
#endif
template<class T>
explicit shuffle_output(T seed) : _rng(seed) { init(); }
explicit shuffle_output(const base_type & rng) : _rng(rng) { init(); }
template<class T>
void seed(T s) { _rng.seed(s); init(); }
result_type operator()() {
// calculating the range every time may seem wasteful. However, this
// makes the information locally available for the optimizer.
result_type range = max()-min()+1;
int j = k*IntType(y-min())/range;
// assert(0 <= j && j < k);
y = v[j];
v[j] = _rng();
return y;
}
result_type min() const { return _rng.min(); }
result_type max() const { return _rng.max(); }
bool validation(result_type x) const { return val == x; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend std::ostream& operator<<(std::ostream& os, const shuffle_output& s)
{
os << s._rng << " " << s.y << " ";
std::copy(s.v, s.v+k, std::ostream_iterator<result_type>(os, " "));
return os;
}
friend std::istream& operator>>(std::istream& is, shuffle_output& s)
{
is >> s._rng >> std::ws >> s.y >> std::ws;
for(int i = 0; i < s.buffer_size; ++i)
is >> s.v[i] >> std::ws;
return is;
}
friend bool operator==(const shuffle_output& x, const shuffle_output& y)
{ return x._rng == y._rng && x.y == y.y && std::equal(x.v, x.v+k, y.v); }
#else
// Use a member function; Streamable concept not supported.
bool operator==(const shuffle_output& rhs) const
{ return _rng == rhs._rng && y == rhs.y && std::equal(v, v+k, rhs.v); }
#endif
private:
void init()
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
#endif
result_type range = max()-min();
assert(range > 0); // otherwise there would be little choice
if(IntType(k * range) < IntType(range)) // not a sufficient condition
// likely overflow with bucket number computation
assert(!"overflow will occur");
// we cannot use std::generate, because it uses pass-by-value for _rng
for(result_type * p = v; p != v+k; ++p)
*p = _rng();
y = _rng();
}
base_type _rng;
result_type v[k];
result_type y;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, int k,
class IntType,
#ifndef BOOST_NO_DEPENDENT_TYPES_IN_TEMPLATE_VALUE_PARAMETERS
typename UniformRandomNumberGenerator::result_type
#else
uint32_t
#endif
val>
const bool shuffle_output<UniformRandomNumberGenerator, k, IntType, val>::has_fixed_range;
template<class UniformRandomNumberGenerator, int k,
class IntType,
#ifndef BOOST_NO_DEPENDENT_TYPES_IN_TEMPLATE_VALUE_PARAMETERS
typename UniformRandomNumberGenerator::result_type
#else
uint32_t
#endif
val>
const int shuffle_output<UniformRandomNumberGenerator, k, IntType, val>::buffer_size;
#endif
} // namespace random
// validation by experiment from Harry Erwin's generator.h (private e-mail)
typedef random::shuffle_output<
random::linear_congruential<uint32_t, 1366, 150889, 714025, 0>,
97, uint32_t, 139726> kreutzer1986;
} // namespace boost
#endif // BOOST_RANDOM_SHUFFLE_OUTPUT_HPP

80
include/boost/random/triangle_distribution.hpp
View File

@@ -1,80 +0,0 @@
/* boost random/triangle_distribution.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP
#define BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP
#include <cmath>
#include <cassert>
#include <boost/random/uniform_01.hpp>
namespace boost {
// triangle distribution, with a smallest, b most probable, and c largest
// value.
template<class UniformRandomNumberGenerator, class RealType = double>
class triangle_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
triangle_distribution(base_type & rng, result_type a, result_type b,
result_type c)
: _rng(rng), _a(a), _b(b), _c(c),
d1(_b-_a), d2(_c-_a), d3(_c-_b), q1(d1/d2), p1(d1*d2)
{
#ifndef BOOST_NO_STDC_NAMESPACE
using std::sqrt;
#endif
d3 = sqrt(d3);
p1 = sqrt(p1);
assert(_a <= _b && _b <= _c);
}
// compiler-generated copy ctor is fine
// uniform_01 cannot be assigned, neither can this class
result_type operator()()
{
#ifndef BOOST_NO_STDC_NAMESPACE
using std::sqrt;
#endif
result_type u = _rng();
if( u <= q1 )
return _a + p1*sqrt(u);
else
return _c - d3*sqrt(d2*u-d1);
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const triangle_distribution& x,
const triangle_distribution& y)
{ return x._a == y._a && x._b == y._b && x._c == y._c && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const triangle_distribution& rhs) const
{ return _a == rhs._a && _b == rhs._b && _c == rhs._c && _rng == rhs._rng; }
#endif
private:
uniform_01<base_type, result_type> _rng;
result_type _a, _b, _c;
result_type d1, d2, d3, q1, p1;
};
} // namespace boost
#endif // BOOST_RANDOM_TRIANGLE_DISTRIBUTION_HPP

78
include/boost/random/uniform_01.hpp
View File

@@ -1,78 +0,0 @@
/* boost random/uniform_01.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_UNIFORM_01_HPP
#define BOOST_RANDOM_UNIFORM_01_HPP
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
namespace boost {
// Because it is so commonly used: uniform distribution on the real [0..1)
// range. This allows for specializations to avoid a costly FP division
template<class UniformRandomNumberGenerator, class RealType = double>
class uniform_01
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
explicit uniform_01(base_type & rng) : _rng(rng) {
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
#endif
}
// compiler-generated copy ctor is fine
// assignment is disallowed because there is a reference member
result_type operator()() {
return static_cast<result_type>(_rng() - _rng.min()) /
(static_cast<result_type>(_rng.max()-_rng.min()) +
(std::numeric_limits<base_result>::is_integer ? 1.0 : 0.0));
}
result_type min() const { return 0.0; }
result_type max() const { return 1.0; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const uniform_01& x, const uniform_01& y)
{ return x._rng == y._rng; }
#else
// Use a member function
bool operator==(const uniform_01& rhs) const
{ return _rng == rhs._rng; }
#endif
private:
typedef typename base_type::result_type base_result;
base_type & _rng;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, class RealType>
const bool uniform_01<UniformRandomNumberGenerator, RealType>::has_fixed_range;
#endif
} // namespace boost
#endif // BOOST_RANDOM_UNIFORM_01_HPP

139
include/boost/random/uniform_int.hpp
View File

@@ -1,139 +0,0 @@
/* boost random/uniform_int.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-04-08 added min<max assertion (N. Becker)
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_UNIFORM_INT_HPP
#define BOOST_RANDOM_UNIFORM_INT_HPP
#include <cassert>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/uniform_smallint.hpp>
#include <boost/random/detail/signed_unsigned_compare.hpp>
namespace boost {
// uniform integer distribution on [min, max]
template<class UniformRandomNumberGenerator, class IntType = int>
class uniform_int
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
uniform_int(base_type & rng, IntType min, IntType max)
: _rng(rng), _min(min), _max(max), _range(_max - _min),
_bmin(_rng.min()), _brange(_rng.max() - _bmin)
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer);
#endif
assert(min < max);
if(random::equal_signed_unsigned(_brange, _range))
_range_comparison = 0;
else if(random::lessthan_signed_unsigned(_brange, _range))
_range_comparison = -1;
else
_range_comparison = 1;
}
result_type operator()();
result_type min() const { return _min; }
result_type max() const { return _max; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const uniform_int& x, const uniform_int& y)
{ return x._min == y._min && x._max == y._max && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const uniform_int& rhs) const
{ return _min == rhs._min && _max == rhs._max && _rng == rhs._rng; }
#endif
private:
typedef typename base_type::result_type base_result;
base_type & _rng;
result_type _min, _max, _range;
base_result _bmin, _brange;
int _range_comparison;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, class IntType>
const bool uniform_int<UniformRandomNumberGenerator, IntType>::has_fixed_range;
#endif
template<class UniformRandomNumberGenerator, class IntType>
inline IntType uniform_int<UniformRandomNumberGenerator, IntType>::operator()()
{
if(_range_comparison == 0) {
// this will probably never happen in real life
// basically nothing to do; just take care we don't overflow / underflow
return static_cast<result_type>(_rng() - _bmin) + _min;
} else if(_range_comparison < 0) {
// use rejection method to handle things like 0..3 --> 0..4
for(;;) {
// concatenate several invocations of the base RNG
// take extra care to avoid overflows
result_type limit;
if(_range == std::numeric_limits<result_type>::max()) {
limit = _range/(static_cast<result_type>(_brange)+1);
if(_range % static_cast<result_type>(_brange)+1 == static_cast<result_type>(_brange))
++limit;
} else {
limit = (_range+1)/(static_cast<result_type>(_brange)+1);
}
// we consider "result" as expressed to base (_brange+1)
// for every power of (_brange+1), we determine a random factor
result_type result = 0;
result_type mult = 1;
while(mult <= limit) {
result += (_rng() - _bmin) * mult;
mult *= static_cast<result_type>(_brange)+1;
}
if(mult == limit)
// _range+1 is an integer power of _brange+1: no rejections required
return result;
// _range/mult < _brange+1 -> no endless loop
result += uniform_int<base_type,result_type>(_rng, 0, _range/mult)() * mult;
if(result <= _range)
return result + _min;
}
} else { // brange > range
if(_brange / _range > 4 /* quantization_cutoff */ ) {
// the new range is vastly smaller than the source range,
// so quantization effects are not relevant
return uniform_smallint<base_type,result_type>(_rng, _min, _max)();
} else {
// use rejection method to handle things like 0..5 -> 0..4
for(;;) {
base_result result = _rng() - _bmin;
// result and range are non-negative, and result is possibly larger
// than range, so the cast is safe
if(result <= static_cast<base_result>(_range))
return result + _min;
}
}
}
}
} // namespace boost
#endif // BOOST_RANDOM_UNIFORM_INT_HPP

79
include/boost/random/uniform_on_sphere.hpp
View File

@@ -1,79 +0,0 @@
/* boost random/uniform_on_sphere.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
#define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
#include <vector>
#include <algorithm> // std::transform
#include <functional> // std::bind2nd, std::divides
#include <boost/random/normal_distribution.hpp>
namespace boost {
template<class UniformRandomNumberGenerator, class RealType = double,
class Cont = std::vector<RealType> >
class uniform_on_sphere
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef Cont result_type;
explicit uniform_on_sphere(base_type & rng, int dim = 2)
: _rng(rng), _container(dim), _dim(dim) { }
// compiler-generated copy ctor is fine
// normal_distribution cannot be assigned, neither can this class
const result_type & operator()()
{
RealType sqsum = 0;
for(typename Cont::iterator it = _container.begin();
it != _container.end();
++it) {
RealType val = _rng();
*it = val;
sqsum += val * val;
}
#ifndef BOOST_NO_STDC_NAMESPACE
using std::sqrt;
#endif
// for all i: result[i] /= sqrt(sqsum)
std::transform(_container.begin(), _container.end(), _container.begin(),
std::bind2nd(std::divides<RealType>(), sqrt(sqsum)));
return _container;
}
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const uniform_on_sphere& x,
const uniform_on_sphere& y)
{ return x._dim == y._dim && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const uniform_on_sphere& rhs) const
{ return _dim == rhs._dim && _rng == rhs._rng; }
#endif
private:
normal_distribution<base_type, RealType> _rng;
result_type _container;
const int _dim;
};
} // namespace boost
#endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP

78
include/boost/random/uniform_real.hpp
View File

@@ -1,78 +0,0 @@
/* boost random/uniform_real.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-04-08 added min<max assertion (N. Becker)
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_UNIFORM_REAL_HPP
#define BOOST_RANDOM_UNIFORM_REAL_HPP
#include <cassert>
#include <boost/config.hpp>
#include <boost/limits.hpp>
#include <boost/static_assert.hpp>
#include <boost/random/uniform_01.hpp>
namespace boost {
// uniform distribution on a real range
template<class UniformRandomNumberGenerator, class RealType = double>
class uniform_real
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
uniform_real(base_type & rng, RealType min, RealType max)
: _rng(rng), _min(min), _max(max)
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
#endif
assert(min < max);
}
// compiler-generated copy ctor is fine
// uniform_01 cannot be assigned, neither can this class
result_type operator()() { return _rng() * (_max - _min) + _min; }
result_type min() const { return _min; }
result_type max() const { return _max; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const uniform_real& x, const uniform_real& y)
{ return x._min == y._min && x._max == y._max && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const uniform_real& rhs) const
{ return _min == rhs._min && _max == rhs._max && _rng == rhs._rng; }
#endif
private:
uniform_01<base_type, result_type> _rng;
RealType _min, _max;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, class RealType>
const bool uniform_real<UniformRandomNumberGenerator, RealType>::has_fixed_range;
#endif
} // namespace boost
#endif // BOOST_RANDOM_UNIFORM_REAL_HPP

103
include/boost/random/uniform_smallint.hpp
View File

@@ -1,103 +0,0 @@
/* boost random/uniform_smallint.hpp header file
*
* Copyright Jens Maurer 2000-2001
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* See http://www.boost.org for most recent version including documentation.
*
* $Id$
*
* Revision history
* 2001-04-08 added min<max assertion (N. Becker)
* 2001-02-18 moved to individual header files
*/
#ifndef BOOST_RANDOM_UNIFORM_SMALLINT_HPP
#define BOOST_RANDOM_UNIFORM_SMALLINT_HPP
#include <cassert>
#include <boost/config.hpp>
#include <boost/limits.hpp>
namespace boost {
// uniform integer distribution on a small range [min, max]
template<class UniformRandomNumberGenerator, class IntType = int>
class uniform_smallint
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
uniform_smallint(base_type & rng, IntType min, IntType max);
result_type operator()()
{
// we must not use the low bits here, because LCGs get very bad then
return ((_rng() - _rng.min()) / _factor) % _range + _min;
}
result_type min() const { return _min; }
result_type max() const { return _max; }
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
friend bool operator==(const uniform_smallint& x, const uniform_smallint& y)
{ return x._min == y._min && x._max == y._max && x._rng == y._rng; }
#else
// Use a member function
bool operator==(const uniform_smallint& rhs) const
{ return _min == rhs._min && _max == rhs._max && _rng == rhs._rng; }
#endif
private:
typedef typename base_type::result_type base_result;
base_type & _rng;
IntType _min, _max;
base_result _range;
int _factor;
};
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
template<class UniformRandomNumberGenerator, class IntType>
const bool uniform_smallint<UniformRandomNumberGenerator, IntType>::has_fixed_range;
#endif
template<class UniformRandomNumberGenerator, class IntType>
uniform_smallint<UniformRandomNumberGenerator, IntType>::
uniform_smallint(base_type & rng, IntType min, IntType max)
: _rng(rng), _min(min), _max(max),
_range(static_cast<base_result>(_max-_min)+1), _factor(1)
{
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
BOOST_STATIC_ASSERT(std::numeric_limits<IntType>::is_integer);
#endif
assert(min < max);
// LCGs get bad when only taking the low bits.
// (probably put this logic into a partial template specialization)
// Check how many low bits we can ignore before we get too much
// quantization error.
base_result r_base = _rng.max() - _rng.min();
if(r_base == std::numeric_limits<base_result>::max()) {
_factor = 2;
r_base /= 2;
}
r_base += 1;
if(r_base % _range == 0) {
// No quantization effects, good
_factor = r_base / _range;
} else {
const base_result r = 32*_range*_range;
for(; r_base >= r; _factor *= 2)
r_base /= 2;
}
}
} // namespace boost
#endif // BOOST_RANDOM_UNIFORM_SMALLINT_HPP

143
index.html
View File

@@ -1,143 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost Random Number Library</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<table border="1" bgcolor="#007F7F" cellpadding="2">
<tr>
<td bgcolor="#FFFFFF"><img src="../../c++boost.gif" alt="c++boost.gif (8819 bytes)" width="277" height="86"></td>
<td><a href="../../index.htm"><font face="Arial,Helvetica" color="#FFFFFF"><big>Home</big></font></a></td>
<td><a href="../libraries.htm"><font face="Arial,Helvetica" color="#FFFFFF"><big>Libraries</big></font></a></td>
<td><a href="../../people/people.htm"><font face="Arial,Helvetica" color="#FFFFFF"><big>People</big></font></a></td>
<td><a href="../../more/faq.htm"><font face="Arial,Helvetica" color="#FFFFFF"><big>FAQ</big></font></a></td>
<td><a href="../../more/index.htm"><font face="Arial,Helvetica" color="#FFFFFF"><big>More</big></font></a></td>
</tr>
</table>
<h1>Boost Random Number Library</h1>
Random numbers are useful in a variety of applications. The Boost
Random Number Library (Boost.Random for short) provides a vast variety
of generators and distributions to produce random numbers having
useful properties, such as uniform distribution.
<p>
You should read the
<a href="random-concepts.html">concepts documentation</a>
for an introduction and the definition of the basic concepts. For a
quick start, it may be sufficient to have a look at <a
href="random_demo.cpp">random_demo.cpp</a>.
<h2>Library Organization</h2>
The library is separated into several header files, all within the
<code>boost/random/</code> directory. Additionally, a convenience
header file which includes all other headers in
<code>boost/random/</code> is available as
<code><a href="../../boost/random.hpp">boost/random.hpp</a></code>.
<p>
Several random number generators are available in the following
header files; please read the
<a href="random-generators.html">documentation</a> about these.
<ul>
<li><code><a href="../../boost/random/linear_congruential.hpp">boost/random/linear_congruential.hpp</a></code>
<li><code><a href="../../boost/random/additive_combine.hpp">boost/random/additive_combine.hpp</a></code>
<li><code><a href="../../boost/random/inversive_congruential.hpp">boost/random/inversive_congruential.hpp</a></code>
<li><code><a href="../../boost/random/shuffle_output.hpp">boost/random/shuffle_output.hpp</a></code>
<li><code><a href="../../boost/random/mersenne_twister.hpp">boost/random/mersenne_twister.hpp</a></code>
<li><code><a href="../../boost/random/lagged_fibonacci.hpp">boost/random/lagged_fibonacci.hpp</a></code>
</ul>
Similarly, several random number distributions are available in the
following header files; please read the
<a href="random-distributions.html">documentation</a> about these.
<ul>
<li><code><a href="../../boost/random/uniform_smallint.hpp">boost/random/uniform_smallint.hpp</a></code>
<li><code><a href="../../boost/random/uniform_int.hpp">boost/random/uniform_int.hpp</a></code>
<li><code><a href="../../boost/random/uniform_01.hpp">boost/random/uniform_01.hpp</a></code>
<li><code><a href="../../boost/random/uniform_real.hpp">boost/random/uniform_real.hpp</a></code>
<li><code><a href="../../boost/random/triangle_distribution.hpp">boost/random/triangle_distribution.hpp</a></code>
<li><code><a href="../../boost/random/bernoulli_distribution.hpp">boost/random/bernoulli_distribution.hpp</a></code>
<li><code><a href="../../boost/random/cauchy_distribution.hpp">boost/random/cauchy_distribution.hpp</a></code>
<li><code><a href="../../boost/random/exponential_distribution.hpp">boost/random/exponential_distribution.hpp</a></code>
<li><code><a href="../../boost/random/geometric_distribution.hpp">boost/random/geometric_distribution.hpp</a></code>
<li><code><a href="../../boost/random/normal_distribution.hpp">boost/random/normal_distribution.hpp</a></code>
<li><code><a href="../../boost/random/lognormal_distribution.hpp">boost/random/lognormal_distribution.hpp</a></code>
<li><code><a href="../../boost/random/uniform_on_sphere.hpp">boost/random/uniform_on_sphere.hpp</a></code>
</ul>
Additionally, non-deterministic random number generators are available
in the header
<code><a href="../../boost/nondet_random.hpp">&lt;boost/nondet_random.hpp&gt;</a></code>.
<a href="nondet_random.html">Documentation</a> is also available.
<p>
In order to map the interface of the generators and distribution functions
to other concepts, some <a href="random-misc.html">decorators</a> are available.
<h2>Tests</h2>
An extensive test suite for the pseudo-random number generators and
distributions is available as
<a href="random_test.cpp">random_test.cpp</a>.
<p>
Some <a href="random-performance.html">performance results</a> obtained
using <a href="random_speed.cpp">random_speed.cpp</a> are also available.
<h2>Rationale</h2>
The methods for generating and evaluating deterministic and
non-deterministic random numbers differ radically. Furthermore, due
to the inherent deterministic design of present-day computers, it is
often difficult to implement non-deterministic random number
generation facilities. Thus, the random number library is split into
separate header files, mirroring the two different application
domains.
<h2>History and Acknowledgements</h2>
In November 1999, Jeet Sukumaran proposed a framework based on virtual
functions, and later sketched a template-based approach. Ed Brey
pointed out that Microsoft Visual C++ does not support in-class member
initializations and suggested the <code>enum</code> workaround. Dave
Abrahams highlighted quantization issues.
<p>
The first public release of this random number library materialized in
March 2000 after extensive discussions on the boost mailing list.
Many thanks to Beman Dawes for his original <code>min_rand</code>
class, portability fixes, documentation suggestions, and general
guidance. Harry Erwin sent a header file which provided additional
insight into the requirements. Ed Brey and Beman Dawes wanted an
iterator-like interface.
<p>
Beman Dawes managed the formal review, during which Matthias Troyer,
Csaba Szepesvari, and Thomas Holenstein gave detailed comments. The
reviewed version became an official part of boost on 17 June 2000.
<p>
Gary Powell contributed suggestions for code cleanliness. Dave
Abrahams and Howard Hinnant suggested to move the basic generator
templates from namespace <code>boost::detail</code> to
<code>boost::random</code>.
<p>
Ed Brey asked to remove superfluous warnings and helped with
<code>uint64_t</code> handling. Andreas Scherer tested with MSVC.
Matthias Troyer contributed a lagged Fibonacci generator. Michael
Stevens found a bug in the copy semantics of normal_distribution and
suggested documentation improvements.
<p>
<hr>
<a href="../../people/jens_maurer.htm">Jens Maurer</a>,
2001-08-31
</body>
</html>

84
integrate.hpp
View File

@@ -1,84 +0,0 @@
/* integrate.hpp header file
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* Revision history
* 01 April 2001: Modified to use new <boost/limits.hpp> header. (JMaddock)
*/
#ifndef INTEGRATE_HPP
#define INTEGRATE_HPP
#include <boost/limits.hpp>
template<class UnaryFunction>
inline typename UnaryFunction::result_type
trapezoid(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::argument_type b, int n)
{
typename UnaryFunction::result_type tmp = 0;
for(int i = 1; i <= n-1; ++i)
tmp += f(a+(b-a)/n*i);
return (b-a)/2/n * (f(a) + f(b) + 2*tmp);
}
template<class UnaryFunction>
inline typename UnaryFunction::result_type
simpson(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::argument_type b, int n)
{
typename UnaryFunction::result_type tmp1 = 0;
for(int i = 1; i <= n-1; ++i)
tmp1 += f(a+(b-a)/n*i);
typename UnaryFunction::result_type tmp2 = 0;
for(int i = 1; i <= n ; ++i)
tmp2 += f(a+(b-a)/2/n*(2*i-1));
return (b-a)/6/n * (f(a) + f(b) + 2*tmp1 + 4*tmp2);
}
// compute b so that f(b) = y; assume f is monotone increasing
template<class UnaryFunction>
inline typename UnaryFunction::argument_type
invert_monotone_inc(UnaryFunction f, typename UnaryFunction::result_type y,
typename UnaryFunction::argument_type lower = -1,
typename UnaryFunction::argument_type upper = 1)
{
while(upper-lower > 1e-6) {
double middle = (upper+lower)/2;
if(f(middle) > y)
upper = middle;
else
lower = middle;
}
return (upper+lower)/2;
}
// compute b so that I(f(x), a, b) == y
template<class UnaryFunction>
inline typename UnaryFunction::argument_type
quantil(UnaryFunction f, typename UnaryFunction::argument_type a,
typename UnaryFunction::result_type y,
typename UnaryFunction::argument_type step)
{
typedef typename UnaryFunction::result_type result_type;
if(y >= 1.0)
return std::numeric_limits<result_type>::infinity();
typename UnaryFunction::argument_type b = a;
for(result_type result = 0; result < y; b += step)
result += step*f(b);
return b;
}
#endif /* INTEGRATE_HPP */

121
nondet_random.html
View File

@@ -1,121 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost RNG Library - Non-Deterministic Random Number Generators</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<h1><img src="../../c++boost.gif" alt="c++boost.gif (8819 bytes)"
align="center" width="277" height="86">Header
<a href="../../boost/nondet_random.hpp">&lt;boost/nondet_random.hpp&gt;</a></h1>
<ul>
<li><a href="#synopsis">Synopsis</a>
<li><a href="#random_device">Class <code>random_device</code></a>
<li><a href="#performance">Performance</a>
</ul>
<h2><a name="synopsis">Header</a><code>&lt;boost/nondet_random.hpp&gt;</code>
Synopsis</h2>
<pre>
namespace boost {
class random_device;
} // namespace boost
</pre>
<h2><a name="random_device">Class <code>random_device</code></a></h2>
<h3>Synopsis</h3>
<pre>
class random_device : noncopyable
{
public:
typedef unsigned int result_type;
static const bool has_fixed_range = true;
static const result_type min_value = /* implementation defined */;
static const result_type max_value = /* implementation defined */;
result_type min() const;
result_type max() const;
explicit random_device(const std::string& token = default_token);
~random_device();
unsigned int operator()();
};
</pre>
<h3>Description</h3>
Class <code>random_device</code> models a
<a href="random-concepts.html#nondet-rng">non-deterministic random number
generator</a>.
It uses one or more implementation-defined stochastic processes to
generate a sequence of uniformly distributed non-deterministic random
numbers. For those environments where a non-deterministic random
number generator is not available, class <code>random_device</code>
must not be implemented. See
<blockquote>
"Randomness Recommendations for Security", D. Eastlake, S.
Crocker, J. Schiller, Network Working Group, RFC 1750, December 1994
</blockquote>
for further discussions.
<p>
<em>Note:</em> Some operating systems abstract the computer hardware
enough to make it difficult to non-intrusively monitor stochastic
processes. However, several do provide a special device for exactly
this purpose. It seems to be impossible to emulate the functionality
using Standard C++ only, so users should be aware that this class may
not be available on all platforms.
<h3>Members</h3>
<pre>explicit random_device(const std::string& token = default_token)</pre>
<strong>Effects:</strong> Constructs a <code>random_device</code>,
optionally using the given <code>token</code> as an access
specification (for example, a URL) to some implementation-defined
service for monitoring a stochastic process.
<h3>Implementation Note for Linux</h3>
On the Linux operating system, <code>token</code> is interpreted as a
filesystem path. It is assumed that this path denotes an operating
system pseudo-device which generates a stream of non-deterministic
random numbers. The pseudo-device should never signal an error or
end-of-file. Otherwise, <code>std::ios_base::failure</code> is
thrown. By default, <code>random_device</code> uses the
<code>/dev/urandom</code> pseudo-device to retrieve the random
numbers. Another option would be to specify the
<code>/dev/random</code> pseudo-device, which blocks on reads if the
entropy pool has no more random bits available.
<h2><a name="performance">Performance</a></h2>
The test program <a href="nondet_random_speed.cpp">nondet_random_speed.cpp</a>
measures the execution times of the
<a href="../../boost/nondet_random.hpp">nondet_random.hpp</a> implementation of the above
algorithms in a tight loop. The performance has been evaluated on a
Pentium Pro 200 MHz with gcc 2.95.2, Linux 2.2.13, glibc 2.1.2.
<p>
<table border=1>
<tr><th>class</th><th>time per invocation [usec]</th></tr>
<tr><td>random_device</td><td>92.0</td></tr>
</table>
<p>
The measurement error is estimated at +/- 1 usec.
<p>
<hr>
Jens Maurer, 2000-06-19
</body>
</html>

69
nondet_random_speed.cpp
View File

@@ -1,69 +0,0 @@
/* boost nondet_random_speed.cpp performance test
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#include <iostream>
#include <string>
#include <boost/timer.hpp>
#include <boost/nondet_random.hpp>
// set to your CPU frequency in MHz
static const double cpu_frequency = 200 * 1e6;
static void show_elapsed(double end, int iter, const std::string & name)
{
double usec = end/iter*1e6;
double cycles = usec * cpu_frequency/1e6;
std::cout << name << ": "
<< usec*1e3 << " nsec/loop = "
<< cycles << " CPU cycles"
<< std::endl;
}
template<class Result, class RNG>
static void timing(RNG & rng, int iter, const std::string& name)
{
volatile Result tmp; // make sure we're not optimizing too much
boost::timer t;
for(int i = 0; i < iter; i++)
tmp = rng();
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
void run(int iter, const std::string & name)
{
RNG rng;
timing<long>(rng, iter, name);
}
int main(int argc, char*argv[])
{
if(argc != 2) {
std::cerr << "usage: " << argv[0] << " iterations" << std::endl;
return 1;
}
int iter = std::atoi(argv[1]);
#ifdef __linux__
boost::random_device dev;
timing<unsigned int>(dev, iter, "random_device");
#else
#error The non-deterministic random device is currently available on Linux only.
#endif
return 0;
}

344
random-concepts.html
View File

@@ -1,344 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost Random Number Library Concepts</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<h1>Random Number Generator Library Concepts</h1>
<h2>Introduction</h2>
Random numbers are required in a number of different problem domains,
such as
<ul>
<li>numerics (simulation, Monte-Carlo integration)
<li>games (non-deterministic enemy behavior)
<li>security (key generation)
<li>testing (random coverage in white-box tests)
</ul>
The Boost Random Number Generator Library provides a framework
for random number generators with well-defined properties so that the
generators can be used in the demanding numerics and security domains.
For a general introduction to random numbers in numerics, see
<blockquote>
"Numerical Recipes in C: The art of scientific computing", William
H. Press, Saul A. Teukolsky, William A. Vetterling, Brian P. Flannery,
2nd ed., 1992, pp. 274-328
</blockquote>
Depending on the requirements of the problem domain, different
variations of random number generators are appropriate:
<ul>
<li>non-deterministic random number generator
<li>pseudo-random number generator
<li>quasi-random number generator
</ul>
All variations have some properties in common, these concepts (in the
STL sense) are called NumberGenerator and
UniformRandomNumberGenerator. Each concept will be defined in a
subsequent section.
<p>
The goals for this library are the following:
<ul>
<li>allow easy integration of third-party random-number generators
<li>define a validation interface for the generators
<li>provide easy-to-use front-end classes which model popular
distributions
<li>provide maximum efficiency
<li>allow control on quantization effects in front-end processing
(not yet done)
</ul>
<h2><a name="number_generator">Number Generator</a></h2>
A number generator is a <em>function object</em> (std:20.3
[lib.function.objects]) that takes zero arguments. Each call to
<code>operator()</code> returns a number.
In the following table, <code>X</code> denotes a number generator
class returning objects of type <code>T</code>, and <code>u</code> is
a value of <code>X</code>.
<p>
<table border=1>
<tr><th colspan=3 align=center><code>NumberGenerator</code>
requirements</th></tr>
<tr><td>expression</td><td>return&nbsp;type</td>
<td>pre/post-condition</td></tr>
<tr><td><code>X::result_type</code></td><td>T</td>
<td><code>std::numeric_limits&lt;T&gt;::is_specialized</code> is true,
<code>T</code> is <code>LessThanComparable</code></td></tr>
<tr><td><code>u.operator()()</code></td><td>T</td><td>-</td></tr>
</table>
<p>
<em>Note:</em> The NumberGenerator requirements do not impose any
restrictions on the characteristics of the returned numbers.
<h2><a name="uniform-rng">Uniform Random Number Generator</a></h2>
A uniform random number generator is a NumberGenerator that provides a
sequence of random numbers uniformly distributed on a given range.
The range can be compile-time fixed or available (only) after run-time
construction of the object.
<p>
The <em>tight lower bound</em> of some (finite) set S is the (unique)
member l in S, so that for all v in S, l <= v holds. Likewise, the
<em>tight upper bound</em> of some (finite) set S is the (unique)
member u in S, so that for all v in S, v <= u holds.
<p>
In the following table, <code>X</code> denotes a number generator
class returning objects of type <code>T</code>, and <code>v</code> is
a const value of <code>X</code>.
<p>
<table border=1>
<tr><th colspan=3 align=center><code>UniformRandomNumberGenerator</code>
requirements</th></tr>
<tr><td>expression</td><td>return&nbsp;type</td>
<td>pre/post-condition</td></tr>
<tr><td><code>X::has_fixed_range</code></td><td><code>bool</code></td>
<td>compile-time constant; if <code>true</code>, the range on which
the random numbers are uniformly distributed is known at compile-time
and members <code>min_value</code> and <code>max_value</code>
exist. <em>Note:</em> This flag may also be <code>false</code> due to
compiler limitations.</td></tr>
<tr><td><code>X::min_value</code></td><td><code>T</code></td>
<td>compile-time constant; <code>min_value</code> is equal to
<code>v.min()</code></td></tr>
<tr><td><code>X::max_value</code></td><td><code>T</code></td>
<td>compile-time constant; <code>max_value</code> is equal to
<code>v.max()</code></td></tr>
<tr><td><code>v.min()</code></td><td><code>T</code></td>
<td>tight lower bound on the set of all values returned by
<code>operator()</code>. The return value of this function shall not
change during the lifetime of the object.</td></tr>
<tr><td><code>v.max()</code></td><td><code>T</code></td>
<td>if <code>std::numeric_limits&lt;T&gt;::is_integer</code>, tight
upper bound on the set of all values returned by
<code>operator()</code>, otherwise, the smallest representable number
larger than the tight upper bound on the set of all values returned by
<code>operator()</code>. In any case, the return value of this
function shall not change during the lifetime of the
object.</code></td></tr>
</table>
<p>
The member functions <code>min</code>, <code>max</code>, and
<code>operator()</code> shall have amortized constant time complexity.
<p>
<em>Note:</em> For integer generators (i.e. integer <code>T</code>),
the generated values <code>x</code> fulfill <code>min() <= x <=
max()</code>, for non-integer generators (i.e. non-integer
<code>T</code>), the generated values <code>x</code> fulfill
<code>min() <= x < max()</code>.
<br>
<em>Rationale:</em> The range description with <code>min</code> and
<code>max</code> serves two purposes. First, it allows scaling of the
values to some canonical range, such as [0..1). Second, it describes
the significant bits of the values, which may be relevant for further
processing.
<br>
The range is a closed interval [min,max] for integers, because the
underlying type may not be able to represent the half-open interval
[min,max+1). It is a half-open interval [min, max) for non-integers,
because this is much more practical for borderline cases of continuous
distributions.
<p>
<em>Note:</em> The UniformRandomNumberGenerator concept does not
require <code>operator()(long)</code> and thus it does not fulfill the
RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle])
requirements. Use the
<a href="random-misc.html#random_number_generator"><code>random_number_generator</code></a>
adapter for that.
<br>
<em>Rationale:</em> <code>operator()(long)</code> is not provided,
because mapping the output of some generator with integer range to a
different integer range is not trivial.
<h2><a name="nondet-rng">Non-deterministic Uniform Random Number
Generator</a></h2>
A non-deterministic uniform random number generator is a
UniformRandomNumberGenerator that is based on some stochastic process.
Thus, it provides a sequence of truly-random numbers. Examples for
such processes are nuclear decay, noise of a Zehner diode, tunneling
of quantum particles, rolling a die, drawing from an urn, and tossing
a coin. Depending on the environment, inter-arrival times of network
packets or keyboard events may be close approximations of stochastic
processes.
<p>
The class
<code><a href="nondet_random.html#random_device">random_device</a></code>
is a model for a non-deterministic random number generator.
<p>
<em>Note:</em> This type of random-number generator is useful for
security applications, where it is important to prevent that an
outside attacker guesses the numbers and thus obtains your
encryption or authentication key. Thus, models of this concept should
be cautious not to leak any information, to the extent possible by the
environment. For example, it might be advisable to explicitly clear
any temporary storage as soon as it is no longer needed.
<h2><a name="pseudo-rng">Pseudo-Random Number Generator</a></h2>
A pseudo-random number generator is a UniformRandomNumberGenerator
which provides a deterministic sequence of pseudo-random numbers,
based on some algorithm and internal state. Linear congruential and
inversive congruential generators are examples of such pseudo-random
number generators. Often, these generators are very sensitive to
their parameters. In order to prevent wrong implementations from
being used, an external testsuite should check that the generated
sequence and the validation value provided do indeed match.
<p>
Donald E. Knuth gives an extensive overview on pseudo-random number
generation in his book "The Art of Computer Programming, Vol. 2, 3rd
edition, Addison-Wesley, 1997". The descriptions for the specific
generators contain additional references.
<p>
<em>Note:</em> Because the state of a pseudo-random number generator
is necessarily finite, the sequence of numbers returned by the
generator will loop eventually.
<p>
In addition to the UniformRandomNumberGenerator requirements, a
pseudo-random number generator has some additional requirements. In
the following table, <code>X</code> denotes a pseudo-random number
generator class returning objects of type <code>T</code>,
<code>x</code> is a value of <code>T</code>, <code>u</code> is a value
of <code>X</code>, and <code>v</code> is a <code>const</code> value of
<code>X</code>.
<p>
<table border=1>
<tr><th colspan=3 align=center><code>PseudoRandomNumberGenerator</code>
requirements</th></tr>
<tr><td>expression</td><td>return&nbsp;type</td>
<td>pre/post-condition</td></tr>
<tr><td><code>X()</code></td><td>-</td>
<td>creates a generator in some implementation-defined
state. <em>Note:</em> Several generators thusly created may possibly
produce dependent or identical sequences of random numbers.</td></tr>
<tr><td><code>explicit X(...)</code></td><td>-</td>
<td>creates a generator with user-provided state; the implementation
shall specify the constructor argument(s)</td></tr>
<tr><td><code>u.seed(...)</code></td><td>void</td>
<td>sets the current state according to the argument(s); at least
functions with the same signature as the non-default
constructor(s) shall be provided.
<tr><td><code>v.validation(x)</code></td><td><code>bool</code></td>
<td>compares the pre-computed and hardcoded 10001th element in the
generator's random number sequence with <code>x</code>. The generator
must have been constructed by its default constructor and
<code>seed</code> must not have been called for the validation to
be meaningful.
</table>
<p>
<em>Note:</em> The <code>seed</code> member function is similar to the
<code>assign</code> member function in STL containers. However, the
naming did not seem appropriate.
<p>
Classes which model a pseudo-random number generator shall also model
EqualityComparable, i.e. implement <code>operator==</code>. Two
pseudo-random number generators are defined to be <em>equivalent</em>
if they both return an identical sequence of numbers starting from a
given state.
<p>
Classes which model a pseudo-random number generator should also model
the Streamable concept, i.e. implement <code>operator&lt;&lt;</code>
and <code>operator&gt;&gt;</code>. If so,
<code>operator&lt;&lt;</code> writes all current state of the
pseudo-random number generator to the given <code>ostream</code> so
that <code>operator&gt;&gt;</code> can restore the state at a later
time. The state shall be written in a platform-independent manner,
but it is assumed that the <code>locale</code>s used for writing and
reading be the same.
The pseudo-random number generator with the restored state and the
original at the just-written state shall be equivalent.
<p>
Classes which model a pseudo-random number generator may also model
the CopyConstructible and Assignable concepts. However, note that the
sequences of the original and the copy are strongly correlated (in
fact, they are identical), which may make them unsuitable for some
problem domains. Thus, copying pseudo-random number generators is
discouraged; they should always be passed by (non-<code>const</code>)
reference.
<p>
The classes
<code><a href="random-generators.html#rand48">rand48</a></code>,
<code><a href="random-generators.html#linear_congruential">minstd_rand</a></code>,
and
<code><a href="random-generators.html#mersenne_twister">mt19937</a></code>
are models for a pseudo-random number generator.
<p>
<em>Note:</em> This type of random-number generator is useful for
numerics, games and testing. The non-zero arguments constructor(s)
and the <code>seed()</code> member function(s) allow for a
user-provided state to be installed in the generator. This is useful
for debugging Monte-Carlo algorithms and analyzing particular test
scenarios. The Streamable concept allows to save/restore the state of
the generator, for example to re-run a test suite at a later time.
<h2><a name="quasi-rng">Quasi-Random Number Generators</a></h2>
A quasi-random number generator is a Number Generator which provides a
deterministic sequence of numbers, based on some algorithm and
internal state. The numbers do not have any statistical properties
(such as uniform distribution or independence of successive values).
<p>
<em>Note:</em> Quasi-random number generators are useful for
Monte-Carlo integrations where specially crafted sequences of random
numbers will make the approximation converge faster.
<p>
<em>[Does anyone have a model?]</em>
<p>
<hr>
Jens Maurer, 2000-02-23
</body>
</html>

781
random-distributions.html
View File

@@ -1,781 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost Random Number Library Distributions</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<h1>Random Number Library Distributions</h1>
<ul>
<li><a href="#intro">Introduction</a>
<li><a href="#synopsis">Synopsis</a>
<li><a href="#uniform_smallint">Class template
<code>uniform_smallint</code></a>
<li><a href="#uniform_int">Class template <code>uniform_int</code></a>
<li><a href="#uniform_01">Class template <code>uniform_01</code></a>
<li><a href="#uniform_real">Class template
<code>uniform_real</code></a>
<li><a href="#bernoulli_distribution">Class template
<code>bernoulli_distribution</code></a>
<li><a href="#geometric_distribution">Class template
<code>geometric_distribution</code></a>
<li><a href="#triangle_distribution">Class template
<code>triangle_distribution</code></a>
<li><a href="#exponential_distribution">Class template
<code>exponential_distribution</code></a>
<li><a href="#normal_distribution">Class template
<code>normal_distribution</code></a>
<li><a href="#lognormal_distribution">Class template
<code>lognormal_distribution</code></a>
<li><a href="#uniform_on_sphere">Class template
<code>uniform_on_sphere</code></a>
</ul>
<h2><a name="intro">Introduction</a></h2>
In addition to the <a href="random-generators.html">random number
generators</a>, this library provides distribution functions which map
one distribution (often a uniform distribution provided by some
generator) to another.
<p>
Usually, there are several possible implementations of any given
mapping. Often, there is a choice between using more space, more
invocations of the underlying source of random numbers, or more
time-consuming arithmetic such as trigonometric functions. This
interface description does not mandate any specific implementation.
However, implementations which cannot reach certain values of the
specified distribution or otherwise do not converge statistically to
it are not acceptable.
<p>
<table border="1">
<tr><th>distribution</th><th>explanation</th><th>example</th></tr>
<tr>
<td><code><a href="#uniform_smallint">uniform_smallint</a></code></td>
<td>discrete uniform distribution on a small set of integers (much
smaller than the range of the underlying generator)</td>
<td>drawing from an urn</td>
</tr>
<tr>
<td><code><a href="#uniform_int">uniform_int</a></code></td>
<td>discrete uniform distribution on a set of integers; the underlying
generator may be called several times to gather enough randomness for
the output</td>
<td>drawing from an urn</td>
</tr>
<tr>
<td><code><a href="#uniform_01">uniform_01</a></code></td>
<td>continuous uniform distribution on the range [0,1); important
basis for other distributions</td>
<td>-</td>
</tr>
<tr>
<td><code><a href="#uniform_real">uniform_real</a></code></td>
<td>continuous uniform distribution on some range [min, max) of real
numbers</td>
<td>for the range [0, 2pi): randomly dropping a stick and measuring
its angle in radiants (assuming the angle is uniformly
distributed)</td>
</tr>
<tr>
<td><code><a href="#bernoulli_distribution">bernoulli_distribution</a></code></td>
<td>Bernoulli experiment: discrete boolean valued distribution with
configurable probability</td>
<td>tossing a coin (p=0.5)</td>
</tr>
<tr>
<td><code><a href="#geometric_distribution">geometric_distribution</a></code></td>
<td>measures distance between outcomes of repeated Bernoulli experiments</td>
<td>throwing a die several times and counting the number of tries
until a "6" appears for the first time</td>
</tr>
<tr>
<td><code><a href="#triangle_distribution">triangle_distribution</a></code></td>
<td>?</td>
<td>?</td>
</tr>
<tr>
<td><code><a href="#exponential_distribution">exponential_distribution</a></code></td>
<td>exponential distribution</td>
<td>measuring the inter-arrival time of alpha particles emitted by
radioactive matter</td>
</tr>
<tr>
<td><code><a href="#normal_distribution">normal_distribution</a></code></td>
<td>counts outcomes of (infinitely) repeated Bernoulli experiments</td>
<td>tossing a coin 10000 times and counting how many front sides are shown</td>
</tr>
<tr>
<td><code><a href="#lognormal_distribution">lognormal_distribution</a></code></td>
<td>lognormal distribution (sometimes used in simulations)</td>
<td>measuring the job completion time of an assembly line worker</td>
</tr>
<tr>
<td><code><a href="#uniform_on_sphere">uniform_on_sphere</a></code></td>
<td>uniform distribution on a unit sphere of arbitrary dimension</td>
<td>choosing a random point on Earth (assumed to be a sphere) where to
spend the next vacations</td>
</tr>
</table>
<p>
The template parameters of the distribution functions are always in
the order
<ul>
<li>Underlying source of random numbers
<li>If applicable, return type, with a default to a reasonable type.
</ul>
<p>
<em>The distribution functions no longer satisfy the input iterator
requirements (std:24.1.1 [lib.input.iterators]), because this is
redundant given the Generator interface and imposes a run-time
overhead on all users. Moreover, a Generator interface appeals to
random number generation as being more "natural". Use an
<a href="../utility/iterator_adaptors.htm">iterator adaptor</a>
if you need to wrap any of the generators in an input iterator
interface.</em>
<p>
All of the distribution functions described below store a non-const
reference to the underlying source of random numbers. Therefore, the
distribution functions are not Assignable. However, they are
CopyConstructible. Copying a distribution function will copy the
parameter values. Furthermore, both the copy and the original will
refer to the same underlying source of random numbers. Therefore,
both the copy and the original will obtain their underlying random
numbers from a single sequence.
<p>
In this description, I have refrained from documenting those members
in detail which are already defined in the
<a href="random-concepts.html">concept documentation</a>.
<h2><a name="synopsis">Synopsis of the distributions</a> available from header
<code>&lt;boost/random.hpp&gt;</code> </h2>
<pre>
namespace boost {
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class uniform_smallint;
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class uniform_int;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class uniform_01;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class uniform_real;
// discrete distributions
template&lt;class UniformRandomNumberGenerator&gt;
class bernoulli_distribution;
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class geometric_distribution;
// continuous distributions
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class triangle_distribution;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class exponential_distribution;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class normal_distribution;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class lognormal_distribution;
template&lt;class UniformRandomNumberGenerator, class RealType = double,
class Cont = std::vector&lt;RealType&gt; &gt;
class uniform_on_sphere;
}
</pre>
<h2><a name="uniform_smallint">Class template
<code>uniform_smallint</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/uniform_smallint.hpp">boost/random/uniform_smallint.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class uniform_smallint
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
static const bool has_fixed_range = false;
uniform_smallint(base_type & rng, IntType min, IntType max);
result_type operator()();
result_type min() const;
result_type max() const;
};
</pre>
<h3>Description</h3>
The distribution function <code>uniform_smallint</code> models a
<a href="random-concepts.html#uniform-rng">uniform random number
generator</a>. On each invocation, it returns a random integer value
uniformly distributed in the set of integer numbers {min, min+1,
min+2, ..., max}. It assumes that the desired range (max-min+1) is
small compared to the range of the underlying source of random
numbers and thus makes no attempt to limit quantization errors.
<p>
Let r<sub>out</sub>=(max-min+1) the desired range of integer numbers,
and let r<sub>base</sub> be the range of the underlying source of
random numbers. Then, for the uniform distribution, the theoretical
probability for any number i in the range r<sub>out</sub> will be
p<sub>out</sub>(i) = 1/r<sub>out</sub>. Likewise, assume a uniform
distribution on r<sub>base</sub> for the underlying source of random
numbers, i.e. p<sub>base</sub>(i) = 1/r<sub>base</sub>. Let
p<sub>out_s</sub>(i) denote the random distribution generated by
<code>uniform_smallint</code>. Then the sum over all i in
r<sub>out</sub> of (p<sub>out_s</sub>(i)/p<sub>out</sub>(i)
-1)<sup>2</sup> shall not exceed
r<sub>out</sub>/r<sub>base</sub><sup>2</sup> (r<sub>base</sub> mod
r<sub>out</sub>)(r<sub>out</sub> - r<sub>base</sub> mod
r<sub>out</sub>).
<p>
The template parameter <code>UniformRandomNumberGenerator</code> shall
denote a class which models a uniform random number generator.
Additionally, <code>UniformRandomNumberGenerator::result_type</code>
shall denote an integral type. The template parameter
<code>IntType</code> shall denote an integer-like value type.
<p>
<em>Note:</em> The property above is the square sum of the relative
differences in probabilities between the desired uniform distribution
p<sub>out</sub>(i) and the generated distribution
p<sub>out_s</sub>(i). The property can be fulfilled with the
calculation (base_rng mod r<sub>out</sub>), as follows: Let r =
r<sub>base</sub> mod r<sub>out</sub>. The base distribution on
r<sub>base</sub> is folded onto the range r<sub>out</sub>. The
numbers i &lt; r have assigned (r<sub>base</sub> div
r<sub>out</sub>)+1 numbers of the base distribution, the rest has only
(r<sub>base</sub> div r<sub>out</sub>). Therefore,
p<sub>out_s</sub>(i) = ((r<sub>base</sub> div r<sub>out</sub>)+1) /
r<sub>base</sub> for i &lt; r and p<sub>out_s</sub>(i) =
(r<sub>base</sub> div r<sub>out</sub>)/r<sub>base</sub> otherwise.
Substituting this in the above sum formula leads to the desired
result.
<p>
<em>Note:</em> The upper bound for (r<sub>base</sub> mod r<sub>out</sub>)(r<sub>out</sub> - r<sub>base</sub>
mod r<sub>out</sub>) is r<sub>out</sub><sup>2</sup>/4. Regarding the upper bound for the square
sum of the relative quantization error of r<sub>out</sub><sup>3</sup>/(4*r<sub>base</sub><sup>2</sup>), it
seems wise to either choose r<sub>base</sub> so that r<sub>base</sub> &gt; 10*r<sub>out</sub><sup>2</sup> or
ensure that r<sub>base</sub> is divisible by r<sub>out</sub>.
<h3>Members</h3>
<pre>uniform_smallint(base_type & rng, IntType min, IntType max)</pre>
<strong>Effects:</strong> Constructs a <code>uniform_smallint</code>
functor with the uniform random number generator <code>rng</code> as
the underlying source of random numbers. <code>min</code> and
<code>max</code> are the lower and upper bounds of the output range,
respectively.
<h2><a name="uniform_int">Class template <code>uniform_int</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/uniform_int.hpp">boost/random/uniform_int.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class uniform_int
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
static const bool has_fixed_range = false;
uniform_int(base_type & rng, IntType min, IntType max);
IntType operator()();
result_type min() const;
result_type max() const;
};
</pre>
<h3>Description</h3>
The distribution function <code>uniform_int</code> models a
<a href="random-concepts.html#uniform-rng">uniform random number
generator</a>. On each invocation, it returns a random integer
value uniformly distributed in the set of integer numbers
{min, min+1, min+2, ..., max}.
<p>
The template parameter <code>IntType</code> shall denote an
integer-like value type.
<h3>Members</h3>
<pre>uniform_int(base_type & rng, IntType min, IntType max)</pre>
<strong>Effects:</strong> Constructs a <code>uniform_int</code> functor
with the uniform random number generator <code>rng</code> as the
underlying source of random numbers. <code>min</code> and
<code>max</code> are the lower and upper bounds of the output range,
respectively.
<p>
<em>Note:</em> Invocations of <code>operator()</code> may call the
underlying generator several times and concatenate the result to build
the required range. Thus, using this distribution with generators
such as linear congruential ones which tend to produce non-random bits
in some positions is strongly discouraged.
<h2><a name="uniform_01">Class template <code>uniform_01</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/uniform_01.hpp">boost/random/uniform_01.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class uniform_01
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
static const bool has_fixed_range = false;
explicit uniform_01(base_type & rng);
result_type operator()();
result_type min() const;
result_type max() const;
};
</pre>
<h3>Description</h3>
The distribution function <code>uniform_01</code> models a
<a href="random-concepts.html#uniform-rng">uniform random number
generator</a>. On each invocation, it returns a random floating-point
value uniformly distributed in the range [0..1).
The value is computed using
<code>std::numeric_limits&lt;RealType&gt;::digits</code> random binary
digits, i.e. the mantissa of the floating-point value is completely
filled with random bits. [<em>Note:</em> Should this be configurable?]
<p>
The template parameter <code>RealType</code> shall denote a float-like
value type with support for binary operators +, -, and /. It must be
large enough to hold floating-point numbers of value
<code>rng.max()-rng.min()+1</code>.
<p>
<code>base_type::result_type</code> must be a number-like value type,
it must support <code>static_cast&lt;&gt;</code> to
<code>RealType</code> and binary operator -.
<p>
<em>Note:</em> The current implementation is buggy, because it may not
fill all of the mantissa with random bits. I'm unsure how to fill a
(to-be-invented) <code>boost::bigfloat</code> class with random bits
efficiently. It's probably time for a traits class.
<h3>Members</h3>
<pre>explicit uniform_01(base_type & rng)</pre>
<strong>Effects:</strong> Constructs a <code>uniform_01</code> functor
with the given uniform random number generator as the underlying
source of random numbers.
<h2><a name="uniform_real">Class template <code>uniform_real</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/uniform_real.hpp">boost/random/uniform_real.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class uniform_real
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
static const bool has_fixed_range = false;
uniform_real(base_type & rng, RealType min, RealType max);
result_type operator()();
result_type min() const;
result_type max() const;
};
</pre>
<h3>Description</h3>
The distribution function <code>uniform_real</code> models a
<a href="random-concepts.html#uniform-rng">uniform random number
generator</a>. On each invocation, it returns a random floating-point
value uniformly distributed in the range [min..max). The value is
computed using
<code>std::numeric_limits&lt;RealType&gt;::digits</code> random binary
digits, i.e. the mantissa of the floating-point value is completely
filled with random bits.
<p>
<em>Note:</em> The current implementation is buggy, because it may not
fill all of the mantissa with random bits.
<h3>Members</h3>
<pre>explicit uniform_real(base_type & rng, RealType min, RealType max)</pre>
<strong>Effects:</strong> Constructs a <code>uniform_real</code>
functor. <code>rng</code> specifies the uniform random number
generator to be used as the underlying source of random numbers,
<code>min</code> and <code>max</code> are the lower and upper bounds of
the output range.
<h2><a name="bernoulli_distribution">Class template
<code>bernoulli_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/bernoulli_distribution.hpp">boost/random/bernoulli_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator&gt;
class bernoulli_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef bool result_type;
bernoulli_distribution(base_type & rng, double q);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>bernoulli_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into a Bernoulli
one.
<h3>Members</h3>
<pre>bernoulli_distribution(base_type & rng, double q)</pre>
<strong>Effects:</strong> Constructs a
<code>bernoulli_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>q</code> is the parameter for the distribution.
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A random boolean value with p(true) = q and
p(false) = 1-q. For example, with q = 1/2 this can be interpreted as
tossing a coin.
<h2><a name="geometric_distribution">Class template
<code>geometric_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/geometric_distribution.hpp">boost/random/geometric_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class IntType = int&gt;
class geometric_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType result_type;
geometric_distribution(base_type& rng, double q);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>geometric_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into a geometric
one.
<h3>Members</h3>
<pre>geometric_distribution(base_type& rng, const result_type& q)</pre>
<strong>Effects:</strong> Constructs a
<code>geometric_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>q</code> is the parameter for the distribution.
<p>
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A random integer value <em>i</em> >= 1 with
p(i) = (1-q) * q<sup>i-1</sup>. For example, with q = 5/6 this can be
interpreted as the number of times one has to roll a die until a given
number shows up.
<h2><a name="triangle_distribution">Class template
<code>triangle_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/triangle_distribution.hpp">boost/random/triangle_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class triangle_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
triangle_distribution(base_type& rng, result_type a, result_type b, result_type c);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>triangle_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into a triangle
one.
<h3>Members</h3>
<pre>triangle_distribution(base_type& rng, result_type a, result_type b, result_type c)</pre>
<strong>Effects:</strong> Constructs a
<code>triangle_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>a, b, c</code> are the parameters for the distribution.
<p>
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A random floating-point value <code>x</code>
where <code>a <= x <= c</code>; <code>x</code> has a triangle
distribution, where <code>b</code> is the most probable value for
<code>x</code>.
<h2><a name="exponential_distribution">Class template
<code>exponential_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/exponential_distribution.hpp">boost/random/exponential_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class exponential_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
exponential_distribution(base_type& rng, const result_type& lambda);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>exponential_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into an
exponential one.
<h3>Members</h3>
<pre>exponential_distribution(base_type& rng, const result_type& lambda)</pre>
<strong>Effects:</strong> Constructs an
<code>exponential_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>lambda</code> is the parameter for the distribution.
<p>
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A random floating-point value <em>x</em>
&gt; 0 with p(x) = <code>lambda</code> * exp(-<code>lambda</code> *
x).
<h2><a name="normal_distribution">Class template
<code>normal_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/normal_distribution.hpp">boost/random/normal_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class normal_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
explicit normal_distribution(base_type& rng, const result_type& mean = 0,
const result_type& sigma = 1);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>normal_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into a
normal (Gaussian) one.
<h3>Members</h3>
<pre>normal_distribution(base_type& rng, const result_type& mean = 0,
const result_type& sigma = 1)</pre>
<strong>Effects:</strong> Constructs a
<code>normal_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>mean</code> and <code>sigma</code> are the parameters for
the distribution.
<p>
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A random floating-point value <em>x</em>
with p(x) = 1/sqrt(2*pi*sigma) * exp(- (x-mean)<sup>2</sup> /
(2*sigma<sup>2</sup>) ).
<h2><a name="lognormal_distribution">Class template
<code>lognormal_distribution</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/lognormal_distribution.hpp">boost/random/lognormal_distribution.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double&gt;
class lognormal_distribution
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef RealType result_type;
explicit lognormal_distribution(base_type& rng, const result_type& mean,
const result_type& sigma);
result_type operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>lognormal_distribution</code>
model a <a href="random-concepts.html#number_generator">number
generator</a>. It transforms a uniform distribution into a
lognormal one.
<h3>Members</h3>
<pre>lognormal_distribution(base_type& rng, const result_type& mean,
const result_type& sigma)</pre>
<strong>Effects:</strong> Constructs a
<code>lognormal_distribution</code> functor with the uniform random
number generator <code>rng</code> as the underlying source of random
numbers. <code>mean</code> and <code>sigma</code> are the mean and
standard deviation of the lognormal distribution.
<p>
<pre>result_type operator()()</pre>
<!-- <strong>Returns:</strong> A random floating-point value <em>x</em>
with p(x) = 1/(x * sigma * sqrt(2*pi)) * exp(-
(log(x)-mean)<sup>2</sup> / (2*sigma<sup>2</sup>) ) for x > 0.
<p> -->
<strong>Returns:</strong> A random floating-point value <em>x</em>
with p(x) = 1/(x * normal_sigma * sqrt(2*pi)) * exp(
-(log(x)-normal_mean)<sup>2</sup> / (2*normal_sigma<sup>2</sup>) )
for x > 0,
where normal_mean = log(mean<sup>2</sup>/sqrt(sigma<sup>2</sup>
+ mean<sup>2</sup>))
and normal_sigma = sqrt(log(1 + sigma<sup>2</sup>/mean<sup>2</sup>)).
<h2><a name="uniform_on_sphere">Class template
<code>uniform_on_sphere</code></a></h2>
<h3>Synopsis</h3>
<pre>
#include &lt;<a href="../../boost/random/uniform_on_sphere.hpp">boost/random/uniform_on_sphere.hpp</a>&gt;
template&lt;class UniformRandomNumberGenerator, class RealType = double,
class Cont = std::vector&lt;RealType&gt; &gt;
class uniform_on_sphere
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef Cont result_type;
explicit uniform_on_sphere(base_type & rng, int dim = 2);
const result_type & operator()();
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>uniform_on_sphere</code> model a
Generator (std:25.2.6 [lib.alg.generate]). It transforms a uniform
distribution into a uniform distribution on the unit sphere of
arbitrary dimension. The <code>Cont</code> template parameter must be
a STL-like container type with <code>begin</code> and <code>end</code>
operations returning non-const ForwardIterators of type
<code>Cont::iterator</code>.
<h3>Members</h3>
<pre>explicit uniform_on_sphere(base_type & rng, int dim = 2)</pre>
<strong>Effects:</strong> Constructs a <code>uniform_on_sphere</code>
functor with the uniform random number generator <code>rng</code> as
the underlying source of random numbers. <code>dim</code> is the
dimension of the sphere.
<p>
<pre>result_type operator()()</pre>
<strong>Returns:</strong> A position on the unit sphere of
<code>dim</code> dimensions in cartesian coordinates. The positions
are uniformly distributed on the unit sphere.
<p>
<strong>Complexity:</strong> Proportional to the number of dimensions.
<p>
<hr>
Jens Maurer, 2001-04-15
</body>
</html>

1042
random-generators.html
View File

File diff suppressed because it is too large Load Diff

93
random-misc.html
View File

@@ -1,93 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost Random Number Generator Library (Miscellaneous)</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<h1>Random Number Generator Library --- Miscellaneous Decorators</h1>
<ul>
<li><a href="#random_number_generator">Class template
<code>random_number_generator</code></a>
<li><a href="#generator_iterator">Class template
<code>generator_iterator</code></a>
</ul>
<h2>Introduction</h2>
These decorator class templates allow adaptation of the random number
generators and distribution functions to concepts found in the C++
Standard Library, in particular the RandomNumberGenerator and the
InputIterator concepts. The latter adaptation is useful, because the
the basic random number generators do not implement the InputIterator
requirements per se, in contrast to the distribution functions.
<h2><a name="synopsis">Synopsis</a> of miscellaneous decorators in
header <code>&lt;boost/random.hpp&gt;</code></h2>
<pre>
namespace boost {
template&lt;class UniformRandomNumberGenerator, class IntType = long&gt;
class random_number_generator;
template&lt;class Generator&gt;
class generator_iterator;
} // namespace boost
</pre>
<h2><a name="random_number_generator">Class template
<code>random_number_generator</code></a></h2>
<h3>Synopsis</h3>
<pre>
template&lt;class UniformRandomNumberGenerator, class IntType = long&gt;
class random_number_generator
{
public:
typedef UniformRandomNumberGenerator base_type;
typedef IntType argument_type;
typedef IntType result_type;
random_number_generator(base_type & rng);
result_type operator()(argument_type n);
};
</pre>
<h3>Description</h3>
Instantiations of class template <code>random_number_generator</code>
model a RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle]).
On each invocation, it returns a uniformly distributed integer in
the range [0..<code>n</code>).
<p>
The template parameter <code>IntType</code> shall denote some
integer-like value type.
<p>
<em>Note:</em> I consider it unfortunate that the C++ Standard uses
the name RandomNumberGenerator for something rather specific.
<h3>Members</h3>
<pre>random_number_generator(base_type & rng)</pre>
<strong>Effects:</strong> Constructs a
<code>random_number_generator</code> functor with the given uniform
random number generator as the underlying source of random numbers.
<pre>result_type operator()(argument_type n)</pre>
<strong>Returns:</strong> The value of
<code>uniform_int&lt;base_type&gt;(rng, 0, n-1)()</code>.
<p>
<hr>
Jens Maurer, 2001-11-19
</body>
</html>

241
random-performance.html
View File

@@ -1,241 +0,0 @@
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<title>Boost Random Number Library Performance</title>
</head>
<body bgcolor="#FFFFFF" text="#000000">
<h1>Random Number Library Performance</h1>
For some people, performance of random number generation is an
important consideration when choosing a random number generator or a
particular distribution function. This page provides numerous
performance tests with the wide variety of generators and
distributions available in the boost library.
<p>
The performance has been evaluated on a Pentium Pro 200 MHz with gcc
2.95.2, Linux 2.2.13, glibc 2.1.2. The speed is reported in million
random numbers per second (M rn/sec), generated in a tight loop.
<p>
<h2>Basic Generators</h2>
<table border="1">
<tr>
<th>generator</th>
<th>M rn/sec</th>
<th>time per random number [usec]</th>
<th>relative speed compared to fastest [percent]</th>
</tr>
<tr>
<td>rand48</td>
<td>5.38</td>
<td>0.183</td>
<td>61%</td>
</tr>
<tr>
<td>rand48 run-time configurable</td>
<td>1.48</td>
<td>0.677</td>
<td>17%</td>
</tr>
<tr>
<td>lrand48 glibc 2.1.2</td>
<td>1.19</td>
<td>0.843</td>
<td>13%</td>
</tr>
<tr>
<td>minstd_rand</td>
<td>2.39</td>
<td>0.318</td>
<td>35%</td>
</tr>
<tr>
<td>ecuyer1988</td>
<td>1.12</td>
<td>0.892</td>
<td>13%</td>
</tr>
<tr>
<td>kreutzer1986</td>
<td>3.87</td>
<td>0.258</td>
<td>43%</td>
</tr>
<tr>
<td>hellekalek1995 (inversive)</td>
<td>0.20</td>
<td>5.12</td>
<td>2%</td>
</tr>
<tr>
<td>mt11213b</td>
<td>6.07</td>
<td>0.165</td>
<td>68%</td>
</tr>
<tr>
<td>mt19937</td>
<td>6.06</td>
<td>0.165</td>
<td>68%</td>
</tr>
<tr>
<td>mt19937 original</td>
<td>5.33</td>
<td>0.188</td>
<td>60%</td>
</tr>
<tr>
<td>lagged_fibonacci607</td>
<td>8.90</td>
<td>0.112</td>
<td>100%</td>
</tr>
<tr>
<td>lagged_fibonacci4423</td>
<td>8.54</td>
<td>0.117</td>
<td>96%</td>
</tr>
<tr>
<td>lagged_fibonacci19937</td>
<td>7.49</td>
<td>0.133</td>
<td>84%</td>
</tr>
<tr>
<td>lagged_fibonacci23209</td>
<td>6.63</td>
<td>0.151</td>
<td>74%</td>
</tr>
<tr>
<td>lagged_fibonacci44497</td>
<td>4.01</td>
<td>0.250</td>
<td>45%</td>
</tr>
</table>
<p>
Note that the lagged Fibonacci generators produce floating-point
numbers, whereas all others produce integers.
<h2>Distributions</h2>
<table border="1">
<tr>
<th>[M rn/sec]</th>
<th>minstd_rand</th>
<th>kreutzer1986</th>
<th>mt19937</th>
<th>lagged_fibonacci607</th>
</tr>
<tr>
<th>uniform_smallint</th>
<td>1.26</td>
<td>1.55</td>
<td>1.93</td>
<td>-</td>
</tr>
<tr>
<th>uniform_01</th>
<td>1.79</td>
<td>1.88</td>
<td>3.03</td>
<td>7.74</td>
</tr>
<tr>
<th>uniform_real</th>
<td>1.74</td>
<td>1.56</td>
<td>2.34</td>
<td>6.62</td>
</tr>
<tr>
<th>geometric</th>
<td>0.593</td>
<td>0.629</td>
<td>0.753</td>
<td>0.916</td>
</tr>
<tr>
<th>triangle</th>
<td>0.97</td>
<td>1.02</td>
<td>1.35</td>
<td>1.31</td>
</tr>
<tr>
<th>exponential</th>
<td>0.849</td>
<td>0.828</td>
<td>0.887</td>
<td>1.53</td>
</tr>
<tr>
<th>normal (polar method)</th>
<td>0.608</td>
<td>0.626</td>
<td>0.738</td>
<td>0.755</td>
</tr>
<tr>
<th>lognormal</th>
<td>0.417</td>
<td>0.442</td>
<td>0.470</td>
<td>0.481</td>
</tr>
<tr>
<th>uniform_on_sphere</th>
<td>0.154</td>
<td>0.155</td>
<td>0.174</td>
<td>0.218</td>
</tr>
</table>
<p>
Note that the lagged Fibonacci generator is at least 2.5 times faster
than the Mersenne twister when generating uniformly distributed
floating-point numbers. For more sophisticated distributions, the
speed improvement is less. Note however that these distributions have
not been optimized for speed, yet.
<p>
<hr>
Jens Maurer, 2001-04-15
</body>
</html>

126
random_demo.cpp
View File

@@ -1,126 +0,0 @@
/* boost random_demo.cpp profane demo
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* A short demo program how to use the random number library.
*/
#include <iostream>
#include <fstream>
#include <ctime> // std::time
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_smallint.hpp>
#include <boost/random/uniform_01.hpp>
// Sun CC doesn't handle boost::iterator_adaptor yet
#if !defined(__SUNPRO_CC) || (__SUNPRO_CC > 0x530)
#include <boost/generator_iterator.hpp>
#endif
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std {
using ::time;
}
#endif
// This is a typedef for a random number generator.
// Try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand
typedef boost::minstd_rand base_generator_type;
// This is a reproducible simulation experiment. See main().
void experiment(base_generator_type & generator)
{
// Define a uniform random number distribution of integer values between
// 1 and 6 inclusive. The random numbers come from "generator".
typedef boost::uniform_smallint<base_generator_type> generator_type;
generator_type die_gen(generator, 1, 6);
#if !defined(__SUNPRO_CC) || (__SUNPRO_CC > 0x530)
// If you want to use an STL iterator interface, use iterator_adaptors.hpp.
// Unfortunately, this doesn't work on SunCC yet.
boost::generator_iterator_generator<generator_type>::type
die = boost::make_generator_iterator(die_gen);
for(int i = 0; i < 10; i++)
std::cout << *die++ << " ";
std::cout << '\n';
#endif
}
int main()
{
// Define a random number generator and initialize it with a reproducible
// seed.
// (The seed is unsigned, otherwise the wrong overload may be selected
// when using mt19937 as the base_generator_type.)
base_generator_type generator(42u);
std::cout << "10 samples of a uniform distribution in [0..1):\n";
// Define a uniform random number distribution which produces "double"
// values between 0 and 1 (0 inclusive, 1 exclusive).
boost::uniform_01<base_generator_type> uni(generator);
std::cout.setf(std::ios::fixed);
// You can now retrieve random numbers from that distribution by means
// of a STL Generator interface, i.e. calling the generator as a zero-
// argument function.
for(int i = 0; i < 10; i++)
std::cout << uni() << '\n';
/*
* Change seed to something else.
*
* Caveat: std::time(0) is not a very good truly-random seed. When
* called in rapid succession, it could return the same values, and
* thus the same random number sequences could ensue. If not the same
* values are returned, the values differ only slightly in the
* lowest bits. A linear congruential generator with a small factor
* wrapped in a uniform_smallint (see experiment) will produce the same
* values for the first few iterations. This is because uniform_smallint
* takes only the highest bits of the generator, and the generator itself
* needs a few iterations to spread the initial entropy from the lowest bits
* to the whole state.
*/
generator.seed(static_cast<unsigned int>(std::time(0)));
std::cout << "\nexperiment: roll a die 10 times:\n";
// You can save a generator's state by copy construction.
base_generator_type saved_generator = generator;
// When calling other functions which take a generator or distribution
// as a parameter, make sure to always call by reference (or pointer).
// Calling by value invokes the copy constructor, which means that the
// sequence of random numbers at the caller is disconnected from the
// sequence at the callee.
experiment(generator);
std::cout << "redo the experiment to verify it:\n";
experiment(saved_generator);
// After that, both generators are equivalent
assert(generator == saved_generator);
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
{
// You can save the generator state for future use. You can read the
// state back in at any later time using operator>>.
std::ofstream file("rng.saved", std::ofstream::trunc);
file << generator;
}
#endif
// Some compilers don't pay attention to std:3.6.1/5 and issue a
// warning here if "return 0;" is omitted.
return 0;
}

117
random_device.cpp
View File

@@ -1,117 +0,0 @@
/* boost random_device.cpp implementation
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#include <boost/nondet_random.hpp>
#include <string>
#include <cassert>
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
// A definition is required even for integral static constants
const bool boost::random_device::has_fixed_range;
const boost::random_device::result_type boost::random_device::min_value;
const boost::random_device::result_type boost::random_device::max_value;
#endif
#ifdef __linux__
// the default is the unlimited capacity device, using some secure hash
// try "/dev/random" for blocking when the entropy pool has drained
const char * const boost::random_device::default_token = "/dev/urandom";
/*
* This uses the POSIX interface for unbuffered reading.
* Using buffered std::istream would consume entropy which may
* not actually be used. Entropy is a precious good we avoid
* wasting.
*/
#if defined(__GNUC__) && defined(_CXXRT_STD_NAME)
// I have severe difficulty to get the POSIX includes to work with
// -fhonor-std and Dietmar Kühl's standard C++ library. Hack around that
// problem for now.
extern "C" {
static const int O_RDONLY = 0;
extern int open(const char *__file, int __oflag, ...);
extern int read(int __fd, __ptr_t __buf, size_t __nbytes);
extern int close(int __fd);
}
#else
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h> // open
#include <unistd.h> // read, close
#endif
#include <errno.h> // errno
#include <string.h> // strerror
#include <stdexcept> // std::invalid_argument
class boost::random_device::impl
{
public:
impl(const std::string & token) : path(token) {
fd = open(token.c_str(), O_RDONLY);
if(fd < 0)
error("cannot open");
}
~impl() { if(close(fd) < 0) error("could not close"); }
unsigned int next() {
unsigned int result;
long sz = read(fd, reinterpret_cast<char *>(&result), sizeof(result));
if(sz == -1)
error("error while reading");
else if(sz != sizeof(result)) {
errno = 0;
error("EOF while reading");
}
return result;
}
private:
void error(const std::string & msg) {
throw std::invalid_argument("boost::random_device: " + msg +
" random-number pseudo-device " + path +
": " + strerror(errno));
}
const std::string path;
int fd;
};
#endif // __linux__
boost::random_device::random_device(const std::string& token)
: pimpl(new impl(token))
{
assert(std::numeric_limits<result_type>::max() == max_value);
}
boost::random_device::~random_device()
{
// the complete class impl is now visible, so we're safe
// (see comment in random.hpp)
delete pimpl;
}
unsigned int boost::random_device::operator()()
{
return pimpl->next();
}

219
random_speed.cpp
View File

@@ -1,219 +0,0 @@
/* boost random_speed.cpp performance measurements
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*/
#include <iostream>
#include <cstdlib>
#include <string>
#include <boost/random.hpp>
#include <boost/progress.hpp>
/*
* Configuration Section
*/
// define if your C library supports the non-standard drand48 family
#undef HAVE_DRAND48
// define if you have the original mt19937int.c (with commented out main())
#undef HAVE_MT19937INT_C
// set to your CPU frequency in MHz
static const double cpu_frequency = 200 * 1e6;
/*
* End of Configuration Section
*/
/*
* General portability note:
* MSVC mis-compiles explicit function template instantiations.
* For example, f<A>() and f<B>() are both compiled to call f<A>().
* BCC is unable to implicitly convert a "const char *" to a std::string
* when using explicit function template instantiations.
*
* Therefore, avoid explicit function template instantiations.
*/
// provides a run-time configurable linear congruential generator, just
// for comparison
template<class IntType>
class linear_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = false;
#else
enum { has_fixed_range = false };
#endif
linear_congruential(IntType x0, IntType a, IntType c, IntType m)
: _x(x0), _a(a), _c(c), _m(m) { }
// compiler-generated copy ctor and assignment operator are fine
void seed(IntType x0, IntType a, IntType c, IntType m)
{ _x = x0; _a = a; _c = c; _m = m; }
void seed(IntType x0) { _x = x0; }
result_type operator()() { _x = (_a*_x+_c) % _m; return _x; }
result_type min() const { return _c == 0 ? 1 : 0; }
result_type max() const { return _m -1; }
private:
IntType _x, _a, _c, _m;
};
void show_elapsed(double end, int iter, const std::string & name)
{
double usec = end/iter*1e6;
double cycles = usec * cpu_frequency/1e6;
std::cout << name << ": "
<< usec*1e3 << " nsec/loop = "
<< cycles << " CPU cycles"
<< std::endl;
}
template<class Result, class RNG>
void timing(RNG & rng, int iter, const std::string& name, const Result&)
{
volatile Result tmp; // make sure we're not optimizing too much
boost::timer t;
for(int i = 0; i < iter; i++)
tmp = rng();
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
void timing_sphere(RNG & rng, int iter, const std::string & name)
{
boost::timer t;
for(int i = 0; i < iter; i++) {
// the special return value convention of uniform_on_sphere saves 20% CPU
const std::vector<double> & tmp = rng();
(void) tmp[0];
}
show_elapsed(t.elapsed(), iter, name);
}
template<class RNG>
void run(int iter, const std::string & name, const RNG &)
{
RNG rng;
std::cout << (RNG::has_fixed_range ? "fixed-range " : "");
// BCC has trouble with string autoconversion for explicit specializations
timing(rng, iter, std::string(name), 0l);
}
#ifdef HAVE_DRAND48
// requires non-standard C library support for srand48/lrand48
void run(int iter, const std::string & name, int)
{
std::srand48(1);
timing(std::lrand48, iter, name, 0l);
}
#endif
#ifdef HAVE_MT19937INT_C // requires the original mt19937int.c
extern "C" void sgenrand(unsigned long);
extern "C" unsigned long genrand();
void run(int iter, const std::string & name, float)
{
sgenrand(4357);
timing(genrand, iter, name, 0u);
}
#endif
template<class Gen>
void distrib(int iter, const std::string & name, const Gen &)
{
Gen gen;
boost::uniform_smallint<Gen> usmallint(gen, -2, 4);
timing(usmallint, iter, name + " uniform_smallint", 0);
boost::uniform_int<Gen> uint(gen, -2, 4);
timing(uint, iter, name + " uniform_int", 0);
boost::uniform_01<Gen> uni(gen);
timing(uni, iter, name + " uniform_01", 0.0);
boost::uniform_real<Gen> ur(gen, -5.3, 4.8);
timing(ur, iter, name + " uniform_real", 0.0);
boost::geometric_distribution<Gen> geo(gen, 0.5);
timing(geo, iter, name + " geometric", 0);
boost::triangle_distribution<Gen> tria(gen, 1, 2, 7);
timing(tria, iter, name + " triangle", 0.0);
boost::exponential_distribution<Gen> ex(gen, 3);
timing(ex, iter, name + " exponential", 0.0);
boost::normal_distribution<Gen,double> no2(gen);
timing(no2, iter, name + " normal polar", 0.0);
boost::lognormal_distribution<Gen,double> lnorm(gen, 1, 1);
timing(lnorm, iter, name + " lognormal", 0.0);
boost::uniform_on_sphere<Gen> usph(gen, 3);
timing_sphere(usph, iter/10, name + " uniform_on_sphere");
}
int main(int argc, char*argv[])
{
if(argc != 2) {
std::cerr << "usage: " << argv[0] << " iterations" << std::endl;
return 1;
}
// okay, it's ugly, but it's only used here
int iter =
#ifndef BOOST_NO_STDC_NAMESPACE
std::
#endif
atoi(argv[1]);
#if !defined(BOOST_NO_INT64_T) && \
!defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)
run(iter, "rand48", boost::rand48());
linear_congruential<boost::uint64_t>
lcg48(boost::uint64_t(1)<<16 | 0x330e,
boost::uint64_t(0xDEECE66DUL) | (boost::uint64_t(0x5) << 32), 0xB,
boost::uint64_t(1)<<48);
timing(lcg48, iter, "lrand48 run-time", 0l);
#endif
#ifdef HAVE_DRAND48
// requires non-standard C library support for srand48/lrand48
run(iter, "lrand48", 0); // coded for lrand48()
#endif
run(iter, "minstd_rand", boost::minstd_rand());
run(iter, "ecuyer combined", boost::ecuyer1988());
run(iter, "kreutzer1986", boost::kreutzer1986());
run(iter, "hellekalek1995 (inversive)", boost::hellekalek1995());
run(iter, "mt11213b", boost::mt11213b());
run(iter, "mt19937", boost::mt19937());
#ifdef HAVE_MT19937INT_C
// requires the original mt19937int.c
run<float>(iter, "mt19937 original"); // coded for sgenrand()/genrand()
#endif
distrib(iter, "minstd_rand", boost::minstd_rand());
distrib(iter, "kreutzer1986", boost::kreutzer1986());
distrib(iter, "mt19937", boost::mt19937());
}

298
random_test.cpp
View File

@@ -1,298 +0,0 @@
/* boost random_test.cpp various tests
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*/
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
#include <iterator>
#include <boost/random.hpp>
#include <boost/config.hpp>
#define BOOST_INCLUDE_MAIN
#include <boost/test/test_tools.hpp>
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std { using ::fabs; }
#endif
/*
* General portability note:
* MSVC mis-compiles explicit function template instantiations.
* For example, f<A>() and f<B>() are both compiled to call f<A>().
* BCC is unable to implicitly convert a "const char *" to a std::string
* when using explicit function template instantiations.
*
* Therefore, avoid explicit function template instantiations.
*/
/*
* Validate correct implementation
*/
template<class PRNG>
void validate(const std::string & name, const PRNG &)
{
std::cout << "validating " << name << ": ";
PRNG rng;
for(int i = 0; i < 9999; i++)
rng();
typename PRNG::result_type val = rng();
// make sure the validation function is a const member
bool result = const_cast<const PRNG&>(rng).validation(val);
// allow for a simple eyeball check for MSVC instantiation brokenness
// (if the numbers for all generators are the same, it's obviously broken)
std::cout << val << std::endl;
BOOST_TEST(result);
}
void validate_all()
{
using namespace boost;
#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
validate("rand48", rand48());
#endif
validate("minstd_rand", minstd_rand());
validate("minstd_rand0", minstd_rand0());
validate("ecuyer combined", ecuyer1988());
validate("mt19937", mt19937());
validate("kreutzer1986", kreutzer1986());
}
/*
* Check function signatures
*/
template<class URNG, class ResultType>
void instantiate_urng(const std::string & s, const URNG &, const ResultType &)
{
std::cout << "Basic tests for " << s << std::endl;
URNG urng;
int a[URNG::has_fixed_range ? 5 : 10]; // compile-time constant
(void) a; // avoid "unused" warning
typename URNG::result_type x1 = urng();
ResultType x2 = x1;
(void) &x2; // avoid "unused" warning
#ifndef BOOST_MSVC // MSVC brokenness
URNG urng2 = urng; // copy constructor
BOOST_TEST(urng == urng2); // operator==
urng();
urng2 = urng; // assignment
BOOST_TEST(urng == urng2);
#endif // BOOST_MSVC
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
// Streamable concept not supported for broken compilers
{
std::ofstream file("rng.tmp", std::ofstream::trunc);
file << urng;
}
// move forward
urng();
{
// restore old state
std::ifstream file("rng.tmp");
file >> urng;
}
#ifndef BOOST_MSVC // MSVC brokenness
BOOST_TEST(urng == urng2);
#endif // BOOST_MSVC
#endif // BOOST_NO_OPERATORS_IN_NAMESPACE
// instantiate various distributions with this URNG
boost::uniform_smallint<URNG> unismall(urng, 0, 11);
unismall();
boost::uniform_int<URNG> uni_int(urng, -200, 20000);
uni_int();
boost::uniform_real<URNG> uni_real(urng, 0, 2.1);
uni_real();
boost::bernoulli_distribution<URNG> ber(urng, 0.2);
ber();
boost::geometric_distribution<URNG> geo(urng, 0.8);
geo();
boost::triangle_distribution<URNG> tria(urng, 1, 1.5, 7);
tria();
boost::exponential_distribution<URNG> ex(urng, 5);
ex();
boost::normal_distribution<URNG> norm(urng);
norm();
boost::lognormal_distribution<URNG> lnorm(urng, 1, 1);
lnorm();
boost::uniform_on_sphere<URNG> usph(urng, 2);
usph();
}
void instantiate_all()
{
using namespace boost;
#if !defined(BOOST_NO_INT64_T) && !defined(BOOST_NO_INTEGRAL_INT64_T)
instantiate_urng("rand48", rand48(), 0);
rand48 rnd(boost::int32_t(5));
rand48 rnd2(boost::uint64_t(0x80000000) * 42);
rnd.seed(boost::int32_t(17));
rnd2.seed(boost::uint64_t(0x80000000) * 49);
#endif
instantiate_urng("minstd_rand0", minstd_rand0(), 0);
instantiate_urng("minstd_rand", minstd_rand(), 0);
minstd_rand mstd(42);
mstd.seed(17);
instantiate_urng("ecuyer1988", ecuyer1988(), 0);
instantiate_urng("kreutzer1986", kreutzer1986(), 0);
instantiate_urng("hellekalek1995", hellekalek1995(), 0);
instantiate_urng("mt11213b", mt11213b(), 0u);
instantiate_urng("mt19937", mt19937(), 0u);
mt19937 mt(boost::uint32_t(17)); // needs to be an exact type match for MSVC
int i = 42;
mt.seed(boost::uint32_t(i));
mt19937 mt2(mstd);
mt2.seed(mstd);
random_number_generator<mt19937> std_rng(mt2);
(void) std_rng(10);
}
/*
* A few equidistribution tests
*/
// yet to come...
template<class Generator>
void check_uniform_int(Generator & gen, int iter)
{
std::cout << "testing uniform_int(" << gen.min() << "," << gen.max()
<< ")" << std::endl;
int range = gen.max()-gen.min()+1;
std::vector<int> bucket(range);
for(int j = 0; j < iter; j++) {
int result = gen();
if(result < gen.min() || result > gen.max())
std::cerr << " ... delivers " << result << std::endl;
else
bucket[result-gen.min()]++;
}
int sum = 0;
// use a different variable name "k", because MSVC has broken "for" scoping
for(int k = 0; k < range; k++)
sum += bucket[k];
double avg = static_cast<double>(sum)/range;
double threshold = 2*avg/std::sqrt(static_cast<double>(iter));
for(int i = 0; i < range; i++) {
if(std::fabs(bucket[i] - avg) > threshold) {
// 95% confidence interval
std::cout << " ... has bucket[" << i << "] = " << bucket[i]
<< " (distance " << (bucket[i] - avg) << ")"
<< std::endl;
}
}
}
template<class Generator>
void test_uniform_int(Generator & gen)
{
typedef boost::uniform_int<Generator, int> int_gen;
// large range => small range (modulo case)
int_gen uint12(gen,1,2);
BOOST_TEST(uint12.min() == 1);
BOOST_TEST(uint12.max() == 2);
check_uniform_int(uint12, 100000);
int_gen uint16(gen,1,6);
check_uniform_int(uint16, 100000);
// test chaining to get all cases in operator()
typedef boost::uniform_int<int_gen, int> intint_gen;
// identity map
intint_gen uint01(uint12, 0, 1);
check_uniform_int(uint01, 100000);
// small range => larger range
intint_gen uint05(uint12, -3, 2);
check_uniform_int(uint05, 100000);
typedef boost::uniform_int<intint_gen, int> intintint_gen;
#if 0
// This takes a lot of time to run and is of questionable net effect:
// avoid for now.
// small => larger range, not power of two
// (for unknown reasons, this has noticeably uneven distribution)
intintint_gen uint1_49(uint05, 1, 49);
check_uniform_int(uint1_49, 500000);
#endif
// larger => small range, rejection case
intintint_gen uint1_4(uint05, 1, 4);
check_uniform_int(uint1_4, 100000);
}
#if defined(BOOST_MSVC) && _MSC_VER <= 1200
// These explicit instantiations are necessary, otherwise MSVC does
// find the <boost/operators.hpp> inline friends.
// We ease the typing with a suitable preprocessor macro.
#define INSTANT(x) \
template class boost::uniform_smallint<x>; \
template class boost::uniform_int<x>; \
template class boost::uniform_real<x>; \
template class boost::bernoulli_distribution<x>; \
template class boost::geometric_distribution<x>; \
template class boost::triangle_distribution<x>; \
template class boost::exponential_distribution<x>; \
template class boost::normal_distribution<x>; \
template class boost::uniform_on_sphere<x>; \
template class boost::lognormal_distribution<x>;
INSTANT(boost::minstd_rand0)
INSTANT(boost::minstd_rand)
INSTANT(boost::ecuyer1988)
INSTANT(boost::kreutzer1986)
INSTANT(boost::hellekalek1995)
INSTANT(boost::mt19937)
INSTANT(boost::mt11213b)
#endif
int test_main(int, char*[])
{
instantiate_all();
validate_all();
boost::mt19937 mt;
test_uniform_int(mt);
// bug report from Ken Mahler: This lead to an endless loop.
boost::minstd_rand r1;
boost::uniform_int<boost::minstd_rand, unsigned int> r2(r1, 0, 0xffffffff);
r2();
r2();
// bug report from Fernando Cacciola
boost::minstd_rand rnd;
boost::uniform_int<boost::minstd_rand> x(rnd,0,8361); // --> This CTOR loops forever.
return 0;
}

669
statistic_tests.cpp
View File

@@ -1,669 +0,0 @@
/* statistic_tests.cpp file
*
* Copyright Jens Maurer 2000, 2002
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
* Revision history
*/
/*
* NOTE: This is not part of the official boost submission. It exists
* only as a collection of ideas.
*/
#include <iostream>
#include <iomanip>
#include <string>
#include <functional>
#include <math.h> // lgamma is not in namespace std
#include <vector>
#include <algorithm>
#include <boost/cstdint.hpp>
#include <boost/random.hpp>
#include "statistic_tests.hpp"
#include "integrate.hpp"
namespace boost {
namespace random {
// Wikramaratna 1989 ACORN
template<class IntType, int k, IntType m, IntType val>
class additive_congruential
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = 0;
static const result_type max_value = m-1;
#else
enum {
has_fixed_range = true,
min_value = 0,
max_value = m-1
};
#endif
template<class InputIterator>
explicit additive_congruential(InputIterator start) { seed(start); }
template<class InputIterator>
void seed(InputIterator start)
{
for(int i = 0; i <= k; ++i, ++start)
values[i] = *start;
}
result_type operator()()
{
for(int i = 1; i <= k; ++i) {
IntType tmp = values[i-1] + values[i];
if(tmp >= m)
tmp -= m;
values[i] = tmp;
}
return values[k];
}
result_type validation() const { return val; }
private:
IntType values[k+1];
};
template<class IntType, int r, int s, IntType m, IntType val>
class lagged_fibonacci_int
{
public:
typedef IntType result_type;
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
static const bool has_fixed_range = true;
static const result_type min_value = 0;
static const result_type max_value = m-1;
#else
enum {
has_fixed_range = true,
min_value = 0,
max_value = m-1
};
#endif
explicit lagged_fibonacci_int(IntType start) { seed(start); }
template<class Generator>
explicit lagged_fibonacci_int(Generator & gen) { seed(gen); }
void seed(IntType start)
{
linear_congruential<uint32_t, 299375077, 0, 0, 0> init;
seed(init);
}
template<class Generator>
void seed(Generator & gen)
{
assert(r > s);
for(int i = 0; i < 607; ++i)
values[i] = gen();
current = 0;
lag = r-s;
}
result_type operator()()
{
result_type tmp = values[current] + values[lag];
if(tmp >= m)
tmp -= m;
values[current] = tmp;
++current;
if(current >= r)
current = 0;
++lag;
if(lag >= r)
lag = 0;
return tmp;
}
result_type validation() const { return val; }
private:
result_type values[r];
int current, lag;
};
} // namespace random
} // namespace boost
// distributions from Haertel's dissertation
// (additional parameterizations of the basic templates)
namespace Haertel {
typedef boost::random::linear_congruential<boost::uint64_t, 45965, 453816691,
(boost::uint64_t(1)<<31), 0> LCG_Af2;
typedef boost::random::linear_congruential<boost::uint64_t, 211936855, 0,
(boost::uint64_t(1)<<29)-3, 0> LCG_Die1;
typedef boost::random::linear_congruential<boost::uint32_t, 2824527309u, 0,
0, 0> LCG_Fis;
typedef boost::random::linear_congruential<boost::uint64_t, 950706376u, 0,
(boost::uint64_t(1)<<31)-1, 0> LCG_FM;
typedef boost::random::linear_congruential<boost::int32_t, 51081, 0,
2147483647, 0> LCG_Hae;
typedef boost::random::linear_congruential<boost::uint32_t, 69069, 1,
0, 0> LCG_VAX;
typedef boost::random::inversive_congruential<boost::int64_t, 240318, 197,
1000081, 0> NLG_Inv1;
typedef boost::random::inversive_congruential<boost::int64_t, 15707262,
13262967, (1<<24)-17, 0> NLG_Inv2;
typedef boost::random::inversive_congruential<boost::int32_t, 1, 1,
2147483647, 0> NLG_Inv4;
typedef boost::random::inversive_congruential<boost::int32_t, 1, 2,
1<<30, 0> NLG_Inv5;
typedef boost::random::additive_congruential<boost::int32_t, 6,
(1<<30)-35, 0> MRG_Acorn7;
typedef boost::random::lagged_fibonacci_int<boost::uint32_t, 607, 273,
0, 0> MRG_Fib2;
template<class Gen, class T>
inline void check_validation(Gen & gen, T value, const std::string & name)
{
for(int i = 0; i < 100000-1; ++i)
gen();
if(value != gen())
std::cout << name << ": validation failed" << std::endl;
}
// we have validation after 100000 steps with Haertel's generators
template<class Gen, class T>
void validate(T value, const std::string & name)
{
Gen gen(1234567);
check_validation(gen, value, name);
}
void validate_all()
{
validate<LCG_Af2>(183269031u, "LCG_Af2");
validate<LCG_Die1>(522319944u, "LCG_Die1");
validate<LCG_Fis>(-2065162233u, "LCG_Fis");
validate<LCG_FM>(581815473u, "LCG_FM");
validate<LCG_Hae>(28931709, "LCG_Hae");
validate<LCG_VAX>(1508154087u, "LCG_VAX");
validate<NLG_Inv2>(6666884, "NLG_Inv2");
validate<NLG_Inv4>(1521640076, "NLG_Inv4");
validate<NLG_Inv5>(641840839, "NLG_Inv5");
static const int acorn7_init[]
= { 1234567, 7654321, 246810, 108642, 13579, 97531, 555555 };
MRG_Acorn7 acorn7(acorn7_init);
check_validation(acorn7, 874294697, "MRG_Acorn7");
validate<MRG_Fib2>(1234567u, "MRG_Fib2");
}
} // namespace Haertel
double normal_density(double x)
{
const double pi = 3.14159265358979323846;
return 1/std::sqrt(2*pi) * std::exp(-x*x/2);
}
namespace std {
#ifdef _CXXRTCF_H__
using _CS_swamp::lgamma;
#elif defined __SGI_STL_PORT
using ::lgamma;
#endif
}
class chi_square_density : public std::unary_function<double, double>
{
public:
chi_square_density(int freedom)
: _exponent( static_cast<double>(freedom)/2-1 ),
_factor(1/(std::pow(2, _exponent+1) * std::exp(lgamma(_exponent+1))))
{ }
double operator()(double x)
{
return _factor*std::pow(x, _exponent)*std::exp(-x/2);
}
private:
double _exponent, _factor;
};
// computes F(x) or F(y) - F(x)
class chi_square_probability : public distribution_function<double>
{
public:
chi_square_probability(int freedom) : dens(freedom) {}
double operator()(double x) { return operator()(0, x); }
double operator()(double x, double y)
{ return trapezoid(dens, x, y, 1000); }
private:
chi_square_density dens;
};
class uniform_distribution : public distribution_function<double>
{
public:
uniform_distribution(double from, double to) : from(from), to(to)
{ assert(from < to); }
double operator()(double x)
{
if(x < from)
return 0;
else if(x > to)
return 1;
else
return (x-from)/(to-from);
}
double operator()(double x, double delta)
{ return operator()(x+delta) - operator()(x); }
private:
double from, to;
};
class test_environment;
class test_base
{
protected:
explicit test_base(test_environment & env) : environment(env) { }
void check(double val) const;
private:
test_environment & environment;
};
class equidistribution_test : test_base
{
public:
equidistribution_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "equidistribution: " << std::flush;
equidistribution_experiment equi(classes);
uniform_smallint<RNG> uint_linear(rng, 0, classes-1);
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi, uint_linear, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi, uint_linear, n1), 2*n2));
std::cout << " 2D: " << std::flush;
equidistribution_2d_experiment equi_2d(classes);
unsigned int root = static_cast<unsigned int>(std::sqrt(double(classes)));
assert(root * root == classes);
uniform_smallint<RNG> uint_square(rng, 0, root-1);
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi_2d, uint_square, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(equi_2d, uint_square, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class ks_equidistribution_test : test_base
{
public:
ks_equidistribution_test(test_environment & env, unsigned int classes)
: test_base(env),
test_distrib_chi_square(kolmogorov_smirnov_probability(5000),
classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "KS: " << std::flush;
// generator_reference_t<RNG> gen_ref(rng);
RNG& gen_ref(rng);
kolmogorov_experiment ks(n1);
uniform_distribution ud(rng.min_value, rng.max_value);
check(run_experiment(test_distrib_chi_square,
ks_experiment_generator(ks, gen_ref, ud), n2));
check(run_experiment(test_distrib_chi_square,
ks_experiment_generator(ks, gen_ref, ud), 2*n2));
}
private:
distribution_experiment test_distrib_chi_square;
};
class runs_test : test_base
{
public:
runs_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "runs: up: " << std::flush;
runs_experiment<true> r_up(classes);
// generator_reference_t<RNG> gen_ref(rng);
RNG& gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_up, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_up, gen_ref, n1), 2*n2));
std::cout << " down: " << std::flush;
runs_experiment<false> r_down(classes);
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_down, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(r_down, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class gap_test : test_base
{
public:
gap_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "gaps: " << std::flush;
gap_experiment gap(classes, 0.2, 0.8);
// generator_reference_t<RNG> gen_ref(rng);
RNG& gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(gap, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(gap, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class poker_test : test_base
{
public:
poker_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "poker: " << std::flush;
poker_experiment poker(8, classes);
uniform_smallint<RNG> usmall(rng, 0, 7);
check(run_experiment(test_distrib_chi_square,
experiment_generator(poker, usmall, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(poker, usmall, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class coupon_collector_test : test_base
{
public:
coupon_collector_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(classes-1), high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "coupon collector: " << std::flush;
coupon_collector_experiment coupon(5, classes);
uniform_smallint<RNG> usmall(rng, 0, 4);
check(run_experiment(test_distrib_chi_square,
experiment_generator(coupon, usmall, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(coupon, usmall, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class permutation_test : test_base
{
public:
permutation_test(test_environment & env, unsigned int classes,
unsigned int high_classes)
: test_base(env), classes(classes),
test_distrib_chi_square(chi_square_probability(fac<int>(classes)-1),
high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "permutation: " << std::flush;
permutation_experiment perm(classes);
// generator_reference_t<RNG> gen_ref(rng);
RNG& gen_ref(rng);
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), 2*n2));
std::cout << std::endl;
}
private:
unsigned int classes;
distribution_experiment test_distrib_chi_square;
};
class maximum_test : test_base
{
public:
maximum_test(test_environment & env, unsigned int high_classes)
: test_base(env),
test_distrib_chi_square(kolmogorov_smirnov_probability(1000),
high_classes)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "maximum-of-t: " << std::flush;
maximum_experiment<RNG> mx(rng, n1, 5);
check(run_experiment(test_distrib_chi_square, mx, n2));
check(run_experiment(test_distrib_chi_square, mx, 2*n2));
std::cout << std::endl;
}
private:
distribution_experiment test_distrib_chi_square;
};
class birthday_test : test_base
{
public:
birthday_test(test_environment & env)
: test_base(env)
{ }
template<class RNG>
void run(RNG & rng, int n1, int n2)
{
using namespace boost;
std::cout << "birthday spacing: " << std::flush;
uniform_int<RNG> uni(rng, 0, (1<<25));
birthday_spacing_experiment bsp(4, 512, (1<<25));
std::cout << run_experiment(bsp, uni, n1);
#if 0
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), n2));
check(run_experiment(test_distrib_chi_square,
experiment_generator(perm, gen_ref, n1), 2*n2));
#endif
std::cout << std::endl;
}
};
class test_environment
{
public:
static const int classes = 20;
explicit test_environment(double confid)
: confidence(confid),
confidence_chi_square_quantil(quantil(chi_square_density(classes-1), 0, confidence, 1e-4)),
test_distrib_chi_square6(chi_square_probability(7-1), classes),
ksequi_test(*this, classes),
equi_test(*this, 100, classes),
rns_test(*this, 7, classes),
gp_test(*this, 7, classes),
pk_test(*this, 5, classes),
cpn_test(*this, 15, classes),
perm_test(*this, 5, classes),
max_test(*this, classes),
bday_test(*this)
{
std::cout << "Confidence level: " << confid
<< "; 1-alpha = " << (1-confid)
<< "; chi_square(" << (classes-1)
<< ", " << confidence_chi_square_quantil
<< ") = "
<< chi_square_probability(classes-1)(0, confidence_chi_square_quantil)
<< std::endl;
}
bool check_confidence(double val, double chi_square_conf) const
{
std::cout << val;
bool result = (val <= chi_square_conf);
if(!result) {
std::cout << "* [";
double prob = (val > 10*chi_square_conf ? 1 :
chi_square_probability(classes-1)(0, val));
std::cout << (1-prob) << "]";
}
std::cout << " " << std::flush;
return result;
}
bool check(double chi_square_value) const
{
return check_confidence(chi_square_value, confidence_chi_square_quantil);
}
template<class RNG>
void run_test(const std::string & name)
{
using namespace boost;
std::cout << "Running tests on " << name << std::endl;
RNG rng(1234567);
typedef boost::uniform_01<RNG> UGen;
#if 1
ksequi_test.run(rng, 5000, 250);
equi_test.run(rng, 5000, 250);
rns_test.run(rng, 100000, 250);
gp_test.run(rng, 10000, 250);
pk_test.run(rng, 5000, 250);
cpn_test.run(rng, 500, 250);
perm_test.run(rng, 1200, 250);
max_test.run(rng, 1000, 250);
#endif
bday_test.run(rng, 1000, 150);
std::cout << std::endl;
}
private:
double confidence;
double confidence_chi_square_quantil;
distribution_experiment test_distrib_chi_square6;
ks_equidistribution_test ksequi_test;
equidistribution_test equi_test;
runs_test rns_test;
gap_test gp_test;
poker_test pk_test;
coupon_collector_test cpn_test;
permutation_test perm_test;
maximum_test max_test;
birthday_test bday_test;
};
void test_base::check(double val) const
{
environment.check(val);
}
void print_ks_table()
{
std::cout.setf(std::ios::fixed);
std::cout.precision(5);
static const double all_p[] = { 0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99 };
for(int n = 0; n <= 10000; (n < 55 ? ++n : n *= 10)) {
std::cout << std::setw(4) << n << " ";
for(unsigned int i = 0; i < sizeof(all_p)/sizeof(all_p[0]); ++i) {
std::cout << std::setw(8)
<< (n == 0 ? all_p[i] :
invert_monotone_inc(kolmogorov_smirnov_probability(n), all_p[i], 0, 10))
<< " ";
}
std::cout << std::endl;
}
}
int main()
{
// Haertel::validate_all();
test_environment env(0.99);
env.run_test<boost::minstd_rand>("minstd_rand");
env.run_test<boost::mt19937>("mt19937");
env.run_test<Haertel::LCG_Af2>("LCG_Af2");
env.run_test<Haertel::LCG_Die1>("LCG_Die1");
env.run_test<Haertel::LCG_Fis>("LCG_Fis");
env.run_test<Haertel::LCG_FM>("LCG_FM");
env.run_test<Haertel::LCG_Hae>("LCG_Hae");
env.run_test<Haertel::LCG_VAX>("LCG_VAX");
env.run_test<Haertel::NLG_Inv1>("NLG_Inv1");
env.run_test<Haertel::NLG_Inv2>("NLG_Inv2");
env.run_test<Haertel::NLG_Inv4>("NLG_Inv4");
env.run_test<Haertel::NLG_Inv5>("NLG_Inv5");
}

641
statistic_tests.hpp
View File

@@ -1,641 +0,0 @@
/* statistic_tests.hpp header file
*
* Copyright Jens Maurer 2000
* Permission to use, copy, modify, sell, and distribute this software
* is hereby granted without fee provided that the above copyright notice
* appears in all copies and that both that copyright notice and this
* permission notice appear in supporting documentation,
*
* Jens Maurer makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* $Id$
*
*/
#ifndef STATISTIC_TESTS_HPP
#define STATISTIC_TESTS_HPP
#include <stdexcept>
#include <iterator>
#include <vector>
#include <boost/limits.hpp>
#include <algorithm>
#include <cmath>
#include <boost/random.hpp>
#include <boost/config.hpp>
#ifdef BOOST_MSVC
namespace std
{
inline double pow(double a, double b) { return ::pow(a,b); }
inline double ceil(double x) { return ::ceil(x); }
} // namespace std
#endif
template<class T>
inline T fac(int k)
{
T result = 1;
for(T i = 2; i <= k; ++i)
result *= i;
return result;
}
template<class T>
T binomial(int n, int k)
{
if(k < n/2)
k = n-k;
T result = 1;
for(int i = k+1; i<= n; ++i)
result *= i;
return result / fac<T>(n-k);
}
template<class T>
T stirling2(int n, int m)
{
T sum = 0;
for(int k = 0; k <= m; ++k)
sum += binomial<T>(m, k) * std::pow(double(k), n) *
( (m-k)%2 == 0 ? 1 : -1);
return sum / fac<T>(m);
}
/*
* Experiments which create an empirical distribution in classes,
* suitable for the chi-square test.
*/
// std::floor(gen() * classes)
class experiment_base
{
public:
experiment_base(int cls) : _classes(cls) { }
unsigned int classes() const { return _classes; }
protected:
unsigned int _classes;
};
class equidistribution_experiment : public experiment_base
{
public:
explicit equidistribution_experiment(unsigned int classes)
: experiment_base(classes) { }
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
assert(f.min() == 0 &&
static_cast<unsigned int>(f.max()) == classes()-1);
for(int i = 0; i < n; ++i)
count(f());
}
double probability(int i) const { return 1.0/classes(); }
};
// two-dimensional equidistribution experiment
class equidistribution_2d_experiment : public equidistribution_experiment
{
public:
explicit equidistribution_2d_experiment(unsigned int classes)
: equidistribution_experiment(classes) { }
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
unsigned int range = f.max()+1;
assert(f.min() == 0 && range*range == classes());
for(int i = 0; i < n; ++i) {
int y1 = f();
int y2 = f();
count(y1 + range * y2);
}
}
};
// distribution experiment: assume a probability density and
// count events so that an equidistribution results.
class distribution_experiment : public equidistribution_experiment
{
public:
template<class UnaryFunction>
distribution_experiment(UnaryFunction probability , unsigned int classes)
: equidistribution_experiment(classes), limit(classes)
{
for(unsigned int i = 0; i < classes-1; ++i)
limit[i] = invert_monotone_inc(probability, (i+1)*0.05, 0, 1000);
limit[classes-1] = std::numeric_limits<double>::infinity();
#if 0
std::cout << __PRETTY_FUNCTION__ << ": ";
for(unsigned int i = 0; i < classes; ++i)
std::cout << limit[i] << " ";
std::cout << std::endl;
#endif
}
template<class NumberGenerator, class Counter>
void run(NumberGenerator f, Counter & count, int n) const
{
for(int i = 0; i < n; ++i) {
limits_type::const_iterator it =
std::lower_bound(limit.begin(), limit.end(), f());
count(it-limit.begin());
}
}
private:
typedef std::vector<double> limits_type;
limits_type limit;
};
// runs-up/runs-down experiment
template<bool up>
class runs_experiment : public experiment_base
{
public:
explicit runs_experiment(unsigned int classes) : experiment_base(classes) { }
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
result_type init = (up ? f.min() : f.max());
result_type previous = init;
unsigned int length = 0;
for(int i = 0; i < n; ++i) {
result_type val = f();
if(up ? previous <= val : previous >= val) {
previous = val;
++length;
} else {
count(std::min(length, classes())-1);
length = 0;
previous = init;
// don't use this value, so that runs are independent
}
}
}
double probability(unsigned int r) const
{
if(r == classes()-1)
return 1.0/fac<double>(classes());
else
return static_cast<double>(r+1)/fac<double>(r+2);
}
};
// gap length experiment
class gap_experiment : public experiment_base
{
public:
gap_experiment(unsigned int classes, double alpha, double beta)
: experiment_base(classes), alpha(alpha), beta(beta) { }
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
double range = f.max() - f.min() + 1.0;
result_type low = static_cast<result_type>(alpha * range);
result_type high = static_cast<result_type>(beta * range);
unsigned int length = 0;
for(int i = 0; i < n; ) {
result_type value = f() - f.min();
if(value < low || value > high)
++length;
else {
count(std::min(length, classes()-1));
length = 0;
++i;
}
}
}
double probability(unsigned int r) const
{
double p = beta-alpha;
if(r == classes()-1)
return std::pow(1-p, static_cast<double>(r));
else
return p * std::pow(1-p, static_cast<double>(r));
}
private:
double alpha, beta;
};
// poker experiment
class poker_experiment : public experiment_base
{
public:
poker_experiment(unsigned int d, unsigned int k)
: experiment_base(k), range(d)
{
assert(range > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(range-1));
std::vector<result_type> v(classes());
for(int i = 0; i < n; ++i) {
for(unsigned int j = 0; j < classes(); ++j)
v[j] = f();
std::sort(v.begin(), v.end());
result_type prev = v[0];
int r = 1; // count different values in v
for(unsigned int i = 1; i < classes(); ++i) {
if(prev != v[i]) {
prev = v[i];
++r;
}
}
count(r-1);
}
}
double probability(unsigned int r) const
{
++r; // transform to 1 <= r <= 5
double result = range;
for(unsigned int i = 1; i < r; ++i)
result *= range-i;
return result / std::pow(range, static_cast<double>(classes())) *
stirling2<double>(classes(), r);
}
private:
unsigned int range;
};
// coupon collector experiment
class coupon_collector_experiment : public experiment_base
{
public:
coupon_collector_experiment(unsigned int d, unsigned int cls)
: experiment_base(cls), d(d)
{
assert(d > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(d-1));
std::vector<bool> occurs(d);
for(int i = 0; i < n; ++i) {
occurs.assign(d, false);
unsigned int r = 0; // length of current sequence
int q = 0; // number of non-duplicates in current set
for(;;) {
result_type val = f();
++r;
if(!occurs[val]) { // new set element
occurs[val] = true;
++q;
if(q == d)
break; // one complete set
}
}
count(std::min(r-d, classes()-1));
}
}
double probability(unsigned int r) const
{
if(r == classes()-1)
return 1-fac<double>(d)/std::pow(d, static_cast<double>(d+classes()-2))*
stirling2<double>(d+classes()-2, d);
else
return fac<double>(d)/std::pow(d, static_cast<double>(d+r)) *
stirling2<double>(d+r-1, d-1);
}
private:
int d;
};
// permutation test
class permutation_experiment : public equidistribution_experiment
{
public:
permutation_experiment(unsigned int t)
: equidistribution_experiment(fac<int>(t)), t(t)
{
assert(t > 1);
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
std::vector<result_type> v(t);
for(int i = 0; i < n; ++i) {
std::generate_n(v.begin(), t, f);
int x = 0;
for(int r = t-1; r > 0; r--) {
std::vector<result_type>::iterator it =
std::max_element(v.begin(), v.begin()+r+1);
x = (r+1)*x + (it-v.begin());
std::iter_swap(it, v.begin()+r);
}
count(x);
}
}
private:
int t;
};
// birthday spacing experiment test
class birthday_spacing_experiment : public experiment_base
{
public:
birthday_spacing_experiment(unsigned int d, int n, int m)
: experiment_base(d), n(n), m(m)
{
}
template<class UniformRandomNumberGenerator, class Counter>
void run(UniformRandomNumberGenerator f, Counter & count, int n_total) const
{
typedef typename UniformRandomNumberGenerator::result_type result_type;
assert(std::numeric_limits<result_type>::is_integer);
assert(f.min() == 0);
assert(f.max() == static_cast<result_type>(m-1));
for(int j = 0; j < n_total; j++) {
std::vector<result_type> v(n);
std::generate_n(v.begin(), n, f);
std::sort(v.begin(), v.end());
std::vector<result_type> spacing(n);
for(int i = 0; i < n-1; i++)
spacing[i] = v[i+1]-v[i];
spacing[n-1] = v[0] + m - v[n-1];
std::sort(spacing.begin(), spacing.end());
unsigned int k = 0;
for(int i = 0; i < n-1; ++i) {
if(spacing[i] == spacing[i+1])
++k;
}
count(std::min(k, classes()-1));
}
}
double probability(unsigned int r) const
{
assert(classes() == 4);
assert(m == (1<<25));
assert(n == 512);
static const double prob[] = { 0.368801577, 0.369035243, 0.183471182,
0.078691997 };
return prob[r];
}
private:
int n, m;
};
/*
* Misc. helper functions.
*/
template<class Float>
struct distribution_function
{
typedef Float result_type;
typedef Float argument_type;
typedef Float first_argument_type;
typedef Float second_argument_type;
};
// computes P(K_n <= t) or P(t1 <= K_n <= t2). See Knuth, 3.3.1
class kolmogorov_smirnov_probability : public distribution_function<double>
{
public:
kolmogorov_smirnov_probability(int n)
: approx(n > 50), n(n), sqrt_n(std::sqrt(double(n)))
{
if(!approx)
n_n = std::pow(static_cast<double>(n), n);
}
double operator()(double t) const
{
if(approx) {
return 1-std::exp(-2*t*t)*(1-2.0/3.0*t/sqrt_n);
} else {
t *= sqrt_n;
double sum = 0;
for(int k = static_cast<int>(std::ceil(t)); k <= n; k++)
sum += binomial<double>(n, k) * std::pow(k-t, k) *
std::pow(t+n-k, n-k-1);
return 1 - t/n_n * sum;
}
}
double operator()(double t1, double t2) const
{ return operator()(t2) - operator()(t1); }
private:
bool approx;
int n;
double sqrt_n;
double n_n;
};
/*
* Experiments for generators with continuous distribution functions
*/
class kolmogorov_experiment
{
public:
kolmogorov_experiment(int n) : n(n), ksp(n) { }
template<class NumberGenerator, class Distribution>
double run(NumberGenerator gen, Distribution distrib) const
{
const int m = n;
typedef std::vector<double> saved_temp;
saved_temp a(m,1.0), b(m,0);
std::vector<int> c(m,0);
for(int i = 0; i < n; ++i) {
double val = gen();
double y = distrib(val);
int k = static_cast<int>(std::floor(m*y));
if(k >= m)
--k; // should not happen
a[k] = std::min(a[k], y);
b[k] = std::max(b[k], y);
++c[k];
}
double kplus = 0, kminus = 0;
int j = 0;
for(int k = 0; k < m; ++k) {
if(c[k] > 0) {
kminus = std::max(kminus, a[k]-j/static_cast<double>(n));
j += c[k];
kplus = std::max(kplus, j/static_cast<double>(n) - b[k]);
}
}
kplus *= std::sqrt(double(n));
kminus *= std::sqrt(double(n));
// std::cout << "k+ " << kplus << " k- " << kminus << std::endl;
return kplus;
}
double probability(double x) const
{
return ksp(x);
}
private:
int n;
kolmogorov_smirnov_probability ksp;
};
// maximum-of-t test (KS-based)
template<class UniformRandomNumberGenerator>
class maximum_experiment
{
public:
typedef UniformRandomNumberGenerator base_type;
maximum_experiment(base_type & f, int n, int t) : f(f), ke(n), t(t)
{ }
double operator()() const
{
double res = ke.run(generator(f, t),
std::bind2nd(std::ptr_fun(static_cast<double (*)(double, double)>(&std::pow)), t));
return res;
}
private:
struct generator {
generator(base_type & f, int t) : f(f), t(t) { }
double operator()()
{
double mx = f();
for(int i = 1; i < t; ++i)
mx = std::max(mx, f());
return mx;
}
private:
boost::uniform_01<base_type> f;
int t;
};
base_type & f;
kolmogorov_experiment ke;
int t;
};
// compute a chi-square value for the distribution approximation error
template<class ForwardIterator, class UnaryFunction>
typename UnaryFunction::result_type
chi_square_value(ForwardIterator first, ForwardIterator last,
UnaryFunction probability)
{
typedef std::iterator_traits<ForwardIterator> iter_traits;
typedef typename iter_traits::value_type counter_type;
typedef typename UnaryFunction::result_type result_type;
unsigned int classes = std::distance(first, last);
result_type sum = 0;
counter_type n = 0;
for(unsigned int i = 0; i < classes; ++first, ++i) {
counter_type count = *first;
n += count;
sum += (count/probability(i)) * count; // avoid overflow
}
#if 0
for(unsigned int i = 0; i < classes; ++i) {
// std::cout << (n*probability(i)) << " ";
if(n * probability(i) < 5)
std::cerr << "Not enough test runs for slot " << i
<< " p=" << probability(i) << ", n=" << n
<< std::endl;
}
#endif
// std::cout << std::endl;
// throw std::invalid_argument("not enough test runs");
return sum/n - n;
}
template<class RandomAccessContainer>
class generic_counter
{
public:
explicit generic_counter(unsigned int classes) : container(classes, 0) { }
void operator()(int i)
{
assert(i >= 0);
assert(static_cast<unsigned int>(i) < container.size());
++container[i];
}
typename RandomAccessContainer::const_iterator begin() const
{ return container.begin(); }
typename RandomAccessContainer::const_iterator end() const
{ return container.end(); }
private:
RandomAccessContainer container;
};
// chi_square test
template<class Experiment, class Generator>
double run_experiment(const Experiment & experiment, Generator gen, int n)
{
generic_counter<std::vector<int> > v(experiment.classes());
experiment.run(gen, v, n);
return chi_square_value(v.begin(), v.end(),
std::bind1st(std::mem_fun_ref(&Experiment::probability),
experiment));
}
// number generator with experiment results (for nesting)
template<class Experiment, class Generator>
class experiment_generator_t
{
public:
experiment_generator_t(const Experiment & exper, Generator & gen, int n)
: experiment(exper), generator(gen), n(n) { }
double operator()() { return run_experiment(experiment, generator, n); }
private:
const Experiment & experiment;
Generator & generator;
int n;
};
template<class Experiment, class Generator>
experiment_generator_t<Experiment, Generator>
experiment_generator(const Experiment & e, Generator & gen, int n)
{
return experiment_generator_t<Experiment, Generator>(e, gen, n);
}
template<class Experiment, class Generator, class Distribution>
class ks_experiment_generator_t
{
public:
ks_experiment_generator_t(const Experiment & exper, Generator & gen,
const Distribution & distrib)
: experiment(exper), generator(gen), distribution(distrib) { }
double operator()() { return experiment.run(generator, distribution); }
private:
const Experiment & experiment;
Generator & generator;
Distribution distribution;
};
template<class Experiment, class Generator, class Distribution>
ks_experiment_generator_t<Experiment, Generator, Distribution>
ks_experiment_generator(const Experiment & e, Generator & gen,
const Distribution & distrib)
{
return ks_experiment_generator_t<Experiment, Generator, Distribution>
(e, gen, distrib);
}
#endif /* STATISTIC_TESTS_HPP */