blitz Version 0.10
blitz/rand-mt.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // $Id: rand-mt.h,v 1.3 2011/03/25 22:41:16 julianc Exp $
00003 
00004 /* A C-program for MT19937: Integer version (1998/4/6)            */
00005 /*  genrand() generates one pseudorandom unsigned integer (32bit) */
00006 /* which is uniformly distributed among 0 to 2^32-1  for each     */
00007 /* call. sgenrand(seed) set initial values to the working area    */
00008 /* of 624 words. Before genrand(), sgenrand(seed) must be         */
00009 /* called once. (seed is any 32-bit integer except for 0).        */
00010 /*   Coded by Takuji Nishimura, considering the suggestions by    */
00011 /* Topher Cooper and Marc Rieffel in July-Aug. 1997.              */
00012 
00013 /* This library is free software; you can redistribute it and/or   */
00014 /* modify it under the terms of the GNU Library General Public     */
00015 /* License as published by the Free Software Foundation; either    */
00016 /* version 2 of the License, or (at your option) any later         */
00017 /* version.                                                        */
00018 /* This library is distributed in the hope that it will be useful, */
00019 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
00020 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
00021 /* See the GNU Library General Public License for more details.    */
00022 /* You should have received a copy of the GNU Library General      */
00023 /* Public License along with this library; if not, write to the    */
00024 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */ 
00025 /* 02111-1307  USA                                                 */
00026 
00027 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
00028 /* When you use this, send an email to: matumoto@math.keio.ac.jp   */
00029 /* with an appropriate reference to your work.                     */
00030 
00031 /* REFERENCE                                                       */
00032 /* M. Matsumoto and T. Nishimura,                                  */
00033 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform  */
00034 /* Pseudo-Random Number Generator",                                */
00035 /* ACM Transactions on Modeling and Computer Simulation,           */
00036 /* Vol. 8, No. 1, January 1998, pp 3--30.                          */
00037 
00038 // See http://www.math.keio.ac.jp/~matumoto/emt.html
00039 
00040 // 1999-01-25 adapted to STL-like idiom
00041 // allan@stokes.ca (Allan Stokes) www.stokes.ca
00042 
00043 #ifndef BZ_RAND_MT
00044 #define BZ_RAND_MT
00045 
00046 #ifndef BZ_BLITZ_H
00047  #include <blitz/blitz.h>
00048 #endif
00049 
00050 #include <vector>
00051 
00052 BZ_NAMESPACE(blitz)
00053 
00054 // decomposition issues:
00055 //   machine representation of integer types
00056 //   output buffer option verses inline post-conditioning
00057 
00058 class MersenneTwister
00059 {
00060 private:
00061   typedef unsigned int twist_int; // must be at least 32 bits
00062                                   // larger might be faster
00063   typedef vector<twist_int> State;
00064   typedef State::iterator Iter;
00065 
00066   struct BitMixer {
00067     enum { K = 0x9908b0df };
00068     BitMixer() : s0(0) {}
00069     inline friend twist_int low_mask (twist_int s1) {
00070       return (s1&1u) ? K : 0u;
00071     }
00072     inline twist_int high_mask (twist_int s1) const {
00073       return ((s0&0x80000000)|(s1&0x7fffffff))>>1;
00074     }
00075     inline twist_int operator() (twist_int s1) {
00076       twist_int r = high_mask(s1) ^ low_mask(s1);
00077       s0 = s1;
00078       return r;
00079     }
00080     twist_int s0;
00081   };
00082 
00083 enum { N = 624, PF = 397, reference_seed = 4357 }; 
00084   
00085 public: 
00086   MersenneTwister () {} // S empty will trigger auto-seed
00087 
00088   void seed (twist_int seed = reference_seed)
00089   {
00090     if (!S.size()) S.resize(N);
00091     enum { Knuth_A = 69069 }; 
00092     twist_int x = seed & 0xFFFFFFFF;
00093     Iter s = S.begin();
00094     twist_int mask = (seed == reference_seed) ? 0 : 0xFFFFFFFF;
00095     for (int j = 0; j < N; ++j) {
00096       // adding j here avoids the risk of all zeros 
00097       // we suppress this term in "compatibility" mode  
00098       *s++ = (x + (mask & j)) & 0xFFFFFFFF; 
00099       x *= Knuth_A;
00100     }
00101     reload();
00102   }
00103 
00104   void reload (void)
00105   {
00106     if (!S.size()) seed (); // auto-seed detection
00107 
00108     Iter p0 = S.begin();
00109     Iter pM = p0 + PF;
00110     BitMixer twist;
00111     twist (S[0]); // prime the pump
00112     for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
00113       *p0 = *pM ^ twist (p0[1]);
00114     pM = S.begin();
00115     for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
00116       *p0 = *pM ^ twist (p0[1]);
00117     *p0 = *pM ^ twist (S[0]);
00118 
00119     I = S.begin();
00120   }
00121 
00122   inline twist_int random (void)
00123   {
00124     if (I >= S.end()) reload();
00125     twist_int y = *I++;
00126     y ^= (y >> 11);
00127     y ^= (y <<  7) & 0x9D2C5680;
00128     y ^= (y << 15) & 0xEFC60000;
00129     y ^= (y >> 18);
00130     return y;
00131   }
00132 
00133 private:
00134   State   S;
00135   Iter    I;
00136 };
00137 
00138 
00139 // This version returns a double in the range [0,1).
00140 
00141 class MersenneTwisterDouble {
00142 
00143 public:
00144   MersenneTwisterDouble()
00145   {
00146       // f = 1/(2^32);
00147       f = (1.0 / 65536) / 65536;
00148   }
00149 
00150   void randomize(unsigned int seed)
00151   {
00152       gen_.seed(seed);
00153   }
00154 
00155   double random()
00156   {
00157       unsigned long y1 = gen_.random();
00158       unsigned long y2 = gen_.random();
00159 
00160       return ((y1 * f) * y2 * f);
00161   }
00162 
00163 private:
00164   MersenneTwister gen_;
00165   double f;
00166 };
00167 
00168 BZ_NAMESPACE_END
00169 
00170 #endif // BZ_RAND_MT
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines