1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/epoc32/include/stdapis/boost/random/mersenne_twister.hpp Tue Mar 16 16:12:26 2010 +0000
1.3 @@ -0,0 +1,302 @@
1.4 +/* boost random/mersenne_twister.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: mersenne_twister.hpp,v 1.20 2005/07/21 22:04:31 jmaurer 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_MERSENNE_TWISTER_HPP
1.20 +#define BOOST_RANDOM_MERSENNE_TWISTER_HPP
1.21 +
1.22 +#include <iostream>
1.23 +#include <algorithm> // std::copy
1.24 +#include <stdexcept>
1.25 +#include <boost/config.hpp>
1.26 +#include <boost/limits.hpp>
1.27 +#include <boost/static_assert.hpp>
1.28 +#include <boost/integer_traits.hpp>
1.29 +#include <boost/cstdint.hpp>
1.30 +#include <boost/random/linear_congruential.hpp>
1.31 +#include <boost/detail/workaround.hpp>
1.32 +#include <boost/random/detail/ptr_helper.hpp>
1.33 +
1.34 +namespace boost {
1.35 +namespace random {
1.36 +
1.37 +// http://www.math.keio.ac.jp/matumoto/emt.html
1.38 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.39 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.40 +class mersenne_twister
1.41 +{
1.42 +public:
1.43 + typedef UIntType result_type;
1.44 + BOOST_STATIC_CONSTANT(int, word_size = w);
1.45 + BOOST_STATIC_CONSTANT(int, state_size = n);
1.46 + BOOST_STATIC_CONSTANT(int, shift_size = m);
1.47 + BOOST_STATIC_CONSTANT(int, mask_bits = r);
1.48 + BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);
1.49 + BOOST_STATIC_CONSTANT(int, output_u = u);
1.50 + BOOST_STATIC_CONSTANT(int, output_s = s);
1.51 + BOOST_STATIC_CONSTANT(UIntType, output_b = b);
1.52 + BOOST_STATIC_CONSTANT(int, output_t = t);
1.53 + BOOST_STATIC_CONSTANT(UIntType, output_c = c);
1.54 + BOOST_STATIC_CONSTANT(int, output_l = l);
1.55 +
1.56 + BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
1.57 +
1.58 + mersenne_twister() { seed(); }
1.59 +
1.60 +#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
1.61 + // Work around overload resolution problem (Gennadiy E. Rozental)
1.62 + explicit mersenne_twister(const UIntType& value)
1.63 +#else
1.64 + explicit mersenne_twister(UIntType value)
1.65 +#endif
1.66 + { seed(value); }
1.67 + template<class It> mersenne_twister(It& first, It last) { seed(first,last); }
1.68 +
1.69 + template<class Generator>
1.70 + explicit mersenne_twister(Generator & gen) { seed(gen); }
1.71 +
1.72 + // compiler-generated copy ctor and assignment operator are fine
1.73 +
1.74 + void seed() { seed(UIntType(5489)); }
1.75 +
1.76 +#if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)
1.77 + // Work around overload resolution problem (Gennadiy E. Rozental)
1.78 + void seed(const UIntType& value)
1.79 +#else
1.80 + void seed(UIntType value)
1.81 +#endif
1.82 + {
1.83 + // New seeding algorithm from
1.84 + // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
1.85 + // In the previous versions, MSBs of the seed affected only MSBs of the
1.86 + // state x[].
1.87 + const UIntType mask = ~0u;
1.88 + x[0] = value & mask;
1.89 + for (i = 1; i < n; i++) {
1.90 + // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106
1.91 + x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;
1.92 + }
1.93 + }
1.94 +
1.95 + // For GCC, moving this function out-of-line prevents inlining, which may
1.96 + // reduce overall object code size. However, MSVC does not grok
1.97 + // out-of-line definitions of member function templates.
1.98 + template<class Generator>
1.99 + void seed(Generator & gen)
1.100 + {
1.101 +#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
1.102 + BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_signed);
1.103 +#endif
1.104 + // I could have used std::generate_n, but it takes "gen" by value
1.105 + for(int j = 0; j < n; j++)
1.106 + x[j] = gen();
1.107 + i = n;
1.108 + }
1.109 +
1.110 + template<class It>
1.111 + void seed(It& first, It last)
1.112 + {
1.113 + int j;
1.114 + for(j = 0; j < n && first != last; ++j, ++first)
1.115 + x[j] = *first;
1.116 + i = n;
1.117 + if(first == last && j < n)
1.118 + throw std::invalid_argument("mersenne_twister::seed");
1.119 + }
1.120 +
1.121 + result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
1.122 + result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
1.123 + {
1.124 + // avoid "left shift count >= with of type" warning
1.125 + result_type res = 0;
1.126 + for(int i = 0; i < w; ++i)
1.127 + res |= (1u << i);
1.128 + return res;
1.129 + }
1.130 +
1.131 + result_type operator()();
1.132 + static bool validation(result_type v) { return val == v; }
1.133 +
1.134 +#ifndef BOOST_NO_OPERATORS_IN_NAMESPACE
1.135 +
1.136 +#ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS
1.137 + template<class CharT, class Traits>
1.138 + friend std::basic_ostream<CharT,Traits>&
1.139 + operator<<(std::basic_ostream<CharT,Traits>& os, const mersenne_twister& mt)
1.140 + {
1.141 + for(int j = 0; j < mt.state_size; ++j)
1.142 + os << mt.compute(j) << " ";
1.143 + return os;
1.144 + }
1.145 +
1.146 + template<class CharT, class Traits>
1.147 + friend std::basic_istream<CharT,Traits>&
1.148 + operator>>(std::basic_istream<CharT,Traits>& is, mersenne_twister& mt)
1.149 + {
1.150 + for(int j = 0; j < mt.state_size; ++j)
1.151 + is >> mt.x[j] >> std::ws;
1.152 + // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template
1.153 + // value parameter "n" available from the class template scope, so use
1.154 + // the static constant with the same value
1.155 + mt.i = mt.state_size;
1.156 + return is;
1.157 + }
1.158 +#endif
1.159 +
1.160 + friend bool operator==(const mersenne_twister& x, const mersenne_twister& y)
1.161 + {
1.162 + for(int j = 0; j < state_size; ++j)
1.163 + if(x.compute(j) != y.compute(j))
1.164 + return false;
1.165 + return true;
1.166 + }
1.167 +
1.168 + friend bool operator!=(const mersenne_twister& x, const mersenne_twister& y)
1.169 + { return !(x == y); }
1.170 +#else
1.171 + // Use a member function; Streamable concept not supported.
1.172 + bool operator==(const mersenne_twister& rhs) const
1.173 + {
1.174 + for(int j = 0; j < state_size; ++j)
1.175 + if(compute(j) != rhs.compute(j))
1.176 + return false;
1.177 + return true;
1.178 + }
1.179 +
1.180 + bool operator!=(const mersenne_twister& rhs) const
1.181 + { return !(*this == rhs); }
1.182 +#endif
1.183 +
1.184 +private:
1.185 + // returns x(i-n+index), where index is in 0..n-1
1.186 + UIntType compute(unsigned int index) const
1.187 + {
1.188 + // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers
1.189 + return x[ (i + n + index) % (2*n) ];
1.190 + }
1.191 + void twist(int block);
1.192 +
1.193 + // state representation: next output is o(x(i))
1.194 + // x[0] ... x[k] x[k+1] ... x[n-1] x[n] ... x[2*n-1] represents
1.195 + // x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]
1.196 + // The goal is to always have x(i-n) ... x(i-1) available for
1.197 + // operator== and save/restore.
1.198 +
1.199 + UIntType x[2*n];
1.200 + int i;
1.201 +};
1.202 +
1.203 +#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
1.204 +// A definition is required even for integral static constants
1.205 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.206 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.207 +const bool mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range;
1.208 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.209 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.210 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size;
1.211 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.212 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.213 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size;
1.214 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.215 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.216 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits;
1.217 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.218 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.219 +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a;
1.220 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.221 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.222 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u;
1.223 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.224 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.225 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s;
1.226 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.227 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.228 +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b;
1.229 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.230 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.231 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t;
1.232 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.233 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.234 +const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c;
1.235 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.236 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.237 +const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l;
1.238 +#endif
1.239 +
1.240 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.241 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.242 +void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)
1.243 +{
1.244 + const UIntType upper_mask = (~0u) << r;
1.245 + const UIntType lower_mask = ~upper_mask;
1.246 +
1.247 + if(block == 0) {
1.248 + for(int j = n; j < 2*n; j++) {
1.249 + UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
1.250 + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
1.251 + }
1.252 + } else if (block == 1) {
1.253 + // split loop to avoid costly modulo operations
1.254 + { // extra scope for MSVC brokenness w.r.t. for scope
1.255 + for(int j = 0; j < n-m; j++) {
1.256 + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
1.257 + x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
1.258 + }
1.259 + }
1.260 +
1.261 + for(int j = n-m; j < n-1; j++) {
1.262 + UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
1.263 + x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
1.264 + }
1.265 + // last iteration
1.266 + UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
1.267 + x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
1.268 + i = 0;
1.269 + }
1.270 +}
1.271 +
1.272 +template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
1.273 + int s, UIntType b, int t, UIntType c, int l, UIntType val>
1.274 +inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type
1.275 +mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()()
1.276 +{
1.277 + if(i == n)
1.278 + twist(0);
1.279 + else if(i >= 2*n)
1.280 + twist(1);
1.281 + // Step 4
1.282 + UIntType z = x[i];
1.283 + ++i;
1.284 + z ^= (z >> u);
1.285 + z ^= ((z << s) & b);
1.286 + z ^= ((z << t) & c);
1.287 + z ^= (z >> l);
1.288 + return z;
1.289 +}
1.290 +
1.291 +} // namespace random
1.292 +
1.293 +
1.294 +typedef random::mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11,
1.295 + 7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b;
1.296 +
1.297 +// validation by experiment from mt19937.c
1.298 +typedef random::mersenne_twister<uint32_t,32,624,397,31,0x9908b0df,11,
1.299 + 7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;
1.300 +
1.301 +} // namespace boost
1.302 +
1.303 +BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)
1.304 +
1.305 +#endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP