/* boost random/lagged_fibonacci.hpp header file | |
* | |
* Copyright Jens Maurer 2000-2001 | |
* Distributed under the Boost Software License, Version 1.0. (See | |
* accompanying file LICENSE_1_0.txt or copy at | |
* http://www.boost.org/LICENSE_1_0.txt) | |
* | |
* See http://www.boost.org for most recent version including documentation. | |
* | |
* $Id: lagged_fibonacci.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ | |
* | |
* Revision history | |
* 2001-02-18 moved to individual header files | |
*/ | |
#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP | |
#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP | |
#include <boost/config/no_tr1/cmath.hpp> | |
#include <iostream> | |
#include <algorithm> // std::max | |
#include <iterator> | |
#include <boost/config/no_tr1/cmath.hpp> // std::pow | |
#include <boost/config.hpp> | |
#include <boost/limits.hpp> | |
#include <boost/cstdint.hpp> | |
#include <boost/detail/workaround.hpp> | |
#include <boost/random/linear_congruential.hpp> | |
#include <boost/random/uniform_01.hpp> | |
#include <boost/random/detail/config.hpp> | |
#include <boost/random/detail/seed.hpp> | |
#include <boost/random/detail/pass_through_engine.hpp> | |
namespace boost { | |
namespace random { | |
#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300 | |
# define BOOST_RANDOM_EXTRACT_LF | |
#endif | |
#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3) | |
# define BOOST_RANDOM_EXTRACT_LF | |
#endif | |
# ifdef BOOST_RANDOM_EXTRACT_LF | |
namespace detail | |
{ | |
template<class IStream, class F, class RealType> | |
IStream& | |
extract_lagged_fibonacci_01( | |
IStream& is | |
, F const& f | |
, unsigned int& i | |
, RealType* x | |
, RealType modulus) | |
{ | |
is >> i >> std::ws; | |
for(unsigned int i = 0; i < f.long_lag; ++i) | |
{ | |
RealType value; | |
is >> value >> std::ws; | |
x[i] = value / modulus; | |
} | |
return is; | |
} | |
template<class IStream, class F, class UIntType> | |
IStream& | |
extract_lagged_fibonacci( | |
IStream& is | |
, F const& f | |
, unsigned int& i | |
, UIntType* x) | |
{ | |
is >> i >> std::ws; | |
for(unsigned int i = 0; i < f.long_lag; ++i) | |
is >> x[i] >> std::ws; | |
return is; | |
} | |
} | |
# endif | |
/** | |
* Instantiations of class template \lagged_fibonacci model a | |
* \pseudo_random_number_generator. It uses a lagged Fibonacci | |
* algorithm with two lags @c p and @c q: | |
* x(i) = x(i-p) + x(i-q) (mod 2<sup>w</sup>) with p > q. | |
*/ | |
template<class UIntType, int w, unsigned int p, unsigned int q, | |
UIntType val = 0> | |
class lagged_fibonacci | |
{ | |
public: | |
typedef UIntType result_type; | |
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); | |
BOOST_STATIC_CONSTANT(int, word_size = w); | |
BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); | |
BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); | |
/** | |
* Returns: the smallest value that the generator can produce | |
*/ | |
result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } | |
/** | |
* Returns: the largest value that the generator can produce | |
*/ | |
result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; } | |
/** | |
* Creates a new @c lagged_fibonacci generator and calls @c seed() | |
*/ | |
lagged_fibonacci() { init_wordmask(); seed(); } | |
/** | |
* Creates a new @c lagged_fibonacci generator and calls @c seed(value) | |
*/ | |
explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); } | |
/** | |
* Creates a new @c lagged_fibonacci generator and calls @c seed(first, last) | |
*/ | |
template<class It> lagged_fibonacci(It& first, It last) | |
{ init_wordmask(); seed(first, last); } | |
// compiler-generated copy ctor and assignment operator are fine | |
private: | |
/// \cond hide_private_members | |
void init_wordmask() | |
{ | |
wordmask = 0; | |
for(int j = 0; j < w; ++j) | |
wordmask |= (1u << j); | |
} | |
/// \endcond | |
public: | |
/** | |
* Sets the state of the generator to values produced by | |
* a \minstd_rand generator. | |
*/ | |
void seed(uint32_t value = 331u) | |
{ | |
minstd_rand0 gen(value); | |
for(unsigned int j = 0; j < long_lag; ++j) | |
x[j] = gen() & wordmask; | |
i = long_lag; | |
} | |
/** | |
* Sets the state of the generator to values from the iterator | |
* range [first, last). If there are not enough elements in the | |
* range [first, last) throws @c std::invalid_argument. | |
*/ | |
template<class It> | |
void seed(It& first, It last) | |
{ | |
// word size could be smaller than the seed values | |
unsigned int j; | |
for(j = 0; j < long_lag && first != last; ++j, ++first) | |
x[j] = *first & wordmask; | |
i = long_lag; | |
if(first == last && j < long_lag) | |
throw std::invalid_argument("lagged_fibonacci::seed"); | |
} | |
/** | |
* Returns: the next value of the generator | |
*/ | |
result_type operator()() | |
{ | |
if(i >= long_lag) | |
fill(); | |
return x[i++]; | |
} | |
static bool validation(result_type x) | |
{ | |
return x == val; | |
} | |
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE | |
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS | |
template<class CharT, class Traits> | |
friend std::basic_ostream<CharT,Traits>& | |
operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f) | |
{ | |
os << f.i << " "; | |
for(unsigned int i = 0; i < f.long_lag; ++i) | |
os << f.x[i] << " "; | |
return os; | |
} | |
template<class CharT, class Traits> | |
friend std::basic_istream<CharT, Traits>& | |
operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f) | |
{ | |
# ifdef BOOST_RANDOM_EXTRACT_LF | |
return detail::extract_lagged_fibonacci(is, f, f.i, f.x); | |
# else | |
is >> f.i >> std::ws; | |
for(unsigned int i = 0; i < f.long_lag; ++i) | |
is >> f.x[i] >> std::ws; | |
return is; | |
# endif | |
} | |
#endif | |
friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y) | |
{ return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } | |
friend bool operator!=(const lagged_fibonacci& x, | |
const lagged_fibonacci& y) | |
{ return !(x == y); } | |
#else | |
// Use a member function; Streamable concept not supported. | |
bool operator==(const lagged_fibonacci& rhs) const | |
{ return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } | |
bool operator!=(const lagged_fibonacci& rhs) const | |
{ return !(*this == rhs); } | |
#endif | |
private: | |
/// \cond hide_private_members | |
void fill(); | |
/// \endcond | |
UIntType wordmask; | |
unsigned int i; | |
UIntType x[long_lag]; | |
}; | |
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION | |
// A definition is required even for integral static constants | |
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val> | |
const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range; | |
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val> | |
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag; | |
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val> | |
const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag; | |
#endif | |
/// \cond hide_private_members | |
template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val> | |
void lagged_fibonacci<UIntType, w, p, q, val>::fill() | |
{ | |
// two loops to avoid costly modulo operations | |
{ // extra scope for MSVC brokenness w.r.t. for scope | |
for(unsigned int j = 0; j < short_lag; ++j) | |
x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask; | |
} | |
for(unsigned int j = short_lag; j < long_lag; ++j) | |
x[j] = (x[j] + x[j-short_lag]) & wordmask; | |
i = 0; | |
} | |
// lagged Fibonacci generator for the range [0..1) | |
// contributed by Matthias Troyer | |
// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958 | |
template<class T, unsigned int p, unsigned int q> | |
struct fibonacci_validation | |
{ | |
BOOST_STATIC_CONSTANT(bool, is_specialized = false); | |
static T value() { return 0; } | |
static T tolerance() { return 0; } | |
}; | |
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION | |
// A definition is required even for integral static constants | |
template<class T, unsigned int p, unsigned int q> | |
const bool fibonacci_validation<T, p, q>::is_specialized; | |
#endif | |
#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \ | |
template<> \ | |
struct fibonacci_validation<T, P, Q> \ | |
{ \ | |
BOOST_STATIC_CONSTANT(bool, is_specialized = true); \ | |
static T value() { return V; } \ | |
static T tolerance() \ | |
{ return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \ | |
}; | |
// (The extra static_cast<T> in the std::max call above is actually | |
// unnecessary except for HP aCC 1.30, which claims that | |
// numeric_limits<double>::epsilon() doesn't actually return a double.) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14) | |
BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14) | |
#undef BOOST_RANDOM_FIBONACCI_VAL | |
/// \endcond | |
/** | |
* Instantiations of class template @c lagged_fibonacci_01 model a | |
* \pseudo_random_number_generator. It uses a lagged Fibonacci | |
* algorithm with two lags @c p and @c q, evaluated in floating-point | |
* arithmetic: x(i) = x(i-p) + x(i-q) (mod 1) with p > q. See | |
* | |
* @blockquote | |
* "Uniform random number generators for supercomputers", Richard Brent, | |
* Proc. of Fifth Australian Supercomputer Conference, Melbourne, | |
* Dec. 1992, pp. 704-706. | |
* @endblockquote | |
* | |
* @xmlnote | |
* The quality of the generator crucially depends on the choice | |
* of the parameters. User code should employ one of the sensibly | |
* parameterized generators such as \lagged_fibonacci607 instead. | |
* @endxmlnote | |
* | |
* The generator requires considerable amounts of memory for the storage | |
* of its state array. For example, \lagged_fibonacci607 requires about | |
* 4856 bytes and \lagged_fibonacci44497 requires about 350 KBytes. | |
*/ | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
class lagged_fibonacci_01 | |
{ | |
public: | |
typedef RealType result_type; | |
BOOST_STATIC_CONSTANT(bool, has_fixed_range = false); | |
BOOST_STATIC_CONSTANT(int, word_size = w); | |
BOOST_STATIC_CONSTANT(unsigned int, long_lag = p); | |
BOOST_STATIC_CONSTANT(unsigned int, short_lag = q); | |
/** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(). */ | |
lagged_fibonacci_01() { init_modulus(); seed(); } | |
/** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(value). */ | |
BOOST_RANDOM_DETAIL_ARITHMETIC_CONSTRUCTOR(lagged_fibonacci_01, uint32_t, value) | |
{ init_modulus(); seed(value); } | |
/** Constructs a @c lagged_fibonacci_01 generator and calls @c seed(gen). */ | |
BOOST_RANDOM_DETAIL_GENERATOR_CONSTRUCTOR(lagged_fibonacci_01, Generator, gen) | |
{ init_modulus(); seed(gen); } | |
template<class It> lagged_fibonacci_01(It& first, It last) | |
{ init_modulus(); seed(first, last); } | |
// compiler-generated copy ctor and assignment operator are fine | |
private: | |
/// \cond hide_private_members | |
void init_modulus() | |
{ | |
#ifndef BOOST_NO_STDC_NAMESPACE | |
// allow for Koenig lookup | |
using std::pow; | |
#endif | |
_modulus = pow(RealType(2), word_size); | |
} | |
/// \endcond | |
public: | |
/** Calls seed(331u). */ | |
void seed() { seed(331u); } | |
/** | |
* Constructs a \minstd_rand0 generator with the constructor parameter | |
* value and calls seed with it. Distinct seeds in the range | |
* [1, 2147483647) will produce generators with different states. Other | |
* seeds will be equivalent to some seed within this range. See | |
* \linear_congruential for details. | |
*/ | |
BOOST_RANDOM_DETAIL_ARITHMETIC_SEED(lagged_fibonacci_01, uint32_t, value) | |
{ | |
minstd_rand0 intgen(value); | |
seed(intgen); | |
} | |
/** | |
* Sets the state of this @c lagged_fibonacci_01 to the values returned | |
* by p invocations of \uniform_01<code>\<RealType\>()(gen)</code>. | |
* | |
* Complexity: Exactly p invocations of gen. | |
*/ | |
BOOST_RANDOM_DETAIL_GENERATOR_SEED(lagged_fibonacci, Generator, gen) | |
{ | |
// use pass-by-reference, but wrap argument in pass_through_engine | |
typedef detail::pass_through_engine<Generator&> ref_gen; | |
uniform_01<ref_gen, RealType> gen01 = | |
uniform_01<ref_gen, RealType>(ref_gen(gen)); | |
// I could have used std::generate_n, but it takes "gen" by value | |
for(unsigned int j = 0; j < long_lag; ++j) | |
x[j] = gen01(); | |
i = long_lag; | |
} | |
template<class It> | |
void seed(It& first, It last) | |
{ | |
#ifndef BOOST_NO_STDC_NAMESPACE | |
// allow for Koenig lookup | |
using std::fmod; | |
using std::pow; | |
#endif | |
unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set | |
RealType two32 = pow(RealType(2), 32); | |
unsigned int j; | |
for(j = 0; j < long_lag && first != last; ++j) { | |
x[j] = RealType(0); | |
for(int k = 0; k < w/32 && first != last; ++k, ++first) | |
x[j] += *first / pow(two32,k+1); | |
if(first != last && mask != 0) { | |
x[j] += fmod((*first & mask) / _modulus, RealType(1)); | |
++first; | |
} | |
} | |
i = long_lag; | |
if(first == last && j < long_lag) | |
throw std::invalid_argument("lagged_fibonacci_01::seed"); | |
} | |
result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); } | |
result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); } | |
result_type operator()() | |
{ | |
if(i >= long_lag) | |
fill(); | |
return x[i++]; | |
} | |
static bool validation(result_type x) | |
{ | |
result_type v = fibonacci_validation<result_type, p, q>::value(); | |
result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance(); | |
// std::abs is a source of trouble: sometimes, it's not overloaded | |
// for double, plus the usual namespace std noncompliance -> avoid it | |
// using std::abs; | |
// return abs(x - v) < 5 * epsilon | |
return x > v - epsilon && x < v + epsilon; | |
} | |
#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE | |
#ifndef BOOST_RANDOM_NO_STREAM_OPERATORS | |
template<class CharT, class Traits> | |
friend std::basic_ostream<CharT,Traits>& | |
operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f) | |
{ | |
#ifndef BOOST_NO_STDC_NAMESPACE | |
// allow for Koenig lookup | |
using std::pow; | |
#endif | |
os << f.i << " "; | |
std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left); | |
for(unsigned int i = 0; i < f.long_lag; ++i) | |
os << f.x[i] * f._modulus << " "; | |
os.flags(oldflags); | |
return os; | |
} | |
template<class CharT, class Traits> | |
friend std::basic_istream<CharT, Traits>& | |
operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f) | |
{ | |
# ifdef BOOST_RANDOM_EXTRACT_LF | |
return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus); | |
# else | |
is >> f.i >> std::ws; | |
for(unsigned int i = 0; i < f.long_lag; ++i) { | |
typename lagged_fibonacci_01::result_type value; | |
is >> value >> std::ws; | |
f.x[i] = value / f._modulus; | |
} | |
return is; | |
# endif | |
} | |
#endif | |
friend bool operator==(const lagged_fibonacci_01& x, | |
const lagged_fibonacci_01& y) | |
{ return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); } | |
friend bool operator!=(const lagged_fibonacci_01& x, | |
const lagged_fibonacci_01& y) | |
{ return !(x == y); } | |
#else | |
// Use a member function; Streamable concept not supported. | |
bool operator==(const lagged_fibonacci_01& rhs) const | |
{ return i == rhs.i && std::equal(x, x+long_lag, rhs.x); } | |
bool operator!=(const lagged_fibonacci_01& rhs) const | |
{ return !(*this == rhs); } | |
#endif | |
private: | |
/// \cond hide_private_members | |
void fill(); | |
/// \endcond | |
unsigned int i; | |
RealType x[long_lag]; | |
RealType _modulus; | |
}; | |
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION | |
// A definition is required even for integral static constants | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range; | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag; | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag; | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
const int lagged_fibonacci_01<RealType,w,p,q>::word_size; | |
#endif | |
/// \cond hide_private_members | |
template<class RealType, int w, unsigned int p, unsigned int q> | |
void lagged_fibonacci_01<RealType, w, p, q>::fill() | |
{ | |
// two loops to avoid costly modulo operations | |
{ // extra scope for MSVC brokenness w.r.t. for scope | |
for(unsigned int j = 0; j < short_lag; ++j) { | |
RealType t = x[j] + x[j+(long_lag-short_lag)]; | |
if(t >= RealType(1)) | |
t -= RealType(1); | |
x[j] = t; | |
} | |
} | |
for(unsigned int j = short_lag; j < long_lag; ++j) { | |
RealType t = x[j] + x[j-short_lag]; | |
if(t >= RealType(1)) | |
t -= RealType(1); | |
x[j] = t; | |
} | |
i = 0; | |
} | |
/// \endcond | |
} // namespace random | |
#ifdef BOOST_RANDOM_DOXYGEN | |
namespace detail { | |
/** | |
* The specializations lagged_fibonacci607 ... lagged_fibonacci44497 | |
* use well tested lags. | |
* | |
* See | |
* | |
* @blockquote | |
* "On the Periods of Generalized Fibonacci Recurrences", Richard P. Brent | |
* Computer Sciences Laboratory Australian National University, December 1992 | |
* @endblockquote | |
* | |
* The lags used here can be found in | |
* | |
* @blockquote | |
* "Uniform random number generators for supercomputers", Richard Brent, | |
* Proc. of Fifth Australian Supercomputer Conference, Melbourne, | |
* Dec. 1992, pp. 704-706. | |
* @endblockquote | |
*/ | |
struct lagged_fibonacci_doc {}; | |
} | |
#endif | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209; | |
/** | |
* @copydoc boost::detail::lagged_fibonacci_doc | |
*/ | |
typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497; | |
// It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<> | |
// to help the compiler generate efficient code. For GCC, this seems useless, | |
// because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now. | |
} // namespace boost | |
#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP |