thijn/SqMod
1
0
Fork 0
mirror of https://github.com/VCMP-SqMod/SqMod.git synced 2025-07-01 06:27:11 +02:00
Code Releases Activity

Dumped the old implementation. Started with a more simple approach.

Browse Source
This commit is contained in:
Sandu Liviu Catalin
2016-02-21 00:25:00 +02:00
parent 96ded94026
commit 06e598acfb
293 changed files with 37439 additions and 92564 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
include/RandomLib/Config.h Normal file
View File

@ -0,0 +1,13 @@
#define RANDOMLIB_VERSION_STRING "Unconfigured"
#define RANDOMLIB_VERSION_MAJOR -1
#define RANDOMLIB_VERSION_MINOR -1
#define RANDOMLIB_VERSION_PATCH -1
// Define HAVE_SSE2 to be 1 if Intel/AMD CPU with SSE2 support
/* #undef HAVE_SSE2 */
// Define HAVE_ALTIVEC to be 1 if Power PC CPU with AltiVec support
/* #undef HAVE_ALTIVEC */
// Undefine HAVE_LONG_DOUBLE if this type is unknown to the compiler
#define HAVE_LONG_DOUBLE 1

432
include/RandomLib/DiscreteNormal.hpp Normal file
View File

@ -0,0 +1,432 @@
/**
* \file DiscreteNormal.hpp
* \brief Header for DiscreteNormal
*
* Sample exactly from the discrete normal distribution.
*
* Copyright (c) Charles Karney (2013) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_DISCRETENORMAL_HPP)
#define RANDOMLIB_DISCRETENORMAL_HPP 1
#include <vector>
#include <limits>
namespace RandomLib {
/**
* \brief The discrete normal distribution.
*
* Sample integers \e i with probability proportional to
* \f[
* \exp\biggl[-\frac12\biggl(\frac{i-\mu}{\sigma}\biggr)^2\biggr],
* \f]
* where &sigma; and &mu; are given as rationals (the ratio of two integers).
* The sampling is exact (provided that the random generator is ideal). For
* example
* \code
#include <iostream>
#include <RandomLib/Random.hpp>
#include <RandomLib/DiscreteNormal.hpp>
int main() {
RandomLib::Random r; // Create r
r.Reseed(); // and give it a unique seed
int sigma_num = 7, sigma_den = 1, mu_num = 1, mu_den = 3;
RandomLib::DiscreteNormal<int> d(sigma_num, sigma_den,
mu_num, mu_den);
for (int i = 0; i < 100; ++i)
std::cout << d(r) << "\n";
}
\endcode
* prints out 100 samples with &sigma; = 7 and &mu; = 1/3.
*
* The algorithm is much the same as for ExactNormal; for details see
* - C. F. F. Karney, <i>Sampling exactly from the normal distribution</i>,
* http://arxiv.org/abs/1303.6257 (Mar. 2013).
* .
* That algorithm samples the integer part of the result \e k, samples \e x
* in [0,1], and (unless rejected) returns <i>s</i>(\e k + \e x), where \e s
* = &plusmn;1. For the discrete case, we sample \e x in [0,1) such that
* \f[
* s(k + x) = (i - \mu)/\sigma,
* \f]
* or
* \f[
* x = s(i - \mu)/\sigma - k
* \f]
* The value of \e i which results in the smallest \e x &ge; 0 is
* \f[
* i_0 = s\lceil k \sigma + s \mu\rceil
* \f]
* so sample
* \f[
* i = i_0 + sj
* \f]
* where \e j is uniformly distributed in [0, &lceil;&sigma;&rceil;). The
* corresponding value of \e x is
* \f[
* \begin{aligned}
* x &= \bigl(si_0 - (k\sigma + s\mu)\bigr)/\sigma + j/\sigma\\
* &= x_0 + j/\sigma,\\
* x_0 &= \bigl(\lceil k \sigma + s \mu\rceil -
* (k \sigma + s \mu)\bigr)/\sigma.
* \end{aligned}
* \f]
* After \e x is sampled in this way, it should be rejected if \e x &ge; 1
* (this is counted with the next larger value of \e k) or if \e x = 0, \e k
* = 0, and \e s = &minus;1 (to avoid double counting the origin). If \e x
* is accepted (in Step 4 of the ExactNormal algorithm), then return \e i.
*
* When &sigma; and &mu; are given as rationals, all the arithmetic outlined
* above can be carried out exactly. The basic rejection techniques used by
* ExactNormal are exact. Thus the result of this discrete form of the
* algorithm is also exact.
*
* RandomLib provides two classes to sample from this distribution:
* - DiscreteNormal which is tuned for speed on a typical general purpose
* computer. This assumes that random samples can be generated relatively
* quickly.
* - DiscreteNormalAlt, which is a prototype for what might be needed on a
* small device used for cryptography which is using a hardware generator
* for obtaining truly random bits. This assumption here is that the
* random bits are relatively expensive to obtain.
* .
* The basic algorithm is the same in the two cases. The main advantages of
* this method are:
* - exact sampling (provided that the source of random numbers is ideal),
* - no need to cut off the tails of the distribution,
* - a short program involving simple integer operations only,
* - no dependence on external libraries (except to generate random bits),
* - no large tables of constants needed,
* - minimal time to set up for a new &sigma; and &mu; (roughly comparable to
* the time it takes to generate one sample),
* - only about 5&ndash;20 times slower than standard routines to sample from
* a normal distribution using plain double-precision arithmetic.
* - DiscreteNormalAlt exhibits ideal scaling for the consumption of random
* bits, namely a constant + log<sub>2</sub>&sigma;, for large &sigma;,
* where the constant is about 31.
* .
* The possible drawbacks of this method are:
* - &sigma; and &mu; are restricted to rational numbers with sufficiently
* small numerators and denominators to avoid overflow (this is unlikely to
* be a severe restriction especially if the template parameter IntType is
* set to <code>long long</code>),
* - the running time is unbounded (but not in any practical sense),
* - the memory consumption is unbounded (but not in any practical sense),
* - the toll, about 30 bits, is considerably worse than that obtained using
* the Knuth-Yao algorithm, for which the toll is no more than 2 (but this
* requires a large table which is expensive to compute and requires a lot
* of memory to store).
*
* This class uses a mutable private vector. So a single DiscreteNormal
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific DiscreteNormal
* object.
*
* Some timing results for IntType = int, &mu; = 0, and 10<sup>8</sup>
* samples (time = time per sample, including setup time, rv = mean number of
* random variables per sample)
* - &sigma; = 10, time = 219 ns, rv = 17.52
* - &sigma; = 32, time = 223 ns, rv = 17.82
* - &sigma; = 1000, time = 225 ns, rv = 17.95
* - &sigma; = 160000, time = 226 ns, rv = 17.95
*
* @tparam IntType the integer type to use (default int).
**********************************************************************/
template<typename IntType = int> class DiscreteNormal {
public:
/**
* Constructor.
*
* @param[in] sigma_num the numerator of &sigma;.
* @param[in] sigma_den the denominator of &sigma; (default 1).
* @param[in] mu_num the numerator of &mu; (default 0).
* @param[in] mu_den the denominator of &mu; (default 1).
*
* The constructor creates a DiscreteNormal objects for sampling with
* specific values of &sigma; and &mu;. This may throw an exception if the
* parameters are such that overflow is possible. Internally &sigma; and
* &mu; are expressed with a common denominator, so it may be possible to
* avoid overflow by picking the fractions of these quantities so that \e
* sigma_den and \e mu_den have many common factors.
**********************************************************************/
DiscreteNormal(IntType sigma_num, IntType sigma_den = 1,
IntType mu_num = 0, IntType mu_den = 1);
/**
* Return a sample.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return discrete normal integer.
**********************************************************************/
template<class Random>
IntType operator()(Random& r) const;
private:
/**
* sigma = _sig / _d, mu = _imu + _mu / _d, _isig = floor(sigma)
**********************************************************************/
IntType _sig, _mu, _d, _isig, _imu;
typedef unsigned short word;
/**
* Holds as much of intermediate uniform deviates as needed.
**********************************************************************/
mutable std::vector<word> _v;
mutable unsigned _m, _l;
/**
* Increment on size of _v.
**********************************************************************/
static const unsigned alloc_incr = 16;
// ceil(n/d) for d > 0
static IntType iceil(IntType n, IntType d);
// abs(n) needed because Visual Studio's std::abs has problems
static IntType iabs(IntType n);
static IntType gcd(IntType u, IntType v);
// After x = LeadingDigit(p), p/_sig = (x + p'/_sig)/b where p and p' are
// in [0, _sig) and b = 1 + max(word).
word LeadingDigit(IntType& p) const;
/**
* Implement outcomes for choosing with prob (\e x + 2\e k) / (2\e k + 2);
* return:
* - 1 (succeed unconditionally) with prob (\e m &minus; 2) / \e m,
* - 0 (succeed with probability x) with prob 1 / \e m,
* - &minus;1 (fail unconditionally) with prob 1 / \e m.
**********************************************************************/
template<class Random> static int Choose(Random& r, int m);
// Compute v' < v. If true set v = v'.
template<class Random> bool less_than(Random& r) const;
// Compute v < (x + p/_sig)/base (updating v)
template<class Random> bool less_than(Random& r, word x, IntType p) const;
// true with prob (x + p/_sig)/base
template<class Random> bool bernoulli(Random& r, word x, IntType p) const;
/**
* Return true with probability exp(&minus;1/2).
**********************************************************************/
template<class Random> bool ExpProbH(Random& r) const;
/**
* Return true with probability exp(&minus;<i>n</i>/2).
**********************************************************************/
template<class Random> bool ExpProb(Random& r, int n) const;
/**
* Return \e n with probability exp(&minus;<i>n</i>/2)
* (1&minus;exp(&minus;1/2)).
**********************************************************************/
template<class Random> int ExpProbN(Random& r) const;
/**
* Algorithm B: true with prob exp(-x * (2*k + x) / (2*k + 2)) where
* x = (x0 + xn / _sig)/b.
**********************************************************************/
template<class Random>
bool B(Random& r, int k, word x0, IntType xn) const;
};
template<typename IntType> DiscreteNormal<IntType>::DiscreteNormal
(IntType sigma_num, IntType sigma_den,
IntType mu_num, IntType mu_den)
: _v(std::vector<word>(alloc_incr)), _m(0), _l(alloc_incr) {
STATIC_ASSERT(std::numeric_limits<IntType>::is_integer,
"DiscreteNormal: invalid integer type IntType");
STATIC_ASSERT(std::numeric_limits<IntType>::is_signed,
"DiscreteNormal: IntType must be a signed type");
STATIC_ASSERT(!std::numeric_limits<word>::is_signed,
"DiscreteNormal: word must be an unsigned type");
STATIC_ASSERT(std::numeric_limits<IntType>::digits + 1 >=
std::numeric_limits<word>::digits,
"DiscreteNormal: IntType must be at least as wide as word");
if (!( sigma_num > 0 && sigma_den > 0 && mu_den > 0 ))
throw RandomErr("DiscreteNormal: need sigma > 0");
_imu = mu_num / mu_den;
if (_imu == (std::numeric_limits<IntType>::min)())
throw RandomErr("DiscreteNormal: abs(mu) too large");
mu_num -= _imu * mu_den;
IntType l;
l = gcd(sigma_num, sigma_den); sigma_num /= l; sigma_den /= l;
l = gcd(mu_num, mu_den); mu_num /= l; mu_den /= l;
_isig = iceil(sigma_num, sigma_den);
l = gcd(sigma_den, mu_den);
_sig = sigma_num * (mu_den / l);
_mu = mu_num * (sigma_den / l);
_d = sigma_den * (mu_den / l);
// The rest of the constructor tests for possible overflow
// Check for overflow in computing member variables
IntType maxint = (std::numeric_limits<IntType>::max)();
if (!( mu_den / l <= maxint / sigma_num &&
mu_num <= maxint / (sigma_den / l) &&
mu_den / l <= maxint / sigma_den ))
throw RandomErr("DiscreteNormal: sigma or mu overflow");
// The probability that k = kmax is about 10^-543.
int kmax = 50;
// Check that max plausible result fits in an IntType, i.e.,
// _isig * (kmax + 1) + abs(_imu) does not lead to overflow.
if (!( kmax + 1 <= maxint / _isig &&
_isig * (kmax + 1) <= maxint - iabs(_imu) ))
throw RandomErr("DiscreteNormal: possible overflow a");
// Check xn0 = _sig * k + s * _mu;
if (!( kmax <= maxint / _sig &&
_sig * kmax <= maxint - iabs(_mu) ))
throw RandomErr("DiscreteNormal: possible overflow b");
// Check for overflow in LeadingDigit
// p << bits, p = _sig - 1, bits = 8
if (!( _sig <= (maxint >> 8) ))
throw RandomErr("DiscreteNormal: overflow in LeadingDigit");
}
template<typename IntType> template<class Random>
IntType DiscreteNormal<IntType>::operator()(Random& r) const {
for (;;) {
int k = ExpProbN(r);
if (!ExpProb(r, k * (k - 1))) continue;
IntType
s = r.Boolean() ? -1 : 1,
xn = _sig * IntType(k) + s * _mu,
i = iceil(xn, _d) + r.template Integer<IntType>(_isig);
xn = i * _d - xn;
if (xn >= _sig || (k == 0 && s < 0 && xn <= 0)) continue;
if (xn > 0) {
word x0 = LeadingDigit(xn); // Find first digit in expansion in words
int h = k + 1; while (h-- && B(r, k, x0, xn));
if (!(h < 0)) continue;
}
return s * i + _imu;
}
}
template<typename IntType>
IntType DiscreteNormal<IntType>::iceil(IntType n, IntType d)
{ IntType k = n / d; return k + (k * d < n ? 1 : 0); }
template<typename IntType> IntType DiscreteNormal<IntType>::iabs(IntType n)
{ return n < 0 ? -n : n; }
template<typename IntType>
IntType DiscreteNormal<IntType>::gcd(IntType u, IntType v) {
// Knuth, TAOCP, vol 2, 4.5.2, Algorithm A
u = iabs(u); v = iabs(v);
while (v > 0) { IntType r = u % v; u = v; v = r; }
return u;
}
template<typename IntType> typename DiscreteNormal<IntType>::word
DiscreteNormal<IntType>::LeadingDigit(IntType& p) const {
static const unsigned bits = 8;
static const unsigned num = std::numeric_limits<word>::digits / bits;
STATIC_ASSERT(bits * num == std::numeric_limits<word>::digits,
"Number of digits in word must be multiple of 8");
word s = 0;
for (unsigned c = num; c--;) {
p <<= bits; s <<= bits;
word d = word(p / _sig);
s += d;
p -= IntType(d) * _sig;
}
return s;
}
template<typename IntType> template<class Random>
int DiscreteNormal<IntType>::Choose(Random& r, int m) {
int k = r.template Integer<int>(m);
return k == 0 ? 0 : (k == 1 ? -1 : 1);
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::less_than(Random& r) const {
for (unsigned j = 0; ; ++j) {
if (j == _m) {
// Need more bits in the old V
if (_l == _m) _v.resize(_l += alloc_incr);
_v[_m++] = r.template Integer<word>();
}
word w = r.template Integer<word>();
if (w > _v[j])
return false; // New V is bigger, so exit
else if (w < _v[j]) {
_v[j] = w; // New V is smaller, update _v
_m = j + 1; // adjusting its size
return true; // and generate the next V
}
// Else w == _v[j] and we need to check the next word
}
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::less_than(Random& r, word x, IntType p) const {
if (_m == 0) _v[_m++] = r.template Integer<word>();
if (_v[0] != x) return _v[0] < x;
for (unsigned j = 1; ; ++j) {
if (p == 0) return false;
if (j == _m) {
if (_l == _m) _v.resize(_l += alloc_incr);
_v[_m++] = r.template Integer<word>();
}
x = LeadingDigit(p);
if (_v[j] != x) return _v[j] < x;
}
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::bernoulli(Random& r, word x, IntType p) const {
word w = r.template Integer<word>();
if (w != x) return w < x;
for (;;) {
if (p == 0) return false;
x = LeadingDigit(p);
w = r.template Integer<word>();
if (w != x) return w < x;
}
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::ExpProbH(Random& r) const {
static const word half = word(1) << (std::numeric_limits<word>::digits - 1);
_m = 0;
if ((_v[_m++] = r.template Integer<word>()) & half) return true;
// Here _v < 1/2. Now loop finding decreasing V. Exit when first
// increasing one is found.
for (unsigned s = 0; ; s ^= 1) { // Parity of loop count
if (!less_than(r)) return s != 0u;
}
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::ExpProb(Random& r, int n) const {
while (n-- > 0) { if (!ExpProbH(r)) return false; }
return true;
}
template<typename IntType> template<class Random>
int DiscreteNormal<IntType>::ExpProbN(Random& r) const {
int n = 0;
while (ExpProbH(r)) ++n;
return n;
}
template<typename IntType> template<class Random>
bool DiscreteNormal<IntType>::B(Random& r, int k, word x0, IntType xn)
const {
int n = 0, h = 2 * k + 2, f;
_m = 0;
for (;; ++n) {
if ( ((f = k ? 0 : Choose(r, h)) < 0) ||
!(n ? less_than(r) : less_than(r, x0, xn)) ||
((f = k ? Choose(r, h) : f) < 0) ||
(f == 0 && !bernoulli(r, x0, xn)) ) break;
}
return (n % 2) == 0;
}
} // namespace RandomLib
#endif // RANDOMLIB_DISCRETENORMAL_HPP

328
include/RandomLib/DiscreteNormalAlt.hpp Normal file
View File

@ -0,0 +1,328 @@
/**
* \file DiscreteNormalAlt.hpp
* \brief Header for DiscreteNormalAlt
*
* Sample exactly from the discrete normal distribution (alternate version).
*
* Copyright (c) Charles Karney (2013) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_DISCRETENORMALALT_HPP)
#define RANDOMLIB_DISCRETENORMALALT_HPP 1
#include <RandomLib/RandomNumber.hpp>
#include <RandomLib/UniformInteger.hpp>
#include <limits>
namespace RandomLib {
/**
* \brief The discrete normal distribution (alternate version).
*
* Sample integers \e i with probability proportional to
* \f[
* \exp\biggl[-\frac12\biggl(\frac{i-\mu}{\sigma}\biggr)^2\biggr],
* \f]
* where &sigma; and &mu; are given as rationals (the ratio of two integers).
* The sampling is exact (provided that the random generator is ideal). The
* results of this class are UniformInteger objects which represent a
* contiguous range which is a power of 2<sup>\e bits</sup>. This can be
* narrowed down to a specific integer as follows
* \code
#include <iostream>
#include <RandomLib/Random.hpp>
#include <RandomLib/UniformInteger.hpp>
#include <RandomLib/DiscreteNormalAlt.hpp>
int main() {
RandomLib::Random r; // Create r
r.Reseed(); // and give it a unique seed
int sigma_num = 7, sigma_den = 1, mu_num = 1, mu_den = 3;
RandomLib::DiscreteNormalAlt<int,1> d(sigma_num, sigma_den,
mu_num, mu_den);
for (int i = 0; i < 100; ++i) {
RandomLib::UniformInteger<int,1> u = d(r);
std::cout << u << " = ";
std::cout << u(r) << "\n";
}
}
\endcode
* prints out 100 samples with &sigma; = 7 and &mu; = 1/3; the samples are
* first given as a range and then <code>u(r)</code> is used to obtain a
* specific integer. In some applications, it might be possible to avoid
* sampling all the additional digits to get to a specific integer. (An
* example would be drawing a sample which satisfies an inequality.) This
* version is slower (by a factor of about 4 for \e bits = 1) than
* DiscreteNormal on a conventional computer, but may be faster when random
* bits are generated by a slow hardware device.
*
* The basic algorithm is the same as for DiscreteNormal. However randomness
* in metered out \e bits random bits at a time. The algorithm uses the
* least random bits (and is slowest!) when \e bits = 1. This rationing of
* random bits also applies to sampling \e j in the expression
* \f[
* x = x_0 + j/\sigma.
* \f]
* Rather than sampling a definite value for \e j, which then might be
* rejected, \e j is sampled using UniformInteger. UniformInteger uses
* Lumbroso's method from sampling an integer uniformly in [0, \e m) using at
* most 2 + log<sub>2</sub>\e m bits on average (for \e bits = 1). However
* when a UniformInteger object is first constructed it samples at most 3
* bits (on average) to obtain a range for \e j which is power of 2. This
* allows the algorithm to achieve ideal scaling in the limit of large
* &sigma; consuming constant + log<sub>2</sub>&sigma; bits on average. In
* addition it can deliver the resuls in the form of a UniformInteger
* consuming a constant number of bits on average (and then about
* log<sub>2</sub>&sigma; additional bits are required to convert the
* UniformInteger into an integer sample). The consumption of random bits
* (for \e bits = 1) is summarized here \verbatim
min mean max
bits for UniformInteger result 30.4 34.3 35.7
bits for integer result - log2(sigma) 28.8 31.2 32.5
toll 26.8 29.1 30.4 \endverbatim
* These results are averaged over many samples. The "min" results apply
* when &sigma; is a power of 2; the "max" results apply when &sigma; is
* slightly more than a power of two; the "mean" results are averaged over
* &sigma;. The toll is the excess number of bits over the entropy of the
* distribution, which is log<sub>2</sub>(2&pi;\e e)/2 +
* log<sub>2</sub>&sigma; (for &sigma; large).
*
* The data for the first two lines in the table above are taken for the blue
* and red lines in the figure below where the abscissa is the fractional
* part of log<sub>2</sub>&sigma;
* \image html
* discrete-bits.png "Random bits to sample from discrete normal distribution"
*
* This class uses a mutable RandomNumber objects. So a single
* DiscreteNormalAlt object cannot safely be used by multiple threads. In a
* multi-processing environment, each thread should use a thread-specific
* DiscreteNormalAlt object.
*
* @tparam IntType the integer type to use (default int).
**********************************************************************/
template<typename IntType = int, int bits = 1> class DiscreteNormalAlt {
public:
/**
* Constructor.
*
* @param[in] sigma_num the numerator of &sigma;.
* @param[in] sigma_den the denominator of &sigma; (default 1).
* @param[in] mu_num the numerator of &mu; (default 0).
* @param[in] mu_den the denominator of &mu; (default 1).
*
* This may throw an exception is the parameters are such that overflow is
* possible.
**********************************************************************/
DiscreteNormalAlt(IntType sigma_num, IntType sigma_den = 1,
IntType mu_num = 0, IntType mu_den = 1);
/**
* Return a sample.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return discrete normal sample in the form of a UniformInteger object.
**********************************************************************/
template<class Random>
UniformInteger<IntType, bits> operator()(Random& r) const;
private:
/**
* sigma = _sig / _d, mu = _imu + _mu / _d, _isig = ceil(sigma)
**********************************************************************/
IntType _sig, _mu, _d, _isig, _imu;
mutable RandomNumber<bits> _p;
mutable RandomNumber<bits> _q;
// ceil(n/d) for d > 0
static IntType iceil(IntType n, IntType d);
// abs(n) needed because Visual Studio's std::abs has problems
static IntType iabs(IntType n);
static IntType gcd(IntType u, IntType v);
/**
* Implement outcomes for choosing with prob (\e x + 2\e k) / (2\e k + 2);
* return:
* - 1 (succeed unconditionally) with prob (\e m &minus; 2) / \e m,
* - 0 (succeed with probability x) with prob 1 / \e m,
* - &minus;1 (fail unconditionally) with prob 1 / \e m.
**********************************************************************/
template<class Random> static int Choose(Random& r, int m);
/**
* Return true with probability exp(&minus;1/2).
**********************************************************************/
template<class Random> bool ExpProbH(Random& r) const;
/**
* Return true with probability exp(&minus;<i>n</i>/2).
**********************************************************************/
template<class Random> bool ExpProb(Random& r, int n) const;
/**
* Return \e n with probability exp(&minus;<i>n</i>/2)
* (1&minus;exp(&minus;1/2)).
**********************************************************************/
template<class Random> int ExpProbN(Random& r) const;
/**
* Algorithm B: true with prob exp(-x * (2*k + x) / (2*k + 2)) where
* x = (xn0 + _d * j) / _sig
**********************************************************************/
template<class Random>
bool B(Random& r, int k, IntType xn0, UniformInteger<IntType, bits>& j)
const;
};
template<typename IntType, int bits>
DiscreteNormalAlt<IntType, bits>::DiscreteNormalAlt
(IntType sigma_num, IntType sigma_den, IntType mu_num, IntType mu_den) {
STATIC_ASSERT(std::numeric_limits<IntType>::is_integer,
"DiscreteNormalAlt: invalid integer type IntType");
STATIC_ASSERT(std::numeric_limits<IntType>::is_signed,
"DiscreteNormalAlt: IntType must be a signed type");
if (!( sigma_num > 0 && sigma_den > 0 && mu_den > 0 ))
throw RandomErr("DiscreteNormalAlt: need sigma > 0");
_imu = mu_num / mu_den;
if (_imu == (std::numeric_limits<IntType>::min)())
throw RandomErr("DiscreteNormalAlt: abs(mu) too large");
mu_num -= _imu * mu_den;
IntType l;
l = gcd(sigma_num, sigma_den); sigma_num /= l; sigma_den /= l;
l = gcd(mu_num, mu_den); mu_num /= l; mu_den /= l;
_isig = iceil(sigma_num, sigma_den);
l = gcd(sigma_den, mu_den);
_sig = sigma_num * (mu_den / l);
_mu = mu_num * (sigma_den / l);
_d = sigma_den * (mu_den / l);
// The rest of the constructor tests for possible overflow
// Check for overflow in computing member variables
IntType maxint = (std::numeric_limits<IntType>::max)();
if (!( mu_den / l <= maxint / sigma_num &&
iabs(mu_num) <= maxint / (sigma_den / l) &&
mu_den / l <= maxint / sigma_den ))
throw RandomErr("DiscreteNormalAlt: sigma or mu overflow");
if (!UniformInteger<IntType, bits>::Check(_isig, _d))
throw RandomErr("DiscreteNormalAlt: overflow in UniformInteger");
// The probability that k = kmax is about 10^-543.
int kmax = 50;
// Check that max plausible result fits in an IntType, i.e.,
// _isig * (kmax + 1) + abs(_imu) does not lead to overflow.
if (!( kmax + 1 <= maxint / _isig &&
_isig * (kmax + 1) <= maxint - iabs(_imu) ))
throw RandomErr("DiscreteNormalAlt: possible overflow a");
// Check xn0 = _sig * k + s * _mu;
if (!( kmax <= maxint / _sig &&
_sig * kmax <= maxint - iabs(_mu) ))
throw RandomErr("DiscreteNormalAlt: possible overflow b");
// Check for overflow in RandomNumber::LessThan(..., UniformInteger& j)
// p0 * b, p0 = arg2 = xn0 = _d - 1
// c *= b, c = arg3 = _d
// digit(D, k) * q, digit(D, k) = b - 1, q = arg4 = _sig
if (!( _d <= maxint >> bits ))
throw std::runtime_error("DiscreteNormalAlt: overflow in RandomNumber a");
if (!( (IntType(1) << bits) - 1 <= maxint / _sig ))
throw std::runtime_error("DiscreteNormalAlt: overflow in RandomNumber b");
}
template<typename IntType, int bits> template<class Random>
UniformInteger<IntType, bits>
DiscreteNormalAlt<IntType, bits>::operator()(Random& r) const {
for (;;) {
int k = ExpProbN(r);
if (!ExpProb(r, k * (k - 1))) continue;
UniformInteger<IntType, bits> j(r, _isig);
IntType
s = r.Boolean() ? -1 : 1,
xn0 = _sig * IntType(k) + s * _mu,
i0 = iceil(xn0, _d); // i = i0 + j
xn0 = i0 * _d - xn0; // xn = xn0 + _d * j
if (!j.LessThan(r, _sig - xn0, _d) ||
(k == 0 && s < 0 && !j.GreaterThan(r, -xn0, _d))) continue;
int h = k + 1; while (h-- && B(r, k, xn0, j));
if (!(h < 0)) continue;
j.Add(i0 + s * _imu);
if (s < 0) j.Negate();
return j;
}
}
template<typename IntType, int bits>
IntType DiscreteNormalAlt<IntType, bits>::iceil(IntType n, IntType d)
{ IntType k = n / d; return k + (k * d < n ? 1 : 0); }
template<typename IntType, int bits>
IntType DiscreteNormalAlt<IntType, bits>::iabs(IntType n)
{ return n < 0 ? -n : n; }
template<typename IntType, int bits>
IntType DiscreteNormalAlt<IntType, bits>::gcd(IntType u, IntType v) {
// Knuth, TAOCP, vol 2, 4.5.2, Algorithm A
u = iabs(u); v = iabs(v);
while (v > 0) { IntType r = u % v; u = v; v = r; }
return u;
}
template<typename IntType, int bits> template<class Random>
int DiscreteNormalAlt<IntType, bits>::Choose(Random& r, int m) {
// Limit base to 2^15 to avoid integer overflow
const int b = bits > 15 ? 15 : bits;
const unsigned mask = (1u << b) - 1;
int n1 = m - 2, n2 = m - 1;
// Evaluate u < n/m where u is a random real number in [0,1). Write u =
// (d + u') / 2^b where d is a random integer in [0,2^b) and u' is in
// [0,1). Then u < n/m becomes u' < n'/m where n' = 2^b * n - d * m and
// exit if n' <= 0 (false) or n' >= m (true).
for (;;) {
int d = (mask & RandomNumber<bits>::RandomDigit(r)) * m;
n1 = (std::max)((n1 << b) - d, 0);
if (n1 >= m) return 1;
n2 = (std::min)((n2 << b) - d, m);
if (n2 <= 0) return -1;
if (n1 == 0 && n2 == m) return 0;
}
}
template<typename IntType, int bits> template<class Random>
bool DiscreteNormalAlt<IntType, bits>::ExpProbH(Random& r) const {
_p.Init();
if (_p.Digit(r, 0) >> (bits - 1)) return true;
for (;;) {
_q.Init(); if (!_q.LessThan(r, _p)) return false;
_p.Init(); if (!_p.LessThan(r, _q)) return true;
}
}
template<typename IntType, int bits> template<class Random>
bool DiscreteNormalAlt<IntType, bits>::ExpProb(Random& r, int n) const {
while (n-- > 0) { if (!ExpProbH(r)) return false; }
return true;
}
template<typename IntType, int bits> template<class Random>
int DiscreteNormalAlt<IntType, bits>::ExpProbN(Random& r) const {
int n = 0;
while (ExpProbH(r)) ++n;
return n;
}
template<typename IntType, int bits> template<class Random> bool
DiscreteNormalAlt<IntType, bits>::B(Random& r, int k, IntType xn0,
UniformInteger<IntType, bits>& j)
const {
int n = 0, m = 2 * k + 2, f;
for (;; ++n) {
if ( ((f = k ? 0 : Choose(r, m)) < 0) ||
(_q.Init(),
!(n ? _q.LessThan(r, _p) : _q.LessThan(r, xn0, _d, _sig, j))) ||
((f = k ? Choose(r, m) : f) < 0) ||
(f == 0 && (_p.Init(), !_p.LessThan(r, xn0, _d, _sig, j))) )
break;
_p.swap(_q); // an efficient way of doing p = q
}
return (n % 2) == 0;
}
} // namespace RandomLib
#endif // RANDOMLIB_DISCRETENORMALALT_HPP

241
include/RandomLib/ExactExponential.hpp Normal file
View File

@ -0,0 +1,241 @@
/**
* \file ExactExponential.hpp
* \brief Header for ExactExponential
*
* Sample exactly from an exponential distribution.
*
* Copyright (c) Charles Karney (2006-2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_EXACTEXPONENTIAL_HPP)
#define RANDOMLIB_EXACTEXPONENTIAL_HPP 1
#include <RandomLib/RandomNumber.hpp>
#if defined(_MSC_VER)
// Squelch warnings about constant conditional expressions
# pragma warning (push)
# pragma warning (disable: 4127)
#endif
namespace RandomLib {
/**
* \brief Sample exactly from an exponential distribution.
*
* Sample \e x &ge; 0 from exp(&minus;\e x). See:
* - J. von Neumann, Various Techniques used in Connection with Random
* Digits, J. Res. Nat. Bur. Stand., Appl. Math. Ser. 12, 36--38
* (1951), reprinted in Collected Works, Vol. 5, 768--770 (Pergammon,
* 1963).
* - M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions
* (National Bureau of Standards, 1964), Sec. 26.8.6.c(2).
* - G. E. Forsythe, Von Neumann's Comparison Method for Random Sampling from
* Normal and Other Distributions, Math. Comp. 26, 817--826 (1972).
* - D. E. Knuth, TAOCP, Vol 2, Sec 3.4.1.C.3.
* - D. E. Knuth and A. C. Yao, The Complexity of Nonuniform Random Number
* Generation, in "Algorithms and Complexity" (Academic Press, 1976),
* pp. 357--428.
* - P. Flajolet and N. Saheb, The Complexity of Generating an
* Exponentially Distributed Variate, J. Algorithms 7, 463--488 (1986).
*
* The following code illustrates the basic method given by von Neumann:
* \code
// Return a random number x >= 0 distributed with probability exp(-x).
double ExpDist(RandomLib::Random& r) {
for (unsigned k = 0; ; ++k) {
double x = r.Fixed(), // executed 1/(1-exp(-1)) times on average
p = x, q;
do {
q = r.Fixed(); // executed exp(x)*cosh(x) times on average
if (!(q < p)) return k + x;
p = r.Fixed(); // executed exp(x)*sinh(x) times on average
} while (p < q);
}
}
\endcode
* This returns a result consuming exp(1)/(1 &minus; exp(-1)) = 4.30 random
* numbers on average. (Von Neumann incorrectly states that the method takes
* (1 + exp(1))/(1 &minus; exp(-1)) = 5.88 random numbers on average.)
* Because of the finite precision of Random::Fixed(), the code snippet above
* only approximates exp(&minus;\e x). Instead, it returns \e x with
* probability \e h(1 &minus; \e h)<sup><i>x</i>/<i>h</i></sup> for \e x = \e
* ih, \e h = 2<sup>&minus;53</sup>, and integer \e i &ge; 0.
*
* The above is precisely von Neumann's method. Abramowitz and Stegun
* consider a variant: sample uniform variants until the first is less than
* the sum of the rest. Forsythe converts the < ranking for the runs to &le;
* which makes the analysis of the discrete case more difficult. He also
* drops the "trick" by which the integer part of the deviate is given by the
* number of rejections when finding the fractional part.
*
* Von Neumann says of his method: "The machine has in effect computed a
* logarithm by performing only discriminations on the relative magnitude of
* numbers in (0,1). It is a sad fact of life, however, that under the
* particular conditions of the Eniac it was slightly quicker to use a
* truncated power series for log(1&minus;\e T) than to carry out all the
* discriminations."
*
* Here the code is modified to make it more efficient:
* \code
// Return a random number x >= 0 distributed with probability exp(-x).
double ExpDist(RandomLib::Random& r) {
for (unsigned k = 0; ; ++k) {
double x = r.Fixed(); // executed 1/(1-exp(-1/2)) times on average
if (x >= 0.5) continue;
double p = x, q;
do {
q = r.Fixed(); // executed exp(x)*cosh(x) times on average
if (!(q < p)) return 0.5 * k + x;
p = r.Fixed(); // executed exp(x)*sinh(x) times on average
} while (p < q);
}
}
\endcode
* In addition, the method is extended to use infinite precision uniform
* deviates implemented by RandomNumber and returning \e exact results for
* the exponential distribution. This is possible because only comparisons
* are done in the method. The template parameter \e bits specifies the
* number of bits in the base used for RandomNumber (i.e., base =
* 2<sup><i>bits</i></sup>).
*
* For example the following samples from an exponential distribution and
* prints various representations of the result.
* \code
#include <RandomLib/RandomNumber.hpp>
#include <RandomLib/ExactExponential.hpp>
RandomLib::Random r;
const int bits = 1;
RandomLib::ExactExponential<bits> edist;
for (size_t i = 0; i < 10; ++i) {
RandomLib::RandomNumber<bits> x = edist(r); // Sample
std::pair<double, double> z = x.Range();
std::cout << x << " = " // Print in binary with ellipsis
<< "(" << z.first << "," << z.second << ")"; // Print range
double v = x.Value<double>(r); // Round exactly to nearest double
std::cout << " = " << v << "\n";
}
\endcode
* Here's a possible result: \verbatim
0.0111... = (0.4375,0.5) = 0.474126
10.000... = (2,2.125) = 2.05196
1.00... = (1,1.25) = 1.05766
0.010... = (0.25,0.375) = 0.318289
10.1... = (2.5,3) = 2.8732
0.0... = (0,0.5) = 0.30753
0.101... = (0.625,0.75) = 0.697654
0.00... = (0,0.25) = 0.0969214
0.0... = (0,0.5) = 0.194053
0.11... = (0.75,1) = 0.867946 \endverbatim
* First number is in binary with ... indicating an infinite sequence of
* random bits. Second number gives the corresponding interval. Third
* number is the result of filling in the missing bits and rounding exactly
* to the nearest representable double.
*
* This class uses some mutable RandomNumber objects. So a single
* ExactExponential object cannot safely be used by multiple threads. In a
* multi-processing environment, each thread should use a thread-specific
* ExactExponential object. In addition, these should be invoked with
* thread-specific random generator objects.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 1> class ExactExponential {
public:
/**
* Return a random deviate with an exponential distribution, exp(&minus;\e
* x) for \e x &ge; 0. Requires 7.232 bits per invocation for \e bits = 1.
* The average number of bits in the fraction = 1.743. The relative
* frequency of the results for the fractional part with \e bits = 1 is
* shown by the histogram
* \image html exphist.png
* The base of each rectangle gives the range represented by the
* corresponding binary number and the area is proportional to its
* frequency. A PDF version of this figure is given
* <a href="exphist.pdf">here</a>. This allows the figure to be magnified
* to show the rectangles for all binary numbers up to 9 bits. Note that
* this histogram was generated using an earlier version of
* ExactExponential (thru version 1.3) that implements von Neumann's
* original method. The histogram generated with the current version of
* ExactExponential is the same as this figure for \e u in [0, 1/2]. The
* histogram for \e u in [1/2, 1] is obtained by shifting and scaling the
* part for \e u in [0, 1/2] to fit under the exponential curve.
*
* Another way of assessing the efficiency of the algorithm is thru the
* mean value of the balance = (number of random bits consumed) &minus;
* (number of bits in the result). If we code the result in mixed Knuth
* and Yao's unary-binary notation (integer is given in unary, followed by
* "0" as a separator, followed by the fraction in binary), then the mean
* balance is 3.906. (Flajolet and Saheb analyzed the algorithm based on
* the original von Neumann method and showed that the balance is 5.680 in
* that case.)
*
* For \e bits large, the mean number of random digits is exp(1/2)/(1
* &minus; exp(&minus;1/2)) = 4.19 (versus 4.30 digits for the original
* method).
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return the random sample.
**********************************************************************/
template<class Random> RandomNumber<bits> operator()(Random& r) const;
private:
/**
* Return true with probability exp(&minus;\e p) for \e p in (0,1/2).
**********************************************************************/
template<class Random> bool
ExpFraction(Random& r, RandomNumber<bits>& p) const;
mutable RandomNumber<bits> _x;
mutable RandomNumber<bits> _v;
mutable RandomNumber<bits> _w;
};
template<int bits> template<class Random> RandomNumber<bits>
ExactExponential<bits>::operator()(Random& r) const {
// A simple rejection method gives the 1/2 fractional part. The number of
// rejections gives the multiples of 1/2.
//
// bits: used fract un-bin balance double
// original stats: 9.31615 2.05429 3.63628 5.67987 61.59456
// new stats: 7.23226 1.74305 3.32500 3.90725 59.82198
//
// The difference between un-bin and fract is exp(1)/(exp(1)-1) = 1.58198
_x.Init();
int k = 0;
while (!ExpFraction(r, _x)) { // Executed 1/(1 - exp(-1/2)) on average
++k;
_x.Init();
}
if (k & 1) _x.RawDigit(0) += 1U << (bits - 1);
_x.AddInteger(k >> 1);
return _x;
}
template<int bits> template<class Random> bool
ExactExponential<bits>::ExpFraction(Random& r, RandomNumber<bits>& p)
const {
// The early bale out
if (p.Digit(r, 0) >> (bits - 1)) return false;
// Implement the von Neumann algorithm
_w.Init();
if (!_w.LessThan(r, p)) // if (w < p)
return true;
while (true) { // Unroll loop to avoid copying RandomNumber
_v.Init(); // v = r.Fixed();
if (!_v.LessThan(r, _w)) // if (v < w)
return false;
_w.Init(); // w = r.Fixed();
if (!_w.LessThan(r, _v)) // if (w < v)
return true;
}
}
} // namespace RandomLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // RANDOMLIB_EXACTEXPONENTIAL_HPP

417
include/RandomLib/ExactNormal.hpp Normal file
View File

@ -0,0 +1,417 @@
/**
* \file ExactNormal.hpp
* \brief Header for ExactNormal
*
* Sample exactly from a normal distribution.
*
* Copyright (c) Charles Karney (2011-2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_EXACTNORMAL_HPP)
#define RANDOMLIB_EXACTNORMAL_HPP 1
#include <RandomLib/RandomNumber.hpp>
#include <algorithm> // for max/min
#if defined(_MSC_VER)
// Squelch warnings about constant conditional expressions
# pragma warning (push)
# pragma warning (disable: 4127)
#endif
namespace RandomLib {
/**
* \brief Sample exactly from a normal distribution.
*
* Sample \e x from exp(&minus;<i>x</i><sup>2</sup>/2) / sqrt(2&pi;). For
* background, see:
* - J. von Neumann, Various Techniques used in Connection with Random
* Digits, J. Res. Nat. Bur. Stand., Appl. Math. Ser. 12, 36--38
* (1951), reprinted in Collected Works, Vol. 5, 768--770 (Pergammon,
* 1963).
* - D. E. Knuth and A. C. Yao, The Complexity of Nonuniform Random Number
* Generation, in "Algorithms and Complexity" (Academic Press, 1976),
* pp. 357--428.
* - P. Flajolet and N. Saheb, The Complexity of Generating an Exponentially
* Distributed Variate, J. Algorithms 7, 463--488 (1986).
*
* The algorithm is given in
* - C. F. F. Karney, <i>Sampling exactly from the normal distribution</i>,
* http://arxiv.org/abs/1303.6257 (Mar. 2013).
* .
* In brief, the algorithm is:
* -# Select an integer \e k &ge; 0 with probability
* exp(&minus;<i>k</i>/2) (1&minus;exp(&minus;1/2)).
* -# Accept with probability
* exp(&minus; \e k (\e k &minus; 1) / 2); otherwise, reject and start
* over at step 1.
* -# Sample a random number \e x uniformly from [0,1).
* -# Accept with probability exp(&minus; \e x (\e x + 2\e k) / 2);
* otherwise, reject and start over at step 1.
* -# Set \e x = \e k + \e x.
* -# With probability 1/2, negate \e x.
* -# Return \e x.
* .
* It is easy to show that this algorithm returns samples from the normal
* distribution with zero mean and unit variance. Futhermore, all these
* steps can be carried out exactly as follows:
* - Step 1:
* - \e k = 0;
* - while (ExpProb(&minus;1/2)) increment \e k by 1.
* - Step 2:
* - \e n = \e k (\e k &minus; 1) / 2;
* - while (\e n > 0)
* { if (!ExpProb(&minus;1/2)) go to step 1; decrement \e n by 1; }
* - Step 4:
* - repeat \e k + 1 times:
* if (!ExpProb(&minus; \e x (\e x + 2\e k) / (2\e k + 2))) go to step 1.
* .
* Here, ExpProb(&minus;\e p) returns true with probability exp(&minus;\e p).
* With \e p = 1/2 (steps 1 and 2), this is implemented with von Neumann's
* rejection technique:
* - Generate a sequence of random numbers <i>U</i><sub><i>i</i></sub> and
* find the greatest \e n such that 1/2 > <i>U</i><sub>1</sub> >
* <i>U</i><sub>2</sub> > . . . > <i>U</i><sub><i>n</i></sub>. (The
* resulting value of \e n may be 0.)
* - If \e n is even, accept and return true; otherwise (\e n odd), reject
* and return false.
* .
* For \e p = \e x (\e x + 2\e k) / (2\e k + 2) (step 4), we generalize von
* Neumann's procedure as follows:
* - Generate two sequences of random numbers <i>U</i><sub><i>i</i></sub>
* and <i>V</i><sub><i>i</i></sub> and find the greatest \e n such that
* both the following conditions hold
* - \e x > <i>U</i><sub>1</sub> > <i>U</i><sub>2</sub> > . . . >
* <i>U</i><sub><i>n</i></sub>;
* - <i>V</i><sub><i>i</i></sub> &lt; (\e x + 2 \e k) / (2 \e k + 2) for
* all \e i in [1, \e n].
* .
* (The resulting value of \e n may be 0.)
* - If \e n is even, accept (return true); otherwise (\e n odd), reject
* (return false).
* .
* Here, instead of testing <i>V</i><sub><i>i</i></sub> &lt; (\e x + 2 \e k)
* / (2 \e k + 2), we carry out the following tests:
* - return true, with probability 2 \e k / (2 \e k + 2);
* - return false, with probability 1 / (2 \e k + 2);
* - otherwise (also with probability 1 / (2 \e k + 2)),
* return \e x > <i>V</i><sub><i>i</i></sub>.
* .
* The resulting method now entails evaluation of simple fractional
* probabilities (e.g., 1 / (2 \e k + 2)), or comparing random numbers (e.g.,
* <i>U</i><sub>1</sub> > <i>U</i><sub>2</sub>). These may be carried out
* exactly with a finite mean running time.
*
* With \e bits = 1, this consumes 30.1 digits on average and the result has
* 1.19 digits in the fraction. It takes about 676 ns to generate a result
* (1460 ns, including the time to round it to a double). With bits = 32, it
* takes 437 ns to generate a result (621 ns, including the time to round it
* to a double). In contrast, NormalDistribution takes about 44 ns to
* generate a double result.
*
* Another way of assessing the efficiency of the algorithm is thru the mean
* value of the balance = (number of random bits consumed) &minus; (number of
* bits in the result). If we code the result in Knuth & Yao's unary-binary
* notation, then the mean balance is 26.6.
*
* For example the following samples from a normal exponential distribution
* and prints various representations of the result.
* \code
#include <RandomLib/RandomNumber.hpp>
#include <RandomLib/ExactNormal.hpp>
RandomLib::Random r;
const int bits = 1;
RandomLib::ExactNormal<bits> ndist;
for (size_t i = 0; i < 10; ++i) {
RandomLib::RandomNumber<bits> x = ndist(r); // Sample
std::pair<double, double> z = x.Range();
std::cout << x << " = " // Print in binary with ellipsis
<< "(" << z.first << "," << z.second << ")"; // Print range
double v = x.Value<double>(r); // Round exactly to nearest double
std::cout << " = " << v << "\n";
}
\endcode
* Here's a possible result: \verbatim
-1.00... = (-1.25,-1) = -1.02142
-0.... = (-1,0) = -0.319708
0.... = (0,1) = 0.618735
-0.0... = (-0.5,0) = -0.396591
0.0... = (0,0.5) = 0.20362
0.0... = (0,0.5) = 0.375662
-1.111... = (-2,-1.875) = -1.88295
-1.10... = (-1.75,-1.5) = -1.68088
-0.... = (-1,0) = -0.577547
-0.... = (-1,0) = -0.890553
\endverbatim
* First number is in binary with ... indicating an infinite sequence of
* random bits. Second number gives the corresponding interval. Third
* number is the result of filling in the missing bits and rounding exactly
* to the nearest representable double.
*
* This class uses some mutable RandomNumber objects. So a single
* ExactNormal object cannot safely be used by multiple threads. In a
* multi-processing environment, each thread should use a thread-specific
* ExactNormal object. In addition, these should be invoked with
* thread-specific random generator objects.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 1> class ExactNormal {
public:
/**
* Return a random deviate with a normal distribution of mean 0 and
* variance 1.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return the random sample.
**********************************************************************/
template<class Random> RandomNumber<bits> operator()(Random& r) const;
private:
/**
* Return true with probability exp(&minus;1/2). For \e bits = 1, this
* consumes, on average, \e t = 2.846 random digits. We have \e t = \e a
* (1&minus;exp(&minus;1/2)) + \e b exp(&minus;1/2), where \e a is the mean
* bit count for false result = 3.786 and \e b is the mean bit count for
* true result = 2.236.
**********************************************************************/
template<class Random> bool ExpProbH(Random& r) const;
/**
* Return true with probability exp(&minus;<i>n</i>/2). For \e bits = 1,
* this consumes, on average, \e t
* (1&minus;exp(&minus;<i>n</i>/2))/(1&minus;exp(&minus;1/2)) random
* digits. A true result uses \e n \e b random digits. A false result
* uses \e a + \e b [exp(&minus;1/2)/(1&minus;exp(&minus;1/2)) &minus;
* <i>n</i> exp(&minus;<i>n</i>/2)/(1&minus;exp(&minus;<i>n</i>/2))] random
* digits.
**********************************************************************/
template<class Random> bool ExpProb(Random& r, unsigned n) const;
/**
* Return \e n with probability exp(&minus;<i>n</i>/2)
* (1&minus;exp(&minus;1/2)). For \e bits = 1, this consumes \e n \e a +
* \e b random digits if the result is \e n. Averaging over \e n this
* becomes (\e b &minus; (\e b &minus; \e a) exp(&minus;1/2))/(1 &minus;
* exp(&minus;1/2)) digits.
**********************************************************************/
template<class Random> unsigned ExpProbN(Random& r) const;
/**
* Return true with probability 1/2. This is similar to r.Boolean() but
* forces all the random results to come thru RandomNumber::RandomDigit.
**********************************************************************/
template<class Random> static bool Boolean(Random& r) {
// A more general implementation which deals with the case where the base
// might be negative is:
//
// const unsigned base = 1u << bits;
// unsigned b;
// do
// b = RandomNumber<bits>::RandomDigit(r);
// while (b == (base / 2) * 2);
// return b & 1u;
return RandomNumber<bits>::RandomDigit(r) & 1u;
}
/**
* Implement outcomes for choosing with prob (\e x + 2\e k) / (2\e k + 2);
* return:
* - 1 (succeed unconditionally) with prob (2\e k) / (2\e k + 2),
* - 0 (succeed with probability x) with prob 1 / (2\e k + 2),
* - &minus;1 (fail unconditionally) with prob 1 / (2\e k + 2).
* .
* This simulates \code
double x = r.Fixed(); // Uniform in [0,1)
x *= (2 * k + 2);
return x < 2 * k ? 1 : (x < 2 * k + 1 ? 0 : -1);
\endcode
**********************************************************************/
template<class Random> static int Choose(Random& r, int k) {
// Limit base to 2^15 to avoid integer overflow
const int b = bits > 15 ? 15 : bits;
const unsigned mask = (1u << b) - 1;
const int m = 2 * k + 2;
int n1 = m - 2, n2 = m - 1;
// Evaluate u < n/m where u is a random real number in [0,1). Write u =
// (d + u') / 2^b where d is a random integer in [0,2^b) and u' is in
// [0,1). Then u < n/m becomes u' < n'/m where n' = 2^b * n - d * m and
// exit if n' <= 0 (false) or n' >= m (true).
while (true) {
int d = (mask & RandomNumber<bits>::RandomDigit(r)) * m;
n1 = (std::max)((n1 << b) - d, 0);
if (n1 >= m) return 1;
n2 = (std::min)((n2 << b) - d, m);
if (n2 <= 0) return -1;
if (n1 == 0 && n2 == m) return 0;
}
}
mutable RandomNumber<bits> _x;
mutable RandomNumber<bits> _p;
mutable RandomNumber<bits> _q;
};
template<int bits> template<class Random>
bool ExactNormal<bits>::ExpProbH(Random& r) const {
// Bit counts
// ExpProbH: 2.846 = 3.786 * (1-exp(-1/2)) + 2.236 * exp(-1/2)
// t = a * (1-exp(-1/2)) + b * exp(-1/2)
// t = mean bit count for result = 2.846
// a = mean bit count for false result = 3.786
// b = mean bit count for true result = 2.236
//
// for bits large
// t = exp(1/2) = 1.6487
// a = exp(1/2)/(2*(1-exp(-1/2))) = 2.0951
// b = exp(1/2)/(2*exp(-1/2)) = 1.3591
//
// Results for Prob(exp(-1)), omitting first test
// total = 5.889, false = 5.347, true = 6.826
//
// Results for Prob(exp(-1)) using ExpProbH(r) && ExpProbH(r),
// total = 4.572 = (1 - exp(-1)) * a + (1 + exp(-1/2)) * exp(-1/2) * b
// false = 4.630 = a + b * exp(-1/2)/(1 + exp(-1/2)),
// true = 4.472 = 2 * b
_p.Init();
if (_p.Digit(r, 0) >> (bits - 1)) return true;
while (true) {
_q.Init(); if (!_q.LessThan(r, _p)) return false;
_p.Init(); if (!_p.LessThan(r, _q)) return true;
}
}
template<int bits> template<class Random>
bool ExactNormal<bits>::ExpProb(Random& r, unsigned n) const {
// Bit counts
// ExpProb(n): t * (1-exp(-n/2))/(1-exp(-1/2))
// ExpProb(n) = true: n * b
// ExpProb(n) = false: a +
// b * (exp(-1/2)/(1-exp(-1/2)) - n*exp(-n/2)/(1-exp(-n/2)))
while (n--) { if (!ExpProbH(r)) return false; }
return true;
}
template<int bits> template<class Random>
unsigned ExactNormal<bits>::ExpProbN(Random& r) const {
// Bit counts
// ExpProbN() = n: n * a + b
unsigned n = 0;
while (ExpProbH(r)) ++n;
return n;
}
template<int bits> template<class Random> RandomNumber<bits>
ExactNormal<bits>::operator()(Random& r) const {
// With bits = 1,
// - mean number of bits used = 30.10434
// - mean number of bits in fraction = 1.18700
// - mean number of bits in result = 3.55257 (unary-binary)
// - mean balance = 30.10434 - 3.55257 = 26.55177
// - mean number of bits to generate a double = 83.33398
// .
// Note
// - unary-binary notation (Knuth + Yao, 1976): write x = n + y, with n =
// integer and y in [0,1). If n >=0, then write (n+1) 1's followed by a
// 0; otherwise (n < 0), write (-n) 0's followed by a 1. Write y as a
// binary fraction.
// - (bits in result) - (bits in fraction) = 2 (for encoding overhead for
// the integer part) + 0.36557, where 0.36557 = (bits used for integer
// part) = sum(k*int(sqrt(2/pi)*exp(-x^2/2), x=k..k+1), k=0..inf)
// - (bits for double) approx (bits used) - (bits in fraction) + 1 (for
// guard bit) + 53.41664 where 53.41664 = (bits in fraction of double) =
// sum((52-l)*int(sqrt(2/pi)*exp(-x^2/2), x=2^l,2^(l+1)), l=-inf..inf)
// This is approximate because it doesn't account for the minimum
// exponent, denormalized numbers, and rounding changing the exponent.
//
while (true) {
// Executed sqrt(2/pi)/(1-exp(-1/2)) = 2.027818889827955 times on
// average.
unsigned k = ExpProbN(r); // the integer part of the result.
if (ExpProb(r, (k - 1) * k)) {
// Probability that this test succeeds is
// (1 - exp(-1/2)) * sum(exp(-k/2) * exp(-(k-1)*k/2), k=0..inf))
// = (1 - exp(-1/2)) * G = 0.689875359564630
// where G = sum(exp(-k^2/2, k=0..inf) = 1.75331414402145
// For k == 0, sample from exp(-x^2/2) for x in [0,1]. This succeeds
// with probability int(exp(-x^2/2),x=0..1).
//
// For general k, substitute x' = x + k in exp(-x'^2/2), and obtain
// exp(-k^2/2) * exp(-x*(x+2*k)/2). So sample from exp(-x*(x+2*k)/2).
// This succeeds with probability int(exp(-x*(x+2*k)/2),x=0..1) =
// int(exp(-x^2/2),x=k..k+1)*exp(k^2/2) =
//
// 0.8556243918921 for k = 0
// 0.5616593588061 for k = 1
// 0.3963669350376 for k = 2
// 0.2974440159655 for k = 3
// 0.2345104783458 for k = 4
// 0.1921445042826 for k = 5
//
// Returns a result with prob sqrt(pi/2) / G = 0.714825772431666;
// otherwise another trip through the outer loop is taken.
_x.Init();
unsigned s = 1;
for (unsigned j = 0; j <= k; ++j) { // execute k + 1 times
bool first;
for (s = 1, first = true; ; s ^= 1, first = false) {
// A simpler algorithm is indicated by ALT, results in
// - mean number of bits used = 29.99968
// - mean number of bits in fraction = 1.55580
// - mean number of bits in result = 3.92137 (unary-binary)
// - mean balance = 29.99968 - 3.92137 = 26.07831
// - mean number of bits to generate a double = 82.86049
// .
// This has a smaller balance (by 0.47 bits). However the number
// of bits in the fraction is larger by 0.37
if (first) { // ALT: if (false) {
// This implements the success prob (x + 2*k) / (2*k + 2).
int y = Choose(r, k);
if (y < 0) break; // the y test fails
_q.Init();
if (y > 0) { // the y test succeeds just test q < x
if (!_q.LessThan(r, _x)) break;
} else { // the y test is ambiguous
// Test max(q, p) < x. List _q before _p since it ends up with
// slightly more digits generated (and these will be used
// subsequently). (_p's digits are immediately thrown away.)
_p.Init(); if (!_x.GreaterPair(r, _q, _p)) break;
}
} else {
// Split off the failure test for k == 0, i.e., factor the prob
// x/2 test into the product: 1/2 (here) times x (in assignment
// of y).
if (k == 0 && Boolean(r)) break;
// ALT: _q.Init(); if (!_q.LessThan(r, first ? _x : _p)) break;
_q.Init(); if (!_q.LessThan(r, _p)) break;
// succeed with prob k == 0 ? x : (x + 2*k) / (2*k + 2)
int y = k == 0 ? 0 : Choose(r, k);
if (y < 0)
break;
else if (y == 0) {
_p.Init(); if (!_p.LessThan(r, _x)) break;
}
}
_p.swap(_q); // a fast way of doing p = q
}
if (s == 0) break;
}
if (s != 0) {
_x.AddInteger(k);
if (Boolean(r)) _x.Negate(); // half of the numbers are negative
return _x;
}
}
}
}
} // namespace RandomLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // RANDOMLIB_EXACTNORMAL_HPP

100
include/RandomLib/ExactPower.hpp Normal file
View File

@ -0,0 +1,100 @@
/**
* \file ExactPower.hpp
* \brief Header for ExactPower
*
* Sample exactly from a power distribution.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_EXACTPOWER_HPP)
#define RANDOMLIB_EXACTPOWER_HPP 1
#include <RandomLib/RandomNumber.hpp>
namespace RandomLib {
/**
* \brief Sample exactly from a power distribution.
*
* Sample exactly from power distribution (<i>n</i> + 1)
* <i>x</i><sup><i>n</i></sup> for \e x in (0,1) and integer \e n &ge; 0 using
* infinite precision. The template parameter \e bits specifies the number
* of bits in the base used for RandomNumber (i.e., base =
* 2<sup><i>bits</i></sup>).
*
* This class uses some mutable RandomNumber objects. So a single ExactPower
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific ExactPower object.
* In addition, these should be invoked with thread-specific random generator
* objects.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 1> class ExactPower {
public:
/**
* Return the random deviate with a power distribution, (<i>n</i> + 1)
* <i>x</i><sup><i>n</i></sup> for \e x in (0,1) and integer \e n &ge; 0.
* Returned result is a RandomNumber with base 2<sup><i>bits</i></sup>.
* For \e bits = 1, the number of random bits in the result and consumed
* are as follows: \verbatim
n random bits
result consumed
0 0 0
1 2 4
2 2.33 6.67
3 2.67 9.24
4 2.96 11.71
5 3.20 14.11
6 3.41 16.45
7 3.59 18.75
8 3.75 21.01
9 3.89 23.25
10 4.02 25.47
\endverbatim
* The relative frequency of the results with \e bits = 1 and \e n = 2 can
* be is shown by the histogram
* \image html powerhist.png
* The base of each rectangle gives the range represented by the
* corresponding binary number and the area is proportional to its
* frequency. A PDF version of this figure
* <a href="powerhist.pdf">here</a>. This allows the figure to be
* magnified to show the rectangles for all binary numbers up to 9 bits.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] n the power.
* @return the random sample.
**********************************************************************/
template<class Random>
RandomNumber<bits> operator()(Random& r, unsigned n) const;
private:
mutable RandomNumber<bits> _x;
mutable RandomNumber<bits> _y;
};
template<int bits> template<class Random> RandomNumber<bits>
ExactPower<bits>::operator()(Random& r, unsigned n) const {
// Return max(u_0, u_1, u_2, ..., u_n). Equivalent to taking the
// (n+1)th root of u_0.
_x.Init();
for (; n--;) {
_y.Init();
// For bits = 1, we can save 1 bit on the first iteration by using a
// technique suggested by Knuth and Yao (1976). When comparing the
// digits of x and y, use 1 bit to determine if the digits are the same.
// If they are, then use another bit for the value of the digit. If they
// are not, then the next bit of the maximum must be 1 (avoiding using a
// second bit). Applying this optimization to subsequent iterations
// (when x already is filled with some bits) gets trickier.
if (_x.LessThan(r, _y))
_x.swap(_y); // x = y;
}
return _x;
}
} // namespace RandomLib
#endif // RANDOMLIB_EXACTPOWER_HPP

69
include/RandomLib/ExponentialDistribution.hpp Normal file
View File

@ -0,0 +1,69 @@
/**
* \file ExponentialDistribution.hpp
* \brief Header for ExponentialDistribution
*
* Sample from an exponential distribution.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_EXPONENTIALDISTRIBUTION_HPP)
#define RANDOMLIB_EXPONENTIALDISTRIBUTION_HPP 1
#include <cmath>
namespace RandomLib {
/**
* \brief The exponential distribution.
*
* Sample from the distribution exp(&minus;<i>x</i>/&mu;) for \e x &ge; 0.
* This uses the logarithm method, see Knuth, TAOCP, Vol 2, Sec 3.4.1.D.
* Example \code
#include <RandomLib/ExponentialDistribution.hpp>
RandomLib::Random r;
std::cout << "Seed set to " << r.SeedString() << "\n";
RandomLib::ExponentialDistribution<double> expdist;
std::cout << "Select from exponential distribution:";
for (size_t i = 0; i < 10; ++i)
std::cout << " " << expdist(r);
std::cout << "\n";
\endcode
*
* @tparam RealType the real type of the results (default double).
**********************************************************************/
template<typename RealType = double> class ExponentialDistribution {
public:
/**
* The type returned by ExponentialDistribution::operator()(Random&)
**********************************************************************/
typedef RealType result_type;
/**
* Return a sample of type RealType from the exponential distribution and
* mean &mu;. This uses Random::FloatU() which avoids taking log(0) and
* allows rare large values to be returned. If &mu; = 1, minimum returned
* value = 0 with prob 1/2<sup><i>p</i></sup>; maximum returned value =
* log(2)(\e p + \e e) with prob 1/2<sup><i>p</i> + <i>e</i></sup>. Here
* \e p is the precision of real type RealType and \e e is the exponent
* range.
*
* @tparam Random the type of RandomCanonical generator.
* @param[in,out] r the RandomCanonical generator.
* @param[in] mu the mean value of the exponential distribution (default 1).
* @return the random sample.
**********************************************************************/
template<class Random>
RealType operator()(Random& r, RealType mu = RealType(1)) const throw();
};
template<typename RealType> template<class Random> inline RealType
ExponentialDistribution<RealType>::operator()(Random& r, RealType mu) const
throw() {
return -mu * std::log(r.template FloatU<RealType>());
}
} // namespace RandomLib
#endif // RANDOMLIB_EXPONENTIALDISTRIBUTION_HPP

166
include/RandomLib/ExponentialProb.hpp Normal file
View File

@ -0,0 +1,166 @@
/**
* \file ExponentialProb.hpp
* \brief Header for ExponentialProb
*
* Return true with probabililty exp(-\e p).
*
* Copyright (c) Charles Karney (2006-2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_EXPONENTIALPROB_HPP)
#define RANDOMLIB_EXPONENTIALPROB_HPP 1
#include <vector>
#include <limits>
#if defined(_MSC_VER)
// Squelch warnings about constant conditional expressions
# pragma warning (push)
# pragma warning (disable: 4127)
#endif
namespace RandomLib {
/**
* \brief The exponential probability.
*
* Return true with probability exp(&minus;\e p). Basic method taken from:\n
* J. von Neumann,\n Various Techniques used in Connection with Random
* Digits,\n J. Res. Nat. Bur. Stand., Appl. Math. Ser. 12, 36--38
* (1951),\n reprinted in Collected Works, Vol. 5, 768--770 (Pergammon,
* 1963).\n See also the references given for the ExactExponential class.
*
* Here the method is extended to be exact by generating sufficient bits in
* the random numbers in the algorithm to allow the unambiguous comparisons
* to be made.
*
* Here's one way of sampling from a normal distribution with zero mean and
* unit variance in the interval [&minus;1,1] with reasonable accuracy:
* \code
#include <RandomLib/Random.hpp>
#include <RandomLib/ExponentialProb.hpp>
double Normal(RandomLib::Random& r) {
double x;
RandomLib::ExponentialProb e;
do
x = r.FloatW();
while ( !e(r, - 0.5 * x * x) );
return x;
}
\endcode
* (Note that the ExactNormal class samples from the normal distribution
* exactly.)
*
* This class uses a mutable private vector. So a single ExponentialProb
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific ExponentialProb
* object.
**********************************************************************/
class ExponentialProb {
private:
typedef unsigned word;
public:
ExponentialProb() : _v(std::vector<word>(alloc_incr)) {}
/**
* Return true with probability exp(&minus;\e p). Returns false if \e p
* &le; 0. For in \e p (0,1], it requires about exp(\e p) random deviates.
* For \e p large, it requires about exp(1)/(1 &minus; exp(&minus;1))
* random deviates.
*
* @tparam RealType the real type of the argument.
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the probability.
* @return true with probability \e p.
**********************************************************************/
template<typename RealType, class Random>
bool operator()(Random& r, RealType p) const;
private:
/**
* Return true with probability exp(&minus;\e p) for \e p in [0,1].
**********************************************************************/
template<typename RealType, class Random>
bool ExpFraction(Random& r, RealType p) const;
/**
* Holds as much of intermediate uniform deviates as needed.
**********************************************************************/
mutable std::vector<word> _v;
/**
* Increment on size of _v.
**********************************************************************/
static const unsigned alloc_incr = 16;
};
template<typename RealType, class Random>
bool ExponentialProb::operator()(Random& r, RealType p) const {
STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer,
"ExponentialProb(): invalid real type RealType");
return p <= 0 || // True if p <=0
// Ensure p - 1 < p. Also deal with IsNaN(p)
( p < RealType(1)/std::numeric_limits<RealType>::epsilon() &&
// exp(a+b) = exp(a) * exp(b)
ExpFraction(r, p < RealType(1) ? p : RealType(1)) &&
( p <= RealType(1) || operator()(r, p - RealType(1)) ) );
}
template<typename RealType, class Random>
bool ExponentialProb::ExpFraction(Random& r, RealType p) const {
// Base of _v is 2^c. Adjust so that word(p) doesn't lose precision.
static const int c = // The Intel compiler needs this to be static??
std::numeric_limits<word>::digits <
std::numeric_limits<RealType>::digits ?
std::numeric_limits<word>::digits :
std::numeric_limits<RealType>::digits;
// m gives number of valid words in _v
unsigned m = 0, l = unsigned(_v.size());
if (p < RealType(1))
while (true) {
if (p <= RealType(0))
return true;
// p in (0, 1)
if (l == m)
_v.resize(l += alloc_incr);
_v[m++] = r.template Integer<word, c>();
p *= std::pow(RealType(2), c); // p in (0, 2^c)
word w = word(p); // w in [0, 2^c)
if (_v[m - 1] > w)
return true;
else if (_v[m - 1] < w)
break;
else // _v[m - 1] == w
p -= RealType(w); // p in [0, 1)
}
// Here _v < p. Now loop finding decreasing V. Exit when first increasing
// one is found.
for (unsigned s = 0; ; s ^= 1) { // Parity of loop count
for (unsigned j = 0; ; ++j) {
if (j == m) {
// Need more bits in the old V
if (l == m)
_v.resize(l += alloc_incr);
_v[m++] = r.template Integer<word, c>();
}
word w = r.template Integer<word, c>();
if (w > _v[j])
return s != 0u; // New V is bigger, so exit
else if (w < _v[j]) {
_v[j] = w; // New V is smaller, update _v
m = j + 1; // adjusting its size
break; // and generate the next V
}
// Else w == _v[j] and we need to check the next c bits
}
}
}
} // namespace RandomLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // RANDOMLIB_EXPONENTIALPROB_HPP

67
include/RandomLib/InverseEProb.hpp Normal file
View File

@ -0,0 +1,67 @@
/**
* \file InverseEProb.hpp
* \brief Header for InverseEProb
*
* Return true with probabililty 1/\e e.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_INVERSEEPROB_HPP)
#define RANDOMLIB_INVERSEEPROB_HPP 1
#include <vector>
#include <RandomLib/Random.hpp>
namespace RandomLib {
/**
* \brief Return true with probability 1/\e e = exp(&minus;1).
*
* InverseEProb p; p(Random& r) returns true with prob 1/\e e using von
* Neumann's rejection method. It consumes 4.572 bits per call on average.
*
* This class illustrates how to return an exact result using coin tosses
* only. A more efficient way of returning an exact result would be to use
* ExponentialProb p; p(r, 1.0f);
**********************************************************************/
class InverseEProb {
private:
mutable std::vector<bool> _p;
template<class Random> bool exph(Random& r) {
// Return true with prob 1/sqrt(e).
if (r.Boolean()) return true;
_p.clear(); // vector of bits in p
_p.push_back(false);
for (bool s = false; ; s = !s) { // s is a parity
for (size_t i = 0; ; ++i) { // Compare bits of p and q
if (i == _p.size())
_p.push_back(r.Boolean()); // Generate next bit of p if necessary
if (r.Boolean()) { // Half the time the bits differ
if (_p[i]) { // p's bit is 1, so q is smaller, update p
_p[i] = false; // Last bit of q 0
if (++i < _p.size()) _p.resize(i); // p = q
break;
} else
return s; // p's bit is 0, so q is bigger, return parity
} // The other half of the time the bits match, so go to next bit
}
}
}
public:
/**
* Return true with probability 1/\e e.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return true with probability 1/\e e.
**********************************************************************/
template<class Random> bool operator()(Random& r)
{ return exph(r) && exph(r); }
};
} // namespace RandomLib
#endif // RANDOMLIB_INVERSEEPROB_HPP

150
include/RandomLib/InversePiProb.hpp Normal file
View File

@ -0,0 +1,150 @@
/**
* \file InversePiProb.hpp
* \brief Header for InversePiProb
*
* Return true with probabililty 1/&pi;.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_INVERSEPIPROB_HPP)
#define RANDOMLIB_INVERSEPIPROB_HPP 1
#include <cstdlib> // for abs(int)
#include <RandomLib/Random.hpp>
namespace RandomLib {
/**
* \brief Return true with probability 1/&pi;.
*
* InversePiProb p; p(Random& r) returns true with prob 1/&pi; using the
* method of Flajolet et al. It consumes 9.6365 bits per call on average.
*
* The method is given in Section 3.3 of
* - P. Flajolet, M. Pelletier, and M. Soria,<br>
* On Buffon Machines and Numbers,<br> Proc. 22nd ACM-SIAM Symposium on
* Discrete Algorithms (SODA), Jan. 2011.<br>
* http://www.siam.org/proceedings/soda/2011/SODA11_015_flajoletp.pdf <br>
* .
* using the identity
* \f[ \frac 1\pi = \sum_{n=0}^\infty
* {{2n}\choose n}^3 \frac{6n+1}{2^{8n+2}} \f]
*
* It is based on the expression for 1/&pi; given by Eq. (28) of<br>
* - S. Ramanujan,<br>
* Modular Equations and Approximations to &pi;,<br>
* Quart. J. Pure App. Math. 45, 350--372 (1914);<br>
* In Collected Papers, edited by G. H. Hardy, P. V. Seshu Aiyar,
* B. M. Wilson (Cambridge Univ. Press, 1927; reprinted AMS, 2000).<br>
* http://books.google.com/books?id=oSioAM4wORMC&pg=PA36 <br>
* .
* \f[\frac4\pi = 1 + \frac74 \biggl(\frac 12 \biggr)^3
* + \frac{13}{4^2} \biggl(\frac {1\cdot3}{2\cdot4} \biggr)^3
* + \frac{19}{4^3} \biggl(\frac {1\cdot3\cdot5}{2\cdot4\cdot6} \biggr)^3
* + \ldots \f]
*
* The following is a description of how to carry out the algorithm "by hand"
* with a real coin, together with a worked example:
* -# Perform three coin tossing experiments in which you toss a coin until
* you get tails, e.g., <tt>HHHHT</tt>; <tt>HHHT</tt>; <tt>HHT</tt>. Let
* <i>h</i><sub>1</sub> = 4, <i>h</i><sub>2</sub> = 3,
* <i>h</i><sub>3</sub> = 2 be the numbers of heads tossed in each
* experiment.
* -# Compute <i>n</i> = &lfloor;<i>h</i><sub>1</sub>/2&rfloor; +
* &lfloor;<i>h</i><sub>2</sub>/2&rfloor; +
* mod(&lfloor;(<i>h</i><sub>3</sub> &minus; 1)/3&rfloor;, 2) = 2 + 1 + 0
* = 3. Here is a table of the 3 contributions to <i>n</i>:\verbatim
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 h
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 floor(h1/2)
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 floor(h2/2)
1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 mod(floor((h3-1)/3), 2)
\endverbatim
* -# Perform three additional coin tossing experiments in each of which you
* toss a coin 2<i>n</i> = 6 times, e.g., <tt>TTHHTH</tt>;
* <tt>HHTHH|H</tt>; <tt>THHHHH</tt>. Are the number of heads and tails
* equal in each experiment? <b>yes</b> and <b>no</b> and <b>no</b> &rarr;
* <b>false</b>. (Here, you can give up at the |.)
* .
* The final result in this example is <b>false</b>. The most common way a
* <b>true</b> result is obtained is with <i>n</i> = 0, in which case the
* last step vacuously returns <b>true</b>.
*
* Proof of the algorithm: Flajolet et al. rearrange Ramanujan's identity as
* \f[ \frac 1\pi = \sum_{n=0}^\infty
* \biggl[{2n\choose n} \frac1{2^{2n}} \biggr]^3
* \frac{6n+1}{2^{2n+2}}. \f]
* Noticing that
* \f[ \sum_{n=0}^\infty
* \frac{6n+1}{2^{2n+2}} = 1, \f]
* the algorithm becomes:
* -# pick <i>n</i> &ge; 0 with prob (6<i>n</i>+1) / 2<sup>2<i>n</i>+2</sup>
* (mean <i>n</i> = 11/9);
* -# return <b>true</b> with prob (binomial(2<i>n</i>, <i>n</i>) /
* 2<sup>2<i>n</i></sup>)<sup>3</sup>.
*
* Implement (1) as
* - geom4(r) + geom4(r) returns <i>n</i> with probability 9(<i>n</i> +
* 1) / 2<sup>2<i>n</i>+4</sup>;
* - geom4(r) + geom4(r) + 1 returns <i>n</i> with probability
* 36<i>n</i> / 2<sup>2<i>n</i>+4</sup>;
* - combine these with probabilities [4/9, 5/9] to yield (6<i>n</i> +
* 1) / 2<sup>2<i>n</i>+2</sup>, as required.
* .
* Implement (2) as the outcome of 3 coin tossing experiments of 2<i>n</i>
* tosses with success defined as equal numbers of heads and tails in each
* trial.
*
* This class illustrates how to return an exact result using coin tosses
* only. A more efficient implementation (which is still exact) would
* replace prob59 by r.Prob(5,9) and geom4 by LeadingZeros z; z(r)/2.
**********************************************************************/
class InversePiProb {
private:
template<class Random> bool prob59(Random& r) {
// true with prob 5/9 = 0.1 000 111 000 111 000 111 ... (binary expansion)
if (r.Boolean()) return true;
for (bool res = false; ; res = !res)
for (int i = 3; i--; ) if (r.Boolean()) return res;
}
template<class Random> int geom4(Random& r) { // Geom(1/4)
int sum = 0;
while (r.Boolean() && r.Boolean()) ++sum;
return sum;
}
template<class Random> bool binom(Random& r, int n) {
// Probability of equal heads and tails on 2*n tosses
// = binomial(2*n, n) / 2^(2*n)
int d = 0;
for (int k = n; k--; ) d += r.Boolean() ? 1 : -1;
for (int k = n; k--; ) {
d += r.Boolean() ? 1 : -1;
// This optimization saves 0.1686 bit per call to operator() on average.
if (std::abs(d) > k) return false;
}
return true;
}
public:
/**
* Return true with probability 1/&pi;.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return true with probability 1/&pi;.
**********************************************************************/
template<class Random> bool operator()(Random& r) {
// Return true with prob 1/pi.
int n = geom4(r) + geom4(r) + (prob59(r) ? 1 : 0);
for (int j = 3; j--; ) if (!binom(r, n)) return false;
return true;
}
};
} // namespace RandomLib
#endif // RANDOMLIB_INVERSEPIPROB_HPP

84
include/RandomLib/LeadingZeros.hpp Normal file
View File

@ -0,0 +1,84 @@
/**
* \file LeadingZeros.hpp
* \brief Header for LeadingZeros
*
* Count the leading zeros in a real number.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_LEADINGZEROS_HPP)
#define RANDOMLIB_LEADINGZEROS_HPP 1
namespace RandomLib {
/**
* \brief Count of leading zeros.
*
* Count of leading zero bits after the binary point in a real number
* uniformly distributed in (0,1). (This is equivalent to the geometric
* distribution with probability 1/2.) For example
* \code
#include <RandomLib/LeadingZeros.hpp>
RandomLib::Random r; // A RandomGenerator works here too
std::cout << "Seed set to " << r.SeedString() << "\n";
LeadingZeros zeros;
std::cout << "Count of leading zeros:";
for (size_t i = 0; i < 20; ++i)
std::cout << " " << zeros(r);
std::cout << "\n";
\endcode
**********************************************************************/
class LeadingZeros {
public:
/**
* Return the number of zero bits after the binary point in a real number
* uniformly distributed in (0,1). Thus \e k is returned with probability
* 1/2<sup><i>k</i>+1</sup>. Because MT19937 is \e not a perfect random
* number generator, this always returns a result in [0, 19937).
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return the random sample.
**********************************************************************/
template<class Random> unsigned operator()(Random& r) const throw();
};
template<class Random>
unsigned LeadingZeros::operator()(Random& r) const throw() {
// It's simpler to count the number of trailing ones in each w-bit block
// stopping when we get to a zero bit.
//
// Process a word in chunks of size m. The algorithm here can deal with
// any m assuming that z is modified accordingly. m = 4 is an approximate
// optimum.
//
// Can also adapt this routine to use RandomNumber::highest_bit_idx
// instead. However the result is considerably slower.
const int m = 4;
STATIC_ASSERT(m <= Random::width, "LeadingZeros: m too large");
// mask with m low bits set
const typename Random::result_type mask = ~(Random::max << m);
// Number of trailing 1 bits in [0, 1<<m). However, correct results are
// also obtained with any permutation of this array. This particular
// permutation is useful since the initial 1/2, 1/4, etc. can be used for
// m-1, m-2, etc. To generate the array for the next higher m, append a
// duplicate of the array and increment the last entry by one.
const unsigned z[1 << m] =
{ 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, };
typename Random::result_type x = r();
for (unsigned b = m, n = 0; b < Random::width; b += m) {
n += z[x & mask]; // count trailing 1s in chunk
if (n < b) // chunk contains a 0
return n;
x >>= m; // shift out the chunk we've processed
}
// x is all ones (prob 1/2^w); process the next word.
return Random::width + operator()(r);
}
} // namespace RandomLib
#endif // RANDOMLIB_LEADINGZEROS_HPP

92
include/RandomLib/MPFRExponential.hpp Normal file
View File

@ -0,0 +1,92 @@
/**
* \file MPFRExponential.hpp
* \brief Header for MPFRExponential
*
* Sampling exactly from the normal distribution for MPFR.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFREXPONENTIAL_HPP)
#define RANDOMLIB_MPFREXPONENTIAL_HPP 1
#include <RandomLib/MPFRRandom.hpp>
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The exponential distribution for MPFR.
*
* This is a transcription of ExactExponential (version 1.4) for use with
* MPFR.
*
* This class uses mutable private objects. So a single MPFRExponential
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRExponential
* object.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 32> class MPFRExponential {
public:
/**
* Initialize the MPFRExponential object.
**********************************************************************/
MPFRExponential() {};
/**
* Sample from the exponential distribution with mean 1 returning a
* MPFRRandom.
*
* @param[out] t the MPFRRandom result.
* @param[in,out] r a GMP random generator.
**********************************************************************/
void operator()(MPFRRandom<bits>& t, gmp_randstate_t r) const
{ Compute(r); _x.swap(t); }
/**
* Sample from the exponential distribution with mean 1.
*
* @param[out] val the sample from the exponential distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn; if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
{ Compute(r); return _x(val, r, round); }
private:
// Disable copy constructor and assignment operator
MPFRExponential(const MPFRExponential&);
MPFRExponential& operator=(const MPFRExponential&);
int ExpFraction(gmp_randstate_t r, MPFRRandom<bits>& p) const {
// The early bale out
if (p.TestHighBit(r)) return 0;
// Implement the von Neumann algorithm
_w.Init();
if (!_w.LessThan(r, p)) return 1;
while (true) {
_v.Init(); if (!_v.LessThan(r, _w)) return 0;
_w.Init(); if (!_w.LessThan(r, _v)) return 1;
}
}
void Compute(gmp_randstate_t r) const {
_x.Init();
unsigned k = 0;
while (!ExpFraction(r, _x)) { ++k; _x.Init(); }
if (k & 1) _x.SetHighBit(r);
_x.AddInteger(k >> 1);
return;
}
mutable MPFRRandom<bits> _x;
mutable MPFRRandom<bits> _v;
mutable MPFRRandom<bits> _w;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFREXPONENTIAL_HPP

123
include/RandomLib/MPFRExponentialL.hpp Normal file
View File

@ -0,0 +1,123 @@
/**
* \file MPFRExponentialL.hpp
* \brief Header for MPFRExponentialL
*
* Sampling exactly from the exponential distribution for MPFR using the
* traditional method.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFREXPONENTIALL_HPP)
#define RANDOMLIB_MPFREXPONENTIALL_HPP 1
#include <cmath> // for log
#include <mpfr.h>
#define HAVE_MPFR (MPFR_VERSION_MAJOR >= 3)
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The exponential distribution for MPFR (the log method).
*
* This class is <b>DEPRECATED</b>. It is included for illustrative purposes
* only. The MPFRExponential class provides a much more efficient method for
* sampling from the exponential distribution.
*
* This is an adaption of ExponentialDistribution to MPFR. The changes are
* - Use MPFR's random number generator
* - Use sufficient precision internally to ensure that a correctly rounded
* result is returned.
*
* This class uses mutable private objects. So a single MPFRExponentialL
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRExponentialL
* object.
**********************************************************************/
class MPFRExponentialL {
private:
// The number of bits of randomness to add at a time.
static const long chunk_ = 32;
public:
/**
* Initialize the MPFRExponentialL object.
**********************************************************************/
MPFRExponentialL() {
mpz_init(_vi);
mpfr_init2(_eps, chunk_);
mpfr_init2(_v1, chunk_);
mpfr_init2(_v2, chunk_);
}
/**
* Destroy the MPFRExponentialL object.
**********************************************************************/
~MPFRExponentialL() {
mpfr_clear(_v2);
mpfr_clear(_v1);
mpfr_clear(_eps);
mpz_clear(_vi);
}
/**
* Sample from the exponential distribution with mean 1.
*
* @param[out] val the sample from the exponential distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn; if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const {
mpfr_prec_t prec0 = mpfr_get_prec (val);
mpfr_prec_t prec = prec0 + 10; // A rough optimum
mpz_urandomb(_vi, r, prec);
mpfr_set_ui_2exp(_eps, 1u, -prec, MPFR_RNDN);
mpfr_set_prec(_v1, prec);
mpfr_set_z_2exp(_v1, _vi, -prec, MPFR_RNDN);
mpfr_set_prec(_v2, prec);
mpfr_add(_v2, _v1, _eps, MPFR_RNDN);
while (true) {
int f2 = mpfr_log(val, _v2, round); // val = log(upper bound)
mpfr_set_prec(_v2, prec0);
int f1 = mpfr_log(_v2, _v1, round); // v2 = log(lower bound)
if (f1 == f2 && mpfr_equal_p(val, _v2)) {
mpfr_neg(val, val, MPFR_RNDN);
return -f1;
}
prec = Refine(r, prec);
}
}
private:
// disable copy constructor and assignment operator
MPFRExponentialL(const MPFRExponentialL&);
MPFRExponentialL& operator=(const MPFRExponentialL&);
// Refine the random interval
mpfr_prec_t Refine(gmp_randstate_t r, mpfr_prec_t prec)
const {
prec += chunk_;
mpfr_div_2ui(_eps, _eps, chunk_, MPFR_RNDN);
mpz_urandomb(_vi, r, chunk_);
mpfr_set_prec(_v2, prec);
mpfr_set_z_2exp(_v2, _vi, -prec, MPFR_RNDN);
mpfr_add(_v2, _v1, _v2, MPFR_RNDN);
mpfr_swap(_v1, _v2); // v1 = v2;
mpfr_set_prec(_v2, prec);
mpfr_add(_v2, _v1, _eps, MPFR_RNDN);
return prec;
}
mutable mpz_t _vi;
mutable mpfr_t _eps;
mutable mpfr_t _v1;
mutable mpfr_t _v2;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFREXPONENTIALL_HPP

144
include/RandomLib/MPFRNormal.hpp Normal file
View File

@ -0,0 +1,144 @@
/**
* \file MPFRNormal.hpp
* \brief Header for MPFRNormal
*
* Sampling exactly from the normal distribution for MPFR.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRNORMAL_HPP)
#define RANDOMLIB_MPFRNORMAL_HPP 1
#include <algorithm> // for max/min
#include <RandomLib/MPFRRandom.hpp>
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The normal distribution for MPFR.
*
* This is a transcription of ExactNormal (version 1.3) for use with MPFR.
*
* This class uses mutable private objects. So a single MPFRNormal object
* cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRNormal object.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 32> class MPFRNormal {
public:
/**
* Initialize the MPFRNormal object.
**********************************************************************/
MPFRNormal() { mpz_init(_tt); }
/**
* Destroy the MPFRNormal object.
**********************************************************************/
~MPFRNormal() { mpz_clear(_tt); }
/**
* Sample from the normal distribution with mean 0 and variance 1 returning
* a MPFRRandom.
*
* @param[out] t the MPFRRandom result.
* @param[in,out] r a GMP random generator.
**********************************************************************/
void operator()(MPFRRandom<bits>& t,gmp_randstate_t r) const
{ Compute(r); return _x.swap(t); }
/**
* Sample from the normal distribution with mean 0 and variance 1.
*
* @param[out] val the sample from the normal distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn;1 if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
{ Compute(r); return _x(val, r, round); }
private:
// Disable copy constructor and assignment operator
MPFRNormal(const MPFRNormal&);
MPFRNormal& operator=(const MPFRNormal&);
// True with prob exp(-1/2)
int ExpProbH(gmp_randstate_t r) const {
_p.Init(); if (_p.TestHighBit(r)) return 1;
// von Neumann rejection
while (true) {
_q.Init(); if (!_q.LessThan(r, _p)) return 0;
_p.Init(); if (!_p.LessThan(r, _q)) return 1;
}
}
// True with prob exp(-n/2)
int ExpProb(gmp_randstate_t r, unsigned n) const {
while (n--) { if (!ExpProbH(r)) return 0; }
return 1;
}
// n with prob (1-exp(-1/2)) * exp(-n/2)
unsigned ExpProbN(gmp_randstate_t r) const {
unsigned n = 0;
while (ExpProbH(r)) ++n;
return n;
}
// Return:
// 1 with prob 2k/(2k + 2)
// 0 with prob 1/(2k + 2)
// -1 with prob 1/(2k + 2)
int Choose(gmp_randstate_t r, int k) const {
const int b = 15; // To avoid integer overflow on multiplication
const int m = 2 * k + 2;
int n1 = m - 2, n2 = m - 1;
while (true) {
mpz_urandomb(_tt, r, b);
int d = int( mpz_get_ui(_tt) ) * m;
n1 = (std::max)((n1 << b) - d, 0);
if (n1 >= m) return 1;
n2 = (std::min)((n2 << b) - d, m);
if (n2 <= 0) return -1;
if (n1 == 0 && n2 == m) return 0;
}
}
void Compute(gmp_randstate_t r) const {
while (true) {
unsigned k = ExpProbN(r); // the integer part of the result.
if (ExpProb(r, (k - 1) * k)) {
_x.Init();
unsigned s = 1;
for (unsigned j = 0; j <= k; ++j) { // execute k + 1 times
bool first;
for (s = 1, first = true; ; s ^= 1, first = false) {
if (k == 0 && _x.Boolean(r)) break;
_q.Init(); if (!_q.LessThan(r, first ? _x : _p)) break;
int y = k == 0 ? 0 : Choose(r, k);
if (y < 0)
break;
else if (y == 0) {
_p.Init(); if (!_p.LessThan(r, _x)) break;
}
_p.swap(_q); // a fast way of doing p = q
}
if (s == 0) break;
}
if (s != 0) {
_x.AddInteger(k);
if (_x.Boolean(r)) _x.Negate();
return;
}
}
}
}
mutable mpz_t _tt; // A temporary
mutable MPFRRandom<bits> _x;
mutable MPFRRandom<bits> _p;
mutable MPFRRandom<bits> _q;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRNORMAL_HPP

138
include/RandomLib/MPFRNormalK.hpp Normal file
View File

@ -0,0 +1,138 @@
/**
* \file MPFRNormalK.hpp
* \brief Header for MPFRNormalK
*
* Sampling exactly from the normal distribution for MPFR.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRNORMALK_HPP)
#define RANDOMLIB_MPFRNORMALK_HPP 1
#include <algorithm> // for max
#include <RandomLib/MPFRRandom.hpp>
#include <RandomLib/MPFRExponential.hpp>
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The normal distribution for MPFR (Kahn algorithm).
*
* This class is <b>DEPRECATED</b>. It is included for illustrative purposes
* only. The MPFRNormal class provides a somewhat more efficient method for
* sampling from the normal distribution.
*
* Refs:
* - H. Kahn, Rand Report RM-1237-AEC, p. 41 (1954).
* - M. Abramowitz and I. A. Stegun, p. 953, Sec. 26.8.6.a(4) (1964).
* .
* N.B. Damien Stehle' drew my attention to this algorithm as a useful way to
* compute normal deviates exactly.
*
* This class uses mutable private objects. So a single MPFRNormalK object
* cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRNormalK object.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 32> class MPFRNormalK {
public:
/**
* Initialize the MPFRNormalK object.
**********************************************************************/
MPFRNormalK()
{ mpfr_init2(_xf, MPFR_PREC_MIN); mpfr_init2(_zf, MPFR_PREC_MIN); }
/**
* Destroy the MPFRNormalK object.
**********************************************************************/
~MPFRNormalK()
{ mpfr_clear(_zf); mpfr_clear(_xf); }
/**
* Sample from the normal distribution with mean 0 and variance 1 returning
* a MPFRRandom.
*
* @param[out] t the MPFRRandom result.
* @param[in,out] r a GMP random generator.
**********************************************************************/
void operator()(MPFRRandom<bits>& t, gmp_randstate_t r) const
{ Compute(r); _x.swap(t); }
/**
* Sample from the normal distribution with mean 0 and variance 1.
*
* @param[out] val the sample from the normal distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn;1 if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
{ Compute(r); return _x(val, r, round); }
private:
// disable copy constructor and assignment operator
MPFRNormalK(const MPFRNormalK&);
MPFRNormalK& operator=(const MPFRNormalK&);
void Compute(gmp_randstate_t r) const {
// The algorithm is sample x and z from the exponential distribution; if
// (x-1)^2 < 2*z, return (random sign)*x; otherwise repeat. Probability
// of acceptance is sqrt(pi/2) * exp(-1/2) = 0.7602.
while (true) {
_edist(_x, r);
_edist(_z, r);
for (mp_size_t k = 1; ; ++k) {
_x.ExpandTo(r, k - 1);
_z.ExpandTo(r, k - 1);
mpfr_prec_t prec = (std::max)(mpfr_prec_t(MPFR_PREC_MIN), k * bits);
mpfr_set_prec(_xf, prec);
mpfr_set_prec(_zf, prec);
// Try for acceptance first; so compute upper limit on (y-1)^2 and
// lower limit on 2*z.
if (_x.UInteger() == 0) {
_x(_xf, MPFR_RNDD);
mpfr_ui_sub(_xf, 1u, _xf, MPFR_RNDU);
} else {
_x(_xf, MPFR_RNDU);
mpfr_sub_ui(_xf, _xf, 1u, MPFR_RNDU);
}
mpfr_sqr(_xf, _xf, MPFR_RNDU);
_z(_zf, MPFR_RNDD);
mpfr_mul_2ui(_zf, _zf, 1u, MPFR_RNDD);
if (mpfr_cmp(_xf, _zf) < 0) { // (y-1)^2 < 2*z, so accept
if (_x.Boolean(r)) _x.Negate(); // include a random sign
return;
}
// Try for rejection; so compute lower limit on (y-1)^2 and upper
// limit on 2*z.
if (_x.UInteger() == 0) {
_x(_xf, MPFR_RNDU);
mpfr_ui_sub(_xf, 1u, _xf, MPFR_RNDD);
} else {
_x(_xf, MPFR_RNDD);
mpfr_sub_ui(_xf, _xf, 1u, MPFR_RNDD);
}
mpfr_sqr(_xf, _xf, MPFR_RNDD);
_z(_zf, MPFR_RNDU);
mpfr_mul_2ui(_zf, _zf, 1u, MPFR_RNDU);
if (mpfr_cmp(_xf, _zf) > 0) // (y-1)^2 > 2*z, so reject
break;
// Otherwise repeat with more precision
}
// Reject and start over with a new y and z
}
}
mutable MPFRRandom<bits> _x;
mutable MPFRRandom<bits> _z;
mutable mpfr_t _xf;
mutable mpfr_t _zf;
const MPFRExponential<bits> _edist;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRNORMALK_HPP

255
include/RandomLib/MPFRNormalR.hpp Normal file
View File

@ -0,0 +1,255 @@
/**
* \file MPFRNormalR.hpp
* \brief Header for MPFRNormalR
*
* Sampling exactly from the normal distribution for MPFR using the ratio
* method.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRNORMALR_HPP)
#define RANDOMLIB_MPFRNORMALR_HPP 1
#include <algorithm> // for max/min
#include <cmath> // for pow
#include <mpfr.h>
#define HAVE_MPFR (MPFR_VERSION_MAJOR >= 3)
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The normal distribution for MPFR (ratio method).
*
* This class is <b>DEPRECATED</b>. It is included for illustrative purposes
* only. The MPFRNormal class provides a much more efficient method for
* sampling from the normal distribution.
*
* This is an adaption of NormalDistribution to MPFR. The changes are
* - Use MPFR's random number generator
* - Use sufficient precision internally to ensure that a correctly rounded
* result is returned.
*
* This class uses a mutable private object. So a single MPFRNormalR
* object cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRNormalR
* object.
**********************************************************************/
class MPFRNormalR {
private:
// The number of bits of randomness to add at a time. Require that Leva's
// bounds "work" at a precision of 2^-chunk and that an unsigned long can
// hold this many bits.
static const long chunk_ = 32;
static const unsigned long m = 3684067834; // ceil(2^chunk*sqrt(2/e))
public:
/**
* Initialize the MPFRNormalR object.
**********************************************************************/
MPFRNormalR() {
mpz_init(_ui);
mpz_init(_vi);
mpfr_init2(_eps, chunk_);
mpfr_init2(_u, chunk_);
mpfr_init2(_v, chunk_);
mpfr_init2(_up, chunk_);
mpfr_init2(_vp, chunk_);
mpfr_init2(_vx, chunk_);
mpfr_init2(_x1, chunk_);
mpfr_init2(_x2, chunk_);
}
/**
* Destroy the MPFRNormalR object.
**********************************************************************/
~MPFRNormalR() {
mpfr_clear(_x2);
mpfr_clear(_x1);
mpfr_clear(_vx);
mpfr_clear(_vp);
mpfr_clear(_up);
mpfr_clear(_v);
mpfr_clear(_u);
mpfr_clear(_eps);
mpz_clear(_vi);
mpz_clear(_ui);
}
/**
* Sample from the normal distribution with mean 0 and variance 1.
*
* @param[out] val the sample from the normal distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn;1 if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const {
const double
s = 0.449871, // Constants from Leva
t = -0.386595,
a = 0.19600 ,
b = 0.25472 ,
r1 = 0.27597 ,
r2 = 0.27846 ,
u1 = 0.606530, // sqrt(1/e) rounded down and up
u2 = 0.606531,
scale = std::pow(2.0, -chunk_); // for turning randoms into doubles
while (true) {
mpz_urandomb(_vi, r, chunk_);
if (mpz_cmp_ui(_vi, m) >= 0) continue; // Very early reject
double vf = (mpz_get_ui(_vi) + 0.5) * scale;
mpz_urandomb(_ui, r, chunk_);
double uf = (mpz_get_ui(_ui) + 0.5) * scale;
double
x = uf - s,
y = vf - t,
Q = x*x + y * (a*y - b*x);
if (Q >= r2) continue; // Early reject
mpfr_set_ui_2exp(_eps, 1u, -chunk_, MPFR_RNDN);
mpfr_prec_t prec = chunk_;
mpfr_set_prec(_u, prec);
mpfr_set_prec(_v, prec);
// (u,v) = sw corner of range
mpfr_set_z_2exp(_u, _ui, -prec, MPFR_RNDN);
mpfr_set_z_2exp(_v, _vi, -prec, MPFR_RNDN);
mpfr_set_prec(_up, prec);
mpfr_set_prec(_vp, prec);
// (up,vp) = ne corner of range
mpfr_add(_up, _u, _eps, MPFR_RNDN);
mpfr_add(_vp, _v, _eps, MPFR_RNDN);
// Estimate how many extra bits will be needed to achieve the desired
// precision.
mpfr_prec_t prec_guard = 3 + chunk_ -
(std::max)(mpz_sizeinbase(_ui, 2), mpz_sizeinbase(_vi, 2));
if (Q > r1) {
int reject;
while (true) {
// Rejection curve v^2 + 4 * u^2 * log(u) < 0 has a peak at u =
// exp(-1/2) = 0.60653066. So treat uf in (0.606530, 0.606531) =
// (u1, u2) specially
// Try for rejection first
if (uf <= u1)
reject = Reject(_u, _vp, prec, MPFR_RNDU);
else if (uf >= u2)
reject = Reject(_up, _vp, prec, MPFR_RNDU);
else { // u in (u1, u2)
mpfr_set_prec(_vx, prec);
mpfr_add(_vx, _vp, _eps, MPFR_RNDN);
reject = Reject(_u, _vx, prec, MPFR_RNDU); // Could use _up too
}
if (reject < 0) break; // tried to reject but failed, so accept
// Try for acceptance
if (uf <= u1)
reject = Reject(_up, _v, prec, MPFR_RNDD);
else if (uf >= u2)
reject = Reject(_u, _v, prec, MPFR_RNDD);
else { // u in (u2, u2)
mpfr_sub(_vx, _v, _eps, MPFR_RNDN);
reject = Reject(_u, _vx, prec, MPFR_RNDD); // Could use _up too
}
if (reject > 0) break; // tried to accept but failed, so reject
prec = Refine(r, prec); // still can't decide, to refine
}
if (reject > 0) continue; // reject, back to outer loop
}
// Now evaluate v/u to the necessary precision
mpfr_prec_t prec0 = mpfr_get_prec (val);
// while (prec < prec0 + prec_guard) prec = Refine(r, prec);
if (prec < prec0 + prec_guard)
prec = Refine(r, prec,
(prec0 + prec_guard - prec + chunk_ - 1) / chunk_);
mpfr_set_prec(_x1, prec0);
mpfr_set_prec(_x2, prec0);
int flag;
while (true) {
int
f1 = mpfr_div(_x1, _v, _up, round), // min slope
f2 = mpfr_div(_x2, _vp, _u, round); // max slope
if (f1 == f2 && mpfr_equal_p(_x1, _x2)) {
flag = f1;
break;
}
prec = Refine(r, prec);
}
mpz_urandomb(_ui, r, 1);
if (mpz_tstbit(_ui, 0)) {
flag = -flag;
mpfr_neg(val, _x1, MPFR_RNDN);
} else
mpfr_set(val, _x1, MPFR_RNDN);
// std::cerr << uf << " " << vf << " " << Q << "\n";
return flag;
}
}
private:
// disable copy constructor and assignment operator
MPFRNormalR(const MPFRNormalR&);
MPFRNormalR& operator=(const MPFRNormalR&);
// Refine the random square
mpfr_prec_t Refine(gmp_randstate_t r, mpfr_prec_t prec, long num = 1)
const {
if (num <= 0) return prec;
// Use _vx as scratch
prec += num * chunk_;
mpfr_div_2ui(_eps, _eps, num * chunk_, MPFR_RNDN);
mpz_urandomb(_ui, r, num * chunk_);
mpfr_set_prec(_up, prec);
mpfr_set_z_2exp(_up, _ui, -prec, MPFR_RNDN);
mpfr_set_prec(_vx, prec);
mpfr_add(_vx, _u, _up, MPFR_RNDN);
mpfr_swap(_u, _vx); // u = vx
mpfr_add(_up, _u, _eps, MPFR_RNDN);
mpz_urandomb(_vi, r, num * chunk_);
mpfr_set_prec(_vp, prec);
mpfr_set_z_2exp(_vp, _vi, -prec, MPFR_RNDN);
mpfr_set_prec(_vx, prec);
mpfr_add(_vx, _v, _vp, MPFR_RNDN);
mpfr_swap(_v, _vx); // v = vx
mpfr_add(_vp, _v, _eps, MPFR_RNDN);
return prec;
}
// Evaluate the sign of the rejection condition v^2 + 4*u^2*log(u)
int Reject(mpfr_t u, mpfr_t v, mpfr_prec_t prec, mpfr_rnd_t round) const {
// Use x1, x2 as scratch
mpfr_set_prec(_x1, prec);
mpfr_log(_x1, u, round);
mpfr_mul(_x1, _x1, u, round); // Important to do the multiplications in
mpfr_mul(_x1, _x1, u, round); // this order so that rounding works right.
mpfr_mul_2ui(_x1, _x1, 2u, round); // 4*u^2*log(u)
mpfr_set_prec(_x2, prec);
mpfr_mul(_x2, v, v, round); // v^2
mpfr_add(_x1, _x1, _x2, round); // v^2 + 4*u^2*log(u)
return mpfr_sgn(_x1);
}
mutable mpz_t _ui;
mutable mpz_t _vi;
mutable mpfr_t _eps;
mutable mpfr_t _u;
mutable mpfr_t _v;
mutable mpfr_t _up;
mutable mpfr_t _vp;
mutable mpfr_t _vx;
mutable mpfr_t _x1;
mutable mpfr_t _x2;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRNORMALR_HPP

383
include/RandomLib/MPFRRandom.hpp Normal file
View File

@ -0,0 +1,383 @@
/**
* \file MPFRRandom.hpp
* \brief Header for MPFRRandom
*
* Utility class for MPFRUniform, MPFRExponential, and MPFRNormal.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRRANDOM_HPP)
#define RANDOMLIB_MPFRRANDOM_HPP 1
#include <algorithm> // for swap
#include <mpfr.h>
#define HAVE_MPFR (MPFR_VERSION_MAJOR >= 3)
#if HAVE_MPFR || defined(DOXYGEN)
/**
* A compile-time assert. Use C++11 static_assert, if available.
**********************************************************************/
#if !defined(STATIC_ASSERT)
# if defined(__GXX_EXPERIMENTAL_CXX0X__)
# define STATIC_ASSERT static_assert
# elif defined(_MSC_VER) && _MSC_VER >= 1600
# define STATIC_ASSERT static_assert
# else
# define STATIC_ASSERT(cond,reason) \
{ enum{ STATIC_ASSERT_ENUM = 1/int(cond) }; }
# endif
#endif
namespace RandomLib {
/**
* \brief Handling random numbers in MPFR.
*
* This class provides roughly the same capabilities as RandomNumber. The
* fraction is represented by a mpz integer \e f and an exponent \e e. We
* have \e e &ge; 0 and 0 &le; \e f < <i>b</i><sup><i>e</i></sup>, and \e b =
* 2<sup><i>bits</i></sup>. This represents the number \e x = \e f
* <i>b</i><sup>&minus;<i>e</i></sup>, with x in [0, 1).
*
* @tparam bits the number of bits in each digit.
*
* \e bits must divide GMP_LIMB_BITS. The default value \e bits = 32 yields
* portable results on all MPFR platforms.
**********************************************************************/
template<int bits = 32> class MPFRRandom {
private:
static const int limb_ = GMP_LIMB_BITS; // How many bits in a limb
static const int loglimb_ = (limb_ == 32 ? 5 :
(limb_ == 64 ? 6 :
(limb_ == 128 ? 7 : -1)));
static const int logbits_ = (bits == 1 ? 0 :
(bits == 2 ? 1 :
(bits == 4 ? 2 :
(bits == 8 ? 3 :
(bits == 16 ? 4 :
(bits == 32 ? 5 :
(bits == 64 ? 6 :
(bits == 128 ? 7 : -1))))))));
static const mp_limb_t mask_ = (bits == limb_ ? ~0UL : // Digit mask
~(~0UL << (bits < limb_ ? bits : 0)));
static const int logw_ = loglimb_ - logbits_; // 2^logw digits per limb
static const unsigned wmask_ = ~(~0U << logw_);
mutable mpz_t _tt; // A temporary
mpz_t _f; // The fraction
mp_size_t _e; // Count of digits
unsigned long _n; // Integer part
int _s; // Sign
void AddDigits(gmp_randstate_t r, long num = 1) { // Add num more digits
if (num <= 0) return;
mpz_mul_2exp(_f, _f, num << logbits_);
mpz_urandomb(_tt, r, num << logbits_);
mpz_add(_f, _f, _tt);
_e += num;
}
// return k'th digit counting k = 0 as most significant
mp_limb_t Digit(gmp_randstate_t r, mp_size_t k) {
ExpandTo(r, k); // Now e > k
k = _e - 1 - k; // Reverse k so k = 0 is least significant
// (k >> logw) is the limb index
// (k & wmask) is the digit position within the limb
return mask_ &
(mpz_getlimbn(_f, k >> logw_) >> ((k & wmask_) << logbits_));
}
// Return index [0..32] of highest bit set. Return 0 if x = 0, 32 is if x
// = ~0. (From Algorithms for programmers by Joerg Arndt.)
static int highest_bit_idx(unsigned long x) throw() {
if (x == 0) return 0;
int r = 1;
// STILL TO DO: handle 64-bit unsigned longs.
if (x & 0xffff0000UL) { x >>= 16; r += 16; }
if (x & 0x0000ff00UL) { x >>= 8; r += 8; }
if (x & 0x000000f0UL) { x >>= 4; r += 4; }
if (x & 0x0000000cUL) { x >>= 2; r += 2; }
if (x & 0x00000002UL) { r += 1; }
return r;
}
public:
/**
* Initialize the MPFRRandom object.
**********************************************************************/
MPFRRandom() : _e(0u), _n(0u), _s(1) {
STATIC_ASSERT(logbits_ >= 0 && loglimb_ >= 0 && logbits_ <= loglimb_,
"MPRFRandom: unsupported value for bits");
mpz_init(_f); mpz_init(_tt);
}
/**
* Initialize the MPFRRandom object from another one.
*
* @param[in] t the MPFRRandom to copy.
**********************************************************************/
MPFRRandom(const MPFRRandom& t) : _e(t._e), _n(t._n), _s(t._s)
{ mpz_init(_f); mpz_set(_f, t._f); mpz_init(_tt); }
/**
* Destroy the MPFRRandom object.
**********************************************************************/
~MPFRRandom() { mpz_clear(_f); mpz_clear(_tt); }
/**
* Assignment operator. (But swapping is typically faster.)
*
* @param[in] t the MPFRRandom to copy.
**********************************************************************/
MPFRRandom& operator=(const MPFRRandom& t) {
_e = t._e;
_n = t._n;
_s = t._s;
mpz_set(_f, t._f); // Don't copy _tt
return *this;
}
/**
* Swap with another MPFRRandom. This is a fast way of doing an
* assignment.
*
* @param[in,out] t the MPFRRandom to swap with.
**********************************************************************/
void swap(MPFRRandom& t) throw() {
if (this != &t) {
std::swap(_e, t._e);
std::swap(_n, t._n);
std::swap(_s, t._s);
mpz_swap(_f, t._f); // Don't swap _tt
}
}
/**
* Reinitialize the MPFRRandom object, setting its value to [0,1].
**********************************************************************/
void Init() { mpz_set_ui(_f, 0u); _e = 0; _n = 0; _s = 1; }
/**
* @return the sign of the MPFRRandom (&plusmn; 1).
**********************************************************************/
int Sign() const throw() { return _s; }
/**
* Change the sign of the MPFRRandom.
**********************************************************************/
void Negate() throw() { _s *= -1; }
/**
* @return the floor of the MPFRRandom
**********************************************************************/
long Floor() const throw() { return _s > 0 ? long(_n) : -1 - long(_n); }
/**
* @return the ceiling of the MPFRRandom
**********************************************************************/
long Ceiling() const throw() { return _s > 0 ? 1 + long(_n) : -long(_n); }
/**
* @return the unsigned integer component of the MPFRRandom.
**********************************************************************/
unsigned long UInteger() const throw() { return _n; }
/**
* @return the number of digits in fraction
**********************************************************************/
unsigned long Size() const throw() { return unsigned(_e); }
/**
* Add integer \e k to the MPRFRandom.
*
* @param[in] k the integer to add.
**********************************************************************/
void AddInteger(long k) {
k += Floor(); // The new floor
int ns = k < 0 ? -1 : 1; // The new sign
if (ns != _s) { // If sign changes, set f = 1 - f
mpz_set_ui(_tt, 1u);
mpz_mul_2exp(_tt, _tt, _e << logbits_);
mpz_sub_ui(_tt, _tt, 1u);
mpz_sub(_f, _tt, _f);
_s = ns;
}
_n = ns > 0 ? k : -(k + 1);
}
/**
* Compare with another MPFRRandom, *this < \e t.
*
* @param[in,out] r a random generator.
* @param[in,out] t a MPFRRandom to compare.
* @return true if *this < \e t.
**********************************************************************/
int LessThan(gmp_randstate_t r, MPFRRandom& t) {
if (this == &t) return false; // same object
if (_s != t._s) return _s < t._s;
if (_n != t._n) return (_s < 0) ^ (_n < t._n);
for (mp_size_t k = 0; ; ++k) {
mp_limb_t x = Digit(r, k);
mp_limb_t y = t.Digit(r, k);
if (x != y) return (_s < 0) ^ (x < y);
}
}
/**
* Set high bit of fraction to 1.
*
* @param[in,out] r a random generator.
**********************************************************************/
void SetHighBit(gmp_randstate_t r) { // Set the msb to 1
ExpandTo(r, 0); // Generate msb if necessary
mpz_setbit(_f, (_e << logbits_) - 1);
}
/**
* Test high bit of fraction.
*
* @param[in,out] r a random generator.
**********************************************************************/
int TestHighBit(gmp_randstate_t r) { // test the msb of f
ExpandTo(r, 0); // Generate msb if necessary
return mpz_tstbit(_f, (_e << logbits_) - 1);
}
/**
* Return the position of the most significant bit in the MPFRRandom.
*
* @param[in,out] r a random generator.
*
* The bit position is numbered such the 1/2 bit is 0, the 1/4 bit is -1,
* etc.
**********************************************************************/
mp_size_t LeadingBit(gmp_randstate_t r) {
if (_n) return highest_bit_idx(_n);
while (true) {
int sgn = mpz_sgn(_f);
if (sgn != 0)
return mp_size_t(mpz_sizeinbase(_f, 2)) - mp_size_t(_e << logbits_);
AddDigits(r);
}
}
/**
* Ensure that the k'th digit of the fraction is computed.
*
* @param[in,out] r a random generator.
* @param[in] k the digit number (0 is the most significant, 1 is the next
* most significant, etc.
**********************************************************************/
void ExpandTo(gmp_randstate_t r, mp_size_t k)
{ if (_e <= k) AddDigits(r, k - _e + 1); }
/**
* Convert to a MPFR number \e without adding more bits.
*
* @param[out] val the value of s * (n + *this).
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn; if val is larger/smaller than
* the exact sample).
*
* If round is MPFR_RNDN, then the rounded midpoint of the interval
* represented by the MPFRRandom is returned. Otherwise it is the rounded
* lower or upper bound of the interval (whichever is appropriate).
**********************************************************************/
int operator()(mpfr_t val, mpfr_rnd_t round)
{ return operator()(val, NULL, round); }
/**
* Convert to a MPFR number.
*
* @param[out] val the value of s * (n + *this).
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn; if val is larger/smaller than
* the exact sample).
*
* If \e r is NULL, then no additional random bits are generated and the
* lower bound, midpoint, or upper bound of the MPFRRandom interval is
* returned, depending on the value of \e round.
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) {
// The value is constructed as a positive quantity, so adjust rounding
// mode to account for this.
switch (round) {
case MPFR_RNDD:
case MPFR_RNDU:
case MPFR_RNDN:
break;
case MPFR_RNDZ:
round = _s < 0 ? MPFR_RNDU : MPFR_RNDD;
break;
case MPFR_RNDA:
round = _s < 0 ? MPFR_RNDD : MPFR_RNDU;
break;
default:
round = MPFR_RNDN; // New rounding modes are variants of N
break;
} // Now round is one of MPFR_RND{D,N,U}
mp_size_t excess;
mpfr_exp_t expt;
if (r == NULL) {
// If r is NULL then all the bits currently generated are considered
// significant. Thus no excess bits need to be squeezed out.
excess = 0;
// And the exponent shift in mpfr_set_z_2exp is just...
expt = -(_e << logbits_);
// However, if rounding to nearest, we need to make room for the
// midpoint bit.
if (round == MPFR_RNDN) {
excess = -1;
--expt;
}
} else { // r is non-NULL
// Generate enough digits, i.e., enough to generate prec significant
// figures for RNDD and RNDU; for RNDN we need to generate an
// additional guard bit.
mp_size_t lead = LeadingBit(r);
mpfr_prec_t prec = mpfr_get_prec (val);
mp_size_t trail = lead - prec; // position one past trailing bit
mp_size_t guard = trail + (round == MPFR_RNDN ? 0 : 1); // guard bit pos
// Generate the bits needed.
if (guard <= 0) ExpandTo(r, (-guard) >> logbits_);
// Unless bits = 1, the generation process will typically have
// generated too many bits. We figure out how many, but leaving room
// for one additional "inexact" bit. The inexact bit is set to 1 in
// order to force MPFR to treat the result as inexact, to break RNDN
// ties, and to get the ternary value set correctly.
//
// expt is the exponent used when forming the number using
// mpfr_set_z_2exp. Without the inexact bit, it's (guard - 1).
// Subtract 1 to account for the inexact bit.
expt = guard - 2;
// The number of excess bits is now the difference between the number
// of bits in the fraction (e << logbits) and -expt. Note that this
// may be -1 (meaning we'll need to shift the number left to
// accommodate the inexact bit).
excess = (_e << logbits_) + expt;
}
mpz_set_ui(_tt, _n); // The integer part
mpz_mul_2exp(_tt, _tt, _e << logbits_); // Shift to allow for fraction
mpz_add(_tt, _tt, _f); // Add fraction
if (excess > 0)
mpz_tdiv_q_2exp(_tt, _tt, excess);
else if (excess < 0)
mpz_mul_2exp(_tt, _tt, -excess);
if (r || round == MPFR_RNDN)
// Set the inexact bit (or compute the midpoint if r is NULL).
mpz_setbit(_tt, 0);
else if (round == MPFR_RNDU)
// If r is NULL, compute the upper bound.
mpz_add_ui(_tt, _tt, 1u);
// Convert to a mpfr number. If r is specified, then there are
// sufficient bits in tt that the result is inexact and that (in the case
// of RNDN) there are no ties.
int flag = mpfr_set_z_2exp(val, _tt, expt, round);
if (_s < 0) {
mpfr_neg (val, val, MPFR_RNDN);
flag = -flag;
}
return flag;
}
/**
* A coin toss. (This should really be a static function. But it uses the
* MPFRRandom temporary variable.)
*
* @param[in,out] r a GMP random generator.
* @return true or false.
**********************************************************************/
int Boolean(gmp_randstate_t r) const {
mpz_urandomb(_tt, r, 1);
return mpz_tstbit(_tt, 0);
}
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRRANDOM_HPP

72
include/RandomLib/MPFRUniform.hpp Normal file
View File

@ -0,0 +1,72 @@
/**
* \file MPFRUniform.hpp
* \brief Header for MPFRUniform
*
* Sampling exactly from a uniform distribution for MPFR.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRUNIFORM_HPP)
#define RANDOMLIB_MPFRUNIFORM_HPP 1
#include <RandomLib/MPFRRandom.hpp>
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The uniform distribution for MPFR.
*
* This is just a thin layer on top of MPFRRandom to provide random numbers
* uniformly distributed in [0,1].
*
* This class uses a mutable private object. So a single MPFRUniform object
* cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRUniform object.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 32> class MPFRUniform {
public:
/**
* Initialize the MPFRUniform object.
**********************************************************************/
MPFRUniform() {};
/**
* Sample from the uniform distribution in [0,1] returning a MPFRRandom.
* This function takes an unused GMP random generator as a parameter, in
* order to parallel the usage of MPFRExponential and MPFRNormal.
*
* @param[out] t the MPFRRandom result.
* @param[in,out] r a GMP random generator (unused).
**********************************************************************/
void operator()(MPFRRandom<bits>& t, gmp_randstate_t r) const
{ Compute(r); _x.swap(t); }
/**
* Sample from the uniform distribution in [0,1].
*
* @param[out] val the sample from the uniform distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (&plusmn; if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
{ Compute(r); return _x(val, r, round); }
private:
// disable copy constructor and assignment operator
MPFRUniform(const MPFRUniform&);
MPFRUniform& operator=(const MPFRUniform&);
void Compute(gmp_randstate_t /* r */) const { _x. Init(); }
mutable MPFRRandom<bits> _x;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRUNIFORM_HPP

114
include/RandomLib/NormalDistribution.hpp Normal file
View File

@ -0,0 +1,114 @@
/**
* \file NormalDistribution.hpp
* \brief Header for NormalDistribution
*
* Compute normal deviates.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_NORMALDISTRIBUTION_HPP)
#define RANDOMLIB_NORMALDISTRIBUTION_HPP 1
#include <cmath> // for std::log
namespace RandomLib {
/**
* \brief Normal deviates
*
* Sample from the normal distribution.
*
* This uses the ratio method; see Knuth, TAOCP, Vol 2, Sec. 3.4.1.C,
* Algorithm R. Unlike the Box-Muller method which generates two normal
* deviates at a time, this method generates just one. This means that this
* class has no state that needs to be saved when checkpointing a
* calculation. Original citation is\n A. J. Kinderman, J. F. Monahan,\n
* Computer Generation of Random Variables Using the Ratio of Uniform
* Deviates,\n ACM TOMS 3, 257--260 (1977).
*
* Improved "quadratic" bounds are given by\n J. L. Leva,\n A Fast Normal
* Random Number Generator,\n ACM TOMS 18, 449--453 and 454--455
* (1992).
*
* The log is evaluated 1.369 times per normal deviate with no bounds, 0.232
* times with Knuth's bounds, and 0.012 times with the quadratic bounds.
* Time is approx 0.3 us per deviate (1GHz machine, optimized, RealType =
* float).
*
* Example
* \code
* #include <RandomLib/NormalDistribution.hpp>
*
* RandomLib::Random r;
* std::cout << "Seed set to " << r.SeedString() << "\n";
* RandomLib::NormalDistribution<double> normdist;
* std::cout << "Select from normal distribution:";
* for (size_t i = 0; i < 10; ++i)
* std::cout << " " << normdist(r);
* std::cout << "\n";
* \endcode
*
* @tparam RealType the real type of the results (default double).
**********************************************************************/
template<typename RealType = double> class NormalDistribution {
public:
/**
* The type returned by NormalDistribution::operator()(Random&)
**********************************************************************/
typedef RealType result_type;
/**
* Return a sample of type RealType from the normal distribution with mean
* &mu; and standard deviation &sigma;.
*
* For &mu; = 0 and &sigma; = 1 (the defaults), the distribution is
* symmetric about zero and is nonzero. The maximum result is less than 2
* sqrt(log(2) \e p) where \e p is the precision of real type RealType.
* The minimum positive value is approximately 1/2<sup><i>p</i>+1</sup>.
* Here \e p is the precision of real type RealType.
*
* @tparam Random the type of RandomCanonical generator.
* @param[in,out] r the RandomCanonical generator.
* @param[in] mu the mean value of the normal distribution (default 0).
* @param[in] sigma the standard deviation of the normal distribution
* (default 1).
* @return the random sample.
**********************************************************************/
template<class Random>
RealType operator()(Random& r, RealType mu = RealType(0),
RealType sigma = RealType(1)) const throw();
};
template<typename RealType> template<class Random> inline RealType
NormalDistribution<RealType>::operator()(Random& r, RealType mu,
RealType sigma) const throw() {
// N.B. These constants can be regarded as "exact", so that the same number
// of significant figures are used in all versions. (They serve to
// "bracket" the real boundary specified by the log expression.)
const RealType
m = RealType( 1.7156 ), // sqrt(8/e) (rounded up)
s = RealType( 0.449871), // Constants from Leva
t = RealType(-0.386595),
a = RealType( 0.19600 ),
b = RealType( 0.25472 ),
r1 = RealType( 0.27597 ),
r2 = RealType( 0.27846 );
RealType u, v, Q;
do { // This loop is executed 1.369 times on average
// Pick point P = (u, v)
u = r.template FixedU<RealType>(); // Sample u in (0,1]
v = m * r.template FixedS<RealType>(); // Sample v in (-m/2, m/2); avoid 0
// Compute quadratic form Q
const RealType x = u - s;
const RealType y = (v < 0 ? -v : v) - t; // Sun has no long double abs!
Q = x*x + y * (a*y - b*x);
} while ( Q >= r1 && // accept P if Q < r1
( Q > r2 || // reject P if Q > r2
v*v > - 4 * u*u * std::log(u) ) ); // accept P if v^2 <= ...
return mu + sigma * (v / u); // return the slope of P (note u != 0)
}
} // namespace RandomLib
#endif // RANDOMLIB_NORMALDISTRIBUTION_HPP

141
include/RandomLib/Random.hpp Normal file
View File

@ -0,0 +1,141 @@
/**
* \file Random.hpp
* \brief Header for Random, RandomGenerator.
*
* This loads up the header for RandomCanonical, RandomEngine, etc., to
* provide access to random integers of various sizes, random reals with
* various precisions, a random probability, etc.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOM_HPP)
#define RANDOMLIB_RANDOM_HPP 1
#include <RandomLib/Config.h>
#if defined(_MSC_VER)
typedef unsigned uint32_t;
typedef unsigned long long uint64_t;
#else
#include <stdint.h>
#endif
/**
* Use table, Power2::power2, for pow2? This isn't necessary with g++ 4.0
* because calls to std::pow are optimized. g++ 4.1 seems to have lost this
* capability though! And it's back in g++ 4.4. So, for simplicity, assume
* that all "current" versions of g++ perform the optimization.
**********************************************************************/
#if !defined(RANDOMLIB_POWERTABLE)
#if defined(__GNUC__)
#define RANDOMLIB_POWERTABLE 0
#else
// otherwise use a lookup table
#define RANDOMLIB_POWERTABLE 1
#endif
#endif
#if !HAVE_LONG_DOUBLE || defined(_MSC_VER)
#define RANDOMLIB_LONGDOUBLEPREC 53
#elif defined(__sparc)
#define RANDOMLIB_LONGDOUBLEPREC 113
#else
/**
* The precision of long doubles, used for sizing Power2::power2. 64 on
* Linux/Intel, 106 on MaxOS/PowerPC
**********************************************************************/
#define RANDOMLIB_LONGDOUBLEPREC __LDBL_MANT_DIG__
#endif
/**
* A compile-time assert. Use C++11 static_assert, if available.
**********************************************************************/
#if !defined(STATIC_ASSERT)
# if __cplusplus >= 201103
# define STATIC_ASSERT static_assert
# elif defined(__GXX_EXPERIMENTAL_CXX0X__)
# define STATIC_ASSERT static_assert
# elif defined(_MSC_VER) && _MSC_VER >= 1600
# define STATIC_ASSERT static_assert
# else
# define STATIC_ASSERT(cond,reason) \
{ enum{ STATIC_ASSERT_ENUM = 1/int(cond) }; }
# endif
#endif
/**
* Are denormalized reals of type RealType supported?
**********************************************************************/
#define RANDOMLIB_HASDENORM(RealType) 1
#if defined(_MSC_VER) && defined(RANDOMLIB_SHARED_LIB) && RANDOMLIB_SHARED_LIB
# if RANDOMLIB_SHARED_LIB > 1
# error RANDOMLIB_SHARED_LIB must be 0 or 1
# elif defined(RandomLib_EXPORTS)
# define RANDOMLIB_EXPORT __declspec(dllexport)
# else
# define RANDOMLIB_EXPORT __declspec(dllimport)
# endif
#else
# define RANDOMLIB_EXPORT
#endif
#include <stdexcept>
/**
* \brief Namespace for %RandomLib
*
* All of %RandomLib is defined within the RandomLib namespace. In addtiion
* all the header files are included via %RandomLib/filename. This minimizes
* the likelihood of conflicts with other packages.
**********************************************************************/
namespace RandomLib {
/**
* \brief Exception handling for %RandomLib
*
* A class to handle exceptions. It's derived from std::runtime_error so it
* can be caught by the usual catch clauses.
**********************************************************************/
class RandomErr : public std::runtime_error {
public:
/**
* Constructor
*
* @param[in] msg a string message, which is accessible in the catch
* clause, via what().
**********************************************************************/
RandomErr(const std::string& msg) : std::runtime_error(msg) {}
};
} // namespace RandomLib
#include <RandomLib/RandomCanonical.hpp>
#if !defined(RANDOMLIB_BUILDING_LIBRARY)
namespace RandomLib {
#if !defined(RANDOMLIB_DEFAULT_GENERATOR)
#define RANDOMLIB_DEFAULT_GENERATOR SRandomGenerator32
#endif
/**
* Point Random to one of a specific MT19937 generators.
**********************************************************************/
typedef RANDOMLIB_DEFAULT_GENERATOR RandomGenerator;
/**
* Hook Random to RandomGenerator
**********************************************************************/
typedef RandomCanonical<RandomGenerator> Random;
} // namespace RandomLib
#endif // !defined(RANDOMLIB_BUILDING_LIBRARY)
#endif // RANDOMLIB_RANDOM_HPP

384
include/RandomLib/RandomAlgorithm.hpp Normal file
View File

@ -0,0 +1,384 @@
/**
* \file RandomAlgorithm.hpp
* \brief Header for MT19937 and SFMT19937.
*
* This provides an interface to the Mersenne Twister
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
* MT19937</a> and SIMD oriented Fast Mersenne Twister
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html">
* SFMT19937</a> random number engines.
*
* Interface routines written by Charles Karney <charles@karney.com> and
* licensed under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMALGORITHM_HPP)
#define RANDOMLIB_RANDOMALGORITHM_HPP 1
#include <RandomLib/RandomType.hpp>
#include <stdexcept>
#include <string>
#if defined(HAVE_SSE2) && HAVE_SSE2
#include <emmintrin.h>
#endif
#if (defined(HAVE_SSE2) && HAVE_SSE2) && (defined(HAVE_ALTIVEC) && HAVE_ALTIVEC)
#error "HAVE_SSE2 and HAVE_ALTIVEC should not both be defined"
#endif
#if defined(_MSC_VER)
// Squelch warnings about casts truncating constants
# pragma warning (push)
# pragma warning (disable: 4310)
#endif
namespace RandomLib {
/**
* \brief The %MT19937 random number engine.
*
* This provides an interface to Mersenne Twister random number engine,
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
* MT19937</a>. See\n Makoto Matsumoto and Takuji Nishimura,\n Mersenne
* Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number
* Generator,\n ACM TOMACS 8, 3--30 (1998)
*
* This is adapted from the 32-bit and 64-bit C versions available at
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html and
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html
*
* The template argument give the type \e RandomType of the "natural" result.
* This incorporates the bit width and the C++ type of the result. Although
* the two versions of MT19937 produce different sequences, the
* implementations here are portable across 32-bit and 64-bit architectures.
*
* The class chiefly supplies the method for advancing the state by
* Transition.
*
* @tparam RandomType the type of the results, either Random_u32 or
* Random_u64.
*
* Interface routines written by Charles Karney <charles@karney.com> and
* licensed under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
template<class RandomType> class RANDOMLIB_EXPORT MT19937 {
public:
/**
* The result RandomType
**********************************************************************/
typedef RandomType engine_t;
/**
* The internal numeric type for MT19337::Transition
**********************************************************************/
typedef typename engine_t::type internal_type;
private:
/**
* The unsigned type of engine_t
**********************************************************************/
typedef typename engine_t::type engine_type;
/**
* The width of the engine_t
**********************************************************************/
static const unsigned width = engine_t::width;
enum {
/**
* The Mersenne prime is 2<sup><i>P</i></sup> &minus; 1
**********************************************************************/
P = 19937,
/**
* The short lag for MT19937
**********************************************************************/
M = width == 32 ? 397 : 156,
/**
* The number of ignored bits in the first word of the state
**********************************************************************/
R = ((P + width - 1)/width) * width - P
};
static const engine_type mask = engine_t::mask;
/**
* Magic matrix for MT19937
**********************************************************************/
static const engine_type magic =
width == 32 ? 0x9908b0dfULL : 0xb5026f5aa96619e9ULL;
/**
* Mask for top \e width &minus; \e R bits of a word
**********************************************************************/
static const engine_type upper = mask << R & mask;
/**
* Mask for low \e R bits of a <i>width</i>-bit word
**********************************************************************/
static const engine_type lower = ~upper & mask;
public:
/**
* A version number "EnMT" or "EnMU" to ensure safety of Save/Load. This
* needs to be unique across RandomAlgorithms.
**********************************************************************/
static const unsigned version = 0x456e4d54UL + (engine_t::width/32 - 1);
enum {
/**
* The size of the state. This is the long lag for MT19937.
**********************************************************************/
N = (P + width - 1)/width
};
/**
* Advance state by \e count batches. For speed all \e N words of state
* are advanced together. If \e count is negative, the state is stepped
* backwards. This is the meat of the MT19937 engine.
*
* @param[in] count how many batches to advance.
* @param[in,out] statev the internal state of the random number generator.
**********************************************************************/
static void Transition(long long count, internal_type statev[]) throw();
/**
* Manipulate a word of the state prior to output.
*
* @param[in] y a word of the state.
* @return the result.
**********************************************************************/
static engine_type Generate(engine_type y) throw();
/**
* Convert an arbitrary state into a legal one. This consists of (a)
* turning on one bit if the state is all zero and (b) making 31 bits of
* the state consistent with the other 19937 bits.
*
* @param[in,out] state the state of the generator.
**********************************************************************/
static void NormalizeState(engine_type state[]) throw();
/**
* Check that the state is legal, throwing an exception if it is not. At
* the same time, accumulate a checksum of the state.
*
* @param[in] state the state of the generator.
* @param[in,out] check an accumulated checksum.
**********************************************************************/
static void CheckState(const engine_type state[], Random_u32::type& check);
/**
* Return the name of the engine
*
* @return the name.
**********************************************************************/
static std::string Name() throw() {
return "MT19937<Random_u" + std::string(width == 32 ? "32" : "64") + ">";
}
};
/// \cond SKIP
template<>
inline Random_u32::type MT19937<Random_u32>::Generate(engine_type y) throw() {
y ^= y >> 11;
y ^= y << 7 & engine_type(0x9d2c5680UL);
y ^= y << 15 & engine_type(0xefc60000UL);
y ^= y >> 18;
return y;
}
template<>
inline Random_u64::type MT19937<Random_u64>::Generate(engine_type y) throw() {
// Specific tempering instantiation for width = 64 given in
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html
y ^= y >> 29 & engine_type(0x5555555555555555ULL);
y ^= y << 17 & engine_type(0x71d67fffeda60000ULL);
y ^= y << 37 & engine_type(0xfff7eee000000000ULL);
y ^= y >> 43;
return y;
}
/// \endcond
/**
* \brief The SFMT random number engine.
*
* This provides an implementation of the SIMD-oriented Fast Mersenne Twister
* random number engine,
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html">
* SFMT</a>. See\n Mutsuo Saito,\n An Application of Finite Field: Design
* and Implementation of 128-bit Instruction-Based Fast Pseudorandom Number
* Generator,\n Master's Thesis, Dept. of Math., Hiroshima University
* (Feb. 2007).\n
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/M062821.pdf
* Mutsuo Saito and Makoto Matsumoto,\n
* SIMD-oriented Fast Mersenne Twister: a 128-bit Pseudorandom Number
* Generator,\n accepted in the proceedings of MCQMC2006\n
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/sfmt.pdf
*
* The template argument gives the type \e RandomType of the "natural"
* result. This incorporates the bit width and the C++ type of the result.
* The 32-bit and 64-bit versions of SFMT19937 produce the same sequences and
* the differing only in whether how the state is represented. The
* implementation includes a version using 128-bit SSE2 instructions. On
* machines without these instructions, portable implementations using
* traditional operations are provided. With the same starting seed,
* SRandom32::Ran64() and SRandom64::Ran64() produces the same sequences.
* Similarly SRandom64::Ran32() produces every other member of the sequence
* produced by SRandom32::Ran32().
*
* The class chiefly supplies the method for advancing the state by
* Transition.
*
* @tparam RandomType the type of the results, either Random_u32 or
* Random_u64.
*
* Written by Charles Karney <charles@karney.com> and licensed under the
* MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
template<class RandomType> class RANDOMLIB_EXPORT SFMT19937 {
public:
/**
* The result RandomType
**********************************************************************/
typedef RandomType engine_t;
#if defined(HAVE_SSE2) && HAVE_SSE2
typedef __m128i internal_type;
#elif defined(HAVE_ALTIVEC) && HAVE_ALTIVEC
typedef vector unsigned internal_type;
#else
/**
* The internal numeric type for SFMT19337::Transition
**********************************************************************/
typedef typename engine_t::type internal_type;
#endif
private:
/**
* The unsigned type of engine_t
**********************************************************************/
typedef typename engine_t::type engine_type;
/**
* The width of the engine_t
**********************************************************************/
static const unsigned width = engine_t::width;
enum {
/**
* The Mersenne prime is 2<sup><i>P</i></sup> &minus; 1
**********************************************************************/
P = 19937,
/**
* The long lag for SFMT19937 in units of 128-bit words
**********************************************************************/
N128 = (P + 128 - 1)/128,
/**
* How many width words per 128-bit word.
**********************************************************************/
R = 128 / width,
/**
* The short lag for SFMT19937 in units of 128-bit words
**********************************************************************/
M128 = 122,
/**
* The short lag for SFMT19937
**********************************************************************/
M = M128 * R
};
#if (defined(HAVE_SSE2) && HAVE_SSE2) || (defined(HAVE_ALTIVEC) && HAVE_ALTIVEC)
static const Random_u32::type magic0 = 0x1fffefUL;
static const Random_u32::type magic1 = 0x1ecb7fUL;
static const Random_u32::type magic2 = 0x1affffUL;
static const Random_u32::type magic3 = 0x1ffff6UL;
#else
/**
* Magic matrix for SFMT19937. Only the low 21 (= 32 &minus; 11) bits need
* to be set. (11 is the right shift applied to the words before masking.
**********************************************************************/
static const engine_type
magic0 = width == 32 ? 0x1fffefULL : 0x1ecb7f001fffefULL;
static const engine_type
magic1 = width == 32 ? 0x1ecb7fULL : 0x1ffff6001affffULL;
static const engine_type
magic2 = width == 32 ? 0x1affffULL : 0ULL;
static const engine_type
magic3 = width == 32 ? 0x1ffff6ULL : 0ULL;
#endif
/**
* Mask for simulating u32 << 18 with 64-bit words
**********************************************************************/
static const engine_type mask18 = engine_type(0xfffc0000fffc0000ULL);
/**
* Magic constants needed by "period certification"
**********************************************************************/
static const engine_type PARITY0 = 1U;
static const engine_type PARITY1 = width == 32 ? 0U : 0x13c9e68400000000ULL;
static const engine_type PARITY2 = 0U;
static const engine_type PARITY3 = width == 32 ? 0x13c9e684UL : 0U;
/**
* Least significant bit of PARITY
**********************************************************************/
static const unsigned PARITY_LSB = 0;
static const engine_type mask = engine_t::mask;
public:
/**
* A version number "EnSM" or "EnSN" to ensure safety of Save/Load. This
* needs to be unique across RandomAlgorithms.
**********************************************************************/
static const unsigned version = 0x456e534dUL + (engine_t::width/32 - 1);
enum {
/**
* The size of the state. The long lag for SFMT19937
**********************************************************************/
N = N128 * R
};
/**
* Advance state by \e count batches. For speed all \e N words of state
* are advanced together. If \e count is negative, the state is stepped
* backwards. This is the meat of the SFMT19937 engine.
*
* @param[in] count how many batches to advance.
* @param[in,out] statev the internal state of the random number generator.
**********************************************************************/
static void Transition(long long count, internal_type statev[])
throw();
/**
* Manipulate a word of the state prior to output. This is a no-op for
* SFMT19937.
*
* @param[in] y a word of the state.
* @return the result.
**********************************************************************/
static engine_type Generate(engine_type y) throw() { return y; }
/**
* Convert an arbitrary state into a legal one. This consists a "period
* certification to ensure that the period of the generator is at least
* 2<sup><i>P</i></sup> &minus; 1.
*
* @param[in,out] state the state of the generator.
**********************************************************************/
static void NormalizeState(engine_type state[]) throw();
/**
* Check that the state is legal, throwing an exception if it is not. This
* merely verifies that the state is not all zero. At the same time,
* accumulate a checksum of the state.
*
* @param[in] state the state of the generator.
* @param[in,out] check an accumulated checksum.
**********************************************************************/
static void CheckState(const engine_type state[], Random_u32::type& check);
/**
* Return the name of the engine
*
* @return the name.
**********************************************************************/
static std::string Name() throw() {
return "SFMT19937<Random_u" +
std::string(width == 32 ? "32" : "64") + ">";
}
};
} // namespace RandomLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // RANDOMLIB_RANDOMALGORITHM_HPP

