From 2ec49b3e5a2c262d2f8af051d003b44bb2f7c080 Mon Sep 17 00:00:00 2001 From: nobody Date: Sat, 9 Sep 2000 10:20:25 +0000 Subject: [PATCH] This commit was manufactured by cvs2svn to create branch 'boost-graph-library'. [SVN r7698] --- .gitattributes | 96 ++ histogram.cpp | 147 +++ include/boost/nondet_random.hpp | 72 ++ include/boost/random.hpp | 1493 +++++++++++++++++++++++++++++++ index.html | 112 +++ integrate.hpp | 83 ++ nondet_random.html | 121 +++ nondet_random_speed.cpp | 67 ++ random-concepts.html | 344 +++++++ random-distributions.html | 738 +++++++++++++++ random-generators.html | 785 ++++++++++++++++ random-misc.html | 159 ++++ random_demo.cpp | 80 ++ random_device.cpp | 108 +++ random_speed.cpp | 217 +++++ random_test.cpp | 309 +++++++ statistic_tests.cpp | 662 ++++++++++++++ statistic_tests.hpp | 628 +++++++++++++ 18 files changed, 6221 insertions(+) create mode 100644 .gitattributes create mode 100644 histogram.cpp create mode 100644 include/boost/nondet_random.hpp create mode 100644 include/boost/random.hpp create mode 100644 index.html create mode 100644 integrate.hpp create mode 100644 nondet_random.html create mode 100644 nondet_random_speed.cpp create mode 100644 random-concepts.html create mode 100644 random-distributions.html create mode 100644 random-generators.html create mode 100644 random-misc.html create mode 100644 random_demo.cpp create mode 100644 random_device.cpp create mode 100644 random_speed.cpp create mode 100644 random_test.cpp create mode 100644 statistic_tests.cpp create mode 100644 statistic_tests.hpp diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..3e84d7c --- /dev/null +++ b/.gitattributes @@ -0,0 +1,96 @@ +* text=auto !eol svneol=native#text/plain +*.gitattributes text svneol=native#text/plain + +# Scriptish formats +*.bat text svneol=native#text/plain +*.bsh text svneol=native#text/x-beanshell +*.cgi text svneol=native#text/plain +*.cmd text svneol=native#text/plain +*.js text svneol=native#text/javascript +*.php text svneol=native#text/x-php +*.pl text svneol=native#text/x-perl +*.pm text svneol=native#text/x-perl +*.py text svneol=native#text/x-python +*.sh eol=lf svneol=LF#text/x-sh +configure eol=lf svneol=LF#text/x-sh + +# Image formats +*.bmp binary svneol=unset#image/bmp +*.gif binary svneol=unset#image/gif +*.ico binary svneol=unset#image/ico +*.jpeg binary svneol=unset#image/jpeg +*.jpg binary svneol=unset#image/jpeg +*.png binary svneol=unset#image/png +*.tif binary svneol=unset#image/tiff +*.tiff binary svneol=unset#image/tiff +*.svg text svneol=native#image/svg%2Bxml + +# Data formats +*.pdf binary svneol=unset#application/pdf +*.avi binary svneol=unset#video/avi +*.doc binary svneol=unset#application/msword +*.dsp text svneol=crlf#text/plain +*.dsw text svneol=crlf#text/plain +*.eps binary svneol=unset#application/postscript +*.gz binary svneol=unset#application/gzip +*.mov binary svneol=unset#video/quicktime +*.mp3 binary svneol=unset#audio/mpeg +*.ppt binary svneol=unset#application/vnd.ms-powerpoint +*.ps binary svneol=unset#application/postscript +*.psd binary svneol=unset#application/photoshop +*.rdf binary svneol=unset#text/rdf +*.rss text svneol=unset#text/xml +*.rtf binary svneol=unset#text/rtf +*.sln text svneol=native#text/plain +*.swf binary svneol=unset#application/x-shockwave-flash +*.tgz binary svneol=unset#application/gzip +*.vcproj text svneol=native#text/xml +*.vcxproj text svneol=native#text/xml +*.vsprops text svneol=native#text/xml +*.wav binary svneol=unset#audio/wav +*.xls binary svneol=unset#application/vnd.ms-excel +*.zip binary svneol=unset#application/zip + +# Text formats +.htaccess text svneol=native#text/plain +*.bbk text svneol=native#text/xml +*.cmake text svneol=native#text/plain +*.css text svneol=native#text/css +*.dtd text svneol=native#text/xml +*.htm text svneol=native#text/html +*.html text svneol=native#text/html +*.ini text svneol=native#text/plain +*.log text svneol=native#text/plain +*.mak text svneol=native#text/plain +*.qbk text svneol=native#text/plain +*.rst text svneol=native#text/plain +*.sql text svneol=native#text/x-sql +*.txt text svneol=native#text/plain +*.xhtml text svneol=native#text/xhtml%2Bxml +*.xml text svneol=native#text/xml +*.xsd text svneol=native#text/xml +*.xsl text svneol=native#text/xml +*.xslt text svneol=native#text/xml +*.xul text svneol=native#text/xul +*.yml text svneol=native#text/plain +boost-no-inspect text svneol=native#text/plain +CHANGES text svneol=native#text/plain +COPYING text svneol=native#text/plain +INSTALL text svneol=native#text/plain +Jamfile text svneol=native#text/plain +Jamroot text svneol=native#text/plain +Jamfile.v2 text svneol=native#text/plain +Jamrules text svneol=native#text/plain +Makefile* text svneol=native#text/plain +README text svneol=native#text/plain +TODO text svneol=native#text/plain + +# Code formats +*.c text svneol=native#text/plain +*.cpp text svneol=native#text/plain +*.h text svneol=native#text/plain +*.hpp text svneol=native#text/plain +*.ipp text svneol=native#text/plain +*.tpp text svneol=native#text/plain +*.jam text svneol=native#text/plain +*.java text svneol=native#text/plain diff --git a/histogram.cpp b/histogram.cpp new file mode 100644 index 0000000..fc00f75 --- /dev/null +++ b/histogram.cpp @@ -0,0 +1,147 @@ +/* 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 free 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 +#include +#include +#include +#include +#include +#include + + +void plot_histogram(const std::vector& 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 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 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 +void histogram(RNG base, int samples, double from, double to, + const std::string & name) +{ + typedef squaresum_result, double > SRNG; + SRNG gen((sum_result(base))); + const int nSlots = 60; + std::vector 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 +void histograms() +{ + PRNG rng; + using namespace boost; + histogram(uniform_01(rng), 100000, -0.5, 1.5, "uniform_01(0,1)"); + histogram(bernoulli_distribution(rng, 0.2), 100000, -0.5, 1.5, + "bernoulli(0.2)"); + histogram(triangle_distribution(rng, 1, 2, 8), 100000, 0, 10, + "triangle(1,2,8)"); + histogram(geometric_distribution(rng, 5.0/6.0), 100000, 0, 10, + "geometric(5/6)"); + histogram(exponential_distribution(rng, 0.3), 100000, 0, 10, + "exponential(0.3)"); + histogram(cauchy_distribution(rng), 100000, -5, 5, + "cauchy"); + histogram(lognormal_distribution(rng, 3, 2), 100000, 0, 10, + "lognormal"); + histogram(normal_distribution(rng), 100000, -3, 3, + "normal"); + histogram(normal_distribution(rng, 0.5, 0.5), 100000, -3, 3, + "normal(0.5, 0.5)"); +} + + +int main() +{ + histograms(); +} + diff --git a/include/boost/nondet_random.hpp b/include/boost/nondet_random.hpp new file mode 100644 index 0000000..67a337d --- /dev/null +++ b/include/boost/nondet_random.hpp @@ -0,0 +1,72 @@ +/* boost nondet_random.hpp header file + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 + * 2000-02-18 Portability fixes (thanks to Beman Dawes) + */ + +#ifndef BOOST_NONDET_RANDOM_HPP +#define BOOST_NONDET_RANDOM_HPP + +#include // std::abs +#include // std::min +#include +#include +#include // noncopyable +#include // compile-time integral limits + +namespace boost { + +// use some OS service to generate non-deterministic random numbers +class random_device : noncopyable +{ +public: + typedef unsigned int result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = true; + static const result_type min_value = integer_traits::const_min; + static const result_type max_value = integer_traits::const_max; +#else + enum { + has_fixed_range = true, + min_value = integer_traits::const_min, + max_value = integer_traits::const_max + }; +#endif + result_type min() const { return integer_traits::min(); } + result_type max() const { return integer_traits::max(); } + explicit random_device(const std::string& token = default_token); + ~random_device(); + unsigned int operator()(); + +private: + static const char * const default_token; + + /* + * std:5.3.5/5 [expr.delete]: "If the object being deleted has incomplete + * class type at the point of deletion and the complete class has a + * non-trivial destructor [...], the behavior is undefined". + * This disallows the use of scoped_ptr<> with pimpl-like classes + * having a non-trivial destructor. + */ + class impl; + impl * pimpl; +}; + + +// TODO: put Schneier's Yarrow-160 algorithm here. + +} // namespace boost + +#endif /* BOOST_NONDET_RANDOM_HPP */ diff --git a/include/boost/random.hpp b/include/boost/random.hpp new file mode 100644 index 0000000..65e34bb --- /dev/null +++ b/include/boost/random.hpp @@ -0,0 +1,1493 @@ +/* boost random.hpp header file + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 + * 2000-02-18 Portability fixes (thanks to Beman Dawes) + * 2000-02-21 shuffle_output, inversive_congruential_schrage, + * generator_iterator, uniform_smallint + * 2000-02-23 generic modulus arithmetic helper, removed *_schrage classes, + * implemented Streamable and EqualityComparable concepts for + * generators, added Bernoulli distribution and Box-Muller + * transform + * 2000-03-01 cauchy, lognormal, triangle distributions; fixed + * uniform_smallint; renamed gaussian to normal distribution + * 2000-03-05 implemented iterator syntax for distribution functions + * 2000-04-21 removed some optimizations for better BCC/MSVC compatibility + * 2000-05-10 adapted to BCC and MSVC + * 2000-06-13 incorporated review results + * 2000-07-06 moved basic templates from namespace detail to random + */ + +#ifndef BOOST_RANDOM_HPP +#define BOOST_RANDOM_HPP + +#include +#include +#include // for std::log etc. +#include // std::ptrdiff_t +#include +#include +#include +#include // std::transform, std::copy +#include // std::iterator_category + +#include +#include // exact bit-sized integers +#include +#include + + +/* + * Hack around various namespace challenged compilers + */ +#ifdef BOOST_NO_STDC_NAMESPACE +namespace std { +#ifndef BOOST_MSVC + using ::exp; + using ::log; + using ::sqrt; + using ::sin; + using ::cos; + using ::tan; +#else + // Even the "using" workaround fails with MSVC. We reimplement the + // functions. This will get into trouble if we need these functions in + // another file. Future plan: Move this stuff into a separate header file. + inline double exp(double x) { return ::exp(x); } + inline double log(double x) { return ::log(x); } + inline double sqrt(double x) { return ::sqrt(x); } + inline double sin(double x) { return ::sin(x); } + inline double cos(double x) { return ::cos(x); } + inline double tan(double x) { return ::tan(x); } +#endif +} // namespace std +#endif + +namespace boost { + +#ifdef __GNUC__ +// Special gcc workaround: gcc does not yet support using-declarations +// in template classes (confirmed by gcc author Martin v. Loewis) + using std::sqrt; + using std::exp; + using std::log; + using std::sin; + using std::cos; + using std::tan; +#endif + +/* + * End of hacks for various namespace challenged compilers + */ + + +/* + * Some random number generators. + */ + +namespace random { + +/* + * Some random number generators require modular arithmetic. Put + * everything we need here. + * IntType must be an integral type. + */ + +template +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 if(traits::is_signed) + return add_signed(x, c); + else { + // difficult + assert(!"const_mod::add with c too large"); + return 0; + } + } + + 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 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 add_signed(IntType x, IntType c) + { + x += (c-m); + if(x < 0) + 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 + assert(c > 0); + assert(m > 0); + assert(boost::integer_traits::is_signed); + 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 +{ + 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(); +}; + +// the modulus is some power of 2: rely partly on machine overflow handling +// we only specialize for rand48 at the moment +#ifdef BOOST_STDINT_H_HAS_UINT64_T +template<> +class const_mod +{ + 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_STDINT_H_HAS_UINT64_T */ + + +// compile-time configurable linear congruential generator +template +class linear_congruential +{ +public: + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const IntType multiplier = a; + static const IntType increment = c; + static const IntType modulus = m; + 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 + enum { + multiplier = a, + increment = c, + modulus = m, + has_fixed_range = false + }; +#endif + 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::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, 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==(linear_congruential x, linear_congruential y) + { return x._x == y._x; } +#else + // Use a member function; Streamable concept not supported. + bool operator==(linear_congruential rhs) const + { return _x == rhs._x; } +#endif +private: + IntType _x; +}; + +} // namespace random + +// validation values from the publications +typedef random::linear_congruential minstd_rand0; +typedef random::linear_congruential minstd_rand; + + +#ifdef BOOST_STDINT_H_HAS_UINT64_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::const_max; +#else + enum { has_fixed_range = false }; +#endif + int32_t min() const { return 0; } + int32_t max() const { return std::numeric_limits::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 lcf; + static uint64_t cnv(int32_t x) + { return (static_cast(x) << 16) | 0x330e; } +}; +#endif /* BOOST_STDINT_H_HAS_UINT64_T */ + + +namespace random { + +// L'Ecuyer 1988 +template +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, + random::linear_congruential, + /* unknown */ 0> ecuyer1988; + + +namespace random { + +// Carter Bays and S.D. Durham 1979 +template +class shuffle_output +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef typename base_type::result_type result_type; + +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; + static const int buffer_size = k; +#else + enum { + has_fixed_range = false, + buffer_size = k + }; +#endif + 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 + explicit shuffle_output(T seed) : _rng(seed) { init(); } + explicit shuffle_output(const base_type & rng) : _rng(rng) { init(); } + template + 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(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() { + // TODO: should be a compile-time assert + assert(std::numeric_limits::is_integer); + 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; +}; + +} // namespace random + +// validation by experiment from Harry Erwin's generator.h (private e-mail) +typedef random::shuffle_output< + random::linear_congruential, + 97, uint32_t, 139726> kreutzer1986; + + +namespace random { + +// Eichenauer and Lehn 1986 +template +class inversive_congruential +{ +public: + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; + static const result_type min_value = (b == 0 ? 1 : 0); + static const result_type max_value = p-1; + static const result_type multiplier = a; + static const result_type increment = b; + static const result_type modulus = p; +#else + enum { + has_fixed_range = false, + multiplier = a, + increment = b, + modulus = p + }; +#endif + result_type min() const { return b == 0 ? 1 : 0; } + result_type max() const { return p-1; } + + explicit inversive_congruential(IntType y0 = 1) : value(y0) { + if(b == 0) + assert(y0 > 0); + // TODO: should be compile-time asserts + assert(b >= 0); + assert(p > 1); + assert(a >= 1); + } + void seed(IntType y0) { value = y0; if(b == 0) assert(y0 > 0); } + IntType operator()() { + typedef const_mod 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; +}; + +} // namespace random + +typedef random::inversive_congruential hellekalek1995; + +namespace random { + +// http://www.math.keio.ac.jp/matumoto/emt.html +template +class mersenne_twister +{ +public: + typedef DataType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const int state_size = n; + static const int shift_size = m; + static const int mask_bits = r; + static const DataType parameter_a = a; + static const int output_u = u; + static const int output_s = s; + static const DataType output_b = b; + static const int output_t = t; + static const DataType output_c = c; + static const int output_l = l; + static const bool has_fixed_range = true; + static const result_type min_value = integer_traits::const_min; + static const result_type max_value = integer_traits::const_max; +#else + enum { + state_size = n, shift_size = m, mask_bits = r, + parameter_a = a, output_u = u, output_s = s, output_b = b, + output_t = t, output_c = c, output_l = l, has_fixed_range = false }; +#endif + result_type min() const { return std::numeric_limits::min(); } + result_type max() const { return std::numeric_limits::max(); } + + mersenne_twister() { seed(); } + explicit mersenne_twister(DataType value) { seed(value); } + template + explicit mersenne_twister(Generator & gen) { seed(gen); } + // compiler-generated copy ctor and assignment operator are fine + void seed() { seed(4357u); } + void seed(DataType value) { + random::linear_congruential + gen(value); + seed(gen); + } + template + void seed(Generator & gen) { + // TODO: should be a compile-time assert + assert(!std::numeric_limits::is_signed); + // 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 x) const { return val == x; } + +#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(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]; +}; + +template +void mersenne_twister::twist() +{ + const int upper_mask = static_cast(-1) << r; + const int 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 +inline typename mersenne_twister::result_type +mersenne_twister::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 mt11213b; + +// validation by experiment from mt19937.c +typedef random::mersenne_twister mt19937; + + + +/* + * Some decorators providing additional interfaces and misc. functionality + */ + +namespace detail { + +// TODO: check out if this is really useful. + +template +class generator_reference_t +{ +public: + typedef typename NumberGenerator::result_type result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + // doing this properly requires partial specialization, avoid for now + static const bool has_fixed_range = false; +#else + enum { has_fixed_range = false }; +#endif + result_type min() const { return generator->min(); } + result_type max() const { return generator->max(); } + explicit generator_reference_t(NumberGenerator & gen) : generator(&gen) { } + // compiler-generated copy ctor and assignment operator are fine + + result_type operator()() { return (*generator)(); } +private: + NumberGenerator * generator; // no ownership, is never 0 +}; + +template +generator_reference_t generator_reference(Generator & gen) +{ + return generator_reference_t(gen); +} + +} // namespace detail + + +template +class generator_iterator + : equality_comparable >, + incrementable >, + dereferenceable, + typename Generator::result_type> +{ +public: + // Note: deriving from std::iterator<> is useless + typedef typename Generator::result_type value_type; + typedef std::ptrdiff_t difference_type; + typedef const typename Generator::result_type * pointer; + typedef const typename Generator::result_type & reference; + typedef std::input_iterator_tag iterator_category; + + explicit generator_iterator(Generator & g) : gen(g), value(gen()) { } + generator_iterator& operator++() { value = gen(); return *this; } + reference operator*() const { return value; } + + friend bool operator==(const generator_iterator& x, + const generator_iterator& y) + { return x.gen == y.gen; } +private: + Generator & gen; + value_type value; +}; + + + +/* + * Some distribution functions + */ + +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. + */ + +template +int equal_signed_unsigned(T1 x, T2 y) +{ + typedef std::numeric_limits x_traits; + typedef std::numeric_limits y_traits; + if(x_traits::is_signed == y_traits::is_signed) { + // no difference in signedness, cast to the larger type + if(x_traits::digits > y_traits::digits) + return x == static_cast(y); + else + return static_cast(x) == y; + } else if(x_traits::is_signed) { + // y must be unsigned, i.e. non-negative + if(x < 0) + return false; + else + return static_cast(x) == y; + } else if(y_traits::is_signed) { + if(y < 0) + return false; + else + return x == static_cast(y); + } +} + +template +int lessthan_signed_unsigned(T1 x, T2 y) +{ + typedef std::numeric_limits x_traits; + typedef std::numeric_limits y_traits; + if(x_traits::is_signed == y_traits::is_signed) { + // no difference in signedness, everything ok + if(x_traits::digits > y_traits::digits) + return x < static_cast(y); + else + return static_cast(x) < y; + } else if(x_traits::is_signed) { + // y must be unsigned, i.e. non-negative + if(x < 0) + return true; + else + return static_cast(x) < y; + } else if(y_traits::is_signed) { + // x must be unsigned, i.e. non-negative + if(y < 0) + return false; + else + return x < static_cast(y); + } +} + +} // namespace random + +// must be in boost namespace, otherwise the inline friend trick fails +template +class generator_iterator_mixin_adapter + : incrementable, decrementable, + equality_comparable +{ +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(*this); } + value_type v; +}; + + +// uniform integer distribution on a small range [min, max] +template +class uniform_smallint + : public generator_iterator_mixin_adapter< + uniform_smallint, IntType > +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; +#else + enum { has_fixed_range = false }; +#endif + 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; +}; + +template +uniform_smallint:: +uniform_smallint(base_type & rng, IntType min, IntType max) + : _rng(rng), _min(min), _max(max), + _range(static_cast(_max-_min)+1), _factor(1) +{ + // TODO: should be a compile-time assert + assert(std::numeric_limits::is_integer); + + // 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::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; + } + iterator_init(); // initialize iterator interface +} + + +// uniform integer distribution on [min, max] +template +class uniform_int + : public generator_iterator_mixin_adapter< + uniform_int, IntType > +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef IntType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; +#else + enum { has_fixed_range = false }; +#endif + 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) { + // TODO: should be a compile-time assert + assert(std::numeric_limits::is_integer); + iterator_init(); + } + 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; +}; + +template +inline IntType uniform_int::operator()() +{ + if(random::equal_signed_unsigned(_range, _brange)) { + // this will probably never happen in real life + // basically nothing to do; just take care we don't overflow / underflow + return static_cast(_rng() - _bmin) + _min; + } else if(random::lessthan_signed_unsigned(_brange, _range)) { + // use rejection method to handle things like 0..3 --> 0..4 + // note: this still does not have perfect efficiency + for(;;) { + // we have to concatenate several invocations of the base RNG + result_type result = 0; + for(result_type mult = 1; + mult-1 <= _range; + mult *= static_cast(_brange)+1) { + result += (_rng() - _bmin) * 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(_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(_range)) + return result + _min; + } + } + } +} + + +// a model for RandomNumberGenerator std:25.2.11 [lib.alg.random.shuffle] +template +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) { + // TODO: should be compile-time asserts + // MSVC requires the typedef workaround + typedef typename base_type::result_type base_result_type; + assert(std::numeric_limits::is_integer); + assert(std::numeric_limits::is_integer); + } + // compiler-generated copy ctor is fine + // assignment is disallowed because there is a reference member + + result_type operator()(argument_type n) { + return uniform_int(_rng, 0, n-1)(); + } + +private: + base_type & _rng; +}; + + +// 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 uniform_01 + : public generator_iterator_mixin_adapter< + uniform_01, RealType > +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef RealType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; +#else + enum { has_fixed_range = false }; +#endif + explicit uniform_01(base_type & rng) : _rng(rng) { + // TODO: should be a compile-time assert + assert(!std::numeric_limits::is_integer); + iterator_init(); + } + // compiler-generated copy ctor is fine + // assignment is disallowed because there is a reference member + + result_type operator()() { + return static_cast(_rng() - _rng.min()) / + (static_cast(_rng.max()-_rng.min()) + + (std::numeric_limits::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; +}; + + +// uniform distribution on a real range +template +class uniform_real + : public generator_iterator_mixin_adapter< + uniform_real, RealType> +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef RealType result_type; +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION + static const bool has_fixed_range = false; +#else + enum { has_fixed_range = false }; +#endif + uniform_real(base_type & rng, RealType min, RealType max) + : _rng(rng), _min(min), _max(max) { + // TODO: should be a compile-time assert + assert(!std::numeric_limits::is_integer); + iterator_init(); + } + // 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 _rng; + RealType _min, _max; +}; + + +// Bernoulli distribution: p(true) = p, p(false) = 1-p (boolean) +template +class bernoulli_distribution + : public generator_iterator_mixin_adapter< + bernoulli_distribution, bool> +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef bool result_type; + bernoulli_distribution(base_type & rng, double p) + : _rng(rng), + _threshold(static_cast + (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; +}; + +// geometric distribution: p(i) = (1-p) * pow(p, i-1) (integer) +template +class geometric_distribution + : public generator_iterator_mixin_adapter< + geometric_distribution, IntType> +{ +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); + using std::log; + _log_p = log(p); + iterator_init(); + } + // compiler-generated copy ctor is fine + // uniform_01 cannot be assigned, neither can this class + + result_type operator()() { + using std::log; + 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 _rng; + typename uniform_01::result_type _log_p; +}; + + +// triangle distribution, with a smallest, b most probable, and c largest +// value. +template +class triangle_distribution + : public generator_iterator_mixin_adapter< + triangle_distribution, RealType> +{ +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) + { + d3 = sqrt(d3); + p1 = sqrt(p1); + assert(_a <= _b && _b <= _c); + iterator_init(); + } + // compiler-generated copy ctor is fine + // uniform_01 cannot be assigned, neither can this class + result_type operator()() { + using std::sqrt; + 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 _rng; + result_type _a, _b, _c; + result_type d1, d2, d3, q1, p1; +}; + + +// exponential distribution: p(x) = lambda * exp(-lambda * x) +template +class exponential_distribution + : public generator_iterator_mixin_adapter< + exponential_distribution, RealType> +{ +public: + typedef UniformRandomNumberGenerator base_type; + typedef RealType result_type; + + exponential_distribution(base_type& rng, result_type lambda) + : _rng(rng), _lambda(lambda) { assert(lambda > 0); iterator_init(); } + // compiler-generated copy ctor is fine + // uniform_01 cannot be assigned, neither can this class + result_type operator()() { + using std::log; + 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 _rng; + const result_type _lambda; +}; + + +// Cauchy distribution: p(x) = sigma/(pi*(sigma**2 + (x-median)**2)) +template +class cauchy_distribution + : public generator_iterator_mixin_adapter< + cauchy_distribution, RealType> +{ +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) { iterator_init(); } + // compiler-generated copy constructor is fine + // uniform_01 cannot be assigned, neither can this class + result_type operator()() + { + const double pi = 3.14159265358979323846; + using std::tan; + 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 _rng; + result_type _median, _sigma; +}; + +// deterministic polar method, uses trigonometric functions +template +class normal_distribution + : public generator_iterator_mixin_adapter< + normal_distribution, RealType> +{ +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); iterator_init(); } + // compiler-generated copy constructor is fine + // uniform_01 cannot be assigned, neither can this class + result_type operator()() + { + // allow for Koenig lookup + using std::sqrt; using std::log; using std::sin; using std::cos; + 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 _rng; + const result_type _mean, _sigma; + result_type _r1, _r2, _cached_rho; + bool _valid; +}; + +template +class lognormal_distribution + : public generator_iterator_mixin_adapter< + lognormal_distribution, RealType> +{ +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); + iterator_init(); + } + // compiler-generated copy constructor is fine + // normal_distribution cannot be assigned, neither can this class + result_type operator()() { + // allow for Koenig lookup + using std::exp; + 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 _rng; +}; + +template > +class uniform_on_sphere + : public generator_iterator_mixin_adapter< + uniform_on_sphere, + std::vector > +{ +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) { iterator_init(); } + // 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; + } + using std::sqrt; + // for all i: result[i] /= sqrt(sqsum) + std::transform(_container.begin(), _container.end(), _container.begin(), + std::bind2nd(std::divides(), 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 _rng; + result_type _container; + const int _dim; +}; + +} // namespace boost + +#endif // BOOST_RANDOM_HPP diff --git a/index.html b/index.html new file mode 100644 index 0000000..db86c74 --- /dev/null +++ b/index.html @@ -0,0 +1,112 @@ + + + + + +Boost Random Number Library + + + + + + + + + + + + + +
c++boost.gif (8819 bytes)HomeLibrariesPeopleFAQMore
+ +

