blitz Version 0.10
|
00001 // -*- C++ -*- 00002 /*************************************************************************** 00003 * blitz/rand-tt800.h Matsumoto and Kurita's TT800 uniform random 00004 * number generator. 00005 * 00006 * $Id: rand-tt800.h,v 1.3 2011/03/25 22:41:16 julianc Exp $ 00007 * 00008 *************************************************************************** 00009 * 00010 * The class TT800 encapsulates Makoto Matsumoto and Yoshiharu Kurita's 00011 * TT800 twisted generalized feedback shift register (TGFSR) random number 00012 * generator. The generator has period 2^800 - 1. 00013 * 00014 * Contact: M. Matsumoto <matumoto@math.keio.ac.jp> 00015 * 00016 * See: M. Matsumoto and Y. Kurita, Twisted GFSR Generators II, 00017 * ACM Transactions on Modelling and Computer Simulation, 00018 * Vol. 4, No. 3, 1994, pages 254-266. 00019 * 00020 * (c) 1994 Association for Computing Machinery. 00021 * 00022 * Distributed with consent of the authors. 00023 * 00024 ***************************************************************************/ 00025 00026 #ifndef BZ_RAND_TT800_H 00027 #define BZ_RAND_TT800_H 00028 00029 #ifndef BZ_RANDOM_H 00030 #include <blitz/random.h> 00031 #endif 00032 00033 BZ_NAMESPACE(blitz) 00034 00035 class TT800 { 00036 00037 public: 00038 typedef double T_numtype; 00039 00040 TT800(double low = 0.0, double high = 1.0, double = 0.0) 00041 : low_(low), length_(high-low) 00042 { 00043 // Initial 25 seeds 00044 x[0] = 0x95f24dab; x[1] = 0x0b685215; x[2] = 0xe76ccae7; 00045 x[3] = 0xaf3ec239; x[4] = 0x715fad23; x[5] = 0x24a590ad; 00046 x[6] = 0x69e4b5ef; x[7] = 0xbf456141; x[8] = 0x96bc1b7b; 00047 x[9] = 0xa7bdf825; x[10] = 0xc1de75b7; x[11] = 0x8858a9c9; 00048 x[12] = 0x2da87693; x[13] = 0xb657f9dd; x[14] = 0xffdc8a9f; 00049 x[15] = 0x8121da71; x[16] = 0x8b823ecb; x[17] = 0x885d05f5; 00050 x[18] = 0x4e20cd47; x[19] = 0x5a9ad5d9; x[20] = 0x512c0c03; 00051 x[21] = 0xea857ccd; x[22] = 0x4cc1d30f; x[23] = 0x8891a8a1; 00052 x[24] = 0xa6b7aadb; 00053 00054 // Magic vector 'a', don't change 00055 mag01[0] = 0; 00056 mag01[1] = 0x8ebfd028; 00057 00058 k = 0; 00059 00060 // f = 1/(2^32); 00061 00062 f = (1.0 / 65536) / 65536; 00063 } 00064 00065 void randomize() 00066 { 00067 BZ_NOT_IMPLEMENTED(); // NEEDS_WORK 00068 } 00069 00070 unsigned long randomUint32() 00071 { 00072 if (k==N) 00073 generate(); 00074 00075 unsigned long y = x[k]; 00076 y ^= (y << 7) & 0x2b5b2500; /* s and b, magic vectors */ 00077 y ^= (y << 15) & 0xdb8b0000; /* t and c, magic vectors */ 00078 y &= 0xffffffff; /* you may delete this line if word size = 32 */ 00079 00080 // the following line was added by Makoto Matsumoto in the 1996 version 00081 // to improve lower bit's corellation. 00082 // Delete this line to use the code published in 1994. 00083 00084 y ^= (y >> 16); /* added to the 1994 version */ 00085 k++; 00086 } 00087 00088 double random() 00089 { 00090 unsigned long y1 = randomUint32(); 00091 unsigned long y2 = randomUint32(); 00092 00093 return low_ + length_ * ((y1 * f) * y2 * f); 00094 } 00095 00096 protected: 00097 void generate() 00098 { 00099 /* generate N words at one time */ 00100 int kk; 00101 for (kk=0;kk<N-M;kk++) { 00102 x[kk] = x[kk+M] ^ (x[kk] >> 1) ^ mag01[x[kk] % 2]; 00103 } 00104 for (; kk<N;kk++) { 00105 x[kk] = x[kk+(M-N)] ^ (x[kk] >> 1) ^ mag01[x[kk] % 2]; 00106 } 00107 k=0; 00108 } 00109 00110 private: 00111 enum { N = 25, M = 7 }; 00112 00113 double low_, length_; 00114 double f; 00115 int k; 00116 unsigned long x[N]; 00117 unsigned long mag01[2]; 00118 }; 00119 00120 BZ_NAMESPACE_END 00121 00122 #endif // BZ_RAND_TT800_H 00123