1301
include/RandomLib/RandomCanonical.hpp Normal file
View File

File diff suppressed because it is too large Load Diff

640
include/RandomLib/RandomEngine.hpp Normal file
View File

@ -0,0 +1,640 @@
/**
* \file RandomEngine.hpp
* \brief Header for RandomEngine.
*
* Copyright (c) Charles Karney (2006-2012) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMENGINE_HPP)
#define RANDOMLIB_RANDOMENGINE_HPP 1
#include <RandomLib/RandomSeed.hpp>
#include <RandomLib/RandomAlgorithm.hpp>
#include <RandomLib/RandomMixer.hpp>
#include <limits>
#include <string>
#include <algorithm>
#if defined(HAVE_SSE2) && HAVE_SSE2 && defined(_MSC_VER) && !defined(_WIN64)
#include <new>
#endif
#if !defined(RANDOMLIB_BUILDING_LIBRARY) && \
defined(HAVE_BOOST_SERIALIZATION) && HAVE_BOOST_SERIALIZATION
#include <boost/serialization/nvp.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/vector.hpp>
#endif
namespace RandomLib {
/**
* \brief Uniform random number generator.
*
* This implements a generic random number generator. Such a generator
* requires two data holders RandomSeed, to hold the seed, and RandomEngine,
* to hold the state. In addition we need two piece of machinery, a "Mixer"
* to convert the seed into an initial state and an "Algorithm" to advance the
* state.
*
* @tparam Algorithm the random number algorithm.
* @tparam Mixer the way seeds are turned into state.
*
* RandomSeed is responsible for setting and reporting the seed.
*
* Mixer has no state and implements only static methods. It needs to have
* the following public interface
* - typedef mixer_t: a RandomType giving the output type
* - unsigned version: an identifying version number
* - static std::string Name(): an identifying name for the mixer
* - static method SeedToState: converts a seed into n words of state.
*
* Algorithm has no state and implements only static methods. It needs to
* have the following public interface
* - typedef engine_t: a RandomType giving the output type
* - typedef internal_type: a integer type used by Transition. This is
* usually the same as engine_t::type. However it allows the use of
* vector instructions on some platforms. We require that engine_t::type
* and internal_type line up properly in a union so that there is no need
* to convert the data explicitly between internal_type and
* engine_t::type.
* - unsigned version: an identifying version number
* - static std::string Name(): an identifying name for the mixer
* - enum N: the size of the state in units of engine_t.
* - static method Transition: steps the generator forwards or backwards.
* - static method Generate: tempers the state immediately prior to output
* - static method NormalizeState: force the initial state (the result of
* the Mixer) into a legal state.
* - static method CheckState accumulates the checksum for the state into
* check. In addition it throws an exception if the state is bad.
*
* RandomEngine is the glue that holds everything together. It repacks
* the mixer_t data from Mixer into engine_t if necessary. It deals with
* delivering individual random results, stepping the state forwards and
* backwards, leapfrogging the generator, I/O of the generator, etc.
*
* Written by Charles Karney <charles@karney.com> and licensed under the
* MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
template<class Algorithm, class Mixer>
class RANDOMLIB_EXPORT RandomEngine : public RandomSeed {
private:
/**
* The result RandomType (carried over from the \e Algorithm).
**********************************************************************/
typedef typename Algorithm::engine_t result_t;
/**
* The RandomType used by the \e Mixer.
**********************************************************************/
typedef typename Mixer::mixer_t mixer_t;
/**
* The internal_type used by the Algorithm::Transition().
**********************************************************************/
typedef typename Algorithm::internal_type engine_type;
public:
/**
* The number of random bits produced by Ran().
**********************************************************************/
enum {
width = result_t::width
};
/**
* A type large enough to hold \e width bits. This is used for the
* internal state of the generator and the result returned by Ran().
**********************************************************************/
typedef typename result_t::type result_type;
/**
* The minimum result returned by Ran() = 0.
**********************************************************************/
static const result_type min = result_t::min;
/**
* The maximum result returned by Ran() = 2<sup><i>w</i></sup> &minus; 1.
**********************************************************************/
static const result_type max = result_t::max;
protected:
/**
* The mask for the result_t.
**********************************************************************/
static const result_type mask = result_t::mask;
private:
/**
* A version number "RandLib0" to ensure safety of Save/Load. The first 7
* bytes can be regarded as a "signature" and the 8th byte a version
* number.
**********************************************************************/
static const u64::type version = 0x52616e644c696230ULL; // 'RandLib0'
/**
* Marker for uninitialized object
**********************************************************************/
static const unsigned UNINIT = 0xffffffffU;
enum {
/**
* The size of the state in units of result_type
**********************************************************************/
N = Algorithm::N,
/**
* The size of the state in units of mixer_t::type
**********************************************************************/
NU = (N * width + mixer_t::width - 1) / mixer_t::width,
/**
* The size of the state in units of engine_type.
**********************************************************************/
NV = N * sizeof(result_type) / sizeof(engine_type)
};
/**
* \brief Union for the state.
*
* A union to hold the state in the result_type, mixer_t::type, and
* engine_type representations.
**********************************************************************/
union {
/**
* the result_type representation returned by Ran()
**********************************************************************/
result_type _state[N];
/**
* the mixer_t::type representation returned by Mixer::SeedToState.
**********************************************************************/
typename mixer_t::type _stateu[NU];
/**
* the engine_type representation returned by Algorithm::Transition.
**********************************************************************/
engine_type _statev[NV];
};
/**
* The index for the next random value
**********************************************************************/
unsigned _ptr;
/**
* How many times has Transition() been called
**********************************************************************/
long long _rounds;
/**
* Stride for leapfrogging
**********************************************************************/
unsigned _stride;
public:
/**
* \name Constructors
**********************************************************************/
///@{
/**
* Initialize from a vector. Only the low \e 32 bits of each element are
* used.
*
* @tparam IntType the integral type of the elements of the vector.
* @param[in] v the vector of elements.
**********************************************************************/
template<typename IntType>
explicit RandomEngine(const std::vector<IntType>& v) { Reseed(v); }
/**
* Initialize from a pair of iterators setting seed to [\e a, \e b). The
* iterator must produce results which can be converted into seed_type.
* Only the low \e 32 bits of each element are used.
*
* @tparam InputIterator the type of the iterator.
* @param[in] a the beginning iterator.
* @param[in] b the ending iterator.
**********************************************************************/
template<typename InputIterator>
RandomEngine(InputIterator a, InputIterator b) { Reseed(a, b); }
/**
* Initialize with seed [\e n]. Only the low \e width bits of \e n are
* used.
*
* @param[in] n the new seed to use.
**********************************************************************/
explicit RandomEngine(seed_type n) { Reseed(n); }
/**
* Initialize with seed []. This can be followed by a call to Reseed() to
* select a unique seed.
**********************************************************************/
RandomEngine() { unsigned long s[1]; Reseed(s, s); }
/**
* Initialize from a string. See Reseed(const std::string& s)
*
* @param[in] s the string to be decoded into a seed.
**********************************************************************/
explicit RandomEngine(const std::string& s) { Reseed(s); }
///@}
/**
* \name Functions for returning random data
**********************************************************************/
///@{
/**
* Return \e width bits of randomness. This is the natural unit of random
* data produced random number generator.
*
* @return the next random number of width \e width.
**********************************************************************/
result_type Ran() throw() {
if (_ptr >= N)
Next();
result_type y = _state[_ptr];
_ptr += _stride;
return Algorithm::Generate(y);
}
/**
* Return 32 bits of randomness.
*
* @return a 32-bit random number.
**********************************************************************/
u32::type Ran32() throw() {
// return width > 32 ? u32::cast(Ran()) : Ran();
return u32::cast(Ran());
}
/**
* Return 64 bits of randomness.
*
* @return a 64-bit random number.
**********************************************************************/
u64::type Ran64() throw() {
const u64::type x = Ran();
return width > 32 ? x : u64::cast(Ran()) << (64 - width) | x;
}
/**
* Return \e width bits of randomness. Result is in [0,
* 2<sup><i>w</i></sup>). (This just calls Ran().)
*
* @return the next random number of width \e width.
**********************************************************************/
result_type operator()() throw() { return Ran(); }
///@}
#if defined(HAVE_SSE2) && HAVE_SSE2 && defined(_MSC_VER) && !defined(_WIN64)
/**
* new operator with alignment (needed for Visual Studio)
**********************************************************************/
void* operator new(size_t n) {
void* p = _aligned_malloc(n, __alignof(RandomEngine));
if (p == 0) throw std::bad_alloc();
return p;
}
/**
* delete operator with alignment (needed for Visual Studio)
**********************************************************************/
void operator delete(void* p) { _aligned_free(p); }
/**
* new[] operator with alignment (needed for Visual Studio)
**********************************************************************/
void* operator new[](size_t n) {
void* p = _aligned_malloc(n, __alignof(RandomEngine));
if (p == 0) throw std::bad_alloc();
return p;
}
/**
* delete[] operator with alignment (needed for Visual Studio)
**********************************************************************/
void operator delete[](void* p) { _aligned_free(p); }
#endif
/**
* \name Comparing Random objects
**********************************************************************/
///@{
/**
* Test equality of two Random objects. This test that the seeds match and
* that they have produced the same number of random numbers.
*
* @param[in] r the RandomEngine object to compare.
* @return true if the RandomEngine objects produce the same results.
**********************************************************************/
bool operator==(const RandomEngine& r) const throw()
// Ensure that the two Random objects behave the same way. Note however
// that the internal states may still be different, e.g., the following all
// result in Random objects which are == (with Count() == 0) but which all
// have different internal states:
//
// Random r(0); _ptr == UNINIT
// r.StepCount( 1); r.StepCount(-1); _ptr == 0, _rounds == 0
// r.StepCount(-1); r.StepCount( 1); _ptr == N, _rounds == -1
{ return Count() == r.Count() && _seed == r._seed &&
_stride == r._stride; }
/**
* Test inequality of two Random objects. See Random::operator==
*
* @param[in] r the RandomEngine object to compare.
* @return true if the RandomEngine objects produce different results.
**********************************************************************/
bool operator!=(const RandomEngine& r) const throw()
{ return !operator==(r); }
///@}
/**
* \name Interchanging Random objects
**********************************************************************/
///@{
/**
* Swap with another Random object.
*
* @param[in,out] t the RandomEngine object to swap with.
**********************************************************************/
void swap(RandomEngine& t) throw() {
_seed.swap(t._seed);
std::swap(_ptr, t._ptr);
std::swap(_stride, t._stride);
std::swap(_rounds, t._rounds);
std::swap_ranges(_state, _state + N, t._state);
}
///@}
/**
* \name Writing to and reading from a stream
**********************************************************************/
///@{
/**
* Save the state of the Random object to an output stream. Format is a
* sequence of unsigned 32-bit integers written either in decimal (\e bin
* false, text format) or in network order with most significant byte first
* (\e bin true, binary format). Data consists of:
*
* - RandomLib magic string + version (2 words)
* - Algorithm version (1 word)
* - Mixer version (1 word)
* - _seed.size() (1 word)
* - _seed data (_seed.size() words)
* - _ptr (1 word)
* - _stride (1 word)
* - if _ptr != UNINIT, _rounds (2 words)
* - if _ptr != UNINIT, _state (N words or 2 N words)
* - checksum
*
* Shortest possible saved result consists of 8 words. This corresponds to
* RandomSeed() = [] and Count() = 0.
*
* @param[in,out] os the output stream.
* @param[in] bin if true (the default) save in binary mode.
**********************************************************************/
void Save(std::ostream& os, bool bin = true) const;
/**
* Restore the state of the Random object from an input stream. If \e bin,
* read in binary, else use text format. See documentation of
* RandomEngine::Save for the format. Include error checking on data to
* make sure the input has not been corrupted. If an error occurs while
* reading, the Random object is unchanged.
*
* @param[in,out] is the input stream.
* @param[in] bin if true (the default) load in binary mode.
* @exception RandomErr if the state read from \e is is illegal.
**********************************************************************/
void Load(std::istream& is, bool bin = true) {
// Read state into temporary so as not to change object on error.
RandomEngine t(is, bin);
_seed.reserve(t._seed.size());
*this = t;
}
///@}
/**
* \name Basic I/O
**********************************************************************/
///@{
/**
* Write the state of a generator to stream \e os as text
*
* @param[in,out] os the output stream.
* @param[in] r the RandomEngine object to be saved.
**********************************************************************/
friend std::ostream& operator<<(std::ostream& os, const RandomEngine& r) {
r.Save(os, false);
return os;
}
/**
* Read the state of a generator from stream \e is as text
*
* @param[in,out] is the output stream.
* @param[in] r the RandomEngine object to be loaded.
* @exception RandomErr if the state read from \e is is illegal.
**********************************************************************/
friend std::istream& operator>>(std::istream& is, RandomEngine& r) {
r.Load(is, false);
return is;
}
///@}
/**
* \name Examining and advancing the Random generator
**********************************************************************/
///@{
/**
* Return the number of random numbers used. This needs to return a long
* long result since it can reasonably exceed 2<sup>31</sup>. (On a 1GHz
* machine, it takes about a minute to produce 2<sup>32</sup> random
* numbers.) More precisely this is the (zero-based) index of the next
* random number to be produced. (This distinction is important when
* leapfrogging is in effect.)
*
* @return the count of random numbers used.
**********************************************************************/
long long Count() const throw()
{ return _ptr == UNINIT ? 0 : _rounds * N + _ptr; }
/**
* Step the generator forwards or backwards so that the value returned
* by Count() is \e n
*
* @param[in] n the new count.
**********************************************************************/
void SetCount(long long n) throw() { StepCount(n - Count()); }
/**
* Step the generator forward \e n steps. \e n can be negative.
*
* @param[in] n how much to step the generator forward.
**********************************************************************/
void StepCount(long long n) throw();
/**
* Resets the sequence. Equivalent to SetCount(0), but works by
* reinitializing the Random object from its seed, rather than by stepping
* the sequence backwards. In addition, this undoes leapfrogging.
**********************************************************************/
void Reset() throw() { _ptr = UNINIT; _stride = 1; }
///@}
/**
* \name Leapfrogging
**********************************************************************/
///@{
/**
* Set leapfrogging stride to a positive number \e n and increment Count()
* by \e k < \e n. If the current Count() is \e i, then normally the next
* 3 random numbers would have (zero-based) indices \e i, \e i + 1, \e i +
* 2, and the new Count() is \e i + 2. However, after SetStride(\e n, \e
* k) the next 3 random numbers have indices \e i + \e k, \e i + \e k + \e
* n, \e i + \e k + 2\e n, and the new Count() is \e i + \e k + 3\e n.
* With leapfrogging in effect, the time to produce raw random numbers is
* roughly proportional to 1 + (\e n &minus; 1)/3. Reseed(...) and Reset()
* both reset the stride back to 1. See \ref parallel for a description of
* how to use this facility.
*
* @param[in] n the stride (default 1).
* @param[in] k the initial increment (default 0).
* @exception RandomErr if \e n is 0 or too large or if \e k is not less
* than \e n.
**********************************************************************/
void SetStride(unsigned n = 1, unsigned k = 0) {
// Limit stride to UNINIT/2. This catches negative numbers that have
// been cast into unsigned. In reality the stride should be no more than
// 10-100.
if (n == 0 || n > UNINIT/2)
throw RandomErr("RandomEngine: Invalid stride");
if (k >= n)
throw RandomErr("RandomEngine: Invalid offset");
_stride = n;
StepCount(k);
}
/**
* Return leapfrogging stride.
*
* @return the stride.
**********************************************************************/
unsigned GetStride() const throw() { return _stride; }
///@}
/**
* Tests basic engine.
*
* @exception RandomErr if any of the tests fail.
**********************************************************************/
static void SelfTest();
/**
* Return the name of the generator. This incorporates the names of the \e
* Algorithm and \e Mixer.
*
* @return the name of the generator.
**********************************************************************/
static std::string Name() {
return "RandomEngine<" + Algorithm::Name() + "," + Mixer::Name() + ">";
}
private:
/**
* Compute initial state from seed
**********************************************************************/
void Init() throw();
/**
* The interface to Transition used by Ran().
**********************************************************************/
void Next() throw() {
if (_ptr == UNINIT)
Init();
_rounds += _ptr/N;
Algorithm::Transition(_ptr/N, _statev);
_ptr %= N;
}
u32::type Check(u64::type v, u32::type e, u32::type m) const;
static result_type SelfTestResult(unsigned) throw() { return 0; }
/**
* Read from an input stream. Potentially corrupts object. This private
* constructor is used by RandomEngine::Load so that it can avoid
* corrupting its state on bad input.
**********************************************************************/
explicit RandomEngine(std::istream& is, bool bin);
#if !defined(RANDOMLIB_BUILDING_LIBRARY) && \
defined(HAVE_BOOST_SERIALIZATION) && HAVE_BOOST_SERIALIZATION
friend class boost::serialization::access;
/**
* Save to a boost archive. Boost versioning isn't very robust. (It
* allows a RandomGenerator32 to be read back in as a RandomGenerator64.
* It doesn't interact well with templates.) So we do our own versioning
* and supplement this with a checksum.
**********************************************************************/
template<class Archive> void save(Archive& ar, const unsigned int) const {
u64::type _version = version;
u32::type _eversion = Algorithm::version,
_mversion = Mixer::version,
_checksum = Check(_version, _eversion, _mversion);
ar & boost::serialization::make_nvp("version" , _version )
& boost::serialization::make_nvp("eversion", _eversion)
& boost::serialization::make_nvp("mversion", _mversion)
& boost::serialization::make_nvp("seed" , _seed )
& boost::serialization::make_nvp("ptr" , _ptr )
& boost::serialization::make_nvp("stride" , _stride );
if (_ptr != UNINIT)
ar & boost::serialization::make_nvp("rounds", _rounds )
& boost::serialization::make_nvp("state" , _state );
ar & boost::serialization::make_nvp("checksum", _checksum);
}
/**
* Load from a boost archive. Do this safely so that the current object is
* not corrupted if the archive is bogus.
**********************************************************************/
template<class Archive> void load(Archive& ar, const unsigned int) {
u64::type _version;
u32::type _eversion, _mversion, _checksum;
ar & boost::serialization::make_nvp("version" , _version )
& boost::serialization::make_nvp("eversion", _eversion )
& boost::serialization::make_nvp("mversion", _mversion );
RandomEngine<Algorithm, Mixer> t(std::vector<seed_type>(0));
ar & boost::serialization::make_nvp("seed" , t._seed )
& boost::serialization::make_nvp("ptr" , t._ptr )
& boost::serialization::make_nvp("stride" , t._stride );
if (t._ptr != UNINIT)
ar & boost::serialization::make_nvp("rounds", t._rounds )
& boost::serialization::make_nvp("state" , t._state );
ar & boost::serialization::make_nvp("checksum", _checksum );
if (t.Check(_version, _eversion, _mversion) != _checksum)
throw RandomErr("RandomEngine: Checksum failure");
_seed.reserve(t._seed.size());
*this = t;
}
/**
* Glue the boost save and load functionality together---a bit of boost
* magic.
**********************************************************************/
template<class Archive>
void serialize(Archive &ar, const unsigned int file_version)
{ boost::serialization::split_member(ar, *this, file_version); }
#endif // HAVE_BOOST_SERIALIZATION
};
typedef RandomEngine<MT19937 <Random_u32>, MixerSFMT> MRandomGenerator32;
typedef RandomEngine<MT19937 <Random_u64>, MixerSFMT> MRandomGenerator64;
typedef RandomEngine<SFMT19937<Random_u32>, MixerSFMT> SRandomGenerator32;
typedef RandomEngine<SFMT19937<Random_u64>, MixerSFMT> SRandomGenerator64;
} // namespace RandomLib
namespace std {
/**
* Swap two RandomEngines. This is about 3x faster than the default swap.
*
* @tparam Algorithm the algorithm for the RandomEngine.
* @tparam Mixer the mixer for the RandomEngine.
* @param[in,out] r the first RandomEngine to swap.
* @param[in,out] s the second RandomEngine to swap.
**********************************************************************/
template<class Algorithm, class Mixer>
void swap(RandomLib::RandomEngine<Algorithm, Mixer>& r,
RandomLib::RandomEngine<Algorithm, Mixer>& s) throw() {
r.swap(s);
}
} // namespace std
#endif // RANDOMLIB_RANDOMENGINE_HPP