Boost Random Number Library

+ +By Jens Maurer + +

+You should read the +concepts documentation before +moving on. + +

+For a quick start, it may be sufficient to have a look at +random_demo.cpp. + +

+ + + + + + + + + + + + + + + +
HeaderContents
<boost/random.hpp> +

documentation for generators
+documentation for distributions
+documentation for miscellaneous classes

Pseudo-random number generators (linear congruential, inversive +congruential, and others), distribution functions (uniform, normal, +exponential, and others) and miscellaneous decorator classes such as +random_number_generator for use in simulations, +games, and testing.
<boost/nondet_random.hpp> +

documentation

Non-deterministic random number generators.
+ +

+ +An extensive test suite for the pseudo-random number generators and +distributions is available as +random_test.cpp. +

+All of the above, this documentation, and additional examples are +available in the random.zip archive for easy +download. + +

Rationale

+ +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. + + +

History and Acknowledgements

+ +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 enum workaround. Dave +Abrahams highlighted quantization issues. +

+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 min_rand +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. +

+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. +

+Gary Powell contributed suggestions for code cleanliness. Dave +Abrahams and Howard Hinnant suggested to move the basic generator +templates from namespace boost::detail to +boost::random. + +

+

+Jens Maurer, 2000-07-06 + + + diff --git a/integrate.hpp b/integrate.hpp new file mode 100644 index 0000000..2e894ec --- /dev/null +++ b/integrate.hpp @@ -0,0 +1,83 @@ +/* integrate.hpp header file + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 + */ + +#ifndef INTEGRATE_HPP +#define INTEGRATE_HPP + +#include + +template +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 +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 +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 +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::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 */ diff --git a/nondet_random.html b/nondet_random.html new file mode 100644 index 0000000..8759296 --- /dev/null +++ b/nondet_random.html @@ -0,0 +1,121 @@ + + + + + +Boost RNG Library - Non-Deterministic Random Number Generators + + + + +

c++boost.gif (8819 bytes)Header +<boost/nondet_random.hpp>

+ + + +

Header<boost/nondet_random.hpp> +Synopsis

+ +
+namespace boost {
+  class random_device;
+} // namespace boost
+
+ + +

Class random_device

+ +

Synopsis

+ +
+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()();
+};
+
+ +

Description

+ +Class random_device models a +non-deterministic random number +generator. +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 random_device +must not be implemented. See +
+"Randomness Recommendations for Security", D. Eastlake, S. +Crocker, J. Schiller, Network Working Group, RFC 1750, December 1994 +
+for further discussions. + +

+Note: 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. + +

Members

+ +
explicit random_device(const std::string& token = default_token)
+ +Effects: Constructs a random_device, +optionally using the given token as an access +specification (for example, a URL) to some implementation-defined +service for monitoring a stochastic process. + +

Implementation Note for Linux

+ +On the Linux operating system, token 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, std::ios_base::failure is +thrown. By default, random_device uses the +/dev/urandom pseudo-device to retrieve the random +numbers. Another option would be to specify the +/dev/random pseudo-device, which blocks on reads if the +entropy pool has no more random bits available. + + +

Performance

+ +The test program nondet_random_speed.cpp +measures the execution times of the +nondet_random.hpp 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. + +

+ + + +
classtime per invocation [usec]
random_device92.0
+ +

+The measurement error is estimated at +/- 1 usec. + +

+

+Jens Maurer, 2000-06-19 + + + + \ No newline at end of file diff --git a/nondet_random_speed.cpp b/nondet_random_speed.cpp new file mode 100644 index 0000000..74efc3e --- /dev/null +++ b/nondet_random_speed.cpp @@ -0,0 +1,67 @@ +/* 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 free 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 +#include +#include +#include + +// 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 +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 +void run(int iter, const std::string & name) +{ + RNG rng; + timing(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(dev, iter, "random_device"); +#endif + + return 0; +} diff --git a/random-concepts.html b/random-concepts.html new file mode 100644 index 0000000..2a8d0d4 --- /dev/null +++ b/random-concepts.html @@ -0,0 +1,344 @@ + + + + + +Boost Random Number Library Concepts + + + + +

Random Number Generator Library Concepts

+ + +

Introduction

+ +Random numbers are required in a number of different problem domains, +such as +
    +
  • numerics (simulation, Monte-Carlo integration) +
  • games (non-deterministic enemy behavior) +
  • security (key generation) +
  • testing (random coverage in white-box tests) +
+ +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 +
+"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 +
+ +Depending on the requirements of the problem domain, different +variations of random number generators are appropriate: + +
    +
  • non-deterministic random number generator +
  • pseudo-random number generator +
  • quasi-random number generator +
+ +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. + +

+ +The goals for this library are the following: +

    +
  • allow easy integration of third-party random-number generators +
  • define a validation interface for the generators +
  • provide easy-to-use front-end classes which model popular +distributions +
  • provide maximum efficiency +
  • allow control on quantization effects in front-end processing +(not yet done) +
+ + +

Number Generator

+ +A number generator is a function object (std:20.3 +[lib.function.objects]) that takes zero arguments. Each call to +operator() returns a number. + + +In the following table, X denotes a number generator +class returning objects of type T, and u is +a value of X. + +

+ + + + + + + + +
NumberGenerator +requirements
expressionreturn typepre/post-condition
X::result_typeTstd::numeric_limits<T>::is_specialized is true, +T is LessThanComparable
u.operator()()T-
+ +

+ +Note: The NumberGenerator requirements do not impose any +restrictions on the characteristics of the returned numbers. + + +

Uniform Random Number Generator

+ +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. + +

+The tight lower bound of some (finite) set S is the (unique) +member l in S, so that for all v in S, l <= v holds. Likewise, the +tight upper bound of some (finite) set S is the (unique) +member u in S, so that for all v in S, v <= u holds. + +

+In the following table, X denotes a number generator +class returning objects of type T, and v is +a const value of X. + +

+ + + + + + + + + + + + + + + +
UniformRandomNumberGenerator +requirements
expressionreturn typepre/post-condition
X::has_fixed_rangeboolcompile-time constant; if true, the range on which +the random numbers are uniformly distributed is known at compile-time +and members min_value and max_value +exist. Note: This flag may also be false due to +compiler limitations.
X::min_valueTcompile-time constant; min_value is equal to +v.min()
X::max_valueTcompile-time constant; max_value is equal to +v.max()
v.min()Ttight lower bound on the set of all values returned by +operator(). The return value of this function shall not +change during the lifetime of the object.
v.max()Tif std::numeric_limits<T>::is_integer, tight +upper bound on the set of all values returned by +operator(), otherwise, the smallest representable number +larger than the tight upper bound on the set of all values returned by +operator(). In any case, the return value of this +function shall not change during the lifetime of the +object.
+ +

+The member functions min, max, and +operator() shall have amortized constant time complexity. + +

+Note: For integer generators (i.e. integer T), +the generated values x fulfill min() <= x <= +max(), for non-integer generators (i.e. non-integer +T), the generated values x fulfill +min() <= x < max(). +
+Rationale: The range description with min and +max 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. +
+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. +

+ +Note: The UniformRandomNumberGenerator concept does not +require operator()(long) and thus it does not fulfill the +RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle]) +requirements. Use the +random_number_generator +adapter for that. +
+ +Rationale: operator()(long) is not provided, +because mapping the output of some generator with integer range to a +different integer range is not trivial. + + +

Non-deterministic Uniform Random Number +Generator

+ +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. +

+ +The class +random_device +is a model for a non-deterministic random number generator. + +

+ +Note: 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. + + +

Pseudo-Random Number Generator

+ +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. +

+ +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. +

+ +Note: Because the state of a pseudo-random number generator +is necessarily finite, the sequence of numbers returned by the +generator will loop eventually. + +

+In addition to the UniformRandomNumberGenerator requirements, a +pseudo-random number generator has some additional requirements. In +the following table, X denotes a pseudo-random number +generator class returning objects of type T, +x is a value of T, u is a value +of X, and v is a const value of +X. + +

+ + + + + + + + + + + +
PseudoRandomNumberGenerator +requirements
expressionreturn typepre/post-condition
X()-creates a generator in some implementation-defined +state. Note: Several generators thusly created may possibly +produce dependent or identical sequences of random numbers.
explicit X(...)-creates a generator with user-provided state; the implementation +shall specify the the constructor argument(s)
u.seed(...)voidsets the current state according to the argument(s); at least +functions with the same signature as the non-default +constructor(s) shall be provided. +
v.validation(x)boolcompares the pre-computed and hardcoded 10001th element in the +generator's random number sequence with x. The generator +must have been constructed by its default constructor and +seed must not have been called for the validation to +be meaningful. +
+ +

+ +Note: The seed member function is similar to the +assign member function in STL containers. However, the +naming did not seem appropriate. + +

+ +Classes which model a pseudo-random number generator shall also model +EqualityComparable, i.e. implement operator==. Two +pseudo-random number generators are defined to be equivalent +if they both return an identical sequence of numbers starting from a +given state. +

+Classes which model a pseudo-random number generator should also model +the Streamable concept, i.e. implement operator<< +and operator>>. If so, +operator<< writes all current state of the +pseudo-random number generator to the given ostream so +that operator>> can restore the state at a later +time. The state shall be written in a platform-independent manner, +but it is assumed that the locales 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. +

+ +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-const) +reference. + +

+ +The classes +rand48, +minstd_rand, +and +mt19937 +are models for a pseudo-random number generator. + +

+Note: This type of random-number generator is useful for +numerics, games and testing. The none-zero arguments constructor(s) +and the seed() 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. + + +

Quasi-Random Number Generators

+ +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). + +

+ +Note: Quasi-random number generators are useful for +Monte-Carlo integrations where specially crafted sequences of random +numbers will make the approximation converge faster. + +

+ +[Does anyone have a model?] + +

+

+Jens Maurer, 2000-02-23 + + + diff --git a/random-distributions.html b/random-distributions.html new file mode 100644 index 0000000..16541f6 --- /dev/null +++ b/random-distributions.html @@ -0,0 +1,738 @@ + + + + + + +Boost Random Number Library Distributions + + + + +

Random Number Library Distributions

+ + + +

Introduction

+ +In addition to the random number +generators, this library provides distribution functions which map +one distribution (often a uniform distribution provided by some +generator) to another. + +

+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. + +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
distributionexplanationexample
uniform_smallintdiscrete uniform distribution on a small set of integers (much +smaller than the range of the underlying generator)drawing from an urn
uniform_intdiscrete uniform distribution on a set of integers; the underlying +generator may be called several times to gather enough randomness for +the outputdrawing from an urn
uniform_01continuous uniform distribution on the range [0,1); important +basis for other distributions-
uniform_realcontinuous uniform distribution on some range [min, max) of real +numbersfor the range [0, 2pi): randomly dropping a stick and measuring +its angle in radiants (assuming the angle is uniformly +distributed)
bernoulli_distributionBernoulli experiment: discrete boolean valued distribution with +configurable probabilitytossing a coin (p=0.5)
geometric_distributionmeasures distance between outcomes of repeated Bernoulli experimentsthrowing a die several times and counting the number of tries +until a "6" appears for the first time
triangle_distribution??
exponential_distributionexponential distributionmeasuring the inter-arrival time of alpha particles emitted by +radioactive matter
normal_distributioncounts outcomes of (infinitely) repeated Bernoulli experimentstossing a coin 10000 times and counting how many front sides are shown
lognormal_distributionlognormal distribution (sometimes used in simulations)measuring the job completion time of an assembly line worker
uniform_on_sphereuniform distribution on a unit sphere of arbitrary dimensionchoosing a random point on Earth (assumed to be a sphere) where to +spend the next vacations
+ +

+ +The template parameters of the distribution functions are always in +the order +

    +
  • Underlying source of random numbers +
  • If applicable, return type, with a default to a reasonable type. +
+ +

+All distribution functions satisfy the input iterator requirements +(std:24.1.1 [lib.input.iterators]) in addition to the NumberGenerator +requirements. After an invocation of operator(), the +effects of invocations of operator* are undefined until +the next call to operator++. + +

+In this description, I have refrained from documenting those members +in detail which are already defined in the +concept documentation. + + +

Synopsis of the distributions in header +<boost/random.hpp>

+ +
+namespace boost {
+  template<class UniformRandomNumberGenerator, class IntType = int>
+  class uniform_smallint;
+  template<class UniformRandomNumberGenerator, class IntType = int>
+  class uniform_int;
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class uniform_01;
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class uniform_real;
+
+  // discrete distributions
+  template<class UniformRandomNumberGenerator>
+  class bernoulli_distribution;
+  template<class UniformRandomNumberGenerator, class IntType = int>
+  class geometric_distribution;
+
+  // continuous distributions
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class triangle_distribution;
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class exponential_distribution;
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class normal_distribution;
+  template<class UniformRandomNumberGenerator, class RealType = double>
+  class lognormal_distribution;
+  template<class UniformRandomNumberGenerator, class RealType = double,
+    class Cont = std::vector<RealType> >
+  class uniform_on_sphere;
+}
+
+ +

Class template +uniform_smallint

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, class IntType = int>
+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;
+};
+
+ +

Description

+ +The distribution function uniform_smallint models a +uniform random number +generator. 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. +

+Let rout=(max-min+1) the desired range of integer numbers, +and let rbase 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 rout will be +pout(i) = 1/rout. Likewise, assume a uniform +distribution on rbase for the underlying source of random +numbers, i.e. pbase(i) = 1/rbase. Let +pout_s(i) denote the random distribution generated by +uniform_smallint. Then the sum over all i in +rout of (pout_s(i)/pout(i) +-1)2 shall not exceed +rout/rbase2 (rbase mod +rout)(rout - rbase mod +rout). +

+ +The template parameter UniformRandomNumberGenerator shall +denote a class which models a uniform random number generator. +Additionally, UniformRandomNumberGenerator::result_type +shall denote an integral type. The template parameter +IntType shall denote an integer-like value type. + +

+Note: The property above is the square sum of the relative +differences in probabilities between the desired uniform distribution +pout(i) and the generated distribution +pout_s(i). The property can be fulfilled with the +calculation (base_rng mod rout), as follows: Let r = +rbase mod rout. The base distribution on +rbase is folded onto the range rout. The +numbers i < r have assigned (rbase div +rout)+1 numbers of the base distribution, the rest has only +(rbase div rout). Therefore, +pout_s(i) = ((rbase div rout)+1) / +rbase for i < r and pout_s(i) = +(rbase div rout)/rbase otherwise. +Substituting this in the above sum formula leads to the desired +result. +

+Note: The upper bound for (rbase mod rout)(rout - rbase +mod rout) is rout2/4. Regarding the upper bound for the square +sum of the relative quantization error of rout3/(4*rbase2), it +seems wise to either choose rbase so that rbase > 10*rout2 or +ensure that rbase is divisible by rout. + + +

Members

+ +
uniform_smallint(base_type & rng, IntType min, IntType max)
+ +Effects: Constructs a uniform_smallint +functor with the uniform random number generator rng as +the underlying source of random numbers. min and +max are the lower and upper bounds of the output range, +respectively. + + +

Class template uniform_int

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, class IntType = int>
+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;
+};
+
+ +

Description

+ +The distribution function uniform_int models a +uniform random number +generator. On each invocation, it returns a random integer +value uniformly distributed in the set of integer numbers +{min, min+1, min+2, ..., max}. +

+ +The template parameter IntType shall denote an +integer-like value type. + +

Members

+ +
uniform_int(base_type & rng, IntType min, IntType max)
+ +Effects: Constructs a uniform_int functor +with the uniform random number generator rng as the +underlying source of random numbers. min and +max are the lower and upper bounds of the output range, +respectively. + +

+Note: Invocations of operator() 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. + +

Class template uniform_01

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, class RealType = double>
+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;
+};
+
+ +

Description

+ +The distribution function uniform_01 models a +uniform random number +generator. On each invocation, it returns a random floating-point +value uniformly distributed in the range [0..1). + +The value is computed using +std::numeric_limits<RealType>::digits random binary +digits, i.e. the mantissa of the floating-point value is completely +filled with random bits. [Note: Should this be configurable?] + +

+The template parameter RealType 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 +rng.max()-rng.min()+1. +

+base_type::result_type must be a number-like value type, +it must support static_cast<> to +RealType and binary operator -. + +

+ +Note: 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) boost::bigfloat class with random bits +efficiently. It's probably time for a traits class. + +

Members

+ +
explicit uniform_01(base_type & rng)
+ +Effects: Constructs a uniform_01 functor +with the given uniform random number generator as the underlying +source of random numbers. + + +

Class template uniform_real

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, class RealType = double>
+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;
+};
+
+ +

Description

+ +The distribution function uniform_real models a +uniform random number +generator. On each invocation, it returns a random floating-point +value uniformly distributed in the range [min..max). The value is +computed using +std::numeric_limits<RealType>::digits random binary +digits, i.e. the mantissa of the floating-point value is completely +filled with random bits. +

+ +Note: The current implementation is buggy, because it may not +fill all of the mantissa with random bits. + + +

Members

+ +
explicit uniform_real(base_type & rng, RealType min, RealType max)
+ +Effects: Constructs a uniform_real +functor. rng specifies the uniform random number +generator to be used as the underlying source of random numbers, +min and max are the lower and upper bounds of +the output range. + + +

Class template +bernoulli_distribution

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator>
+class bernoulli_distribution
+{
+public:
+  typedef UniformRandomNumberGenerator base_type;
+  typedef bool result_type;
+  bernoulli_distribution(base_type & rng, double q);
+  result_type operator()();
+};
+
+ +

Description

+ +Instantiations of class template bernoulli_distribution +model a number +generator. It transforms a uniform distribution into a Bernoulli +one. + +

Members

+ +
bernoulli_distribution(base_type & rng, double q)
+ +Effects: Constructs a +bernoulli_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. q is the parameter for the distribution. + +
result_type operator()()
+Returns: 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. + + +

Class template +geometric_distribution

+ +

Synopsis

+
+template<class UniformRandomNumberGenerator, class IntType = int>
+class geometric_distribution
+{
+public:
+  typedef UniformRandomNumberGenerator base_type;
+  typedef IntType result_type;
+  geometric_distribution(base_type& rng, double q);
+  result_type operator()();
+};
+
+ + +

Description

+ +Instantiations of class template geometric_distribution +model a number +generator. It transforms a uniform distribution into a geometric +one. + + +

Members

+ +
geometric_distribution(base_type& rng, const result_type& q)
+ +Effects: Constructs a +geometric_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. q is the parameter for the distribution. +

+ +

result_type operator()()
+ +Returns: A random integer value i >= 1 with +p(i) = (1-q) * qi-1. 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. + + +

Class template +triangle_distribution

+ +

Synopsis

+
+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);
+  result_type operator()();
+};
+
+ + +

Description

+ +Instantiations of class template triangle_distribution +model a number +generator. It transforms a uniform distribution into a triangle +one. + + +

Members

+ +
triangle_distribution(base_type& rng, result_type a, result_type b, result_type c)
+ +Effects: Constructs a +triangle_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. a, b, c are the parameters for the distribution. +

+ +

result_type operator()()
+ +Returns: A random floating-point value x +where a <= x <= c; x has a triangle +distribution, where b is the most probable value for +x. + + +

Class template +exponential_distribution

+ +

Synopsis

+
+template<class UniformRandomNumberGenerator, class RealType = double>
+class exponential_distribution
+{
+public:
+  typedef UniformRandomNumberGenerator base_type;
+  typedef RealType result_type;
+  exponential_distribution(base_type& rng, const result_type& lambda);
+  result_type operator()();
+};
+
+ +

Description

+ +Instantiations of class template exponential_distribution +model a number +generator. It transforms a uniform distribution into an +exponential one. + +

Members

+ +
exponential_distribution(base_type& rng, const result_type& lambda)
+ +Effects: Constructs an +exponential_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. lambda is the parameter for the distribution. +

+ +

result_type operator()()
+ +Returns: A random floating-point value x +> 0 with p(x) = lambda * exp(-lambda * +x). + + +

Class template +normal_distribution

+ +

Synopsis

+ +
+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);
+  result_type operator()();
+};
+
+ +

Description

+ +Instantiations of class template normal_distribution +model a number +generator. It transforms a uniform distribution into a +normal (Gaussian) one. + +

Members

+ +
normal_distribution(base_type& rng, const result_type& mean = 0,
+		     const result_type& sigma = 1)
+ +Effects: Constructs a +normal_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. mean and sigma are the parameters for +the distribution. +

+ +

result_type operator()()
+ +Returns: A random floating-point value x +with p(x) = 1/sqrt(2*pi*sigma) * exp(- (x-mean)2 / +(2*sigma2) ). + + +

Class template +lognormal_distribution

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, class RealType = double>
+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()();
+};
+
+ +

Description

+ +Instantiations of class template lognormal_distribution +model a number +generator. It transforms a uniform distribution into a +lognormal one. + +

Members

+ +
lognormal_distribution(base_type& rng, const result_type& mean = 0,
+	   	       const result_type& sigma = 1)
+ +Effects: Constructs a +normal_distribution functor with the uniform random +number generator rng as the underlying source of random +numbers. mean and sigma are the parameters for +the distribution. +

+ +

result_type operator()()
+ +Returns: A random floating-point value x +with p(x) = 1/(x * sigma * sqrt(2*pi*sigma)) * exp(- +(log(x)-mean)2 / (2*sigma2) ) for x > 0. + + +

Class template +uniform_on_sphere

+ +

Synopsis

+ +
+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);
+  const result_type & operator()();
+};
+
+ +

Description

+ +Instantiations of class template uniform_on_sphere 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 Cont template parameter must be +a STL-like container type with begin and end +operations returning non-const ForwardIterators of type +Cont::iterator. + +

Members

+ +
explicit uniform_on_sphere(base_type & rng, int dim = 2)
+ +Effects: Constructs a uniform_on_sphere +functor with the uniform random number generator rng as +the underlying source of random numbers. dim is the +dimension of the sphere. +

+ +

result_type operator()()
+ +Returns: A position on the unit sphere of +dim dimensions in cartesian coordinates. The positions +are uniformly distributed on the unit sphere. +

+Complexity: Proportional to the number of dimensions. + +

+

+Jens Maurer, 2000-06-14 + + + diff --git a/random-generators.html b/random-generators.html new file mode 100644 index 0000000..88cac3b --- /dev/null +++ b/random-generators.html @@ -0,0 +1,785 @@ + + + + + + +Boost Random Number Library Generators + + + + +

Random Number Library Generators

+ + + +

Introduction

+ +This library provides several pseudo-random number generators. The +quality of a pseudo-random number generator crucially depends on both +the algorithm and its parameters. This library implements the +algorithms as class templates with template value parameters, hidden +in namespace boost::random. Any particular choice of +parameters is represented as the appropriately specializing +typedef in namespace boost. + +

+ +The following table gives an overview of some characteristics of the +generators. The cycle length is a rough estimate of the quality of +the generator; the approximate relative speed is a performance +measure, higher numbers mean faster random number generation. + +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
generatorlength of cyclememory requirementsapproximate relative speedcomment
minstd_rand2**31-2sizeof(int32_t)0.4-
rand482**48-1sizeof(uint64_t)0.8-
lrand48 (C library)2**48-1-0.2global state
ecuyer1988approx. 2**612*sizeof(int32_t)0.2-
kreutzer1986?1368*sizeof(uint32_t)0.6-
hellekalek19952**31-1sizeof(int32_t)0.03good uniform distribution in several dimensions
mt11213b +2**11213-1352*sizeof(uint32_t)1good uniform distribution in up to 350 dimensions
mt19937 +2**19937-1625*sizeof(uint32_t)1good uniform distribution in up to 623 dimensions
+ +

+As observable from the table, there is generally a +quality/performance/memory trade-off to be decided upon when choosing +a random-number generator. The multitude of generators provided in +this library allows the application programmer to optimize the +trade-off with regard to his application domain. Additionally, +employing several fundamentally different random number generators for +a given application of Monte Carlo simulation will improve the +confidence in the results. + +

+Note: These random number generators are not indended for use +in applications where non-deterministic random numbers are required. +See nondet_random.html for a choice +of (hopefully) non-deterministic random number generators. + +

+In this description, I have refrained from documenting those members +in detail which are already defined in the +concept documentation. + + +

Synopsis of the generators in header +<boost/random.hpp>

+ +
+namespace boost {
+  namespace random {
+    template<class IntType, IntType m>
+    class const_mod;
+    template<class IntType, IntType a, IntType c, IntType m, IntType val>
+    class linear_congruential;
+  }
+  class rand48;
+  typedef random::linear_congruential< /* ... */ > minstd_rand0;
+  typedef random::linear_congruential< /* ... */ > minstd_rand;
+
+  namespace random {
+    template<class DataType, int n, int m, int r, DataType a, int u,
+        int s, DataType b, int t, DataType c, int l, IntType val>
+  class mersenne_twister;
+  }
+  typedef random::mersenne_twister< /* ... */ > mt11213b;
+  typedef random::mersenne_twister< /* ... */ > mt19937;
+} // namespace boost
+
+ + +

