Random.cpp Source File

00001 /**
00002  * \file Random.cpp
00003  * \brief Implementation code for RandomLib.
00004  *
00005  * Written by <a href="http://charles.karney.info/">Charles Karney</a>
00006  * <charles@karney.com> and licensed under the LGPL.  For more
00007  * information, see http://charles.karney.info/random/
00008  **********************************************************************/
00009 
00010 /**
00011  * \file RandomMixer.cpp
00012  * \brief Code for MixerMT0, MixerMT1, MixerSFMT.
00013  *
00014  * MixerMT0 is adapted from MT19937 (init_by_array) and MT19937_64
00015  * (init_by_array64) by Makoto Matsumoto and Takuji Nishimura.  See
00016  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c and
00017  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
00018  *
00019  * MixerMT1 contains modifications to MixerMT0 by Charles Karney to
00020  * correct defects in MixerMT0.  This is described in W. E. Brown,
00021  * M. Fischler, J. Kowalkowski, M. Paterno, Random Number Generation in C++0X:
00022  * A Comprehensive Proposal, version 3, Sept 2006, Sec. 26.4.7.1,
00023  * http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2006/n2079.pdf
00024  * This has been replaced in the C++0X proposal by MixerSFMT.
00025  *
00026  * RandonMixer2 is adapted from SFMT19937's init_by_array Mutsuo Saito given in
00027  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.2.tar.gz and
00028  * is part of the C++0X proposal; see P. Becker, Working Draft, Standard for
00029  * Programming Language C++, Oct. 2007, Sec. 26.4.7.1,
00030  * http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf
00031  *
00032  * The adaption to the C++ is by <a href="http://charles.karney.info/">
00033  * Charles Karney</a> <charles@karney.com> and licensed under the LGPL.  For
00034  * more information, see http://charles.karney.info/random/
00035  **********************************************************************/
00036 
00037 /**
00038  * \file RandomMixer.cpp
00039  * \brief Code for MT19937<T> and SFMT19937<T>.
00040  *
00041  * MT19937<T> is adapted from MT19937 and MT19937_64 by Makoto Matsumoto and
00042  * Takuji Nishimura.  See
00043  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c and
00044  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
00045  *
00046  * The code for stepping MT19937 backwards is adapted (and simplified) from
00047  * revrand() by Katsumi Hagita.  See
00048  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/REVmt19937b.f
00049  *
00050  * SFMT19937<T> is adapted from SFMT19937 Mutsuo Saito given in
00051  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/M062821.pdf and
00052  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.2.tar.gz
00053  *
00054  * The code for stepping SFMT19937 backwards is by Charles Karney.
00055  *
00056  * The adaption to the C++ is by <a href="http://charles.karney.info/">
00057  * Charles Karney</a> <charles@karney.com> and licensed under the LGPL.  For
00058  * more information, see http://charles.karney.info/random/
00059  **********************************************************************/
00060 
00061 #include "RandomLib/Random.hpp"
00062 
00063 #include <fstream>              // For SeedWord reading /dev/urandom
00064 #include <ctime>                // For SeedWord calling time()
00065 #include <sstream>              // For formatting in Write32/Read32
00066 #include <iomanip>              // For formatting in Write32/Read32
00067 #if !WINDOWS
00068 #include <sys/time.h>           // For SeedWord calling gettimeofday
00069 #include <unistd.h>             // For SeedWord calling getpid(), gethostid()
00070 #else
00071 #include <windows.h>            // For SeedWord calling high prec timer
00072 #include <winbase.h>
00073 #include <process.h>            // For SeedWord calling getpid()
00074 #define getpid _getpid
00075 #define gmtime_r(t,g) gmtime_s(g,t)
00076 #endif
00077 
00078 #if WINDOWS || defined(__CYGWIN__)
00079 #define strtoull strtoul
00080 #endif
00081 
00082 #define RANDOM_CPP "$Id: Random.cpp 6489 2008-11-10 21:37:47Z ckarney $"
00083 
00084 RCSID_DECL(RANDOM_CPP);
00085 RCSID_DECL(RANDOMTYPE_HPP);
00086 RCSID_DECL(RANDOMSEED_HPP);
00087 RCSID_DECL(RANDOMENGINE_HPP);
00088 RCSID_DECL(RANDOMMIXER_HPP);
00089 RCSID_DECL(RANDOMALGORITHM_HPP);
00090 RCSID_DECL(RANDOMPOWER2_HPP);
00091 RCSID_DECL(RANDOMCANONICAL_HPP);
00092 
00093 namespace RandomLib {
00094 
00095   // RandomType implementation
00096 
00097   template<>
00098   void Random_u32::Write32(std::ostream& os, bool bin, int& cnt,
00099                            Random_u32::type x)
00100     throw(std::ios::failure) {
00101     if (bin) {
00102       unsigned char buf[4];
00103       // Use network order -- most significant byte first
00104       buf[3] = (unsigned char)(x);
00105       buf[2] = (unsigned char)(x >>= 8);
00106       buf[1] = (unsigned char)(x >>= 8);
00107       buf[0] = (unsigned char)(x >>= 8);
00108       os.write(reinterpret_cast<const char *>(buf), 4);
00109     } else {
00110       const int longsperline = 72/9;
00111       // Use hexadecimal to minimize storage together with stringstream to
00112       // isolate the effect of changing the base.
00113       std::ostringstream str;
00114       // No spacing before or after
00115       if (cnt > 0)
00116         // Newline every longsperline longs
00117         str << (cnt % longsperline ? ' ' : '\n');
00118       str << std::hex << x;
00119       os << str.str();
00120       ++cnt;
00121     }
00122   }
00123 
00124   template<>
00125   void Random_u32::Read32(std::istream& is, bool bin, Random_u32::type& x)
00126     throw(std::ios::failure) {
00127     if (bin) {
00128       unsigned char buf[4];
00129       is.read(reinterpret_cast<char *>(buf), 4);
00130       // Use network order -- most significant byte first
00131       x = Random_u32::type(buf[0]) << 24 | Random_u32::type(buf[1]) << 16 |
00132         Random_u32::type(buf[2]) << 8 | Random_u32::type(buf[3]);
00133     } else {
00134       std::string s;
00135       is >> std::ws >> s;
00136       // Use hexadecimal to minimize storage together with stringstream to
00137       // isolate the effect of changing the base.
00138       std::istringstream str(s);
00139       str >> std::hex >> x;
00140     }
00141     x &= Random_u32::mask;
00142   }
00143 
00144   template<>
00145   void Random_u64::Write32(std::ostream& os, bool bin, int& cnt,
00146                            Random_u64::type x)
00147     throw(std::ios::failure) {
00148     Random_u32::Write32(os, bin, cnt, Random_u32::cast(x >> 32));
00149     Random_u32::Write32(os, bin, cnt, Random_u32::cast(x      ));
00150   }
00151 
00152   template<>
00153   void Random_u64::Read32(std::istream& is, bool bin, Random_u64::type& x)
00154     throw(std::ios::failure) {
00155     Random_u32::type t;
00156     Random_u32::Read32(is, bin, t);
00157     x = Random_u64::type(t) << 32;
00158     Random_u32::Read32(is, bin, t);
00159     x |= Random_u64::type(t);
00160   }
00161 
00162   // RandomSeed implementation
00163 
00164   RandomSeed::seed_type RandomSeed::SeedWord() {
00165     // Check that the assumptions made about the capabilities of the number
00166     // system are valid.
00167     STATIC_ASSERT(std::numeric_limits<seed_type>::radix == 2 &&
00168                   !std::numeric_limits<seed_type>::is_signed &&
00169                   std::numeric_limits<seed_type>::digits >= 32,
00170                   "seed_type is a bad type");
00171     u32::type t = 0;
00172     // Linux has /dev/urandom to initialize the seed randomly.  (Use
00173     // /dev/urandom instead of /dev/random because it does not block.)
00174     {
00175       std::ifstream f("/dev/urandom", std::ios::binary | std::ios::in);
00176       if (f.good()) {
00177         // Read 32 bits from /dev/urandom
00178         f.read(reinterpret_cast<char *>(&t), sizeof(t));
00179       }
00180     }
00181     std::vector<seed_type> v = SeedVector();
00182     for (size_t i = v.size(); i--;)
00183       u32::CheckSum(v[i], t);
00184     return seed_t::cast(t);
00185   }
00186 
00187   std::vector<RandomSeed::seed_type> RandomSeed::SeedVector() {
00188     std::vector<seed_type> v;
00189     {
00190       // fine-grained timer
00191 #if !WINDOWS
00192       timeval tv;
00193       if (gettimeofday(&tv, 0) == 0)
00194         v.push_back(seed_t::cast(tv.tv_usec));
00195 #else
00196       LARGE_INTEGER taux;
00197       if (QueryPerformanceCounter((LARGE_INTEGER *)&taux)) {
00198         v.push_back(seed_t::cast(taux.LowPart));
00199         v.push_back(seed_t::cast(taux.HighPart));
00200       }
00201 #endif
00202     }
00203     // seconds
00204     const time_t tim = std::time(0);
00205     v.push_back(seed_t::cast(seed_type(tim)));
00206     // PID
00207     v.push_back(seed_t::cast(getpid()));
00208 #if !WINDOWS
00209     // host ID
00210     v.push_back(seed_t::cast(gethostid()));
00211 #endif
00212     {
00213       // year
00214       tm gt;
00215       gmtime_r(&tim, &gt);
00216       v.push_back((seed_type(1900) + seed_t::cast(gt.tm_year)));
00217     }
00218     std::transform(v.begin(), v.end(), v.begin(), seed_t::cast<seed_type>);
00219     return v;
00220   }
00221 
00222   std::vector<RandomSeed::seed_type>
00223   RandomSeed::StringToVector(const std::string& s)
00224     throw(std::bad_alloc) {
00225     std::vector<seed_type> v(0);
00226     const char* c = s.c_str();
00227     char* q;
00228     std::string::size_type p = 0;
00229     while (true) {
00230       p = s.find_first_of("0123456789", p);
00231       if (p == std::string::npos)
00232         break;
00233       v.push_back(seed_t::cast(std::strtoull(c + p, &q, 0)));
00234       p = q - c;
00235     }
00236     return v;
00237   }
00238 
00239   // RandomEngine implementation
00240 
00241   template<class Algorithm, class Mixer>
00242   void RandomEngine<Algorithm, Mixer>::Init() throw() {
00243     // On exit we have _ptr == N.
00244 
00245     STATIC_ASSERT(std::numeric_limits<typename mixer_t::type>::radix == 2 &&
00246                   !std::numeric_limits<typename mixer_t::type>::is_signed &&
00247                   std::numeric_limits<typename mixer_t::type>::digits >=
00248                   int(mixer_t::width),
00249                   "mixer_type is a bad type");
00250 
00251     STATIC_ASSERT(std::numeric_limits<result_type>::radix == 2 &&
00252                   !std::numeric_limits<result_type>::is_signed &&
00253                   std::numeric_limits<result_type>::digits >= width,
00254                   "engine_type is a bad type");
00255 
00256     STATIC_ASSERT(mixer_t::width == 32 || mixer_t::width == 64,
00257                   "Mixer width must be 32 or 64");
00258 
00259     STATIC_ASSERT(width == 32 || width == 64,
00260                   "Algorithm width must be 32 or 64");
00261 
00262     // If the bit-widths are the same then the data sizes must be the same.
00263     STATIC_ASSERT(!(mixer_t::width == width) ||
00264                   sizeof(_stateu) == sizeof(_state),
00265                   "Same bit-widths but different storage");
00266 
00267     // Repacking assumes that narrower data type is at least as wasteful than
00268     // the broader one.
00269     STATIC_ASSERT(!(mixer_t::width < width) ||
00270                   sizeof(_stateu) >= sizeof(_state),
00271                   "Narrow data type uses less storage");
00272 
00273     STATIC_ASSERT(!(mixer_t::width > width) ||
00274                   sizeof(_stateu) <= sizeof(_state),
00275                   "Narrow data type uses less storage");
00276 
00277     // Require that _statev and _state are aligned since no repacking is done
00278     // when calling Transition
00279     STATIC_ASSERT(sizeof(_statev) == sizeof(_state),
00280                   "Storage mismatch with internal engine data type");
00281 
00282     // Convert the seed into state
00283     Mixer::SeedToState(_seed, _stateu, NU);
00284 
00285     // Pack into _state
00286     if (mixer_t::width < width) {
00287       for (size_t i = 0; i < N; ++i)
00288         // Assume 2:1 LSB packing
00289         _state[i] = result_type(_stateu[2*i]) |
00290           result_type(_stateu[2*i + 1]) <<
00291           (mixer_t::width < width ? mixer_t::width : 0);
00292     } else if (mixer_t::width > width) {
00293       for (size_t i = N; i--;)
00294         // Assume 1:2 LSB packing
00295         _state[i] = result_t::cast(_stateu[i>>1] >> width * (i&1u));
00296     } // Otherwise the union takes care of it
00297 
00298     Algorithm::NormalizeState(_state);
00299 
00300     _rounds = -1;
00301     _ptr = N;
00302   }
00303 
00304   template<class Algorithm, class Mixer> Random_u32::type
00305   RandomEngine<Algorithm, Mixer>::Check(u64::type v, u32::type e,
00306                                          u32::type m) const
00307     throw(std::out_of_range) {
00308     if (v != version)
00309       throw std::out_of_range(Name() + ": Unknown version");
00310     if (e != Algorithm::version)
00311       throw std::out_of_range(Name() + ": Algorithm mismatch");
00312     if (m != Mixer::version)
00313       throw std::out_of_range(Name() + ": Mixer mismatch");
00314     u32::type check = 0;
00315     u64::CheckSum(v, check);
00316     u32::CheckSum(e, check);
00317     u32::CheckSum(m, check);
00318     u32::CheckSum(u32::type(_seed.size()), check);
00319     for (std::vector<seed_type>::const_iterator n = _seed.begin();
00320          n != _seed.end(); ++n) {
00321       if (*n != seed_t::cast(*n))
00322         throw std::out_of_range(Name() + ": Illegal seed value");
00323       u32::CheckSum(u32::type(*n), check);
00324     }
00325     u32::CheckSum(_ptr, check);
00326     if (_stride == 0 || _stride > UNINIT/2)
00327       throw std::out_of_range(Name() + ": Invalid stride");
00328     u32::CheckSum(_stride, check);
00329     if (_ptr != UNINIT) {
00330       if (_ptr >= N + _stride)
00331         throw std::out_of_range(Name() + ": Invalid pointer");
00332       u64::CheckSum(_rounds, check);
00333       Algorithm::CheckState(_state, check);
00334     }
00335     return check;
00336   }
00337 
00338   template<typename Algorithm, typename Mixer>
00339   RandomEngine<Algorithm, Mixer>::RandomEngine(std::istream& is, bool bin)
00340     throw(std::ios::failure, std::out_of_range, std::bad_alloc) {
00341     u64::type versionr;
00342     u32::type versione, versionm, t;
00343     u64::Read32(is, bin, versionr);
00344     u32::Read32(is, bin, versione);
00345     u32::Read32(is, bin, versionm);
00346     u32::Read32(is, bin, t);
00347     _seed.resize(size_t(t));
00348     for (std::vector<seed_type>::iterator n = _seed.begin();
00349          n != _seed.end(); ++n) {
00350       u32::Read32(is, bin, t);
00351       *n = seed_type(t);
00352     }
00353     u32::Read32(is, bin, t);
00354     // Don't need to worry about sign extension because _ptr is unsigned.
00355     _ptr = unsigned(t);
00356     u32::Read32(is, bin, t);
00357     _stride = unsigned(t);
00358     if (_ptr != UNINIT) {
00359       u64::type p;
00360       u64::Read32(is, bin, p);
00361       _rounds = (long long)(p);
00362       // Sign extension in case long long is bigger than 64 bits.
00363       _rounds <<= 63 - std::numeric_limits<long long>::digits;
00364       _rounds >>= 63 - std::numeric_limits<long long>::digits;
00365       for (unsigned i = 0; i < N; ++i)
00366         result_t::Read32(is, bin, _state[i]);
00367     }
00368     u32::Read32(is, bin, t);
00369     if (t != Check(versionr, versione, versionm))
00370       throw std::out_of_range(Name() + ": Checksum failure");
00371   }
00372 
00373   template<typename Algorithm, typename Mixer>
00374   void RandomEngine<Algorithm, Mixer>::Save(std::ostream& os,
00375                                             bool bin) const
00376     throw(std::ios::failure) {
00377     u32::type check = Check(version, Algorithm::version, Mixer::version);
00378     int c = 0;
00379     u64::Write32(os, bin, c, version);
00380     u32::Write32(os, bin, c, Algorithm::version);
00381     u32::Write32(os, bin, c, Mixer::version);
00382     u32::Write32(os, bin, c, u32::type(_seed.size()));
00383     for (std::vector<seed_type>::const_iterator n = _seed.begin();
00384          n != _seed.end(); ++n)
00385       u32::Write32(os, bin, c, u32::type(*n));
00386     u32::Write32(os, bin, c, _ptr);
00387     u32::Write32(os, bin, c, _stride);
00388     if (_ptr != UNINIT) {
00389       u64::Write32(os, bin, c, u64::type(_rounds));
00390       for (unsigned i = 0; i < N; ++i)
00391         result_t::Write32(os, bin, c, _state[i]);
00392     }
00393     u32::Write32(os, bin, c, check);
00394   }
00395 
00396   template<typename Algorithm, typename Mixer>
00397   void RandomEngine<Algorithm, Mixer>::StepCount(long long n) throw() {
00398     // On exit we have 0 <= _ptr <= N.
00399     if (_ptr == UNINIT)
00400       Init();
00401     const long long ncount = n + Count(); // new Count()
00402     long long nrounds = ncount / N;
00403     int nptr = int(ncount - nrounds * N);
00404     // We pick _ptr = N or _ptr = 0 depending on which choice involves the
00405     // least work.  We thus avoid doing one (potentially unneeded) call to
00406     // Transition.
00407     if (nptr < 0) {
00408       --nrounds;
00409       nptr += N;
00410     } else if (nptr == 0 && nrounds > _rounds) {
00411       nptr = N;
00412       --nrounds;
00413     }
00414     if (nrounds != _rounds)
00415       Algorithm::Transition(nrounds - _rounds, _statev);
00416     _rounds = nrounds;
00417     _ptr = nptr;
00418   }
00419 
00420   template<typename Algorithm, typename Mixer>
00421   void RandomEngine<Algorithm, Mixer>::SelfTest() {
00422     RandomEngine g(std::vector<seed_type>(0));
00423     g.SetCount(10000-1);
00424     result_type x = g();
00425     if (SelfTestResult(0) && x != SelfTestResult(1))
00426       throw std::out_of_range(Name() + ": Incorrect result with seed " +
00427                               g.SeedString());
00428     seed_type s[] = {0x1234U, 0x5678U, 0x9abcU, 0xdef0U};
00429     //    seed_type s[] = {1, 2, 3, 4};
00430     g.Reseed(s, s+4);
00431     g.StepCount(-20000);
00432     std::string save;
00433     {
00434       std::ostringstream stream;
00435       stream << g << std::endl;
00436       save = stream.str();
00437     }
00438     g.Reset();
00439     {
00440       std::istringstream stream(save);
00441       stream >> g;
00442     }
00443     g.SetCount(10000);
00444     {
00445       std::ostringstream stream;
00446       g.Save(stream, true);
00447       save = stream.str();
00448     }
00449     {
00450       std::istringstream stream(save);
00451       RandomEngine h(std::vector<seed_type>(0));
00452       h.Load(stream, true);
00453       h.SetCount(1000000-1);
00454       x = h();
00455       if (SelfTestResult(0) && x != SelfTestResult(2))
00456         throw std::out_of_range(Name() + ": Incorrect result with seed " +
00457                                 h.SeedString());
00458       g.SetCount(1000000);
00459       if (h != g)
00460         throw std::out_of_range(Name() + ": Comparison failure");
00461     }
00462   }
00463 
00464   template<> Random_u32::type
00465   RandomEngine<MT19937<Random_u32>, MixerMT0<Random_u32> >::
00466   SelfTestResult(unsigned i)
00467     throw() {
00468     return i == 0 ? 1 :
00469       i == 1 ? 4123659995UL : 3016432305UL;
00470   }
00471 
00472   template<> Random_u64::type
00473   RandomEngine<MT19937<Random_u64>, MixerMT0<Random_u64> >::
00474   SelfTestResult(unsigned i)
00475     throw() {
00476     return i == 0 ? 1 :
00477       i == 1 ? 9981545732273789042ULL : 1384037754719008581ULL;
00478   }
00479 
00480   template<> Random_u32::type
00481   RandomEngine<MT19937<Random_u32>, MixerMT1<Random_u32> >::
00482   SelfTestResult(unsigned i)
00483     throw() {
00484     return i == 0 ? 1 :
00485       i == 1 ? 4123659995UL : 2924523180UL;
00486   }
00487 
00488   template<> Random_u64::type
00489   RandomEngine<MT19937<Random_u64>, MixerMT1<Random_u64> >::
00490   SelfTestResult(unsigned i)
00491     throw() {
00492     return i == 0 ? 1 :
00493       i == 1 ? 9981545732273789042ULL : 5481486777409645478ULL;
00494   }
00495 
00496   template<> Random_u32::type
00497   RandomEngine<MT19937<Random_u32>, MixerSFMT>::
00498   SelfTestResult(unsigned i)
00499     throw() {
00500     return i == 0 ? 1 :
00501       i == 1 ? 666528879UL : 2183745132UL;
00502   }
00503 
00504   template<> Random_u64::type
00505   RandomEngine<MT19937<Random_u64>, MixerSFMT>::
00506   SelfTestResult(unsigned i)
00507     throw() {
00508     return i == 0 ? 1 :
00509       i == 1 ? 12176471137395770412ULL : 66914054428611861ULL;
00510   }
00511 
00512   template<> Random_u32::type
00513   RandomEngine<SFMT19937<Random_u32>, MixerSFMT>::
00514   SelfTestResult(unsigned i)
00515     throw() {
00516     return i == 0 ? 1 :
00517       i == 1 ? 2695024307UL : 782200760UL;
00518   }
00519 
00520   template<> Random_u64::type
00521   RandomEngine<SFMT19937<Random_u64>, MixerSFMT>::
00522   SelfTestResult(unsigned i)
00523     throw() {
00524     return i == 0 ? 1 :
00525       i == 1 ? 1464461649847485149ULL : 5050640804923595109ULL;
00526   }
00527 
00528   // RandomMixer implementation
00529 
00530   template<class RandomType> void MixerMT0<RandomType>::
00531   SeedToState(const std::vector<RandomSeed::seed_type>& seed,
00532               mixer_type state[], unsigned n) throw() {
00533     // Adapted from
00534     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
00535     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
00536     const unsigned s = unsigned(seed.size());
00537     const unsigned w = mixer_t::width;
00538 
00539     mixer_type r = s ? a1 : a0;
00540     state[0] = r;
00541     for (unsigned k = 1; k < n; ++k) {
00542       r = b * (r ^ r >> w - 2) + k;
00543       r &= mask;
00544       state[k] = r;
00545     }
00546     if (s > 0) {
00547       const unsigned m = mixer_t::width / 32,
00548         s2 = (s + m - 1)/m;
00549       unsigned i1 = 1;
00550       r = state[0];
00551       for (unsigned k = (n > s2 ? n : s2), j = 0;
00552            k; --k, i1 = i1 == n - 1 ? 1 : i1 + 1, // i1 = i1 + 1 mod n - 1
00553              j = j == s2 - 1 ? 0 : j + 1 ) { // j = j+1 mod s2
00554         r = state[i1] ^ c * (r ^ r >> w - 2);
00555         r += j + mixer_type(seed[m * j]) +
00556           (m == 1 || 2 * j + 1 == s ? mixer_type(0) :
00557            mixer_type(seed[m * j + 1]) << w - 32);
00558         r &= mask;
00559         state[i1] = r;
00560       }
00561       for (unsigned k = n - 1; k; --k,
00562              i1 = i1 == n - 1 ? 1 : i1 + 1) { // i1 = i1 + 1 mod n - 1
00563         r = state[i1] ^ d * (r ^ r >> w - 2);
00564         r -= i1;
00565         r &= mask;
00566         state[i1] = r;
00567       }
00568       state[0] = typename mixer_t::type(1) << w - 1;
00569     }
00570   }
00571 
00572   template<class RandomType> void MixerMT1<RandomType>::
00573   SeedToState(const std::vector<RandomSeed::seed_type>& seed,
00574               mixer_type state[], unsigned n) throw() {
00575     // This is the algorithm given in the seed_seq class described in
00576     // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2006/n2079.pdf It is
00577     // a modification of
00578     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
00579     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
00580     const unsigned s = unsigned(seed.size());
00581     const unsigned w = mixer_t::width;
00582 
00583     mixer_type r = (a + s) & mask;
00584     state[0] = r;
00585     for (unsigned k = 1; k < n; ++k) {
00586       r = b * (r ^ r >> w - 2) + k;
00587       r &= mask;
00588       state[k] = r;
00589     }
00590     if (s > 0) {
00591       const unsigned m = mixer_t::width / 32,
00592         s2 = (s + m - 1)/m;
00593       unsigned i1 = 0;
00594       for (unsigned k = (n > s2 ? n : s2), j = 0;
00595            k; --k, i1 = i1 == n - 1 ? 0 : i1 + 1, // i1 = i1 + 1 mod n
00596              j = j == s2 - 1 ? 0 : j + 1 ) { // j = j+1 mod s2
00597         r = state[i1] ^ c * (r ^ r >> w - 2);
00598         r += j + mixer_type(seed[m * j]) +
00599           (m == 1 || 2 * j + 1 == s ? mixer_type(0) :
00600            mixer_type(seed[m * j + 1]) << w - 32);
00601         r &= mask;
00602         state[i1] = r;
00603       }
00604       for (unsigned k = n; k; --k,
00605              i1 = i1 == n - 1 ? 0 : i1 + 1) { // i1 = i1 + 1 mod n
00606         r = state[i1] ^ d * (r ^ r >> w - 2);
00607         r -= i1;
00608         r &= mask;
00609         state[i1] = r;
00610       }
00611     }
00612   }
00613 
00614   void MixerSFMT::SeedToState(const std::vector<RandomSeed::seed_type>& seed,
00615                                  mixer_type state[], unsigned n)
00616     throw() {
00617     // This is adapted from the routine init_by_array by Mutsuo Saito given in
00618     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/SFMT-src-1.2.tar.gz
00619 
00620     if (n == 0)
00621       return;                   // Nothing to do
00622 
00623     const unsigned s = unsigned(seed.size()),
00624       // Add treatment of small n with lag = (n - 1)/2 for n <= 7.  In
00625       // particular, the first operation (xor or plus) in each for loop
00626       // involves three distinct indices for n > 2.
00627       lag = n >= 623 ? 11 : (n >= 68 ? 7 : (n >= 39 ? 5 :
00628                                             (n >= 7 ? 3 : (n - 1)/2))),
00629       // count = max( s + 1, n )
00630       count = s + 1 > n ? s + 1 : n;
00631 
00632     std::fill(state, state + n, mixer_type(a));
00633     const unsigned w = mixer_t::width;
00634 
00635     unsigned i = 0, k = (n - lag) / 2, l = k + lag;
00636     mixer_type r = state[n - 1];
00637     for (unsigned j = 0; j < count; ++j,
00638            i = i == n - 1 ? 0 : i + 1,
00639            k = k == n - 1 ? 0 : k + 1,
00640            l = l == n - 1 ? 0 : l + 1) {
00641       // Here r = state[(j - 1) mod n]
00642       //      i = j mod n
00643       //      k = (j + (n - lag)/2) mod n
00644       //      l = (j + (n - lag)/2 + lag) mod n
00645       r ^= state[i] ^ state[k];
00646       r &= mask;
00647       r = b * (r ^ r >> (w - 5));
00648       state[k] += r;
00649       r += i + (j > s ? 0 : (j ? mixer_type(seed[j - 1]) : s));
00650       state[l] += r;
00651       state[i] = r;
00652     }
00653 
00654     for (unsigned j = n; j; --j,
00655            i = i == n - 1 ? 0 : i + 1,
00656            k = k == n - 1 ? 0 : k + 1,
00657            l = l == n - 1 ? 0 : l + 1) {
00658       // Here r = state[(i - 1) mod n]
00659       //      k = (i + (n - lag)/2) mod n
00660       //      l = (i + (n - lag)/2 + lag) mod n
00661       r += state[i] + state[k];
00662       r &= mask;
00663       r = c * (r ^ r >> (w - 5));
00664       r &= mask;
00665       state[k] ^= r;
00666       r -= i;
00667       r &= mask;
00668       state[l] ^= r;
00669       state[i] = r;
00670     }
00671   }
00672 
00673   // RandomAlgorithm implementation
00674 
00675   // Here, input is I, J = I + 1, K = I + M; output is I = I + N (mod N)
00676 
00677 #define MT19937_STEP(I, J, K) statev[I] = statev[K] ^           \
00678     (statev[J] & engine_type(1) ? magic : engine_type(0)) ^     \
00679     ((statev[I] & upper) | (statev[J] & lower)) >> 1
00680 
00681   // The code is cleaned up a little from Hagita's Fortran version by getting
00682   // rid of the unnecessary masking by YMASK and by using simpler logic to
00683   // restore the correct value of _state[0].
00684   //
00685   // Here input is J = I + N - 1, K = I + M - 1, and p = y[I] (only the high
00686   // bits are used); output _state[I] and p = y[I - 1].
00687 
00688 #define MT19937_REVSTEP(I, J, K) {                                      \
00689     engine_type q = statev[J] ^ statev[K], s = q >> (width - 1);        \
00690     q = (q ^ (s ? magic : engine_type(0))) << 1 | s;                    \
00691     statev[I] = (p & upper) | (q & lower);                              \
00692     p = q;                                                              \
00693   }
00694 
00695   template<class RandomType>
00696   void MT19937<RandomType>::Transition(long long count, internal_type statev[])
00697     throw() {
00698     if (count > 0)
00699       for (; count; --count) {
00700         // This ONLY uses high bit of statev[0]
00701         unsigned i = 0;
00702         for (; i < N - M; ++i) MT19937_STEP(i, i + 1, i + M    );
00703         for (; i < N - 1; ++i) MT19937_STEP(i, i + 1, i + M - N);
00704         MT19937_STEP(N - 1, 0, M - 1); // i = N - 1
00705       }
00706     else if (count < 0)
00707       for (; count; ++count) {
00708         // This ONLY uses high bit of statev[0]
00709         engine_type p = statev[0];
00710         // Fix low bits of statev[0] and compute y[-1]
00711         MT19937_REVSTEP(0, N - 1, M - 1); // i = N
00712         unsigned i = N - 1;
00713         for (; i > N - M; --i) MT19937_REVSTEP(i, i - 1, i + M - 1 - N);
00714         for (; i        ; --i) MT19937_REVSTEP(i, i - 1, i + M - 1    );
00715         MT19937_REVSTEP(0, N - 1, M - 1); // i = 0
00716       }
00717   }
00718 
00719 #undef MT19937_STEP
00720 #undef MT19937_REVSTEP
00721 
00722   template<class RandomType>
00723   void MT19937<RandomType>::NormalizeState(engine_type state[]) throw () {
00724 
00725     // Perform the MT-specific sanity check on the resulting state ensuring
00726     // that the significant 19937 bits are not all zero.
00727     state[0] &= upper;  // Mask out unused bits
00728     unsigned i = 0;
00729     while (i < N && state[i] == 0)
00730       ++i;
00731     if (i >= N)
00732       state[0] = engine_type(1) << (width - 1); // with prob 2^-19937
00733 
00734     // This sets the low R bits of _state[0] consistent with the rest of the
00735     // state.  Needed to handle SetCount(-N); Ran32(); immediately following
00736     // reseeding.  This wasn't required in the original code because a
00737     // Transition was always done first.
00738     engine_type q = state[N - 1] ^ state[M - 1], s = q >> (width - 1);
00739     q = (q ^ (s ? magic : engine_type(0))) << 1 | s;
00740     state[0] = (state[0] & upper) | (q & lower);
00741   }
00742 
00743   template<class RandomType>
00744   void MT19937<RandomType>::CheckState(const engine_type state[],
00745                                        Random_u32::type& check)
00746     throw(std::out_of_range) {
00747     engine_type x = 0;
00748     Random_u32::type c = check;
00749     for (unsigned i = 0; i < N; ++i) {
00750       engine_t::CheckSum(state[i], c);
00751       x |= state[i];
00752     }
00753     if (x == 0)
00754       throw std::out_of_range("MT19937: All-zero state");
00755 
00756     // There are only width*(N-1) + 1 = 19937 independent bits of state.  Thus
00757     // the low width-1 bits of _state[0] are derivable from the other bits in
00758     // state.  Verify that the redundant bits bits are consistent.
00759     engine_type q = state[N - 1] ^ state[M - 1], s = q >> (width - 1);
00760     q = (q ^ (s ? magic : engine_type(0))) << 1 | s;
00761     if ((q ^ state[0]) & lower)
00762       throw std::out_of_range("MT19937: Invalid state");
00763 
00764     check = c;
00765   }
00766 
00767 #if HAVE_SSE2
00768 
00769   // Transition is from Saito's Master's Thesis
00770   // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/M062821.pdf
00771   //
00772   // This implements
00773   //
00774   //     w_{i+N} = w_i A xor w_M B xor w_{i+N-2} C xor w_{i+N-1} D
00775   //
00776   // where w_i is a 128-bit word and
00777   //
00778   //     w A = (w << 8)_128 xor w
00779   //     w B = (w >> 11)_32 & MSK
00780   //     w C = (w >> 8)_128
00781   //     w D = (w << 18)_32
00782   //
00783   // Here the _128 means shift is of whole 128-bit word.  _32 means the shifts
00784   // are independently done on each 32-bit word.
00785   //
00786   // In SFMT19937_STEP32 and SFMT19937_STEP64 input is I, J = I + M and output
00787   // is I = I + N (mod N).  On input, s and r give state for I + N - 2 and I +
00788   // N - 1; on output s and r give state for I + N - 1 and I + N.  The
00789   // implementation of 128-bit operations is open-coded in a portable fashion
00790   // (with LSB ordering).
00791   //
00792   // N.B. Here N and M are the lags in units of BitWidth words and so are 4
00793   // (for u32 implementation) or 2 (for u64 implementation) times bigger than
00794   // defined in Saito's thesis.
00795 
00796   // This is adapted from SFMT-sse.c in the SFMT 1.2 distribution.
00797   // The order of instructions has been rearranged to increase the
00798   // speed slightly
00799 
00800 #define SFMT19937_STEP128(I, J) {                       \
00801     internal_type x = _mm_load_si128(statev + I),       \
00802       y = _mm_srli_epi32(statev[J], 11),                \
00803       z = _mm_srli_si128(s, 1);                         \
00804     s = _mm_slli_epi32(r, 18);                          \
00805     z = _mm_xor_si128(z, x);                            \
00806     x = _mm_slli_si128(x, 1);                           \
00807     z = _mm_xor_si128(z, s);                            \
00808     y = _mm_and_si128(y, m);                            \
00809     z = _mm_xor_si128(z, x);                            \
00810     s = r;                                              \
00811     r = _mm_xor_si128(z, y);                            \
00812     _mm_store_si128(statev + I, r);                     \
00813   }
00814 
00815   // This undoes SFMT19937_STEP.  Trivially, we have
00816   //
00817   //     w_i A = w_{i+N} xor w_{i+M} B xor w_{i+N-2} C xor w_{i+N-1} D
00818   //
00819   // Given w_i A we can determine w_i from the observation that A^16 =
00820   // identity, thus
00821   //
00822   //     w_i = (w_i A) A^15
00823   //
00824   // Because x A^(2^n) = x << (8*2^n) xor x, the operation y = x A^15 can be
00825   // implemented as
00826   //
00827   //     y'   = (x    << 64)_128 xor x    = x    A^8
00828   //     y''  = (y'   << 32)_128 xor y'   = y'   A^4 = x A^12
00829   //     y''' = (y''  << 16)_128 xor y''  = y''  A^2 = x A^14
00830   //     y    = (y''' <<  8)_128 xor y''' = y''' A   = x A^15
00831   //
00832   // Here input is I = I + N, J = I + M, K = I + N - 2, L = I + N -1, and
00833   // output is I = I.
00834   //
00835   // This is about 15-35% times slower than SFMT19937_STEPNN, because (1) there
00836   // doesn't appear to be a straightforward way of saving intermediate results
00837   // across calls as with SFMT19937_STEPNN and (2) w A^15 is slower to compute
00838   // than w A.
00839 
00840 #define SFMT19937_REVSTEP128(I, J, K, L) {              \
00841     internal_type x = _mm_load_si128(statev + I),       \
00842       y = _mm_srli_epi32(statev[J], 11),                \
00843       z = _mm_slli_epi32(statev[L], 18);                \
00844     y = _mm_and_si128(y, m);                            \
00845     x = _mm_xor_si128(x, _mm_srli_si128(statev[K], 1)); \
00846     x = _mm_xor_si128(x, z);                            \
00847     x = _mm_xor_si128(x, y);                            \
00848     x = _mm_xor_si128(_mm_slli_si128(x, 8), x);         \
00849     x = _mm_xor_si128(_mm_slli_si128(x, 4), x);         \
00850     x = _mm_xor_si128(_mm_slli_si128(x, 2), x);         \
00851     x = _mm_xor_si128(_mm_slli_si128(x, 1), x);         \
00852     _mm_store_si128(statev + I, x);                     \
00853   }
00854 
00855   template<class RandomType>
00856   void SFMT19937<RandomType>::Transition(long long count, internal_type statev[])
00857     throw() {
00858     const internal_type m = _mm_set_epi32(magic3, magic2, magic1, magic0);
00859     if (count > 0) {
00860       internal_type s = _mm_load_si128(statev + N128 - 2),
00861         r = _mm_load_si128(statev + N128 - 1);
00862       for (; count; --count) {
00863         unsigned i = 0;
00864         for (; i + M128 < N128; ++i) SFMT19937_STEP128(i, i + M128       );
00865         for (; i < N128       ; ++i) SFMT19937_STEP128(i, i + M128 - N128);
00866       }
00867     } else if (count < 0)
00868       for (; count; ++count) {
00869         unsigned i = N128;
00870         for (; i + M128 > N128;) {
00871           --i; SFMT19937_REVSTEP128(i, i + M128 - N128, i - 2, i - 1);
00872         }
00873         for (; i > 2;) {
00874           --i; SFMT19937_REVSTEP128(i, i + M128, i - 2, i - 1);
00875         }
00876         SFMT19937_REVSTEP128(1, M128 + 1, N128 - 1, 0       ); // i = 1
00877         SFMT19937_REVSTEP128(0, M128    , N128 - 2, N128 - 1); // i = 0
00878       }
00879   }
00880 
00881 #undef SFMT19937_STEP128
00882 #undef SFMT19937_REVSTEP128
00883 
00884 #elif HAVE_ALTIVEC
00885 
00886   // The Altivec versions of SFMT19937_{,REV}STEP128 are simply translated from
00887   // the SSE2 versions.  The only significant differences arise because of the
00888   // MSB ordering of the PowerPC.  This means that teh 32-bit and 64-bit
00889   // versions are no different because 32-bit and 64-bit words don't pack
00890   // together in the same way as on an SSE2 machine (see the two definitions of
00891   // magic).  This also means that the 128-bit byte shifts on an LSB machine
00892   // change into more complicated byte permutations.
00893 
00894 #define ALTIVEC_PERM(X, P) vec_perm(X, P, P)
00895 
00896 #define SFMT19937_STEP128(I, J) {                               \
00897     internal_type x = statev[I],                                \
00898       z = vec_xor(vec_xor(ALTIVEC_PERM(s, right1), x),          \
00899                   vec_sl(r, bitleft));                          \
00900     s = r;                                                      \
00901     r = vec_xor(z,                                              \
00902                 vec_xor(ALTIVEC_PERM(x, left1),                 \
00903                         vec_and(vec_sr(statev[J], bitright),    \
00904                                 magic)));                       \
00905     statev[I] = r;                                              \
00906   }
00907 
00908 #define SFMT19937_REVSTEP128(I, J, K, L) {              \
00909     internal_type x = statev[I],                        \
00910       y = vec_sr(statev[J], bitright),                  \
00911       z = vec_sl(statev[L], bitleft);                   \
00912     y = vec_and(y, magic);                              \
00913     x = vec_xor(x, ALTIVEC_PERM(statev[K], right1));    \
00914     x = vec_xor(x, z);                                  \
00915     x = vec_xor(x, y);                                  \
00916     x = vec_xor(ALTIVEC_PERM(x, left8), x);             \
00917     x = vec_xor(ALTIVEC_PERM(x, left4), x);             \
00918     x = vec_xor(ALTIVEC_PERM(x, left2), x);             \
00919     statev[I] = vec_xor(ALTIVEC_PERM(x, left1), x);     \
00920   }
00921 
00922   template<class RandomType>
00923   void SFMT19937<RandomType>::Transition(long long count, internal_type statev[])
00924     throw() {
00925     const internal_type magic = width == 32 ?
00926       (vector unsigned)(magic0, magic1, magic2, magic3) :
00927       (vector unsigned)(magic1, magic0, magic3, magic2),
00928       bitleft = (vector unsigned)(18, 18, 18, 18),
00929       bitright = (vector unsigned)(11, 11, 11, 11);
00930     // Shift left and right by 1 byte.  Note that vec_perm(X, Y, P) glues X and
00931     // Y together into a 32-byte quantity and then the 16-byte permutation
00932     // vector P specifies which bytes to put into the 16-byte output.  We
00933     // follow here the convention of using Y = P and using the zero entries in
00934     // P to allow zero bytes to be introduces into the shifted output.  The
00935     // following describes how the left1 table (32-bit version) is produced:
00936     //
00937     // Byte layout of original with LSB ordering
00938     // 33 32 31 30  23 22 21 20  13 12 11 10  03 02 01 00
00939     // shift left by 1 byte (z means zeros enter)
00940     // 32 31 30 23  22 21 20 13  12 11 10 03  02 01 00 zz
00941     //
00942     // Rearrange original to LSB order in 4-byte units
00943     // 03 02 01 00  13 12 11 10  23 22 21 20  33 32 31 30
00944     // with sequential MSB byte indices
00945     // 0  1  2  3   4  5  6  7   8  9 10 11  12 13 14 15
00946     //
00947     // Rearrange shift left verion to LSB order in 4-byte units
00948     // 02 01 00 zz  12 11 10 03  22 21 20 13  32 31 30 23
00949     // with corresponding MSB byte indices
00950     // 1  2  3  z   5  6  7  0   9 10 11  4  13 14 15  8
00951     //
00952     // Replace byte index at x by 16 + index of 0 = 16 + 7 = 23 to give
00953     // 1  2  3 23   5  6  7  0   9 10 11  4  13 14 15  8
00954     const vector unsigned char left1 = width == 32 ?
00955       (vector unsigned char)(1,2,3,23, 5,6,7,0, 9,10,11,4, 13,14,15,8) :
00956       (vector unsigned char)(1,2,3,4,5,6,7,31, 9,10,11,12,13,14,15,0),
00957       right1 = width == 32 ?
00958       (vector unsigned char)(7,0,1,2, 11,4,5,6, 15,8,9,10, 17,12,13,14) :
00959       (vector unsigned char)(15,0,1,2,3,4,5,6, 17,8,9,10,11,12,13,14);
00960     if (count > 0) {
00961       internal_type s = statev[N128 - 2],
00962         r = statev[N128 - 1];
00963       for (; count; --count) {
00964         unsigned i = 0;
00965         for (; i + M128 < N128; ++i) SFMT19937_STEP128(i, i + M128       );
00966         for (; i < N128       ; ++i) SFMT19937_STEP128(i, i + M128 - N128);
00967       }
00968     } else if (count < 0) {
00969       // leftN shifts left by N bytes.
00970       const vector unsigned char left2 = width == 32 ?
00971         (vector unsigned char)(2,3,22,22, 6,7,0,1, 10,11,4,5, 14,15,8,9) :
00972         (vector unsigned char)(2,3,4,5,6,7,30,30, 10,11,12,13,14,15,0,1),
00973         left4 = width == 32 ?
00974         (vector unsigned char)(20,20,20,20, 0,1,2,3, 4,5,6,7, 8,9,10,11) :
00975         (vector unsigned char)(4,5,6,7,28,28,28,28, 12,13,14,15,0,1,2,3),
00976         left8 = (vector unsigned char)(24,24,24,24,24,24,24,24,0,1,2,3,4,5,6,7);
00977       for (; count; ++count) {
00978         unsigned i = N128;
00979         for (; i + M128 > N128;) {
00980           --i; SFMT19937_REVSTEP128(i, i + M128 - N128, i - 2, i - 1);
00981         }
00982         for (; i > 2;) {
00983           --i; SFMT19937_REVSTEP128(i, i + M128, i - 2, i - 1);
00984         }
00985         SFMT19937_REVSTEP128(1, M128 + 1, N128 - 1, 0       ); // i = 1
00986         SFMT19937_REVSTEP128(0, M128    , N128 - 2, N128 - 1); // i = 0
00987       }
00988     }
00989   }
00990 
00991 #undef SFMT19937_STEP128
00992 #undef SFMT19937_REVSTEP128
00993 #undef ALTVEC_PERM
00994 
00995 #else  // neither HAVE_SSE2 or HAVE_ALTIVEC
00996 
00997 #define SFMT19937_STEP32(I, J) {                                \
00998     internal_type t = statev[I] ^ statev[I] << 8 ^              \
00999       statev[J] >> 11 & magic0 ^                                \
01000       (s0 >> 8 | s1 << 24) ^ r0 << 18;                          \
01001     s0 = r0; r0 = t & mask;                                     \
01002     t = statev[I + 1] ^                                         \
01003       (statev[I + 1] << 8 | statev[I] >> 24) ^                  \
01004       statev[J + 1] >> 11 & magic1 ^                            \
01005       (s1 >> 8 | s2 << 24) ^ r1 << 18;                          \
01006     s1 = r1; r1 = t & mask;                                     \
01007     t = statev[I + 2] ^                                         \
01008       (statev[I + 2] << 8 | statev[I + 1] >> 24) ^              \
01009       statev[J + 2] >> 11 & magic2 ^                            \
01010       (s2 >> 8 | s3 << 24) ^ r2 << 18;                          \
01011     s2 = r2; r2 = t & mask;                                     \
01012     t = statev[I + 3] ^                                         \
01013       (statev[I + 3] << 8 | statev[I + 2] >> 24) ^              \
01014       statev[J + 3] >> 11 & magic3 ^ s3 >> 8 ^ r3 << 18;        \
01015     s3 = r3; r3 = t & mask;                                     \
01016     statev[I    ] = r0; statev[I + 1] = r1;                     \
01017     statev[I + 2] = r2; statev[I + 3] = r3;                     \
01018   }
01019 
01020 #define SFMT19937_REVSTEP32(I, J, K, L) {                       \
01021     internal_type                                               \
01022       t0 = (statev[I] ^ statev[J] >> 11 & magic0 ^              \
01023             (statev[K] >> 8 | statev[K + 1] << 24) ^            \
01024             statev[L] << 18) & mask,                            \
01025       t1 = (statev[I + 1] ^                                     \
01026             statev[J + 1] >> 11 & magic1 ^                      \
01027             (statev[K + 1] >> 8 | statev[K + 2] << 24) ^        \
01028             statev[L + 1] << 18) & mask,                        \
01029       t2 = (statev[I + 2] ^                                     \
01030             statev[J + 2] >> 11 & magic2 ^                      \
01031             (statev[K + 2] >> 8 | statev[K + 3] << 24) ^        \
01032             statev[L + 2] << 18) & mask,                        \
01033       t3 = (statev[I + 3] ^                                     \
01034             statev[J + 3] >> 11 & magic3 ^                      \
01035             statev[K + 3] >> 8 ^                                \
01036             statev[L + 3] << 18) & mask;                        \
01037     t3 ^= t1; t2 ^= t0; t3 ^= t2; t2 ^= t1; t1 ^= t0;           \
01038     t3 ^= t2 >> 16 | t3 << 16 & mask;                           \
01039     t2 ^= t1 >> 16 | t2 << 16 & mask;                           \
01040     t1 ^= t0 >> 16 | t1 << 16 & mask;                           \
01041     t0 ^=            t0 << 16 & mask;                           \
01042     statev[I    ] = t0 ^             t0 << 8 & mask;            \
01043     statev[I + 1] = t1 ^ (t0 >> 24 | t1 << 8 & mask);           \
01044     statev[I + 2] = t2 ^ (t1 >> 24 | t2 << 8 & mask);           \
01045     statev[I + 3] = t3 ^ (t2 >> 24 | t3 << 8 & mask);           \
01046  }
01047 
01048   template<>
01049   void SFMT19937<Random_u32>::Transition(