258
include/RandomLib/RandomMixer.hpp Normal file
View File

@ -0,0 +1,258 @@
/**
* \file RandomMixer.hpp
* \brief Header for Mixer classes.
*
* Mixer classes convert a seed vector into a random generator state. An
* important property of this method is that "close" seeds should produce
* "widely separated" states. This allows the seeds to be set is some
* systematic fashion to produce a set of uncorrelated random number
* sequences.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMMIXER_HPP)
#define RANDOMLIB_RANDOMMIXER_HPP 1
#include <vector>
#include <string>
#include <RandomLib/RandomSeed.hpp>
namespace RandomLib {
/**
* \brief The original %MT19937 mixing functionality
*
* This implements the functionality of init_by_array in MT19937
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
* and init_by_array64 in MT19937_64
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
* with the following changes:
* - in the case of an zero-length seed array, behave in the same way if
* MT19937 and MT19937_64 are called without initialization in which case,
* e.g., init_genrand(5489UL) is called. (init_by_array does not allow
* calling with a zero-length seed.)
* - init_by_array64 accepts a seed array of 64-bit unsigned ints. Here with
* seed is an array of 32-bit unsigned ints and these are repacked into
* 64-bit quantities internally using a LSB convention. Thus, to mimic the
* MT19937_64 sample invocation with a seed array {0x12345ULL, 0x23456ULL,
* 0x34567ULL, 0x45678ULL}, MixerMT0<Random_u64>::SeedToState needs to
* be invoked with a seed vector [0x12345UL, 0, 0x23456UL, 0, 0x34567UL, 0,
* 0x45678UL, 0]. (Actually the last 0 is unnecessary.)
*
* The template parameter \e RandomType switches between the 32-bit and
* 64-bit versions.
*
* MixerMT0 is specific to the MT19937 generators and should not be used
* for other generators (e.g., SFMT19937). In addition, MixerMT0 has
* known defects and should only be used to check the operation of the
* MT19937 engines against the original implementation. These defects are
* described in the MixerMT1 which is a modification of MixerMT0
* which corrects these defects. For production use MixerMT1 or,
* preferably, MixerSFMT should be used.
*
* @tparam RandomType the type of the results, either Random_u32 or
* Random_u64.
**********************************************************************/
template<class RandomType> class RANDOMLIB_EXPORT MixerMT0 {
public:
/**
* The RandomType controlling the output of MixerMT0::SeedToState
**********************************************************************/
typedef RandomType mixer_t;
/**
* A version number which should be unique to this RandomMixer. This
* prevents RandomEngine::Load from loading a saved generator with a
* different RandomMixer. Here the version is "MxMT" or "MxMU".
**********************************************************************/
static const unsigned version = 0x4d784d54UL + (mixer_t::width == 64);
private:
/**
* The unsigned type corresponding to mixer_t.
**********************************************************************/
typedef typename mixer_t::type mixer_type;
/**
* The mask for mixer_t.
**********************************************************************/
static const mixer_type mask = mixer_t::mask;
public:
/**
* Mix the seed vector, \e seed, into the state array, \e state, of size \e
* n.
*
* @param[in] seed the input seed vector.
* @param[out] state the generator state.
* @param[in] n the size of the state.
**********************************************************************/
static void SeedToState(const std::vector<RandomSeed::seed_type>& seed,
mixer_type state[], unsigned n) throw();
/**
* Return the name of this class.
*
* @return the name.
**********************************************************************/
static std::string Name() {
return "MixerMT0<Random_u" +
std::string(mixer_t::width == 32 ? "32" : "64") + ">";
}
private:
static const mixer_type a0 = 5489ULL;
static const mixer_type a1 = 19650218ULL;
static const mixer_type
b = mixer_t::width == 32 ? 1812433253ULL : 6364136223846793005ULL;
static const mixer_type
c = mixer_t::width == 32 ? 1664525ULL : 3935559000370003845ULL;
static const mixer_type
d = mixer_t::width == 32 ? 1566083941ULL : 2862933555777941757ULL;
};
/**
* \brief The modified %MT19937 mixing functionality
*
* MixerMT0 has two defects
* - The zeroth word of the state is set to a constant (independent of the
* seed). This is a relatively minor defect which halves the accessible
* state space for MT19937 (but the resulting state space is still huge).
* (Actually, for the 64-bit version, it reduces the accessible states by
* 2<sup>33</sup>. On the other hand the 64-bit has better mixing
* properties.)
* - Close seeds, for example, [1] and [1,0], result in the same state. This
* is a potentially serious flaw which might result is identical random
* number sequences being generated instead of independent sequences.
*
* MixerMT1 fixes these defects in a straightforward manner. The
* resulting algorithm was included in one of the proposals for Random Number
* Generation for C++0X, see Brown, et al.,
* http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2006/n2079.pdf
*
* The template parameter \e RandomType switches between the 32-bit and
* 64-bit versions.
*
* MixerMT1 still has a weakness in that it doesn't thoroughly mix the
* state. This is illustrated by an example given to me by Makoto Matsumoto:
* Consider a seed of length \e N and suppose we consider all \e
* W<sup><i>N</i>/2</sup> values for the first half of the seed (here \e W =
* 2<sup><i>width</i></sup>). MixerMT1 has a bottleneck in the way that
* the state is initialized which results in the second half of the state
* only taking on \e W<sup>2</sup> possible values. MixerSFMT mixes the
* seed into the state much more thoroughly.
*
* @tparam RandomType the type of the results, either Random_u32 or
* Random_u64.
**********************************************************************/
template<class RandomType> class RANDOMLIB_EXPORT MixerMT1 {
public:
/**
* The RandomType controlling the output of MixerMT1::SeedToState
**********************************************************************/
typedef RandomType mixer_t;
/**
* A version number which should be unique to this RandomMixer. This
* prevents RandomEngine::Load from loading a saved generator with a
* different RandomMixer. Here the version is "MxMV" or "MxMW".
**********************************************************************/
static const unsigned version = 0x4d784d56UL + (mixer_t::width == 64);
private:
/**
* The unsigned type corresponding to mixer_t.
**********************************************************************/
typedef typename mixer_t::type mixer_type;
/**
* The mask for mixer_t.
**********************************************************************/
static const mixer_type mask = mixer_t::mask;
public:
/**
* Mix the seed vector, \e seed, into the state array, \e state, of size \e
* n.
*
* @param[in] seed the input seed vector.
* @param[out] state the generator state.
* @param[in] n the size of the state.
**********************************************************************/
static void SeedToState(const std::vector<RandomSeed::seed_type>& seed,
mixer_type state[], unsigned n) throw();
/**
* Return the name of this class.
*
* @return the name.
**********************************************************************/
static std::string Name() {
return "MixerMT1<Random_u" +
std::string(mixer_t::width == 32 ? "32" : "64") + ">";
}
private:
static const mixer_type a = 5489ULL;
static const mixer_type
b = mixer_t::width == 32 ? 1812433253ULL : 6364136223846793005ULL;
static const mixer_type
c = mixer_t::width == 32 ? 1664525ULL : 3935559000370003845ULL;
static const mixer_type
d = mixer_t::width == 32 ? 1566083941ULL : 2862933555777941757ULL;
};
/**
* \brief The SFMT mixing functionality
*
* MixerSFMT is adapted from SFMT's init_by_array Mutsuo Saito given in
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.2.tar.gz
* and is part of the C++11 standard; see P. Becker, Working Draft, Standard
* for Programming Language C++, Oct. 2007, Sec. 26.4.7.1,
* http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf
*
* MixerSFMT contains a single change is to allow it to function properly
* when the size of the state is small.
*
* MixerSFMT mixes the seed much more thoroughly than MixerMT1 and, in
* particular, it removes the mixing bottleneck present in MixerMT1.
* Thus it is the recommended mixing scheme for all production work.
**********************************************************************/
class RANDOMLIB_EXPORT MixerSFMT {
public:
/**
* The RandomType controlling the output of MixerSFMT::SeedToState
**********************************************************************/
typedef Random_u32 mixer_t;
/**
* A version number which should be unique to this RandomMixer. This
* prevents RandomEngine::Load from loading a saved generator with a
* different RandomMixer. Here the version is "MxSM".
**********************************************************************/
static const unsigned version = 0x4d78534dUL;
private:
/**
* The unsigned type corresponding to mixer_t.
**********************************************************************/
typedef mixer_t::type mixer_type;
/**
* The mask for mixer_t.
**********************************************************************/
static const mixer_type mask = mixer_t::mask;
public:
/**
* Mix the seed vector, \e seed, into the state array, \e state, of size \e
* n.
*
* @param[in] seed the input seed vector.
* @param[out] state the generator state.
* @param[in] n the size of the state.
**********************************************************************/
static void SeedToState(const std::vector<RandomSeed::seed_type>& seed,
mixer_type state[], unsigned n) throw();
/**
* Return the name of this class.
*
* @return the name.
**********************************************************************/
static std::string Name() { return "MixerSFMT"; }
private:
static const mixer_type a = 0x8b8b8b8bUL;
static const mixer_type b = 1664525UL;
static const mixer_type c = 1566083941UL;
};
} // namespace RandomLib
#endif // RANDOMLIB_RANDOMMIXER_HPP

472
include/RandomLib/RandomNumber.hpp Normal file
View File

@ -0,0 +1,472 @@
/**
* \file RandomNumber.hpp
* \brief Header for RandomNumber
*
* Infinite precision random numbers.
*
* Copyright (c) Charles Karney (2006-2013) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMNUMBER_HPP)
#define RANDOMLIB_RANDOMNUMBER_HPP 1
#include <vector>
#include <iomanip>
#include <limits>
#include <cmath> // for std::pow
#include <RandomLib/UniformInteger.hpp>
namespace RandomLib {
/**
* \brief Infinite precision random numbers.
*
* Implement infinite precision random numbers. Integer part is non-random.
* Fraction part consists of any some number of digits in base
* 2<sup><i>b</i></sup>. If \e m digits have been generated then the
* fraction is uniformly distributed in the open interval
* &sum;<sub><i>k</i>=1</sub><sup><i>m</i></sup>
* <i>f</i><sub><i>k</i>&minus;1</sub>/2<sup><i>kb</i></sup> +
* (0,1)/2<sup><i>mb</i></sup>. When a RandomNumber is first constructed the
* integer part is zero and \e m = 0, and the number represents (0,1). A
* RandomNumber is able to represent all numbers in the symmetric open
* interval (&minus;2<sup>31</sup>, 2<sup>31</sup>). In this implementation,
* \e b must one of 1, 2, 3, 4, 8, 12, 16, 20, 24, 28, or 32. (This
* restriction allows printing in hexadecimal and can easily be relaxed.
* There's also no essential reason why the base should be a power of 2.)
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 1> class RandomNumber {
public:
/**
* Constructor sets number to a random number uniformly distributed in
* (0,1).
**********************************************************************/
RandomNumber() throw() : _n(0), _s(1) {}
/**
* Swap with another RandomNumber. This is a fast way of doing an
* assignment.
*
* @param[in,out] t the RandomNumber to swap with.
**********************************************************************/
void swap(RandomNumber& t) throw() {
if (this != &t) {
std::swap(_n, t._n);
std::swap(_s, t._s);
_f.swap(t._f);
}
}
/**
* Return to initial state, uniformly distributed in (0,1).
**********************************************************************/
void Init() throw() {
STATIC_ASSERT(bits > 0 && bits <= w && (bits < 4 || bits % 4 == 0),
"RandomNumber: unsupported value for bits");
_n = 0;
_s = 1;
_f.clear();
}
/**
* @return the sign of the RandomNumber (&plusmn; 1).
**********************************************************************/
int Sign() const throw() { return _s; }
/**
* Change the sign of the RandomNumber.
**********************************************************************/
void Negate() throw() { _s *= -1; }
/**
* @return the floor of the RandomNumber.
**********************************************************************/
int Floor() const throw() { return _s > 0 ? int(_n) : -1 - int(_n); }
/**
* @return the ceiling of the RandomNumber.
**********************************************************************/
int Ceiling() const throw() { return _s > 0 ? 1 + int(_n) : - int(_n); }
/**
* @return the unsigned integer component of the RandomNumber.
**********************************************************************/
unsigned UInteger() const throw() { return _n; }
/**
* Add integer \e k to the RandomNumber.
*
* @param[in] k the integer to add.
**********************************************************************/
void AddInteger(int k) throw() {
k += Floor(); // The new floor
int ns = k < 0 ? -1 : 1; // The new sign
if (ns != _s) // If sign changes, set f = 1 - f
for (size_t k = 0; k < Size(); ++k)
_f[k] = ~_f[k] & mask;
_n = ns > 0 ? k : -(k + 1);
}
/**
* Compare with another RandomNumber, *this &lt; \e t
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in,out] t a RandomNumber to compare.
* @return true if *this &lt; \e t.
**********************************************************************/
template<class Random> bool LessThan(Random& r, RandomNumber& t) {
if (this == &t) return false; // same object
if (_s != t._s) return _s < t._s;
if (_n != t._n) return (_s < 0) ^ (_n < t._n);
for (unsigned k = 0; ; ++k) {
// Impose an order on the evaluation of the digits.
const unsigned x = Digit(r,k);
const unsigned y = t.Digit(r,k);
if (x != y) return (_s < 0) ^ (x < y);
// Two distinct numbers are never equal
}
}
/**
* Compare RandomNumber with two others, *this &gt; max(\e u, \e v)
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in,out] u first RandomNumber to compare.
* @param[in,out] v second RandomNumber to compare.
* @return true if *this &gt; max(\e u, \e v).
**********************************************************************/
template<class Random> bool GreaterPair(Random& r,
RandomNumber& u, RandomNumber& v) {
// cmps is set to false as soon as u <= *this, and likewise for cmpt.
bool cmpu = this != &u, cmpv = this != &v && &u != &v;
if (!(cmpu || cmpv)) return true;
// Check signs first
if (cmpu) {
if (u._s > _s) return false; // u > *this
if (u._s < _s) cmpu = false;
}
if (cmpv) {
if (v._s > _s) return false; // v > *this
if (v._s < _s) cmpv = false;
}
if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
// Check integer parts
if (cmpu) {
if ((_s < 0) ^ (u._n > _n)) return false; // u > *this
if ((_s < 0) ^ (u._n < _n)) cmpu = false;
}
if (cmpv) {
if ((_s < 0) ^ (v._n > _n)) return false; // v > *this
if ((_s < 0) ^ (v._n < _n)) cmpv = false;
}
if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
// Check fractions
for (unsigned k = 0; ; ++k) {
// Impose an order on the evaluation of the digits. Note that this is
// asymmetric on interchange of u and v; since u is tested first, more
// digits of u are generated than v (on average).
const unsigned x = Digit(r,k);
if (cmpu) {
const unsigned y = u.Digit(r,k);
if ((_s < 0) ^ (y > x)) return false; // u > *this
if ((_s < 0) ^ (y < x)) cmpu = false;
}
if (cmpv) {
const unsigned y = v.Digit(r,k);
if ((_s < 0) ^ (y > x)) return false; // v > *this
if ((_s < 0) ^ (y < x)) cmpv = false;
}
if (!(cmpu || cmpv)) return true; // u <= *this && v <= *this
}
}
/**
* Compare with a fraction, *this &lt; <i>p</i>/<i>q</i>
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the numerator of the fraction.
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @return true if *this &lt; <i>p</i>/<i>q</i>.
**********************************************************************/
template<class Random, typename IntType>
bool LessThan(Random& r, IntType p, IntType q) {
for (int k = 0;; ++k) {
if (p <= 0) return false;
if (p >= q) return true;
// Here p is in [1,q-1]. Need to avoid overflow in computation of
// (q-1)<<bits and (2^bits-1)*q
p = (p << bits) - Digit(r,k) * q;
}
}
/**
* Compare with a paritally sampled fraction
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p0 the starting point for the numerator.
* @param[in] c the stride for the fraction (require \e c &gt; 0).
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @param[in,out] j the increment for the numerator.
* @return true if *this &lt; (<i>p</i><sub>0</sub> + <i>cj</i>)/<i>q</i>.
**********************************************************************/
template<class Random, typename IntType>
bool LessThan(Random& r, IntType p0, IntType c, IntType q,
UniformInteger<IntType, bits>& j) {
for (int k = 0;; ++k) {
if (j. LessThanEqual(r, - p0, c)) return false;
if (j.GreaterThanEqual(r, q - p0, c)) return true;
p0 = (p0 << bits) - IntType(Digit(r,k)) * q;
c <<= bits;
}
}
/**
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] k the index of a digit of the fraction
* @return digit number \e k, generating it if necessary.
**********************************************************************/
template<class Random> unsigned Digit(Random& r, unsigned k) {
ExpandTo(r, k + 1);
return _f[k];
}
/**
* Add one digit to the fraction.
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
**********************************************************************/
template<class Random> void AddDigit(Random& r)
{ _f.push_back(RandomDigit(r)); }
/**
* @param[in] k the index of a digit of the fraction
* @return a const reference to digit number \e k, without generating new
* digits.
* @exception std::out_of_range if the digit hasn't been generated.
**********************************************************************/
const unsigned& RawDigit(unsigned k) const throw()
{ return (const unsigned&)(_f.at(k)); }
/**
* @param[in] k the index of a digit of the fraction
* @return a non-const reference to digit number \e k, without generating
* new digits.
* @exception std::out_of_range if the digit hasn't been generated.
**********************************************************************/
unsigned& RawDigit(unsigned k) throw()
{ return (unsigned&)(_f.at(k)); }
/**
* Return to initial state, uniformly distributed in \e n + (0,1). This is
* similar to Init but also returns the memory used by the object to the
* system. Normally Init should be used.
**********************************************************************/
void Clear() {
std::vector<unsigned> z(0);
_n = 0;
_s = 1;
_f.swap(z);
}
/**
* @return the number of digits in fraction
**********************************************************************/
unsigned Size() const throw() { return unsigned(_f.size()); }
/**
* Return the fraction part of the RandomNumber as a floating point number
* of type RealType rounded to the nearest multiple of
* 1/2<sup><i>p</i></sup>, where \e p =
* std::numeric_limits<RealType>::digits, and, if necessary, creating
* additional digits of the number.
*
* @tparam RealType the floating point type to convert to.
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator for generating the necessary digits.
* @return the fraction of the RandomNumber rounded to a RealType.
**********************************************************************/
template<typename RealType, typename Random> RealType Fraction(Random& r) {
STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer,
"RandomNumber::Fraction: invalid real type RealType");
const int d = std::numeric_limits<RealType>::digits;
const int k = (d + bits - 1)/bits;
const int kg = (d + bits)/bits; // For guard bit
RealType y = 0;
if (Digit(r, kg - 1) & (1U << (kg * bits - d - 1)))
// if guard bit is set, round up.
y += std::pow(RealType(2), -d);
const RealType fact = std::pow(RealType(2), -bits);
RealType mult = RealType(1);
for (int i = 0; i < k; ++i) {
mult *= fact;
y += mult * RealType(i < k - 1 ? RawDigit(i) :
RawDigit(i) & (~0U << (k * bits - d)));
}
return y;
}
/**
* Return the value of the RandomNumber rounded to nearest floating point
* number of type RealType and, if necessary, creating additional digits of
* the number.
*
* @tparam RealType the floating point type to convert to.
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator for generating the necessary digits.
* @return the value of the RandomNumber rounded to a RealType.
**********************************************************************/
template<typename RealType, class Random> RealType Value(Random& r) {
// Ignore the possibility of overflow here (OK because int doesn't
// currently overflow any real type). Assume the real type supports
// denormalized numbers. Need to treat rounding explicitly since the
// missing digits always imply rounding up.
STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer,
"RandomNumber::Value: invalid real type RealType");
const int digits = std::numeric_limits<RealType>::digits,
min_exp = std::numeric_limits<RealType>::min_exponent;
RealType y;
int lead; // Position of leading bit (0.5 = position 0)
if (_n) lead = highest_bit_idx(_n);
else {
int i = 0;
while ( Digit(r, i) == 0 && i < (-min_exp)/bits ) ++i;
lead = highest_bit_idx(RawDigit(i)) - (i + 1) * bits;
// To handle denormalized numbers set lead = max(lead, min_exp)
lead = lead > min_exp ? lead : min_exp;
}
int trail = lead - digits; // Position of guard bit (0.5 = position 0)
if (trail > 0) {
y = RealType(_n & (~0U << trail));
if (_n & (1U << (trail - 1)))
y += std::pow(RealType(2), trail);
} else {
y = RealType(_n);
int k = (-trail)/bits; // Byte with guard bit
if (Digit(r, k) & (1U << ((k + 1) * bits + trail - 1)))
// If guard bit is set, round bit (some subsequent bit will be 1).
y += std::pow(RealType(2), trail);
// Byte with trailing bit (can be negative)
k = (-trail - 1 + bits)/bits - 1;
const RealType fact = std::pow(RealType(2), -bits);
RealType mult = RealType(1);
for (int i = 0; i <= k; ++i) {
mult *= fact;
y += mult *
RealType(i < k ? RawDigit(i) :
RawDigit(i) & (~0U << ((k + 1) * bits + trail)));
}
}
if (_s < 0) y *= -1;
return y;
}
/**
* Return the range of possible values for the RandomNumber as pair of
* doubles. This doesn't create any additional digits of the result and
* doesn't try to control roundoff.
*
* @return a pair denoting the range with first being the lower limit and
* second being the upper limit.
**********************************************************************/
std::pair<double, double> Range() const throw() {
double y = _n;
const double fact = std::pow(double(2), -bits);
double mult = double(1);
for (unsigned i = 0; i < Size(); ++i) {
mult *= fact;
y += mult * RawDigit(i);
}
return std::pair<double, double>(_s > 0 ? y : -(y + mult),
_s > 0 ? (y + mult) : -y);
}
/**
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @return a random digit in [0, 2<sup><i>bits</i></sup>).
**********************************************************************/
template<class Random> static unsigned RandomDigit(Random& r) throw()
{ return unsigned(r.template Integer<bits>()); }
private:
/**
* The integer part
**********************************************************************/
unsigned _n;
/**
* The sign
**********************************************************************/
int _s;
/**
* The fraction part
**********************************************************************/
std::vector<unsigned> _f;
/**
* Fill RandomNumber to \e k digits.
**********************************************************************/
template<class Random> void ExpandTo(Random& r, size_t k) {
size_t l = _f.size();
if (k <= l)
return;
_f.resize(k);
for (size_t i = l; i < k; ++i)
_f[i] = RandomDigit(r);
}
/**
* Return index [0..32] of highest bit set. Return 0 if x = 0, 32 is if x
* = ~0. (From Algorithms for programmers by Joerg Arndt.)
**********************************************************************/
static int highest_bit_idx(unsigned x) throw() {
if (x == 0) return 0;
int r = 1;
if (x & 0xffff0000U) { x >>= 16; r += 16; }
if (x & 0x0000ff00U) { x >>= 8; r += 8; }
if (x & 0x000000f0U) { x >>= 4; r += 4; }
if (x & 0x0000000cU) { x >>= 2; r += 2; }
if (x & 0x00000002U) { r += 1; }
return r;
}
/**
* The number of bits in unsigned.
**********************************************************************/
static const int w = std::numeric_limits<unsigned>::digits;
public:
/**
* A mask for the digits.
**********************************************************************/
static const unsigned mask =
bits == w ? ~0U : ~(~0U << (bits < w ? bits : 0));
};
/**
* \relates RandomNumber
* Print a RandomNumber. Format is n.dddd... where the base for printing is
* 2<sup>max(4,<i>b</i>)</sup>. The ... represents an infinite sequence of
* ungenerated random digits (uniformly distributed). Thus with \e b = 1,
* 0.0... = (0,1/2), 0.00... = (0,1/4), 0.11... = (3/4,1), etc.
**********************************************************************/
template<int bits>
std::ostream& operator<<(std::ostream& os, const RandomNumber<bits>& n) {
const std::ios::fmtflags oldflags = os.flags();
RandomNumber<bits> t = n;
os << (t.Sign() > 0 ? "+" : "-");
unsigned i = t.UInteger();
os << std::hex << std::setfill('0');
if (i == 0)
os << "0";
else {
bool first = true;
const int w = std::numeric_limits<unsigned>::digits;
const unsigned mask = RandomNumber<bits>::mask;
for (int s = ((w + bits - 1)/bits) * bits - bits; s >= 0; s -= bits) {
unsigned d = mask & (i >> s);
if (d || !first) {
if (first) {
os << d;
first = false;
}
else
os << std::setw((bits+3)/4) << d;
}
}
}
os << ".";
unsigned s = t.Size();
for (unsigned i = 0; i < s; ++i)
os << std::setw((bits+3)/4) << t.RawDigit(i);
os << "..." << std::setfill(' ');
os.flags(oldflags);
return os;
}
} // namespace RandomLib
#endif // RANDOMLIB_RANDOMNUMBER_HPP

77
include/RandomLib/RandomPower2.hpp Normal file
View File

@ -0,0 +1,77 @@
/**
* \file RandomPower2.hpp
* \brief Header for RandomPower2.
*
* Return and multiply by powers of two.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMPOWER2_HPP)
#define RANDOMLIB_RANDOMPOWER2_HPP 1
#include <cmath> // For std::pow
namespace RandomLib {
/**
* \brief Return or multiply by powers of 2
*
* With some compilers it's fastest to do a table lookup of powers of
* 2. If RANDOMLIB_POWERTABLE is 1, a lookup table is used. If
* RANDOMLIB_POWERTABLE is 0, then std::pow is used.
**********************************************************************/
class RANDOMLIB_EXPORT RandomPower2 {
public:
/**
* Return powers of 2 (either using a lookup table or std::pow)
*
* @param[in] n the integer power.
* @return 2<sup><i>n</i></sup>.
**********************************************************************/
template<typename RealType> static inline RealType pow2(int n) throw() {
#if RANDOMLIB_POWERTABLE
return RealType(power2[n - minpow]);
#else
return std::pow(RealType(2), n);
#endif
}
/**
* Multiply a real by a power of 2
*
* @tparam RealType the type of \e x.
* @param[in] x the real number.
* @param[in] n the power (positive or negative).
* @return \e x 2<sup><i>n</i></sup>.
**********************************************************************/
template<typename RealType>
static inline RealType shiftf(RealType x, int n) throw()
// std::ldexp(x, n); is equivalent, but slower
{ return x * pow2<RealType>(n); }
// Constants
enum {
/**
* Minimum power in RandomPower2::power2
**********************************************************************/
#if RANDOMLIB_LONGDOUBLEPREC > 64
minpow = -120,
#else
minpow = -64,
#endif
maxpow = 64 /**< Maximum power in RandomPower2::power2. */
};
private:
#if RANDOMLIB_POWERTABLE
/**
* Table of powers of two
**********************************************************************/
static const float power2[maxpow - minpow + 1]; // Powers of two
#endif
};
} // namespace RandomLib
#endif // RANDOMLIB_RANDOMPOWER2_HPP

251
include/RandomLib/RandomSeed.hpp Normal file
View File

@ -0,0 +1,251 @@
/**
* \file RandomSeed.hpp
* \brief Header for RandomSeed
*
* This provides a base class for random generators.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMSEED_HPP)
#define RANDOMLIB_RANDOMSEED_HPP 1
#include <iostream>
#include <stdexcept>
#include <vector>
#include <iterator>
#include <algorithm> // For std::transform
#include <sstream> // For VectorToString
#include <RandomLib/RandomType.hpp>
#if defined(_MSC_VER)
// Squelch warnings about dll vs vector
#pragma warning (push)
#pragma warning (disable: 4251)
#endif
namespace RandomLib {
/**
* \brief A base class for random generators
*
* This provides facilities for managing the seed and for converting the seed
* into random generator state.
*
* The seed is taken to be a vector of unsigned longs of arbitrary length.
* (Only the low 32 bit of each element of the vector are used.) The class
* provides several methods for setting the seed, static functions for
* producing "random" and "unique" seeds, and facilities for converting the
* seed to a string so that it can be printed easily.
*
* The seeding algorithms are those used by
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
* MT19937</a> with some modifications to make all states accessible and to
* minimize the likelihood of different seeds giving the same state.
*
* Finally some low-level routines are provided to facilitate the creation of
* I/O methods for the random generator.
*
* A random generator class can be written based on this class. The
* generator class would use the base class methods for setting the seed and
* for converting the seed into state. It would provide the machinery for
* advancing the state and for producing random data. It is also responsible
* for the routine to save and restore the generator state (including the
* seed).
*
* Written by Charles Karney <charles@karney.com> and licensed under the
* MIT/X11 License. The seeding algorithms are adapted from those of
* <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
* MT19937</a>. For more information, see http://randomlib.sourceforge.net/
**********************************************************************/
class RANDOMLIB_EXPORT RandomSeed {
public:
typedef Random_u32 u32;
typedef Random_u64 u64;
virtual ~RandomSeed() throw() = 0;
/**
* A type large enough to hold the seed words. This is needs to hold 32
* bits and is an unsigned long for portability.
**********************************************************************/
typedef RandomType<32, unsigned long> seed_t;
typedef seed_t::type seed_type;
/**
* \name Resetting the seed
**********************************************************************/
///@{
/**
* Set the seed to a vector \e v. Only the low \e 32 bits of each element
* are used.
*
* @tparam IntType the integral type of the elements of the vector.
* @param[in] v the vector of elements.
**********************************************************************/
template<typename IntType> void Reseed(const std::vector<IntType>& v) {
Reseed(v.begin(), v.end());
}
/**
* Set the seed to [\e a, \e b) from a pair of iterators. The iterator
* must produce results which can be converted into seed_type. Only the
* low 32 bits of each element are used.
*
* @tparam InputIterator the type of the iterator.
* @param[in] a the beginning iterator.
* @param[in] b the ending iterator.
**********************************************************************/
template<typename InputIterator>
void Reseed(InputIterator a, InputIterator b) {
// Read new seed into temporary so as not to change object on error.
std::vector<seed_type> t;
std::transform(a, b, back_inserter(t),
seed_t::cast<typename std::iterator_traits<InputIterator>
::value_type>);
_seed.swap(t);
Reset();
}
/**
* Set the seed to [\e n]. Only the low 32 bits of \e n are used.
*
* @param[in] n the new seed to use.
**********************************************************************/
void Reseed(seed_type n) {
// Reserve space for new seed so as not to change object on error.
_seed.reserve(1);
_seed.resize(1);
_seed[0] = seed_t::cast(n);
Reset();
}
/**
* Set the seed to [SeedVector()]. This is the standard way to reseed with
* a "unique" seed.
**********************************************************************/
void Reseed() { Reseed(SeedVector()); }
/**
* Set the seed from the string \e s using Random::StringToVector.
*
* @param[in] s the string to be decoded into a seed.
**********************************************************************/
void Reseed(const std::string& s) {
// Read new seed into temporary so as not to change object on error.
std::vector<seed_type> t = StringToVector(s);
_seed.swap(t);
Reset();
}
///@}
/**
* \name Examining the seed
**********************************************************************/
///@{
/**
* Return reference to the seed vector (read-only).
*
* @return the seed vector.
**********************************************************************/
const std::vector<seed_type>& Seed() const throw() { return _seed; }
/**
* Format the current seed suitable for printing.
*
* @return the seedd as a string.
**********************************************************************/
std::string SeedString() const { return VectorToString(_seed); }
///@}
/**
* \name Resetting the random seed
**********************************************************************/
///@{
/**
* Resets the sequence to its just-seeded state. This needs to be declared
* virtual here so that the Reseed functions can call it after saving the
* seed.
**********************************************************************/
virtual void Reset() throw() = 0;
///@}
/**
* \name Static functions for seed management
**********************************************************************/
///@{
/**
* Return a 32 bits of data suitable for seeding the random generator. The
* result is obtained by combining data from /dev/urandom, gettimeofday,
* time, and getpid to provide a reasonably "random" word of data.
* Usually, it is safer to seed the random generator with SeedVector()
* instead of SeedWord().
*
* @return a single "more-or-less random" seed_type to be used as a seed.
**********************************************************************/
static seed_type SeedWord();
/**
* Return a vector of unsigned longs suitable for seeding the random
* generator. The vector is almost certainly unique; however, the results
* of successive calls to Random::SeedVector() will be correlated. If
* several Random objects are required within a single program execution,
* call Random::SeedVector once, print it out (!), push_back additional
* data to identify the instance (e.g., loop index, thread ID, etc.), and
* use the result to seed the Random object. The number of elements
* included in the vector may depend on the operating system. Additional
* elements may be added in future versions of this library.
*
* @return a "unique" vector of seed_type to be uses as a seed.
**********************************************************************/
static std::vector<seed_type> SeedVector();
/**
* Convert a vector into a string suitable for printing or as an argument
* for Random::Reseed(const std::string& s).
*
* @tparam IntType the integral type of the elements of the vector.
* @param[in] v the vector to be converted.
* @return the resulting string.
**********************************************************************/
template<typename IntType>
static std::string VectorToString(const std::vector<IntType>& v) {
std::ostringstream os;
os << "[";
for (typename std::vector<IntType>::const_iterator n = v.begin();
n != v.end(); ++n) {
if (n != v.begin())
os << ",";
// Normalize in case this is called by user.
os << seed_t::cast(*n);
}
os << "]";
return os.str();
}
/**
* Convert a string into a vector of seed_type suitable for printing or as
* an argument for Random::Reseed(const std::vector<seed_type>& v). Reads
* consecutive digits in string. Thus "[1,2,3]" => [1,2,3]; "-0.123e-4" =>
* [0,123,4], etc. strtoul understands C's notation for octal and
* hexadecimal, for example "012 10 0xa" => [10,10,10]. Reading of a
* number stops at the first illegal character for the base. Thus
* "2006-04-08" => [2006,4,0,8] (i.e., 08 becomes two numbers). Note that
* input numbers greater than ULONG_MAX overflow to ULONG_MAX, which
* probably will result in the number being interpreted as LONG_MASK.
*
* @param[in] s the string to be converted.
* @return the resulting vector of seed_type.
**********************************************************************/
static std::vector<seed_type> StringToVector(const std::string& s);
///@}
protected:
/**
* The seed vector
**********************************************************************/
std::vector<seed_type> _seed;
};
inline RandomSeed::~RandomSeed() throw() {}
} // namespace RandomLib
#if defined(_MSC_VER)
#pragma warning (pop)
#endif
#endif // RANDOMLIB_RANDOMSEED_HPP