Class template +random::const_mod

+ +

Synopsis

+ +
+template<class IntType, IntType m>
+class random::const_mod
+{
+public:
+  template<IntType c>
+  static IntType add(IntType x);
+
+  template<IntType a>
+  static IntType mult(IntType x);
+
+  template<IntType a, IntType c>
+  static IntType mult_add(IntType x);
+
+  static IntType invert(IntType x);
+private:
+  const_mod();         // don't instantiate
+};
+
+ +

Description

+ +Class template const_mod provides functions performing +modular arithmetic, carefully avoiding overflows. All member +functions are static; there shall be no objects of type +const_mod<>. +

+ +The template parameter IntType shall denote an integral +type, m is the modulus. +

+ +Note: For modulo multiplications with large m, a trick allows +fast computation under certain conditions, see +

+"A more portable FORTRAN random number generator", Linus Schrage, ACM +Transactions on Mathematical Software, Vol. 5, No. 2, June 1979, pp. 132-138 +
+ + +

Member functions

+ +
template<IntType c> static IntType add(IntType x)
+ +Returns: (x+c) mod m + +
template<IntType a> static IntType mult(IntType x)
+ +Returns: (a*x) mod m + +
template<IntType a, IntType c> static IntType
+mult_add(IntType x)
+ +Returns: (a*x+c) mod m + +
static IntType invert(IntType x)
+ +Returns: i so that (a*i) mod m == 1 +
+Precondition: m is prime + + +

Class template +random::linear_congruential

+ +

Synopsis

+ +
+template<class IntType, IntType a, IntType c, IntType m>
+class linear_congruential
+{
+public:
+  typedef IntType result_type;
+  static const IntType multiplier = a;
+  static const IntType increment = c;
+  static const IntType modulus = m;
+  static const bool has_fixed_range = true;
+  static const result_type min_value;
+  static const result_type max_value;
+  explicit linear_congruential_fixed(IntType x0 = 1);
+  // compiler-generated copy constructor and assignment operator are fine
+  void seed(IntType x0);
+  IntType operator()();
+};
+
+typedef random::linear_congruential<long, 16807L, 0, 2147483647L,
+     1043618065L> minstd_rand0;
+typedef random::linear_congruential<long, 48271L, 0, 2147483647L,
+     399268537L> minstd_rand;
+
+ +

Description

