1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/epoc32/include/stdapis/boost/random/subtract_with_carry.hpp Tue Mar 16 16:12:26 2010 +0000
1.3 @@ -0,0 +1,445 @@
1.4 +/* boost random/subtract_with_carry.hpp header file
1.5 + *
1.6 + * Copyright Jens Maurer 2002
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: subtract_with_carry.hpp,v 1.24 2005/05/21 15:57:00 dgregor Exp $
1.14 + *
1.15 + * Revision history
1.16 + * 2002-03-02 created
1.17 + */
1.18 +
1.19 +#ifndef BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
1.20 +#define BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP
1.21 +
1.22 +#include <cmath>
1.23 +#include <iostream>
1.24 +#include <algorithm> // std::equal
1.25 +#include <stdexcept>
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/static_assert.hpp>
1.31 +#include <boost/detail/workaround.hpp>
1.32 +#include <boost/random/linear_congruential.hpp>
1.33 +
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_SWC_01
1.40 +#endif
1.41 +
1.42 +#if defined(__APPLE_CC__) && defined(__GNUC__) && (__GNUC__ == 3) && (__GNUC_MINOR__ <= 3)
1.43 +# define BOOST_RANDOM_EXTRACT_SWC_01
1.44 +#endif
1.45 +
1.46 +# ifdef BOOST_RANDOM_EXTRACT_SWC_01
1.47 +namespace detail
1.48 +{
1.49 + template <class IStream, class SubtractWithCarry, class RealType>
1.50 + void extract_subtract_with_carry_01(
1.51 + IStream& is
1.52 + , SubtractWithCarry& f
1.53 + , RealType& carry
1.54 + , RealType* x
1.55 + , RealType modulus)
1.56 + {
1.57 + RealType value;
1.58 + for(unsigned int j = 0; j < f.long_lag; ++j) {
1.59 + is >> value >> std::ws;
1.60 + x[j] = value / modulus;
1.61 + }
1.62 + is >> value >> std::ws;
1.63 + carry = value / modulus;
1.64 + }
1.65 +}
1.66 +# endif
1.67 +// subtract-with-carry generator
1.68 +// Marsaglia and Zaman
1.69 +
1.70 +template<class IntType, IntType m, unsigned int s, unsigned int r,
1.71 + IntType val>
1.72 +class subtract_with_carry
1.73 +{
1.74 +public:
1.75 + typedef IntType result_type;
1.76 + BOOST_STATIC_CONSTANT(bool, has_fixed_range = true);
1.77 + BOOST_STATIC_CONSTANT(result_type, min_value = 0);
1.78 + BOOST_STATIC_CONSTANT(result_type, max_value = m-1);
1.79 + BOOST_STATIC_CONSTANT(result_type, modulus = m);
1.80 + BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);
1.81 + BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);
1.82 +
1.83 + subtract_with_carry() {
1.84 + // MSVC fails BOOST_STATIC_ASSERT with std::numeric_limits at class scope
1.85 +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
1.86 + BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_signed);
1.87 + BOOST_STATIC_ASSERT(std::numeric_limits<result_type>::is_integer);
1.88 +#endif
1.89 + seed();
1.90 + }
1.91 + explicit subtract_with_carry(uint32_t value) { seed(value); }
1.92 + template<class Generator>
1.93 + explicit subtract_with_carry(Generator & gen) { seed(gen); }
1.94 + template<class It> subtract_with_carry(It& first, It last) { seed(first,last); }
1.95 +
1.96 + // compiler-generated copy ctor and assignment operator are fine
1.97 +
1.98 + void seed(uint32_t value = 19780503u)
1.99 + {
1.100 + random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> intgen(value);
1.101 + seed(intgen);
1.102 + }
1.103 +
1.104 + // For GCC, moving this function out-of-line prevents inlining, which may
1.105 + // reduce overall object code size. However, MSVC does not grok
1.106 + // out-of-line template member functions.
1.107 + template<class Generator>
1.108 + void seed(Generator & gen)
1.109 + {
1.110 + // I could have used std::generate_n, but it takes "gen" by value
1.111 + for(unsigned int j = 0; j < long_lag; ++j)
1.112 + x[j] = gen() % modulus;
1.113 + carry = (x[long_lag-1] == 0);
1.114 + k = 0;
1.115 + }
1.116 +
1.117 + template<class It>
1.118 + void seed(It& first, It last)
1.119 + {
1.120 + unsigned int j;
1.121 + for(j = 0; j < long_lag && first != last; ++j, ++first)
1.122 + x[j] = *first % modulus;
1.123 + if(first == last && j < long_lag)
1.124 + throw std::invalid_argument("subtract_with_carry::seed");
1.125 + carry = (x[long_lag-1] == 0);
1.126 + k = 0;
1.127 + }
1.128 +
1.129 + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return min_value; }
1.130 + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return max_value; }
1.131 +
1.132 + result_type operator()()
1.133 + {
1.134 + int short_index = k - short_lag;
1.135 + if(short_index < 0)
1.136 + short_index += long_lag;
1.137 + IntType delta;
1.138 + if (x[short_index] >= x[k] + carry) {
1.139 + // x(n) >= 0
1.140 + delta = x[short_index] - (x[k] + carry);
1.141 + carry = 0;
1.142 + } else {
1.143 + // x(n) < 0
1.144 + delta = modulus - x[k] - carry + x[short_index];
1.145 + carry = 1;
1.146 + }
1.147 + x[k] = delta;
1.148 + ++k;
1.149 + if(k >= long_lag)
1.150 + k = 0;
1.151 + return delta;
1.152 + }
1.153 +
1.154 +public:
1.155 + static bool validation(result_type x) { return x == val; }
1.156 +
1.157 +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
1.158 +
1.159 +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
1.160 + template<class CharT, class Traits>
1.161 + friend std::basic_ostream<CharT,Traits>&
1.162 + operator<<(std::basic_ostream<CharT,Traits>& os,
1.163 + const subtract_with_carry& f)
1.164 + {
1.165 + for(unsigned int j = 0; j < f.long_lag; ++j)
1.166 + os << f.compute(j) << " ";
1.167 + os << f.carry << " ";
1.168 + return os;
1.169 + }
1.170 +
1.171 + template<class CharT, class Traits>
1.172 + friend std::basic_istream<CharT,Traits>&
1.173 + operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry& f)
1.174 + {
1.175 + for(unsigned int j = 0; j < f.long_lag; ++j)
1.176 + is >> f.x[j] >> std::ws;
1.177 + is >> f.carry >> std::ws;
1.178 + f.k = 0;
1.179 + return is;
1.180 + }
1.181 +#endif
1.182 +
1.183 + friend bool operator==(const subtract_with_carry& x, const subtract_with_carry& y)
1.184 + {
1.185 + for(unsigned int j = 0; j < r; ++j)
1.186 + if(x.compute(j) != y.compute(j))
1.187 + return false;
1.188 + return true;
1.189 + }
1.190 +
1.191 + friend bool operator!=(const subtract_with_carry& x, const subtract_with_carry& y)
1.192 + { return !(x == y); }
1.193 +#else
1.194 + // Use a member function; Streamable concept not supported.
1.195 + bool operator==(const subtract_with_carry& rhs) const
1.196 + {
1.197 + for(unsigned int j = 0; j < r; ++j)
1.198 + if(compute(j) != rhs.compute(j))
1.199 + return false;
1.200 + return true;
1.201 + }
1.202 +
1.203 + bool operator!=(const subtract_with_carry& rhs) const
1.204 + { return !(*this == rhs); }
1.205 +#endif
1.206 +
1.207 +private:
1.208 + // returns x(i-r+index), where index is in 0..r-1
1.209 + IntType compute(unsigned int index) const
1.210 + {
1.211 + return x[(k+index) % long_lag];
1.212 + }
1.213 +
1.214 + // state representation; next output (state) is x(i)
1.215 + // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents
1.216 + // x(i-k) ... x(i) x(i+1) ... x(i-k+long_lag-1)
1.217 + // speed: base: 20-25 nsec
1.218 + // ranlux_4: 230 nsec, ranlux_7: 430 nsec, ranlux_14: 810 nsec
1.219 + // This state representation makes operator== and save/restore more
1.220 + // difficult, because we've already computed "too much" and thus
1.221 + // have to undo some steps to get at x(i-r) etc.
1.222 +
1.223 + // state representation: next output (state) is x(i)
1.224 + // x[0] ... x[k] x[k+1] ... x[long_lag-1] represents
1.225 + // x(i-k) ... x(i) x(i-long_lag+1) ... x(i-k-1)
1.226 + // speed: base 28 nsec
1.227 + // ranlux_4: 370 nsec, ranlux_7: 688 nsec, ranlux_14: 1343 nsec
1.228 + IntType x[long_lag];
1.229 + unsigned int k;
1.230 + int carry;
1.231 +};
1.232 +
1.233 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.234 +// A definition is required even for integral static constants
1.235 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.236 +const bool subtract_with_carry<IntType, m, s, r, val>::has_fixed_range;
1.237 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.238 +const IntType subtract_with_carry<IntType, m, s, r, val>::min_value;
1.239 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.240 +const IntType subtract_with_carry<IntType, m, s, r, val>::max_value;
1.241 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.242 +const IntType subtract_with_carry<IntType, m, s, r, val>::modulus;
1.243 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.244 +const unsigned int subtract_with_carry<IntType, m, s, r, val>::long_lag;
1.245 +template<class IntType, IntType m, unsigned int s, unsigned int r, IntType val>
1.246 +const unsigned int subtract_with_carry<IntType, m, s, r, val>::short_lag;
1.247 +#endif
1.248 +
1.249 +
1.250 +// use a floating-point representation to produce values in [0..1)
1.251 +template<class RealType, int w, unsigned int s, unsigned int r, int val=0>
1.252 +class subtract_with_carry_01
1.253 +{
1.254 +public:
1.255 + typedef RealType result_type;
1.256 + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
1.257 + BOOST_STATIC_CONSTANT(int, word_size = w);
1.258 + BOOST_STATIC_CONSTANT(unsigned int, long_lag = r);
1.259 + BOOST_STATIC_CONSTANT(unsigned int, short_lag = s);
1.260 +
1.261 +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
1.262 + BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_integer);
1.263 +#endif
1.264 +
1.265 + subtract_with_carry_01() { init_modulus(); seed(); }
1.266 + explicit subtract_with_carry_01(uint32_t value)
1.267 + { init_modulus(); seed(value); }
1.268 + template<class It> subtract_with_carry_01(It& first, It last)
1.269 + { init_modulus(); seed(first,last); }
1.270 +
1.271 +private:
1.272 + void init_modulus()
1.273 + {
1.274 +#ifndef BOOST_NO_STDC_NAMESPACE
1.275 + // allow for Koenig lookup
1.276 + using std::pow;
1.277 +#endif
1.278 + _modulus = pow(RealType(2), word_size);
1.279 + }
1.280 +
1.281 +public:
1.282 + // compiler-generated copy ctor and assignment operator are fine
1.283 +
1.284 + void seed(uint32_t value = 19780503u)
1.285 + {
1.286 +#ifndef BOOST_NO_STDC_NAMESPACE
1.287 + // allow for Koenig lookup
1.288 + using std::fmod;
1.289 +#endif
1.290 + random::linear_congruential<int32_t, 40014, 0, 2147483563, 0> gen(value);
1.291 + unsigned long array[(w+31)/32 * long_lag];
1.292 + for(unsigned int j = 0; j < sizeof(array)/sizeof(unsigned long); ++j)
1.293 + array[j] = gen();
1.294 + unsigned long * start = array;
1.295 + seed(start, array + sizeof(array)/sizeof(unsigned long));
1.296 + }
1.297 +
1.298 + template<class It>
1.299 + void seed(It& first, It last)
1.300 + {
1.301 +#ifndef BOOST_NO_STDC_NAMESPACE
1.302 + // allow for Koenig lookup
1.303 + using std::fmod;
1.304 + using std::pow;
1.305 +#endif
1.306 + unsigned long mask = ~((~0u) << (w%32)); // now lowest (w%32) bits set
1.307 + RealType two32 = pow(RealType(2), 32);
1.308 + unsigned int j;
1.309 + for(j = 0; j < long_lag && first != last; ++j) {
1.310 + x[j] = RealType(0);
1.311 + for(int i = 0; i < w/32 && first != last; ++i, ++first)
1.312 + x[j] += *first / pow(two32,i+1);
1.313 + if(first != last && mask != 0) {
1.314 + x[j] += fmod((*first & mask) / _modulus, RealType(1));
1.315 + ++first;
1.316 + }
1.317 + }
1.318 + if(first == last && j < long_lag)
1.319 + throw std::invalid_argument("subtract_with_carry_01::seed");
1.320 + carry = (x[long_lag-1] ? 0 : 1 / _modulus);
1.321 + k = 0;
1.322 + }
1.323 +
1.324 + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(0); }
1.325 + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return result_type(1); }
1.326 +
1.327 + result_type operator()()
1.328 + {
1.329 + int short_index = k - short_lag;
1.330 + if(short_index < 0)
1.331 + short_index += long_lag;
1.332 + RealType delta = x[short_index] - x[k] - carry;
1.333 + if(delta < 0) {
1.334 + delta += RealType(1);
1.335 + carry = RealType(1)/_modulus;
1.336 + } else {
1.337 + carry = 0;
1.338 + }
1.339 + x[k] = delta;
1.340 + ++k;
1.341 + if(k >= long_lag)
1.342 + k = 0;
1.343 + return delta;
1.344 + }
1.345 +
1.346 + static bool validation(result_type x)
1.347 + { return x == val/pow(RealType(2), word_size); }
1.348 +
1.349 +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
1.350 +
1.351 +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
1.352 + template<class CharT, class Traits>
1.353 + friend std::basic_ostream<CharT,Traits>&
1.354 + operator<<(std::basic_ostream<CharT,Traits>& os,
1.355 + const subtract_with_carry_01& f)
1.356 + {
1.357 +#ifndef BOOST_NO_STDC_NAMESPACE
1.358 + // allow for Koenig lookup
1.359 + using std::pow;
1.360 +#endif
1.361 + std::ios_base::fmtflags oldflags = os.flags(os.dec | os.fixed | os.left);
1.362 + for(unsigned int j = 0; j < f.long_lag; ++j)
1.363 + os << (f.compute(j) * f._modulus) << " ";
1.364 + os << (f.carry * f._modulus);
1.365 + os.flags(oldflags);
1.366 + return os;
1.367 + }
1.368 +
1.369 + template<class CharT, class Traits>
1.370 + friend std::basic_istream<CharT,Traits>&
1.371 + operator>>(std::basic_istream<CharT,Traits>& is, subtract_with_carry_01& f)
1.372 + {
1.373 +# ifdef BOOST_RANDOM_EXTRACT_SWC_01
1.374 + detail::extract_subtract_with_carry_01(is, f, f.carry, f.x, f._modulus);
1.375 +# else
1.376 + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template type
1.377 + // parameter "RealType" available from the class template scope, so use
1.378 + // the member typedef
1.379 + typename subtract_with_carry_01::result_type value;
1.380 + for(unsigned int j = 0; j < long_lag; ++j) {
1.381 + is >> value >> std::ws;
1.382 + f.x[j] = value / f._modulus;
1.383 + }
1.384 + is >> value >> std::ws;
1.385 + f.carry = value / f._modulus;
1.386 +# endif
1.387 + f.k = 0;
1.388 + return is;
1.389 + }
1.390 +#endif
1.391 +
1.392 + friend bool operator==(const subtract_with_carry_01& x,
1.393 + const subtract_with_carry_01& y)
1.394 + {
1.395 + for(unsigned int j = 0; j < r; ++j)
1.396 + if(x.compute(j) != y.compute(j))
1.397 + return false;
1.398 + return true;
1.399 + }
1.400 +
1.401 + friend bool operator!=(const subtract_with_carry_01& x,
1.402 + const subtract_with_carry_01& y)
1.403 + { return !(x == y); }
1.404 +#else
1.405 + // Use a member function; Streamable concept not supported.
1.406 + bool operator==(const subtract_with_carry_01& rhs) const
1.407 + {
1.408 + for(unsigned int j = 0; j < r; ++j)
1.409 + if(compute(j) != rhs.compute(j))
1.410 + return false;
1.411 + return true;
1.412 + }
1.413 +
1.414 + bool operator!=(const subtract_with_carry_01& rhs) const
1.415 + { return !(*this == rhs); }
1.416 +#endif
1.417 +
1.418 +private:
1.419 + RealType compute(unsigned int index) const;
1.420 + unsigned int k;
1.421 + RealType carry;
1.422 + RealType x[long_lag];
1.423 + RealType _modulus;
1.424 +};
1.425 +
1.426 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.427 +// A definition is required even for integral static constants
1.428 +template<class RealType, int w, unsigned int s, unsigned int r, int val>
1.429 +const bool subtract_with_carry_01<RealType, w, s, r, val>::has_fixed_range;
1.430 +template<class RealType, int w, unsigned int s, unsigned int r, int val>
1.431 +const int subtract_with_carry_01<RealType, w, s, r, val>::word_size;
1.432 +template<class RealType, int w, unsigned int s, unsigned int r, int val>
1.433 +const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::long_lag;
1.434 +template<class RealType, int w, unsigned int s, unsigned int r, int val>
1.435 +const unsigned int subtract_with_carry_01<RealType, w, s, r, val>::short_lag;
1.436 +#endif
1.437 +
1.438 +template<class RealType, int w, unsigned int s, unsigned int r, int val>
1.439 +RealType subtract_with_carry_01<RealType, w, s, r, val>::compute(unsigned int index) const
1.440 +{
1.441 + return x[(k+index) % long_lag];
1.442 +}
1.443 +
1.444 +
1.445 +} // namespace random
1.446 +} // namespace boost
1.447 +
1.448 +#endif // BOOST_RANDOM_SUBTRACT_WITH_CARRY_HPP