335
include/RandomLib/RandomSelect.hpp Normal file
View File

@ -0,0 +1,335 @@
/**
* \file RandomSelect.hpp
* \brief Header for RandomSelect.
*
* An implementation of the Walker algorithm for selecting from a finite set.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMSELECT_HPP)
#define RANDOMLIB_RANDOMSELECT_HPP 1
#include <vector>
#include <limits>
#include <stdexcept>
#if defined(_MSC_VER)
// Squelch warnings about constant conditional expressions
# pragma warning (push)
# pragma warning (disable: 4127)
#endif
namespace RandomLib {
/**
* \brief Random selection from a discrete set.
*
* An implementation of Walker algorithm for selecting from a finite set
* (following Knuth, TAOCP, Vol 2, Sec 3.4.1.A). This provides a rapid way
* of selecting one of several choices depending on a discrete set weights.
* Original citation is\n A. J. Walker,\n An Efficient Method for Generating
* Discrete Random Variables and General Distributions,\n ACM TOMS 3,
* 253--256 (1977).
*
* There are two changes here in the setup algorithm as given by Knuth:
*
* - The probabilities aren't sorted at the beginning of the setup; nor are
* they maintained in a sorted order. Instead they are just partitioned on
* the mean. This improves the setup time from O(\e k<sup>2</sup>) to O(\e
* k).
*
* - The internal calculations are carried out with type \e NumericType. If
* the input weights are of integer type, then choosing an integer type for
* \e NumericType yields an exact solution for the returned distribution
* (assuming that the underlying random generator is exact.)
*
* Example:
* \code
#include <RandomLib/RandomSelect.hpp>
// Weights for throwing a pair of dice
unsigned w[] = { 0, 0, 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1 };
// Initialize selection
RandomLib::RandomSelect<unsigned> sel(w, w + 13);
RandomLib::Random r; // Initialize random numbers
std::cout << "Seed set to " << r.SeedString() << "\n";
std::cout << "Throw a pair of dice 100 times:";
for (unsigned i = 0; i < 100; ++i)
std::cout << " " << sel(r);
std::cout << "\n";
\endcode
*
* @tparam NumericType the numeric type to use (default double).
**********************************************************************/
template<typename NumericType = double> class RandomSelect {
public:
/**
* Initialize in a cleared state (equivalent to having a single
* choice).
**********************************************************************/
RandomSelect() : _k(0), _wsum(0), _wmax(0) {}
/**
* Initialize with a weight vector \e w of elements of type \e WeightType.
* Internal calculations are carried out with type \e NumericType. \e
* NumericType needs to allow Choices() * MaxWeight() to be represented.
* Sensible combinations are:
* - \e WeightType integer, \e NumericType integer with
* digits(\e NumericType) &ge; digits(\e WeightType)
* - \e WeightType integer, \e NumericType real
* - \e WeightType real, \e NumericType real with digits(\e NumericType)
* &ge; digits(\e WeightType)
*
* @tparam WeightType the type of the weights.
* @param[in] w the vector of weights.
* @exception RandomErr if any of the weights are negative or if the total
* weight is not positive.
**********************************************************************/
template<typename WeightType>
RandomSelect(const std::vector<WeightType>& w) { Init(w.begin(), w.end()); }
/**
* Initialize with a weight given by a pair of iterators [\e a, \e b).
*
* @tparam InputIterator the type of the iterator.
* @param[in] a the beginning iterator.
* @param[in] b the ending iterator.
* @exception RandomErr if any of the weights are negative or if the total
* weight is not positive.
**********************************************************************/
template<typename InputIterator>
RandomSelect(InputIterator a, InputIterator b);
/**
* Clear the state (equivalent to having a single choice).
**********************************************************************/
void Init() throw()
{ _k = 0; _wsum = 0; _wmax = 0; _Q.clear(); _Y.clear(); }
/**
* Re-initialize with a weight vector \e w. Leave state unaltered in the
* case of an error.
*
* @tparam WeightType the type of the weights.
* @param[in] w the vector of weights.
**********************************************************************/
template<typename WeightType>
void Init(const std::vector<WeightType>& w) { Init(w.begin(), w.end()); }
/**
* Re-initialize with a weight given as a pair of iterators [\e a, \e b).
* Leave state unaltered in the case of an error.
*
* @tparam InputIterator the type of the iterator.
* @param[in] a the beginning iterator.
* @param[in] b the ending iterator.
**********************************************************************/
template<typename InputIterator>
void Init(InputIterator a, InputIterator b) {
RandomSelect<NumericType> t(a, b);
_Q.reserve(t._k);
_Y.reserve(t._k);
*this = t;
}
/**
* Return an index into the weight vector with probability proportional to
* the weight.
*
* @tparam Random the type of RandomCanonical generator.
* @param[in,out] r the RandomCanonical generator.
* @return the random index into the weight vector.
**********************************************************************/
template<class Random>
unsigned operator()(Random& r) const throw() {
if (_k <= 1)
return 0; // Special cases
const unsigned K = r.template Integer<unsigned>(_k);
// redundant casts to type NumericType to prevent warning from MS Project
return (std::numeric_limits<NumericType>::is_integer ?
r.template Prob<NumericType>(NumericType(_Q[K]),
NumericType(_wsum)) :
r.template Prob<NumericType>(NumericType(_Q[K]))) ?
K : _Y[K];
}
/**
* @return the sum of the weights.
**********************************************************************/
NumericType TotalWeight() const throw() { return _wsum; }
/**
* @return the maximum weight.
**********************************************************************/
NumericType MaxWeight() const throw() { return _wmax; }
/**
* @param[in] i the index in to the weight vector.
* @return the weight for sample \e i. Weight(i) / TotalWeight() gives the
* probability of sample \e i.
**********************************************************************/
NumericType Weight(unsigned i) const throw() {
if (i >= _k)
return NumericType(0);
else if (_k == 1)
return _wsum;
const NumericType n = std::numeric_limits<NumericType>::is_integer ?
_wsum : NumericType(1);
NumericType p = _Q[i];
for (unsigned j = _k; j;)
if (_Y[--j] == i)
p += n - _Q[j];
// If NumericType is integral, then p % _k == 0.
// assert(!std::numeric_limits<NumericType>::is_integer || p % _k == 0);
return (p / NumericType(_k)) * (_wsum / n);
}
/**
* @return the number of choices, i.e., the length of the weight vector.
**********************************************************************/
unsigned Choices() const throw() { return _k; }
private:
/**
* Size of weight vector
**********************************************************************/
unsigned _k;
/**
* Vector of cutoffs
**********************************************************************/
std::vector<NumericType> _Q;
/**
* Vector of aliases
**********************************************************************/
std::vector<unsigned> _Y;
/**
* The sum of the weights
**********************************************************************/
NumericType _wsum;
/**
* The maximum weight
**********************************************************************/
NumericType _wmax;
};
template<typename NumericType> template<typename InputIterator>
RandomSelect<NumericType>::RandomSelect(InputIterator a, InputIterator b) {
typedef typename std::iterator_traits<InputIterator>::value_type
WeightType;
// Disallow WeightType = real, NumericType = integer
STATIC_ASSERT(std::numeric_limits<WeightType>::is_integer ||
!std::numeric_limits<NumericType>::is_integer,
"RandomSelect: inconsistent WeightType and NumericType");
// If WeightType and NumericType are the same type, NumericType as precise
// as WeightType
STATIC_ASSERT(std::numeric_limits<WeightType>::is_integer !=
std::numeric_limits<NumericType>::is_integer ||
std::numeric_limits<NumericType>::digits >=
std::numeric_limits<WeightType>::digits,
"RandomSelect: NumericType insufficiently precise");
_wsum = 0;
_wmax = 0;
std::vector<NumericType> p;
for (InputIterator wptr = a; wptr != b; ++wptr) {
// Test *wptr < 0 without triggering compiler warning when *wptr =
// unsigned
if (!(*wptr > 0 || *wptr == 0))
// This also catches NaNs
throw RandomErr("RandomSelect: Illegal weight");
NumericType w = NumericType(*wptr);
if (w > (std::numeric_limits<NumericType>::max)() - _wsum)
throw RandomErr("RandomSelect: Overflow");
_wsum += w;
_wmax = w > _wmax ? w : _wmax;
p.push_back(w);
}
_k = unsigned(p.size());
if (_wsum <= 0)
throw RandomErr("RandomSelect: Zero total weight");
if (_k <= 1) {
// We treat k <= 1 as a special case in operator()
_Q.clear();
_Y.clear();
return;
}
if ((std::numeric_limits<NumericType>::max)()/NumericType(_k) <
NumericType(_wmax))
throw RandomErr("RandomSelect: Overflow");
std::vector<unsigned> j(_k);
_Q.resize(_k);
_Y.resize(_k);
// Pointers to the next empty low and high slots
unsigned u = 0;
unsigned v = _k - 1;
// Scale input and store in p and setup index array j. Note _wsum =
// mean(p). We could scale out _wsum here, but the following is exact when
// w[i] are low integers.
for (unsigned i = 0; i < _k; ++i) {
p[i] *= NumericType(_k);
j[p[i] > _wsum ? v-- : u++] = i;
}
// Pointers to the next low and high slots to use. Work towards the
// middle. This simplifies the loop exit test to u == v.
u = 0;
v = _k - 1;
// For integer NumericType, store the unnormalized probability in _Q and
// select using the exact Prob(_Q[k], _wsum). For real NumericType, store
// the normalized probability and select using Prob(_Q[k]). There will be
// a round off error in performing the division; but there is also the
// potential for round off errors in performing the arithmetic on p. There
// is therefore no point in simulating the division exactly using the
// slower Prob(real, real).
const NumericType n = std::numeric_limits<NumericType>::is_integer ?
NumericType(1) : _wsum;
while (true) {
// A loop invariant here is mean(p[j[u..v]]) == _wsum
_Q[j[u]] = p[j[u]] / n;
// If all arithmetic were exact this assignment could be:
// if (p[j[u]] < _wsum) _Y[j[u]] = j[v];
// But the following is safer:
_Y[j[u]] = j[p[j[u]] < _wsum ? v : u];
if (u == v) {
// The following assertion may fail because of roundoff errors
// assert( p[j[u]] == _wsum );
break;
}
// Update p, u, and v maintaining the loop invariant
p[j[v]] = p[j[v]] - (_wsum - p[j[u]]);
if (p[j[v]] > _wsum)
++u;
else
j[u] = j[v--];
}
return;
}
} // namespace RandomLib
#if defined(_MSC_VER)
# pragma warning (pop)
#endif
#endif // RANDOMLIB_RANDOMSELECT_HPP

137
include/RandomLib/RandomType.hpp Normal file
View File

@ -0,0 +1,137 @@
/**
* \file RandomType.hpp
* \brief Class to hold bit-width and unsigned type
*
* This provides a simple class to couple a bit-width and an unsigned type
* capable of holding all the bits. In addition is offers static methods for
* I/O and checksumming.
*
* Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_RANDOMTYPE_HPP)
#define RANDOMLIB_RANDOMTYPE_HPP 1
#include <limits>
#include <string>
#include <iostream>
namespace RandomLib {
/**
* \brief Class to hold bit-width and unsigned type
*
* This provides a simple class to couple a bit-width and an unsigned type
* capable of holding all the bits. In addition is offers static methods for
* I/O and checksumming.
*
* @tparam bits the number of significant bits.
* @tparam UIntType the C++ unsigned integer type capable of holding the bits.
**********************************************************************/
template<int bits, typename UIntType>
class RANDOMLIB_EXPORT RandomType {
public:
/**
* The unsigned C++ type
**********************************************************************/
typedef UIntType type;
/**
* The number of significant bits
**********************************************************************/
static const unsigned width = bits;
/**
* A mask for the significant bits.
**********************************************************************/
static const type mask =
~type(0) >> (std::numeric_limits<type>::digits - width);
/**
* The minimum representable value
**********************************************************************/
static const type min = type(0);
/**
* The maximum representable value
**********************************************************************/
static const type max = mask;
/**
* A combined masking and casting operation
*
* @tparam IntType the integer type of the \e x.
* @param[in] x the input integer.
* @return the masked and casted result.
**********************************************************************/
template<typename IntType> static type cast(IntType x) throw()
{ return type(x) & mask; }
/**
* Read a data value from a stream of 32-bit quantities (binary or text)
*
* @param[in,out] is the input stream.
* @param[in] bin true if the stream is binary.
* @param[out] x the data value read from the stream.
**********************************************************************/
static void Read32(std::istream& is, bool bin, type& x);
/**
* Write the data value to a stream of 32-bit quantities (binary or text)
*
* @param[in,out] os the output stream.
* @param[in] bin true if the stream is binary.
* @param[in,out] cnt controls the use of spaces and newlines for text
* output.
* @param[in] x the data value to be written to the stream.
*
* \e cnt should be zero on the first invocation of a series of writes.
* This function increments it by one on each call.
**********************************************************************/
static void Write32(std::ostream& os, bool bin, int& cnt, type x);
/**
* Accumulate a checksum of a integer into a 32-bit check. This implements
* a very simple checksum and is intended to avoid accidental corruption
* only.
*
* @param[in] n the number to be included in the checksum.
* @param[in,out] check the running checksum.
**********************************************************************/
static void CheckSum(type n, uint32_t& check) throw();
};
/**
* The standard unit for 32-bit quantities
**********************************************************************/
typedef RandomType<32, uint32_t> Random_u32;
/**
* The standard unit for 64-bit quantities
**********************************************************************/
typedef RandomType<64, uint64_t> Random_u64;
/**
* The integer type of constructing bitsets. This used to be unsigned long.
* C++11 has made this unsigned long long.
**********************************************************************/
#if defined(_MSC_VER) && _MSC_VER >= 1600
typedef unsigned long long bitset_uint_t;
#else
typedef unsigned long bitset_uint_t;
#endif
/// \cond SKIP
// Accumulate a checksum of a 32-bit quantity into check
template<>
inline void Random_u32::CheckSum(Random_u32::type n, Random_u32::type& check)
throw() {
// Circular shift left by one bit and add new word.
check = (check << 1 | (check >> 31 & Random_u32::type(1))) + n;
check &= Random_u32::mask;
}
// Accumulate a checksum of a 64-bit quantity into check
template<>
inline void Random_u64::CheckSum(Random_u64::type n, Random_u32::type& check)
throw() {
Random_u32::CheckSum(Random_u32::cast(n >> 32), check);
Random_u32::CheckSum(Random_u32::cast(n ), check);
}
/// \endcond
} // namespace RandomLib
#endif // RANDOMLIB_RANDOMTYPE_HPP

253
include/RandomLib/UniformInteger.hpp Normal file
View File

@ -0,0 +1,253 @@
/**
* \file UniformInteger.hpp
* \brief Header for UniformInteger
*
* Partially sample a uniform integer distribution.
*
* Copyright (c) Charles Karney (2013) <charles@karney.com> and licensed
* under the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_UNIFORMINTEGER_HPP)
#define RANDOMLIB_UNIFORMINTEGER_HPP 1
#include <limits>
namespace RandomLib {
/**
* \brief The partial uniform integer distribution.
*
* A class to sample in [0, \e m). For background, see:
* - D. E. Knuth and A. C. Yao, The Complexity of Nonuniform Random Number
* Generation, in "Algorithms and Complexity" (Academic Press, 1976),
* pp. 357--428.
* - J. Lumbroso, Optimal Discrete Uniform Generation from Coin Flips,
* and Applications, http://arxiv.org/abs/1304.1916 (2013)
* .
* Lumbroso's algorithm is a realization of the Knuth-Yao method for the case
* of uniform probabilities. This class generalizes the method to accept
* random digits in a base, \e b = 2<sup>\e bits</sup>. An important
* additional feature is that only sufficient random digits are drawn to
* narrow the allowed range to a power of b. Thus after
* <code>UniformInteger<int,1> u(r,5)</code>, \e u represents \verbatim
range prob
[0,4) 8/15
[0,2) 2/15
[2,4) 2/15
4 1/5 \endverbatim
* <code>u.Min()</code> and <code>u.Max()</code> give the extent of the
* closed range. The number of additional random digits needed to fix the
* value is given by <code>u.Entropy()</code>. The comparison operations may
* require additional digits to be drawn and so the range might be narrowed
* down. If you need a definite value then use <code>u(r)</code>.
*
* The DiscreteNormalAlt class uses UniformInteger to achieve an
* asymptotically ideal scaling wherein the number of random bits required
* per sample is constant + log<sub>2</sub>&sigma;. If Lumbroso's algorithm
* for sampling in [0,\e m) were used the log<sub>2</sub>&sigma; term would
* be multiplied by about 1.4.
*
* It is instructive to look at the Knuth-Yao discrete distribution
* generating (DDG) tree for the case \e m = 5 (the binary expansion of 1/5
* is 0.00110011...); Lumbroso's algorithm implements this tree.
* \image html ky-5.png "Knuth-Yao for \e m = 5"
*
* UniformInteger collapses all of the full subtrees above to their parent
* nodes to yield this tree where now some of the outcomes are ranges.
* \image html ky-5-collapse.png "Collapsed Knuth-Yao for \e m = 5"
*
* Averaging over many samples, the maximum number of digits required to
* construct a UniformInteger, i.e., invoking
* <code>UniformInteger(r,m)</code>, is (2\e b &minus; 1)/(\e b &minus; 1).
* (Note that this does not increase as \e m increases.) The maximum number
* of digits required to sample specific integers, i.e., invoking
* <code>UniformInteger(r,m)(r)</code>, is <i>b</i>/(\e b &minus; 1) +
* log<sub>\e b</sub>\e m. The worst cases are when \e m is slightly more
* than a power of \e b.
*
* The number of random bits required for sampling is shown as a function of
* the fractional part of log<sub>2</sub>\e m below. The red line shows what
* Lumbroso calls the "toll", the number of bits in excess of the entropy
* that are required for sampling.
* \image html
* uniform-bits.png "Random bits to sample in [0,\e m) for \e b = 2"
*
* @tparam IntType the type of the integer (must be signed).
* @tparam bits the number of bits in each digit used for sampling;
* the base for sampling is \e b = 2<sup>\e bits</sup>.
**********************************************************************/
template<typename IntType = int, int bits = 1> class UniformInteger {
public:
/**
* Constructor creating a partially sampled integer in [0, \e m)
*
* @param[in] r random object.
* @param[in] m constructed object represents an integer in [0, \e m).
* @param[in] flip if true, rearrange the ranges so that the widest ones
* are at near the upper end of [0, \e m) (default false).
*
* The samples enough random digits to obtain a uniform range whose size is
* a power of the base. The range can subsequently be narrowed by sampling
* additional digits.
**********************************************************************/
template<class Random>
UniformInteger(Random& r, IntType m, bool flip = false);
/**
* @return the minimum of the current range.
**********************************************************************/
IntType Min() const { return _a; }
/**
* @return the maximum of the current range.
**********************************************************************/
IntType Max() const { return _a + (IntType(1) << (_l * bits)) - 1; }
/**
* @return the entropy of the current range (in units of random digits).
*
* Max() + 1 - Min() = 2<sup>Entropy() * \e bits</sup>.
**********************************************************************/
IntType Entropy() const { return _l; }
/**
* Sample until the entropy vanishes, i.e., Min() = Max().
*
* @return the resulting integer sample.
**********************************************************************/
template<class Random> IntType operator()(Random& r)
{ while (_l) Refine(r); return _a; }
/**
* Negate the range, [Min(), Max()] &rarr; [&minus;Max(), &minus;Min()].
**********************************************************************/
void Negate() { _a = -Max(); }
/**
* Add a constant to the range
*
* @param[in] c the constant to be added.
*
* [Min(), Max()] &rarr; [Min() + \e c, Max() + \e c].
**********************************************************************/
void Add(IntType c) { _a += c; }
/**
* Compare with a fraction, *this &lt; <i>p</i>/<i>q</i>
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the numerator of the fraction.
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @return true if *this &lt; <i>p</i>/<i>q</i>.
**********************************************************************/
// test j < p/q (require q > 0)
template<class Random> bool LessThan(Random& r, IntType p, IntType q) {
for (;;) {
if ( (q * Max() < p)) return true;
if (!(q * Min() < p)) return false;
Refine(r);
}
}
/**
* Compare with a fraction, *this &le; <i>p</i>/<i>q</i>
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the numerator of the fraction.
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @return true if *this &le; <i>p</i>/<i>q</i>.
**********************************************************************/
template<class Random>
bool LessThanEqual(Random& r, IntType p, IntType q)
{ return LessThan(r, p + 1, q); }
/**
* Compare with a fraction, *this &gt; <i>p</i>/<i>q</i>
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the numerator of the fraction.
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @return true if *this &gt; <i>p</i>/<i>q</i>.
**********************************************************************/
template<class Random>
bool GreaterThan(Random& r, IntType p, IntType q)
{ return !LessThanEqual(r, p, q); }
/**
* Compare with a fraction, *this &ge; <i>p</i>/<i>q</i>
*
* @tparam Random the type of the random generator.
* @param[in,out] r a random generator.
* @param[in] p the numerator of the fraction.
* @param[in] q the denominator of the fraction (require \e q &gt; 0).
* @return true if *this &ge; <i>p</i>/<i>q</i>.
**********************************************************************/
template<class Random>
bool GreaterThanEqual(Random& r, IntType p, IntType q)
{ return !LessThan(r, p, q); }
/**
* Check that overflow will not happen.
*
* @param[in] mmax the largest \e m in the constructor.
* @param[in] qmax the largest \e q in LessThan().
* @return true if overflow will not happen.
*
* It is important that this check be carried out. If overflow occurs,
* incorrect results are obtained and the constructor may never terminate.
**********************************************************************/
static bool Check(IntType mmax, IntType qmax) {
return ( mmax - 1 <= ((std::numeric_limits<IntType>::max)() >> bits) &&
mmax - 1 <= (std::numeric_limits<IntType>::max)() / qmax );
}
private:
IntType _a, _l; // current range is _a + [0, 2^(bits*_l)).
template<class Random> static unsigned RandomDigit(Random& r) throw()
{ return unsigned(r.template Integer<bits>()); }
template<class Random> void Refine(Random& r) // only gets called if _l > 0.
{ _a += IntType(RandomDigit(r) << (bits * --_l)); }
};
template<typename IntType, int bits> template<class Random>
UniformInteger<IntType, bits>::UniformInteger(Random& r, IntType m, bool flip)
{
STATIC_ASSERT(std::numeric_limits<IntType>::is_integer,
"UniformInteger: invalid integer type IntType");
STATIC_ASSERT(std::numeric_limits<IntType>::is_signed,
"UniformInteger: IntType must be a signed type");
STATIC_ASSERT(bits > 0 && bits < std::numeric_limits<IntType>::digits &&
bits <= std::numeric_limits<unsigned>::digits,
"UniformInteger: bits out of range");
m = m < 1 ? 1 : m;
for (IntType v = 1, c = 0;;) {
_l = 0; _a = c;
for (IntType w = v, a = c, d = 1;;) {
// play out Lumbroso's algorithm without drawing random digits with w
// playing the role of v and c represented by the range [a, a + d).
// Return if both ends of range qualify as return values at the same
// time. Otherwise, fail and draw another random digit.
if (w >= m) {
IntType j = (a / m) * m; a -= j; w -= j;
if (w >= m) {
if (a + d <= m) { _a = !flip ? a : m - a - d; return; }
break;
}
}
w <<= bits; a <<= bits; d <<= bits; ++_l;
}
IntType j = (v / m) * m; v -= j; c -= j;
v <<= bits; c <<= bits; c += IntType(RandomDigit(r));
}
}
/**
* \relates UniformInteger
* Print a UniformInteger. Format is [\e min,\e max] unless the entropy is
* zero, in which case it's \e val.
**********************************************************************/
template<typename IntType, int bits>
std::ostream& operator<<(std::ostream& os,
const UniformInteger<IntType, bits>& u) {
if (u.Entropy())
os << "[" << u.Min() << "," << u.Max() << "]";
else
os << u.Min();
return os;
}
} // namespace RandomLib
#endif // RANDOMLIB_UNIFORMINTEGER_HPP