+ +Instantiations of class template linear_congruential +model a pseudo-random number +generator. Linear congruential pseudo-random number generators +are described in: +
+"Mathematical methods in large-scale computing units", D. H. Lehmer, +Proc. 2nd Symposium on Large-Scale Digital Calculating Machines, +Harvard University Press, 1951, pp. 141-146 +
+ +Let x(n) denote the sequence of numbers returned by +some pseudo-random number generator. Then for the linear congruential +generator, x(n+1) := (a * x(n) + c) mod m. Parameters for the +generator are x(0), a, c, m. + +The template parameter IntType shall denote an +integral type. It must be large enough to hold values a, c, and m. +The template parameters a and c must be smaller than m. + +

+ +Note: The quality of the generator crucially depends on the +choice of the parameters. User code should use one of the sensibly +parameterized generators such as minstd_rand instead. +
+For each choice of the parameters a, c, m, some distinct type is +defined, so that the static members do not interfere with +regard to the one definition rule. + +

Members

+ +
explicit linear_congruential(IntType x0 = 1)
+ +Effects: Constructs a +linear_congruential generator with x(0) := +x0. + +
void seed(IntType x0)
+ +Effects: Changes the current value x(n) of the +generator to x0. + +

Specializations

+ +The specialization minstd_rand0 was originally suggested +in +
+A pseudo-random number generator for the System/360, P.A. Lewis, +A.S. Goodman, J.M. Miller, IBM Systems Journal, Vol. 8, No. 2, 1969, +pp. 136-146 +
+ +It is examined more closely together with minstd_rand in +
+"Random Number Generators: Good ones are hard to find", Stephen +K. Park and Keith W. Miller, Communications of the ACM, Vol. 31, +No. 10, October 1988, pp. 1192-1201 +
+ + +

Class rand48

+ +

Synopsis

+
+class rand48 
+{
+public:
+  typedef int32_t result_type;
+  static const bool has_fixed_range = true;
+  static const int32_t min_value = 0;
+  static const int32_t max_value = 0x7fffffff;
+  
+  explicit rand48(int32_t x0 = 1);
+  explicit rand48(uint64_t x0);
+  // compiler-generated copy ctor and assignment operator are fine
+  void seed(int32_t x0);
+  void seed(uint64_t x0);
+  int32_t operator()();
+};
+
+ +

Description

+ +Class rand48 models a +pseudo-random number +generator. It uses the linear congruential algorithm with the +parameters a = 0x5DEECE66D, c = 0xB, m = 2**48. It delivers identical +results to the lrand48() function available on some +systems (assuming lcong48 has not been called). +

+It is only available on systems where uint64_t is provided. + +

Constructors

+ +
rand48(int32_t x0)
+ +Effects: Constructs a rand48 generator +with x(0) := (x0 << 16) | 0x330e. + +
rand48(uint64_t x0)
+ +Effects: Constructs a rand48 generator +with x(0) := x0. + +

Seeding

+
void seed(int32_t x0)
+ +Effects: Changes the current value x(n) of the +generator to (x0 << 16) | 0x330e. + +
void seed(uint64_t x0)
+ +Effects: Changes the current value x(n) of the +generator to x0. + + +

Class template +random::additive_combine

+ +

Synopsis

+ +
+template<class MLCG1, class MLCG2, typename MLCG1::result_type val>
+class random::additive_combine
+{
+public:
+  typedef MLCG1 first_base;
+  typedef MLCG2 second_base;
+  typedef typename MLCG1::result_type result_type;
+  static const bool has_fixed_range = true;
+  static const result_type min_value = 1;
+  static const result_type max_value = MLCG1::max_value-1;
+  additive_combine();
+  additive_combine(typename MLCG1::result_type seed1, 
+		   typename MLCG2::result_type seed2);
+  result_type operator()();
+  bool validation(result_type x) const;
+};
+
+typedef random::additive_combine<
+    random::linear_congruential<int32_t, 40014, 0, 2147483563, 0>,
+    random::linear_congruential<int32_t, 40692, 0, 2147483399, 0>,
+  /* unknown */ 0> ecuyer1988;
+
+
+ +

Description

+ +Instatiations of class template additive_combine model a +pseudo-random number +generator. It combines two multiplicative linear congruential +number generators, i.e. those with c = 0. It is described in +
+"Efficient and Portable Combined Random Number Generators", Pierre +L'Ecuyer, Communications of the ACM, Vol. 31, No. 6, June 1988, +pp. 742-749, 774 +
+ +The template parameters MLCG1 and MLCG2 +shall denote two different linear congruential number generators, each +with c = 0. Each invocation returns a random number X(n) := (MLCG1(n) +- MLCG2(n)) mod (m1 - 1), where m1 denotes the modulus of +MLCG1. + +

+The template parameter val is the validation value +checked by validation. + + +

Members

+ +
additive_combine()
+ +Effects: Constructs an additive_combine +generator using the default constructors of the two base generators. + +
additive_combine(typename MLCG1::result_type seed1, 
+ 		 typename MLCG2::result_type seed2)
+ +Effects: Constructs an additive_combine +generator, using seed1 and seed2 as the +constructor argument to the first and second base generator, +respectively. + + +

Specialization

+ +The specialization ecuyer1988 was suggested in the above +paper. + + +

Class template +random::shuffle_output

+ +

Synopsis

+ +
+template<class UniformRandomNumberGenerator, int k, 
+  class IntType = typename UniformRandomNumberGenerator::result_type,
+  typename UniformRandomNumberGenerator::result_type val = 0>
+class random::shuffle_output
+{
+public:
+  typedef UniformRandomNumberGenerator base_type;
+  typedef typename base_type::result_type result_type;
+  static const bool has_fixed_range = false;
+
+  shuffle_output();
+  template<class T> explicit shuffle_output(T seed);
+  explicit shuffle_output(const base_type & rng);
+  template<class T> void seed(T s);
+
+  result_type operator()();
+  result_type min() const;
+  result_type max() const;
+  bool validation(result_type) const;
+};
+
+ +

Description

+ +Instatiations of class template shuffle_output model a +pseudo-random number +generator. It mixes the output of some (usually linear +congruential) uniform random number generator to get better +statistical properties. According to Donald E. Knuth, "The Art of +Computer Programming, Vol. 2", the algorithm is described in +
+"Improving a poor random number generator", Carter Bays and +S.D. Durham, ACM Transactions on Mathematical Software, Vol. 2, 1979, +pp. 59-64. +
+The output of the base generator is buffered in an array of length +k. Every output X(n) has a second role: It gives an index into the +array where X(n+1) will be retrieved. Used array elements are +replaced with fresh output from the base generator. + +

+ +Template parameters are the base generator and the array length k, +which should be around 100. As an implementation detail, the template +parameter IntType shall denote an integer-like type which +is large enough to hold integer numbers of value k * +base_type::max(). The template parameter +val is the validation value checked by +validation. + + +

Members

+ +
shuffle_output()
+ +Effects: Constructs a shuffle_output +generator by invoking the default constructor of the base generator. +

+Complexity: Exactly k+1 invocations of the base +generator. + +

template<class T> explicit shuffle_output(T seed)
+ +Effects: Constructs a shuffle_output +generator by invoking the one-argument constructor of the base +generator with the parameter seed. +

+Complexity: Exactly k+1 invocations of the base +generator. + +

explicit shuffle_output(const base_type & rng)
+ +Precondition: The template argument +UniformRandomNumberGenerator shall denote a +CopyConstructible type. +

+Effects: Constructs a shuffle_output +generator by using a copy of the provided generator. +

+Complexity: Exactly k+1 invocations of the base +generator. + +

template<class T> void seed(T s)
+ +Effects: Invokes the one-argument seed +method of the base generator with the parameter seed and +re-initializes the internal buffer array. +

+Complexity: Exactly k+1 invocations of the base +generator. + + +

Specializations

+ +According to Harry Erwin (private e-mail), the specialization +kreutzer1986 was suggested in: +
+"System Simulation: programming Styles and Languages (International +Computer Science Series)", Wolfgang Kreutzer, Addison-Wesley, December +1986. +
+ + +

Class template +random::inversive_congruential

+ +

Synopsis

+ +
+template<class IntType, IntType a, IntType b, IntType p>
+class random::inversive_congruential
+{
+public:
+  typedef IntType result_type;
+  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;
+  static const result_type multiplier = a;
+  static const result_type increment = b;
+  static const result_type modulus = p;
+  explicit inversive_congruential(IntType y0 = 1);
+  void seed(IntType y0);
+  IntType operator()();
+};
+
+typedef random::inversive_congruential<int32_t, 9102, 2147483647-36884165, 2147483647> hellekalek1995;
+
+ +

Description

+ +Instantiations of class template inversive_congruential model a +pseudo-random number +generator. It uses the inversive congruential algorithm (ICG) +described in +
+"Inversive pseudorandom number generators: concepts, results and +links", Peter Hellekalek, In: "Proceedings of the 1995 Winter +Simulation Conference", C. Alexopoulos, K. Kang, W.R. Lilegdon, and +D. Goldsman (editors), 1995, pp. 255-262. +ftp://random.mat.sbg.ac.at/pub/data/wsc95.ps +
+ +The output sequence is defined by x(n+1) = (a*inv(x(n)) - b) (mod p), +where x(0), a, b, and the prime number p are parameters of the +generator. The expression inv(k) denotes the multiplicative inverse +of k in the field of integer numbers modulo p, with inv(0) := 0. + +

+ +The template parameter IntType shall denote a signed +integral type large enough to hold p; a, b, and p are the parameters +of the generators. +

+Note: The implementation currently uses the Euclidian +Algorithm to compute the multiplicative inverse. Therefore, the +inversive generators are about 10-20 times slower than the others (see +section"performance"). However, the paper +talks of only 3x slowdown, so the Euclidian Algorithm is probably not +optimal for calculating the multiplicative inverse. + + +

Members

+ +
inversive_congruential(IntType y0 = 1)
+ +Effects: Constructs an +inversive_congruential generator with +y0 as the initial state. + +
void seed(IntType y0)
+ +Effects: +Changes the current state to y0. + + +

Specialization

+ +The specialization hellekalek1995 was suggested in the +above paper. + + +

Class template +random::mersenne_twister

+ +

Synopsis

+ +
+template<class DataType, int n, int m, int r, DataType a, int u,
+int s, DataType b, int t, DataType c, int l, IntType val>
+class random::mersenne_twister
+{
+public:
+  typedef DataType result_type;
+  static const bool has_fixed_range = true;
+  static const result_type min_value;
+  static const result_type max_value;
+  mersenne_twister();
+  explicit mersenne_twister(DataType value);
+  template<class Generator> explicit mersenne_twister(Generator & gen);
+  // compiler-generated copy ctor and assignment operator are fine
+  void seed();
+  void seed(DataType value);
+  template<class Generator> void seed(Generator & gen);
+  result_type operator()();
+  bool validation(result_type) const;
+};
+
+typedef mersenne_twister<uint32_t,351,175,19,0xccab8ee7,11,7,0x31b6ab00,15,0xffe50000,17, /* unknown */ 0> mt11213b;
+typedef mersenne_twister<uint32_t,624,397,31,0x9908b0df,11,7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
+
+ +

Description

+ +Instantiations of class template mersenne_twister model a +pseudo-random number +generator. It uses the algorithm described in + +
+"Mersenne Twister: A 623-dimensionally equidistributed uniform +pseudo-random number generator", Makoto Matsumoto and Takuji Nishimura, +ACM Transactions on Modeling and Computer Simulation: Special Issue +on Uniform Random Number Generation, Vol. 8, No. 1, January 1998, +pp. 3-30. +http://www.math.keio.ac.jp/matumoto/emt.html +
+ +Note: The boost variant has been implemented from scratch +and does not derive from or use mt19937.c provided on the above WWW +site. However, it was verified that both produce identical output. +
+The quality of the generator crucially depends on the choice of the +parameters. User code should employ one of the sensibly parameterized +generators such as mt19937 instead. +
+The generator requires considerable amounts of memory for the storage +of its state array. For example, mt11213b requires about +1408 bytes and mt19937 requires about 2496 bytes. + +

Constructors

+ +
mersenne_twister()
+ +Effects: Constructs a mersenne_twister +and calls seed(). + +
explicit mersenne_twister(DataType value)
+ +Effects: Constructs a mersenne_twister +and calls seed(value). + +
template<class Generator> explicit mersenne_twister(Generator & gen)
+ +Effects: Constructs a mersenne_twister +and calls seed(gen). + +

Seeding

+ +
void seed()
+ +Effects: Calls seed(4357). + +
void seed(DataType value)
+ +Effects: Constructs a +linear_congruential<uint32_t, 69069, 0, 0, 0> +generator with the constructor parameter +value and calls seed with it. + +
template<class Generator> void seed(Generator & gen)
+ +Effects: Sets the state of this +mersenne_twister to the values returned by n +invocations of gen. + +

+ +Complexity: Exactly n invocations of gen. + +

Specializations

+ +The specializations mt11213b and mt19937 are +from the paper cited above. + + +

Performance

+ +The test program random_speed.cpp +measures the execution times of the +random.hpp 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. + +

+ + + + + + + + + + + + + +
classtime per invocation [usec]
rand480.096
rand48 run-time configurable0.697
lrand48 glibc 2.1.20.844
minstd_rand0.174
ecuyer19880.445
kreutzer19860.249
hellekalek1995 (inversive)4.895
mt11213b0.165
mt199370.165
mt19937 original0.185
+ +

+The measurement error is estimated at +/- 10 nsec. + +

+

+Jens Maurer, 2000-06-13 + + + diff --git a/random-misc.html b/random-misc.html new file mode 100644 index 0000000..d7d100a --- /dev/null +++ b/random-misc.html @@ -0,0 +1,159 @@ + + + + + + +Boost Random Number Generator Library (Miscellaneous) + + + + +

Random Number Generator Library --- Miscellaneous Decorators

+ + + +

Introduction

+ +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. + + +

Synopsis of miscellaneous decorators in +header <boost/random.hpp>

+ +
+namespace boost {
+  template<class UniformRandomNumberGenerator, class IntType = long>
+  class random_number_generator;
+  template<class Generator>
+  class generator_iterator;
+} // namespace boost
+
+ + +

Class template +random_number_generator

+ +

Synopsis

+
+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);
+  result_type operator()(argument_type n);
+};
+
+ +

Description

+ +Instantiations of class template random_number_generator +model a RandomNumberGenerator (std:25.2.11 [lib.alg.random.shuffle]). +On each invocation, it returns a uniformly distributed integer in +the range [0..n). +

+The template parameter IntType shall denote some +integer-like value type. +

+ +Note: I consider it unfortunate that the C++ Standard uses +the name RandomNumberGenerator for something rather specific. + +

Members

+ +
random_number_generator(base_type & rng)
+ +Effects: Constructs a +random_number_generator functor with the given uniform +random number generator as the underlying source of random numbers. + +
result_type operator()(argument_type n)
+ +Returns: The value of +uniform_int<base_type>(rng, 0, n-1)(). + + +

Class template +generator_iterator

+ +

Synopsis

+
+template<class Generator>
+class generator_iterator
+  : equality_comparable<generator_iterator<Generator> >,
+  incrementable<generator_iterator<Generator> >,
+  dereferenceable<generator_iterator<Generator>,
+      typename Generator::result_type>
+{
+public:
+  typedef typename Generator::result_type value_type;
+  typedef std::ptrdiff_t difference_type;
+  typedef const typename Generator::result_type * pointer;
+  typedef const typename Generator::result_type & reference;
+  typedef std::input_iterator_tag iterator_category;
+
+  explicit generator_iterator(Generator & g);
+  generator_iterator& operator++();
+  reference operator*() const;
+  friend bool operator==(const generator_iterator<Generator>& x, 
+			 const generator_iterator<Generator>& y);
+};
+
+ +

Description

+ +Instantiations of class template generator_iterator +satisfy the input iterator requirements (std:24.1.1 +[lib.input.iterators]). It allows iterator-like access to a +generator, e.g. a NumberGenerator. Note that all distribution +functions now satisfy the input iterator requirements as-is. However, +the base generators do not. + +

Members

+ +
explicit generator_iterator(Generator & g)
+ +Effects: Constructs a generator_iterator +with g as the underlying generator. Invokes the +underlying generator functor. + +
generator_iterator& operator++()
+ +Effects: Invokes the underlying generator functor. +

+Returns: *this + +

reference operator*() const
+ +Returns: The value of the last invocation of the +underlying generator functor. + +

Overloaded global operators

+ +
bool operator==(const generator_iterator<Generator>& x, 
+	         const generator_iterator<Generator>& y)
+ +Returns: true iff the x and +y have been initialized with a reference to the same +generator functor and *x == *y. + + +

+

+Jens Maurer, 2000-03-06 + + + diff --git a/random_demo.cpp b/random_demo.cpp new file mode 100644 index 0000000..d807ebc --- /dev/null +++ b/random_demo.cpp @@ -0,0 +1,80 @@ +/* boost random_demo.cpp profane demo + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include +#include // std::time +#include + + +// try boost::mt19937 or boost::ecuyer1988 instead of boost::minstd_rand +typedef boost::minstd_rand base_generator_type; + +// This is a reproducible simulation experiment. +void experiment(base_generator_type & generator) +{ + boost::uniform_smallint die(generator, 1, 6); + + // you can use a STL Iterator interface + for(int i = 0; i < 10; i++) + std::cout << *die++ << " "; + std::cout << '\n'; +} + +int main() +{ + // initialize by reproducible seed + base_generator_type generator(42); + + std::cout << "10 samples of a uniform distribution in [0..1):\n"; + boost::uniform_01 uni(generator); + std::cout.setf(std::ios::fixed); + + // you can also use a STL Generator interface + for(int i = 0; i < 10; i++) + std::cout << uni() << '\n'; + + // change seed to something else + // Note: this is not the preferred way of hacking around missing std:: + generator.seed( +#ifndef BOOST_NO_STDC_NAMESPACE + std:: +#endif + time(0)); + + std::cout << "\nexperiment: roll a die 10 times:\n"; + base_generator_type saved_generator = generator; + 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 + { + // save the generator state for future use, + // can be read again at any time via 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; +} diff --git a/random_device.cpp b/random_device.cpp new file mode 100644 index 0000000..6b7e745 --- /dev/null +++ b/random_device.cpp @@ -0,0 +1,108 @@ +/* boost random_device.cpp implementation + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include + + +#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 +#include +#include // open +#include // read, close +#endif + +#include // errno +#include // strerror +#include // 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(&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::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(); +} diff --git a/random_speed.cpp b/random_speed.cpp new file mode 100644 index 0000000..68b9608 --- /dev/null +++ b/random_speed.cpp @@ -0,0 +1,217 @@ +/* boost random_speed.cpp performance measurements + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include +#include +#include +#include + +/* + * 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() and f() are both compiled to call f(). + * 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 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; +}; + +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 +static 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 +static 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 & tmp = rng(); + (void) tmp[0]; + } + show_elapsed(t.elapsed(), iter, name); +} + +template +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 +void distrib(int iter, const std::string & name, const Gen &) +{ + Gen gen; + + boost::uniform_smallint usmallint(gen, -2, 4); + timing(usmallint, iter, name + " uniform_smallint", 0); + + boost::uniform_int uint(gen, -2, 4); + timing(uint, iter, name + " uniform_int", 0); + + boost::uniform_01 uni(gen); + timing(uni, iter, name + " uniform_01", 0.0); + + boost::uniform_real ur(gen, -5.3, 4.8); + timing(ur, iter, name + " uniform_real", 0.0); + + boost::geometric_distribution geo(gen, 0.5); + timing(geo, iter, name + " geometric", 0); + + boost::triangle_distribution tria(gen, 1, 2, 7); + timing(tria, iter, name + " triangle", 0.0); + + boost::exponential_distribution ex(gen, 3); + timing(ex, iter, name + " exponential", 0.0); + + boost::normal_distribution no2(gen); + timing(no2, iter, name + " normal polar", 0.0); + + boost::lognormal_distribution lnorm(gen, 1, 1); + timing(lnorm, iter, name + " lognormal", 0.0); + + boost::uniform_on_sphere 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]); + +#ifdef BOOST_STDINT_H_HAS_UINT64_T + run(iter, "rand48", boost::rand48); + linear_congruential lcg48(1ULL<<16 | 0x330e, + 0x5DEECE66DULL, 0xB, + (1ULL<<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(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()); +}; diff --git a/random_test.cpp b/random_test.cpp new file mode 100644 index 0000000..d71e826 --- /dev/null +++ b/random_test.cpp @@ -0,0 +1,309 @@ +/* boost random_test.cpp various tests + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include +#include +#include +#include + +/* + * General portability note: + * MSVC mis-compiles explicit function template instantiations. + * For example, f() and f() are both compiled to call f(). + * 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 +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(); + bool result = const_cast(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 << " "; + if(!result) + std::cout << "validation failed or not possible.\n"; + else + std::cout << "ok\n"; +} + +void validate_all() +{ + using namespace boost; +#ifdef BOOST_STDINT_H_HAS_UINT64_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 +void instantiate_iterator_interface(Generator & gen) +{ + ++gen; + gen++; + typename Generator::result_type res = *gen; + typename Generator::pointer p = &res; + (void) p; + typename Generator::reference ref = res; + (void) ref; + typename Generator::difference_type diff = 0; + (void) diff; + std::input_iterator_tag it(typename Generator::iterator_category()); + (void) it; + + assert(res == *gen++); + assert(gen == gen); + Generator gen2 = gen; + assert(gen == gen2); // must be equal to a copy +#if 0 + ++gen2; + // this is only correct for elementary generators, others have ref members + assert(gen != gen2); +#endif +} + +template +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 + assert(urng == urng2); // operator== + urng(); + urng2 = urng; // assignment + assert(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 + assert(urng == urng2); +#endif // BOOST_MSVC +#endif // BOOST_NO_OPERATORS_IN_NAMESPACE + + // instantiate various distributions with this URNG + boost::uniform_smallint unismall(urng, 0, 11); + instantiate_iterator_interface(unismall); + unismall(); + boost::uniform_int uni_int(urng, -200, 20000); + instantiate_iterator_interface(uni_int); + uni_int(); + boost::uniform_real uni_real(urng, 0, 2.1); + instantiate_iterator_interface(uni_real); + uni_real(); + + boost::bernoulli_distribution ber(urng, 0.2); + instantiate_iterator_interface(ber); + ber(); + boost::geometric_distribution geo(urng, 0.8); + instantiate_iterator_interface(geo); + geo(); + boost::triangle_distribution tria(urng, 1, 1.5, 7); + instantiate_iterator_interface(tria); + tria(); + boost::exponential_distribution ex(urng, 5); + instantiate_iterator_interface(ex); + ex(); + boost::normal_distribution norm(urng); + instantiate_iterator_interface(norm); + norm(); + boost::lognormal_distribution lnorm(urng, 1, 1); + instantiate_iterator_interface(lnorm); + lnorm(); + boost::uniform_on_sphere usph(urng, 2); + instantiate_iterator_interface(usph); + usph(); +} + +void instantiate_all() +{ + using namespace boost; + +#ifdef BOOST_STDINT_H_HAS_UINT64_T + instantiate_urng("rand48", rand48, 0); + rand48 rnd(5); + rand48 rnd2(uint64_t(0x80000000) * 42); + rnd.seed(17); + rnd2.seed(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(17u); // needs to be an exact type match for MSVC + mt.seed(42u); + mt19937 mt2(mstd); + mt2.seed(mstd); + + + random_number_generator std_rng(mt2); + (void) std_rng(10); +} + +/* + * A few equidistribution tests + */ + +// yet to come... + +template +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 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]; + int avg = sum/range; + for(int i = 0; i < range; i++) { + if(abs(bucket[i] - avg) > 2*avg/sqrt(iter)) { // 95% confidence interval + std::cout << " ... has bucket[" << i << "] = " << bucket[i] + << " (distance " << (bucket[i] - avg) << ")" + << std::endl; + } + } +} + +template +void test_uniform_int(Generator & gen) +{ + typedef boost::uniform_int int_gen; + + // large range => small range (modulo case) + int_gen uint12(gen,1,2); + assert(uint12.min() == 1); + assert(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 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 intintint_gen; + + // 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); + + // 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; \ +template class boost::uniform_int; \ +template class boost::uniform_real; \ +template class boost::bernoulli_distribution; \ +template class boost::geometric_distribution; \ +template class boost::triangle_distribution; \ +template class boost::exponential_distribution; \ +template class boost::normal_distribution; \ +template class boost::uniform_on_sphere; \ +template class boost::lognormal_distribution; + +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 main() +{ + instantiate_all(); + validate_all(); + boost::mt19937 mt; + test_uniform_int(mt); + // Some compilers don't pay attention to std:3.6.1/5 and issue a + // warning here if "return 0;" is omitted. + return 0; +} diff --git a/statistic_tests.cpp b/statistic_tests.cpp new file mode 100644 index 0000000..170e018 --- /dev/null +++ b/statistic_tests.cpp @@ -0,0 +1,662 @@ +/* statistic_tests.cpp file + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "statistic_tests.hpp" +#include "integrate.hpp" + + +namespace boost { +namespace random { + +// Wikramaratna 1989 ACORN +template +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 + explicit additive_congruential(InputIterator start) { seed(start); } + template + 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 lagged_fibonacci +{ +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(IntType start) { seed(start); } + template + explicit lagged_fibonacci(Generator & gen) { seed(gen); } + void seed(IntType start) + { + linear_congruential init; + seed(init); + } + template + 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 LCG_Af2; + typedef boost::random::linear_congruential LCG_Die1; + typedef boost::random::linear_congruential LCG_Fis; + typedef boost::random::linear_congruential LCG_FM; + typedef boost::random::linear_congruential LCG_Hae; + typedef boost::random::linear_congruential LCG_VAX; + typedef boost::random::inversive_congruential NLG_Inv1; + typedef boost::random::inversive_congruential NLG_Inv2; + typedef boost::random::inversive_congruential NLG_Inv4; + typedef boost::random::inversive_congruential NLG_Inv5; + typedef boost::random::additive_congruential MRG_Acorn7; + typedef boost::random::lagged_fibonacci MRG_Fib2; + + template + 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 + void validate(T value, const std::string & name) + { + Gen gen(1234567); + check_validation(gen, value, name); + } + + void validate_all() + { + validate(183269031u, "LCG_Af2"); + validate(522319944u, "LCG_Die1"); + validate(-2065162233u, "LCG_Fis"); + validate(581815473u, "LCG_FM"); + validate(28931709, "LCG_Hae"); + validate(1508154087u, "LCG_VAX"); + validate(6666884, "NLG_Inv2"); + validate(1521640076, "NLG_Inv4"); + validate(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(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); +} + +#ifdef _CXXRTCF_H__ +namespace std { + using _CS_swamp::lgamma; +} +#endif + +class chi_square_density : public std::unary_function +{ +public: + chi_square_density(int freedom) + : _exponent( static_cast(freedom)/2-1 ), + _factor(1/(std::pow(2, _exponent+1) * std::exp(std::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 +{ +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 +{ +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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "equidistribution: " << std::flush; + equidistribution_experiment equi(classes); + uniform_smallint 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(std::sqrt(classes)); + assert(root * root == classes); + uniform_smallint 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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "KS: " << std::flush; + generator_reference_t 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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "runs: up: " << std::flush; + runs_experiment r_up(classes); + generator_reference_t 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 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 + 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 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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "poker: " << std::flush; + poker_experiment poker(8, classes); + uniform_smallint 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 + 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 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(classes)-1), + high_classes) + { } + + template + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "permutation: " << std::flush; + permutation_experiment perm(classes); + + generator_reference_t 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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "maximum-of-t: " << std::flush; + maximum_experiment 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 + void run(RNG & rng, int n1, int n2) + { + using namespace boost; + std::cout << "birthday spacing: " << std::flush; + uniform_int 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 + 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 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("minstd_rand"); + env.run_test("mt19937"); + env.run_test("LCG_Af2"); + env.run_test("LCG_Die1"); + env.run_test("LCG_Fis"); + env.run_test("LCG_FM"); + env.run_test("LCG_Hae"); + env.run_test("LCG_VAX"); + env.run_test("NLG_Inv1"); + env.run_test("NLG_Inv2"); + env.run_test("NLG_Inv4"); + env.run_test("NLG_Inv5"); +} diff --git a/statistic_tests.hpp b/statistic_tests.hpp new file mode 100644 index 0000000..f66c6d3 --- /dev/null +++ b/statistic_tests.hpp @@ -0,0 +1,628 @@ +/* statistic_tests.hpp header file + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free 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 +#include +#include +#include +#include +#include + + +template +inline T fac(int k) +{ + T result = 1; + for(T i = 2; i <= k; ++i) + result *= i; + return result; +} + +template +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(n-k); +} + +template +T stirling2(int n, int m) +{ + T sum = 0; + for(int k = 0; k <= m; ++k) + sum += binomial(m, k) * std::pow(k, n) * ( (m-k)%2 == 0 ? 1 : -1); + return sum / fac(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 + void run(NumberGenerator f, Counter & count, int n) const + { + assert(f.min() == 0 && + static_cast(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 + 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 + 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::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 + 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 limits_type; + limits_type limit; +}; + +// runs-up/runs-down experiment +template +class runs_experiment : public experiment_base +{ +public: + explicit runs_experiment(unsigned int classes) : experiment_base(classes) { } + + template + 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(classes()); + else + return static_cast(r+1)/fac(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 + 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(alpha * range); + result_type high = static_cast(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(r)); + else + return p * std::pow(1-p, static_cast(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 + void run(UniformRandomNumberGenerator f, Counter & count, int n) const + { + typedef typename UniformRandomNumberGenerator::result_type result_type; + assert(std::numeric_limits::is_integer); + assert(f.min() == 0); + assert(f.max() == static_cast(range-1)); + std::vector 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(classes())) * + stirling2(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 + void run(UniformRandomNumberGenerator f, Counter & count, int n) const + { + typedef typename UniformRandomNumberGenerator::result_type result_type; + assert(std::numeric_limits::is_integer); + assert(f.min() == 0); + assert(f.max() == static_cast(d-1)); + std::vector 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(d)/std::pow(d, static_cast(d+classes()-2))* + stirling2(d+classes()-2, d); + else + return fac(d)/std::pow(d, static_cast(d+r)) * + stirling2(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(t)), t(t) + { + assert(t > 1); + } + + template + void run(UniformRandomNumberGenerator f, Counter & count, int n) const + { + typedef typename UniformRandomNumberGenerator::result_type result_type; + std::vector 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::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 + void run(UniformRandomNumberGenerator f, Counter & count, int n_total) const + { + typedef typename UniformRandomNumberGenerator::result_type result_type; + assert(std::numeric_limits::is_integer); + assert(f.min() == 0); + assert(f.max() == static_cast(m-1)); + + for(int j = 0; j < n_total; j++) { + std::vector v(n); + std::generate_n(v.begin(), n, f); + std::sort(v.begin(), v.end()); + std::vector 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 +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 +{ +public: + kolmogorov_smirnov_probability(int n) + : approx(n > 50), n(n), sqrt_n(std::sqrt(n)) + { + if(!approx) + n_n = std::pow(static_cast(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(std::ceil(t)); k <= n; k++) + sum += binomial(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 + double run(NumberGenerator gen, Distribution distrib) const + { + const int m = n; + typedef std::vector saved_temp; + saved_temp a(m,1.0), b(m,0); + std::vector c(m,0); + for(int i = 0; i < n; ++i) { + double val = gen(); + double y = distrib(val); + int k = static_cast(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(n)); + j += c[k]; + kplus = std::max(kplus, j/static_cast(n) - b[k]); + } + } + kplus *= std::sqrt(n); + kminus *= std::sqrt(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 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(&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 f; + int t; + }; + base_type & f; + kolmogorov_experiment ke; + int t; +}; + +// compute a chi-square value for the distribution approximation error +template +typename UnaryFunction::result_type +chi_square_value(ForwardIterator first, ForwardIterator last, + UnaryFunction probability) +{ + typedef std::iterator_traits 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 generic_counter +{ +public: + explicit generic_counter(unsigned int classes) : container(classes, 0) { } + void operator()(int i) + { + assert(i >= 0); + assert(static_cast(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 +double run_experiment(const Experiment & experiment, Generator gen, int n) +{ + generic_counter > 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_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 +experiment_generator_t +experiment_generator(const Experiment & e, Generator & gen, int n) +{ + return experiment_generator_t(e, gen, n); +} + + +template +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 +ks_experiment_generator_t +ks_experiment_generator(const Experiment & e, Generator & gen, + const Distribution & distrib) +{ + return ks_experiment_generator_t + (e, gen, distrib); +} + + +#endif /* STATISTIC_TESTS_HPP */ +