1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/epoc32/include/stdapis/boost/random/lagged_fibonacci.hpp Tue Mar 16 16:12:26 2010 +0000
1.3 @@ -0,0 +1,464 @@
1.4 +/* boost random/lagged_fibonacci.hpp header file
1.5 + *
1.6 + * Copyright Jens Maurer 2000-2001
1.7 + * Distributed under the Boost Software License, Version 1.0. (See
1.8 + * accompanying file LICENSE_1_0.txt or copy at
1.9 + * http://www.boost.org/LICENSE_1_0.txt)
1.10 + *
1.11 + * See http://www.boost.org for most recent version including documentation.
1.12 + *
1.13 + * $Id: lagged_fibonacci.hpp,v 1.28 2005/05/21 15:57:00 dgregor Exp $
1.14 + *
1.15 + * Revision history
1.16 + * 2001-02-18 moved to individual header files
1.17 + */
1.18 +
1.19 +#ifndef BOOST_RANDOM_LAGGED_FIBONACCI_HPP
1.20 +#define BOOST_RANDOM_LAGGED_FIBONACCI_HPP
1.21 +
1.22 +#include <cmath>
1.23 +#include <iostream>
1.24 +#include <algorithm> // std::max
1.25 +#include <iterator>
1.26 +#include <cmath> // std::pow
1.27 +#include <boost/config.hpp>
1.28 +#include <boost/limits.hpp>
1.29 +#include <boost/cstdint.hpp>
1.30 +#include <boost/detail/workaround.hpp>
1.31 +#include <boost/random/linear_congruential.hpp>
1.32 +#include <boost/random/uniform_01.hpp>
1.33 +#include <boost/random/detail/pass_through_engine.hpp>
1.34 +
1.35 +namespace boost {
1.36 +namespace random {
1.37 +
1.38 +#if BOOST_WORKAROUND(_MSC_FULL_VER, BOOST_TESTED_AT(13102292)) && BOOST_MSVC > 1300
1.39 +# define BOOST_RANDOM_EXTRACT_LF
1.40 +#endif
1.41 +
1.42 +#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
1.43 +# define BOOST_RANDOM_EXTRACT_LF
1.44 +#endif
1.45 +
1.46 +# ifdef BOOST_RANDOM_EXTRACT_LF
1.47 +namespace detail
1.48 +{
1.49 + template<class IStream, class F, class RealType>
1.50 + IStream&
1.51 + extract_lagged_fibonacci_01(
1.52 + IStream& is
1.53 + , F const& f
1.54 + , unsigned int& i
1.55 + , RealType* x
1.56 + , RealType modulus)
1.57 + {
1.58 + is >> i >> std::ws;
1.59 + for(unsigned int i = 0; i < f.long_lag; ++i)
1.60 + {
1.61 + RealType value;
1.62 + is >> value >> std::ws;
1.63 + x[i] = value / modulus;
1.64 + }
1.65 + return is;
1.66 + }
1.67 +
1.68 + template<class IStream, class F, class UIntType>
1.69 + IStream&
1.70 + extract_lagged_fibonacci(
1.71 + IStream& is
1.72 + , F const& f
1.73 + , unsigned int& i
1.74 + , UIntType* x)
1.75 + {
1.76 + is >> i >> std::ws;
1.77 + for(unsigned int i = 0; i < f.long_lag; ++i)
1.78 + is >> x[i] >> std::ws;
1.79 + return is;
1.80 + }
1.81 +}
1.82 +# endif
1.83 +
1.84 +template<class UIntType, int w, unsigned int p, unsigned int q,
1.85 + UIntType val = 0>
1.86 +class lagged_fibonacci
1.87 +{
1.88 +public:
1.89 + typedef UIntType result_type;
1.90 + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
1.91 + BOOST_STATIC_CONSTANT(int, word_size = w);
1.92 + BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
1.93 + BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
1.94 +
1.95 + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
1.96 + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return wordmask; }
1.97 +
1.98 + lagged_fibonacci() { init_wordmask(); seed(); }
1.99 + explicit lagged_fibonacci(uint32_t value) { init_wordmask(); seed(value); }
1.100 + template<class It> lagged_fibonacci(It& first, It last)
1.101 + { init_wordmask(); seed(first, last); }
1.102 + // compiler-generated copy ctor and assignment operator are fine
1.103 +
1.104 +private:
1.105 + void init_wordmask()
1.106 + {
1.107 + wordmask = 0;
1.108 + for(int i = 0; i < w; ++i)
1.109 + wordmask |= (1u << i);
1.110 + }
1.111 +
1.112 +public:
1.113 + void seed(uint32_t value = 331u)
1.114 + {
1.115 + minstd_rand0 gen(value);
1.116 + for(unsigned int j = 0; j < long_lag; ++j)
1.117 + x[j] = gen() & wordmask;
1.118 + i = long_lag;
1.119 + }
1.120 +
1.121 + template<class It>
1.122 + void seed(It& first, It last)
1.123 + {
1.124 + // word size could be smaller than the seed values
1.125 + unsigned int j;
1.126 + for(j = 0; j < long_lag && first != last; ++j, ++first)
1.127 + x[j] = *first & wordmask;
1.128 + i = long_lag;
1.129 + if(first == last && j < long_lag)
1.130 + throw std::invalid_argument("lagged_fibonacci::seed");
1.131 + }
1.132 +
1.133 + result_type operator()()
1.134 + {
1.135 + if(i >= long_lag)
1.136 + fill();
1.137 + return x[i++];
1.138 + }
1.139 +
1.140 + static bool validation(result_type x)
1.141 + {
1.142 + return x == val;
1.143 + }
1.144 +
1.145 +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
1.146 +
1.147 +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
1.148 + template<class CharT, class Traits>
1.149 + friend std::basic_ostream<CharT,Traits>&
1.150 + operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci& f)
1.151 + {
1.152 + os << f.i << " ";
1.153 + for(unsigned int i = 0; i < f.long_lag; ++i)
1.154 + os << f.x[i] << " ";
1.155 + return os;
1.156 + }
1.157 +
1.158 + template<class CharT, class Traits>
1.159 + friend std::basic_istream<CharT, Traits>&
1.160 + operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci& f)
1.161 + {
1.162 +# ifdef BOOST_RANDOM_EXTRACT_LF
1.163 + return detail::extract_lagged_fibonacci(is, f, f.i, f.x);
1.164 +# else
1.165 + is >> f.i >> std::ws;
1.166 + for(unsigned int i = 0; i < f.long_lag; ++i)
1.167 + is >> f.x[i] >> std::ws;
1.168 + return is;
1.169 +# endif
1.170 + }
1.171 +#endif
1.172 +
1.173 + friend bool operator==(const lagged_fibonacci& x, const lagged_fibonacci& y)
1.174 + { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
1.175 + friend bool operator!=(const lagged_fibonacci& x,
1.176 + const lagged_fibonacci& y)
1.177 + { return !(x == y); }
1.178 +#else
1.179 + // Use a member function; Streamable concept not supported.
1.180 + bool operator==(const lagged_fibonacci& rhs) const
1.181 + { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
1.182 + bool operator!=(const lagged_fibonacci& rhs) const
1.183 + { return !(*this == rhs); }
1.184 +#endif
1.185 +
1.186 +private:
1.187 + void fill();
1.188 + UIntType wordmask;
1.189 + unsigned int i;
1.190 + UIntType x[long_lag];
1.191 +};
1.192 +
1.193 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.194 +// A definition is required even for integral static constants
1.195 +template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
1.196 +const bool lagged_fibonacci<UIntType, w, p, q, val>::has_fixed_range;
1.197 +template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
1.198 +const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::long_lag;
1.199 +template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
1.200 +const unsigned int lagged_fibonacci<UIntType, w, p, q, val>::short_lag;
1.201 +#endif
1.202 +
1.203 +template<class UIntType, int w, unsigned int p, unsigned int q, UIntType val>
1.204 +void lagged_fibonacci<UIntType, w, p, q, val>::fill()
1.205 +{
1.206 + // two loops to avoid costly modulo operations
1.207 + { // extra scope for MSVC brokenness w.r.t. for scope
1.208 + for(unsigned int j = 0; j < short_lag; ++j)
1.209 + x[j] = (x[j] + x[j+(long_lag-short_lag)]) & wordmask;
1.210 + }
1.211 + for(unsigned int j = short_lag; j < long_lag; ++j)
1.212 + x[j] = (x[j] + x[j-short_lag]) & wordmask;
1.213 + i = 0;
1.214 +}
1.215 +
1.216 +
1.217 +
1.218 +// lagged Fibonacci generator for the range [0..1)
1.219 +// contributed by Matthias Troyer
1.220 +// for p=55, q=24 originally by G. J. Mitchell and D. P. Moore 1958
1.221 +
1.222 +template<class T, unsigned int p, unsigned int q>
1.223 +struct fibonacci_validation
1.224 +{
1.225 + BOOST_STATIC_CONSTANT(bool, is_specialized = false);
1.226 + static T value() { return 0; }
1.227 + static T tolerance() { return 0; }
1.228 +};
1.229 +
1.230 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.231 +// A definition is required even for integral static constants
1.232 +template<class T, unsigned int p, unsigned int q>
1.233 +const bool fibonacci_validation<T, p, q>::is_specialized;
1.234 +#endif
1.235 +
1.236 +#define BOOST_RANDOM_FIBONACCI_VAL(T,P,Q,V,E) \
1.237 +template<> \
1.238 +struct fibonacci_validation<T, P, Q> \
1.239 +{ \
1.240 + BOOST_STATIC_CONSTANT(bool, is_specialized = true); \
1.241 + static T value() { return V; } \
1.242 + static T tolerance() \
1.243 +{ return (std::max)(E, static_cast<T>(5*std::numeric_limits<T>::epsilon())); } \
1.244 +};
1.245 +// (The extra static_cast<T> in the std::max call above is actually
1.246 +// unnecessary except for HP aCC 1.30, which claims that
1.247 +// numeric_limits<double>::epsilon() doesn't actually return a double.)
1.248 +
1.249 +BOOST_RANDOM_FIBONACCI_VAL(double, 607, 273, 0.4293817707235914, 1e-14)
1.250 +BOOST_RANDOM_FIBONACCI_VAL(double, 1279, 418, 0.9421630240437659, 1e-14)
1.251 +BOOST_RANDOM_FIBONACCI_VAL(double, 2281, 1252, 0.1768114046909004, 1e-14)
1.252 +BOOST_RANDOM_FIBONACCI_VAL(double, 3217, 576, 0.1956232694868209, 1e-14)
1.253 +BOOST_RANDOM_FIBONACCI_VAL(double, 4423, 2098, 0.9499762202147172, 1e-14)
1.254 +BOOST_RANDOM_FIBONACCI_VAL(double, 9689, 5502, 0.05737836943695162, 1e-14)
1.255 +BOOST_RANDOM_FIBONACCI_VAL(double, 19937, 9842, 0.5076528587449834, 1e-14)
1.256 +BOOST_RANDOM_FIBONACCI_VAL(double, 23209, 13470, 0.5414473810619185, 1e-14)
1.257 +BOOST_RANDOM_FIBONACCI_VAL(double, 44497,21034, 0.254135073399297, 1e-14)
1.258 +
1.259 +#undef BOOST_RANDOM_FIBONACCI_VAL
1.260 +
1.261 +template<class RealType, int w, unsigned int p, unsigned int q>
1.262 +class lagged_fibonacci_01
1.263 +{
1.264 +public:
1.265 + typedef RealType result_type;
1.266 + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
1.267 + BOOST_STATIC_CONSTANT(int, word_size = w);
1.268 + BOOST_STATIC_CONSTANT(unsigned int, long_lag = p);
1.269 + BOOST_STATIC_CONSTANT(unsigned int, short_lag = q);
1.270 +
1.271 + lagged_fibonacci_01() { init_modulus(); seed(); }
1.272 + explicit lagged_fibonacci_01(uint32_t value) { init_modulus(); seed(value); }
1.273 + template<class Generator>
1.274 + explicit lagged_fibonacci_01(Generator & gen) { init_modulus(); seed(gen); }
1.275 + template<class It> lagged_fibonacci_01(It& first, It last)
1.276 + { init_modulus(); seed(first, last); }
1.277 + // compiler-generated copy ctor and assignment operator are fine
1.278 +
1.279 +private:
1.280 + void init_modulus()
1.281 + {
1.282 +#ifndef BOOST_NO_STDC_NAMESPACE
1.283 + // allow for Koenig lookup
1.284 + using std::pow;
1.285 +#endif
1.286 + _modulus = pow(RealType(2), word_size);
1.287 + }
1.288 +
1.289 +public:
1.290 + void seed(uint32_t value = 331u)
1.291 + {
1.292 + minstd_rand0 intgen(value);
1.293 + seed(intgen);
1.294 + }
1.295 +
1.296 + // For GCC, moving this function out-of-line prevents inlining, which may
1.297 + // reduce overall object code size. However, MSVC does not grok
1.298 + // out-of-line template member functions.
1.299 + template<class Generator>
1.300 + void seed(Generator & gen)
1.301 + {
1.302 + // use pass-by-reference, but wrap argument in pass_through_engine
1.303 + typedef detail::pass_through_engine<Generator&> ref_gen;
1.304 + uniform_01<ref_gen, RealType> gen01 =
1.305 + uniform_01<ref_gen, RealType>(ref_gen(gen));
1.306 + // I could have used std::generate_n, but it takes "gen" by value
1.307 + for(unsigned int j = 0; j < long_lag; ++j)
1.308 + x[j] = gen01();
1.309 + i = long_lag;
1.310 + }
1.311 +
1.312 + template<class It>
1.313 + void seed(It& first, It last)
1.314 + {
1.315 +#ifndef BOOST_NO_STDC_NAMESPACE
1.316 + // allow for Koenig lookup
1.317 + using std::fmod;
1.318 + using std::pow;
1.319 +#endif
1.320 + unsigned long mask = ~((~0u) << (w%32)); // now lowest w bits set
1.321 + RealType two32 = pow(RealType(2), 32);
1.322 + unsigned int j;
1.323 + for(j = 0; j < long_lag && first != last; ++j, ++first) {
1.324 + x[j] = RealType(0);
1.325 + for(int i = 0; i < w/32 && first != last; ++i, ++first)
1.326 + x[j] += *first / pow(two32,i+1);
1.327 + if(first != last && mask != 0)
1.328 + x[j] += fmod((*first & mask) / _modulus, RealType(1));
1.329 + }
1.330 + i = long_lag;
1.331 + if(first == last && j < long_lag)
1.332 + throw std::invalid_argument("lagged_fibonacci_01::seed");
1.333 + }
1.334 +
1.335 + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
1.336 + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
1.337 +
1.338 + result_type operator()()
1.339 + {
1.340 + if(i >= long_lag)
1.341 + fill();
1.342 + return x[i++];
1.343 + }
1.344 +
1.345 + static bool validation(result_type x)
1.346 + {
1.347 + result_type v = fibonacci_validation<result_type, p, q>::value();
1.348 + result_type epsilon = fibonacci_validation<result_type, p, q>::tolerance();
1.349 + // std::abs is a source of trouble: sometimes, it's not overloaded
1.350 + // for double, plus the usual namespace std noncompliance -> avoid it
1.351 + // using std::abs;
1.352 + // return abs(x - v) < 5 * epsilon
1.353 + return x > v - epsilon && x < v + epsilon;
1.354 + }
1.355 +
1.356 +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
1.357 +
1.358 +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
1.359 + template<class CharT, class Traits>
1.360 + friend std::basic_ostream<CharT,Traits>&
1.361 + operator<<(std::basic_ostream<CharT,Traits>& os, const lagged_fibonacci_01&f)
1.362 + {
1.363 +#ifndef BOOST_NO_STDC_NAMESPACE
1.364 + // allow for Koenig lookup
1.365 + using std::pow;
1.366 +#endif
1.367 + os << f.i << " ";
1.368 + std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
1.369 + for(unsigned int i = 0; i < f.long_lag; ++i)
1.370 + os << f.x[i] * f._modulus << " ";
1.371 + os.flags(oldflags);
1.372 + return os;
1.373 + }
1.374 +
1.375 + template<class CharT, class Traits>
1.376 + friend std::basic_istream<CharT, Traits>&
1.377 + operator>>(std::basic_istream<CharT, Traits>& is, lagged_fibonacci_01& f)
1.378 + {
1.379 +# ifdef BOOST_RANDOM_EXTRACT_LF
1.380 + return detail::extract_lagged_fibonacci_01(is, f, f.i, f.x, f._modulus);
1.381 +# else
1.382 + is >> f.i >> std::ws;
1.383 + for(unsigned int i = 0; i < f.long_lag; ++i) {
1.384 + typename lagged_fibonacci_01::result_type value;
1.385 + is >> value >> std::ws;
1.386 + f.x[i] = value / f._modulus;
1.387 + }
1.388 + return is;
1.389 +# endif
1.390 + }
1.391 +#endif
1.392 +
1.393 + friend bool operator==(const lagged_fibonacci_01& x,
1.394 + const lagged_fibonacci_01& y)
1.395 + { return x.i == y.i && std::equal(x.x, x.x+long_lag, y.x); }
1.396 + friend bool operator!=(const lagged_fibonacci_01& x,
1.397 + const lagged_fibonacci_01& y)
1.398 + { return !(x == y); }
1.399 +#else
1.400 + // Use a member function; Streamable concept not supported.
1.401 + bool operator==(const lagged_fibonacci_01& rhs) const
1.402 + { return i == rhs.i && std::equal(x, x+long_lag, rhs.x); }
1.403 + bool operator!=(const lagged_fibonacci_01& rhs) const
1.404 + { return !(*this == rhs); }
1.405 +#endif
1.406 +
1.407 +private:
1.408 + void fill();
1.409 + unsigned int i;
1.410 + RealType x[long_lag];
1.411 + RealType _modulus;
1.412 +};
1.413 +
1.414 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.415 +// A definition is required even for integral static constants
1.416 +template<class RealType, int w, unsigned int p, unsigned int q>
1.417 +const bool lagged_fibonacci_01<RealType, w, p, q>::has_fixed_range;
1.418 +template<class RealType, int w, unsigned int p, unsigned int q>
1.419 +const unsigned int lagged_fibonacci_01<RealType, w, p, q>::long_lag;
1.420 +template<class RealType, int w, unsigned int p, unsigned int q>
1.421 +const unsigned int lagged_fibonacci_01<RealType, w, p, q>::short_lag;
1.422 +template<class RealType, int w, unsigned int p, unsigned int q>
1.423 +const int lagged_fibonacci_01<RealType,w,p,q>::word_size;
1.424 +
1.425 +#endif
1.426 +
1.427 +template<class RealType, int w, unsigned int p, unsigned int q>
1.428 +void lagged_fibonacci_01<RealType, w, p, q>::fill()
1.429 +{
1.430 + // two loops to avoid costly modulo operations
1.431 + { // extra scope for MSVC brokenness w.r.t. for scope
1.432 + for(unsigned int j = 0; j < short_lag; ++j) {
1.433 + RealType t = x[j] + x[j+(long_lag-short_lag)];
1.434 + if(t >= RealType(1))
1.435 + t -= RealType(1);
1.436 + x[j] = t;
1.437 + }
1.438 + }
1.439 + for(unsigned int j = short_lag; j < long_lag; ++j) {
1.440 + RealType t = x[j] + x[j-short_lag];
1.441 + if(t >= RealType(1))
1.442 + t -= RealType(1);
1.443 + x[j] = t;
1.444 + }
1.445 + i = 0;
1.446 +}
1.447 +
1.448 +} // namespace random
1.449 +
1.450 +typedef random::lagged_fibonacci_01<double, 48, 607, 273> lagged_fibonacci607;
1.451 +typedef random::lagged_fibonacci_01<double, 48, 1279, 418> lagged_fibonacci1279;
1.452 +typedef random::lagged_fibonacci_01<double, 48, 2281, 1252> lagged_fibonacci2281;
1.453 +typedef random::lagged_fibonacci_01<double, 48, 3217, 576> lagged_fibonacci3217;
1.454 +typedef random::lagged_fibonacci_01<double, 48, 4423, 2098> lagged_fibonacci4423;
1.455 +typedef random::lagged_fibonacci_01<double, 48, 9689, 5502> lagged_fibonacci9689;
1.456 +typedef random::lagged_fibonacci_01<double, 48, 19937, 9842> lagged_fibonacci19937;
1.457 +typedef random::lagged_fibonacci_01<double, 48, 23209, 13470> lagged_fibonacci23209;
1.458 +typedef random::lagged_fibonacci_01<double, 48, 44497, 21034> lagged_fibonacci44497;
1.459 +
1.460 +
1.461 +// It is possible to partially specialize uniform_01<> on lagged_fibonacci_01<>
1.462 +// to help the compiler generate efficient code. For GCC, this seems useless,
1.463 +// because GCC optimizes (x-0)/(1-0) to (x-0). This is good enough for now.
1.464 +
1.465 +} // namespace boost
1.466 +
1.467 +#endif // BOOST_RANDOM_LAGGED_FIBONACCI_HPP