blitz/rand-mt.h

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

Generated on Mon Dec 3 09:21:54 2007 for blitz by  doxygen 1.5.1