00001 00002 // This java version of the Mersenne Twister ``MT19937'' is translated 00003 // from the C version by Takuji Nishimura and Makoto Matsumoto 00004 // distributed in the file // mt19937ar-cok.c. They describe 00005 // this as ``a faster version by taking Shawn Cokus's optimization, 00006 // Matthe Bellew's simplification, Isaku Wada's real version.'' 00007 // Hopefully my translation to java has not introduced any new errors. 00008 // One new feature I have added is the capacity to save or restore the 00009 // state of the random number generator, which allows us, for example, 00010 // to rerun the simulation at new parameter values but the same 00011 // random draws, without having to save the entire vector of random draws 00012 // used in the simulation. I also added a method to simulate a Gaussian 00013 // draw. 00014 00015 // Copyright (c) Philip H. Dybvig 2002 00016 // Philip H. Dybvig <http://phildybvig.com> 00017 00018 // One warning: this implementation is not serialized (so it may become 00019 // confused if different threads call the same instantiation of the 00020 // class at the same time). In principle, this should be faster than 00021 // a serialized version but some care is required in the use, for 00022 // example, by using different instantiations in different threads that 00023 // may run at the same time. I would also be wary of code like 00024 00025 // MTwister mt = new MTwister(274983L); 00026 // double x; 00027 // x = mt.genrand_real1() + mt.genrand_real1(); 00028 00029 // because some java engine may try to execute the two calls simultaneously 00030 // (good for MP machines perhaps) in different threads. Just to be safe, 00031 // I would do something like 00032 00033 // MTwister mt = new MTwister(274983L); 00034 // double x; 00035 // x = mt.genrand_real1(); 00036 // x += mt.genrand_real1(); 00037 00038 // instead. 00039 00040 // Notes: (1) I applied the C pre-processor to implement the #define 00041 // lines in the original C program. The resulting code looks ugly 00042 // but should be fast because it moves everything in-line (as in the 00043 // C program). 00044 // (2) I removed technical comments about the algorithm from the original 00045 // program. For reference, I have posted the version of the Matsumoto- 00046 // Nishimura C program I am working from in the directory 00047 // <http://phildybvig.com/misc/MTwister>. 00048 // (3) While I did spot checks of the code against the distributed output 00049 // from the C program, I do not guarantee conformance to that code and 00050 // in general I offer the same ``caveat emptor'' license as originally 00051 // given by Matsumoto and Nishimura below. 00052 // (4) I have made all the unsigned variables into long's. Making them 00053 // int's should be possible since int's in java are 32 bit numbers. 00054 // However, they would have to be converted to longs (using the sign 00055 // bit as high order bit) when using them. I am not sure there would 00056 // be any net savings, and doing it this way requires less debugging. 00057 00058 // Here is the original copyright and contact information for Makoto 00059 // Matsumoto and Takuji Nishimura. 00060 00061 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 00062 // All rights reserved. 00063 00064 // Redistribution and use in source and binary forms, with or without 00065 // modification, are permitted provided that the following conditions 00066 // are met: 00067 00068 // 1. Redistributions of source code must retain the above copyright 00069 // notice, this list of conditions and the following disclaimer. 00070 00071 // 2. Redistributions in binary form must reproduce the above copyright 00072 // notice, this list of conditions and the following disclaimer in the 00073 // documentation and/or other materials provided with the distribution. 00074 00075 // 3. The names of its contributors may not be used to endorse or promote 00076 // products derived from this software without specific prior written 00077 // permission. 00078 00079 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00080 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00081 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00082 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 00083 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00084 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00085 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00086 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00087 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00088 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00089 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00090 00091 00092 // Any feedback is very welcome. 00093 // http://www.math.keio.ac.jp/matumoto/emt.html 00094 // email: matumoto@math.keio.ac.jp 00095 package org.sci2s.eamhco; 00096 00097 public class MTwister { 00098 long[] state=new long[624]; 00099 int left = 1; 00100 int initf = 0; 00101 int inext; 00102 00103 public MTwister() {} 00104 00105 public MTwister(long s) { 00106 init_genrand(s); 00107 } 00108 00109 public MTwister(long[] init_key) { 00110 init_by_array(init_key); 00111 } 00112 00113 public void init_genrand(long s) { 00114 int j; 00115 state[0]= s & 0xffffffffL; 00116 for (j=1; j<624; j++) { 00117 state[j] = (1812433253L * (state[j-1] ^ (state[j-1] >> 30)) + j); 00118 state[j] &= 0xffffffffL; 00119 } 00120 left = 1; 00121 initf = 1; 00122 } 00123 00124 void init_by_array(long[] init_key) { 00125 int i, j, k; 00126 int key_length; 00127 key_length = init_key.length; 00128 init_genrand(19650218L); 00129 i=1; j=0; 00130 k = (624>key_length ? 624 : key_length); 00131 for (; k>0; k--) { 00132 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525L))+ init_key[j] + j; 00133 state[i] &= 0xffffffffL; 00134 i++; j++; 00135 if (i>=624) { 00136 state[0] = state[624 -1]; 00137 i=1; 00138 } 00139 if (j>=key_length) 00140 j=0; 00141 } 00142 for (k=624 -1; k>0; k--) { 00143 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941L)) - i; 00144 state[i] &= 0xffffffffL; 00145 i++; 00146 if (i>=624) { 00147 state[0] = state[624 -1]; i=1; 00148 } 00149 } 00150 state[0] = 0x80000000L; 00151 left = 1; 00152 initf = 1; 00153 } 00154 00155 void next_state() { 00156 int ip=0; 00157 int j; 00158 if (initf==0) init_genrand(5489L); 00159 left = 624; 00160 inext = 0; 00161 for (j=624 -397 +1; (--j)>0; ip++) 00162 state[ip] = state[ip+397] ^ ((( ((state[ip+0]) & 0x80000000L) | ((state[ip+1]) & 0x7fffffffL) ) >> 1) ^ (((state[ip+1]) & 1L) != 0L ? 0x9908b0dfL : 0L)); 00163 for (j=397; (--j)>0; ip++) 00164 state[ip] = state[ip+397 -624] ^ ((( ((state[ip+0]) & 0x80000000L) | ((state[ip+1]) & 0x7fffffffL) ) >> 1) ^ (((state[ip+1]) & 1L) != 0L ? 0x9908b0dfL : 0L)); 00165 state[ip] = state[ip+397 -624] ^ ((( ((state[ip+0]) & 0x80000000L) | ((state[0]) & 0x7fffffffL) ) >> 1) ^ (((state[0]) & 1L) != 0L ? 0x9908b0dfL : 0L)); 00166 } 00167 00168 // generates a random number on [0,0xffffffff]-interval 00169 long genrand_int32() { 00170 long y; 00171 if (--left == 0) next_state(); 00172 y = state[inext++]; 00173 y ^= (y >> 11); 00174 y ^= (y << 7) & 0x9d2c5680L; 00175 y ^= (y << 15) & 0xefc60000L; 00176 y ^= (y >> 18); 00177 return y; 00178 } 00179 00180 // generates a random number on [0,0x7fffffff]-interval 00181 long genrand_int31() { 00182 long y; 00183 if (--left == 0) next_state(); 00184 y = state[inext++]; 00185 y ^= (y >> 11); 00186 y ^= (y << 7) & 0x9d2c5680L; 00187 y ^= (y << 15) & 0xefc60000L; 00188 y ^= (y >> 18); 00189 return (long)(y>>1); 00190 } 00191 00192 // generates a random number on [0,1]-real-interval 00193 public double genrand_real1() { 00194 long y; 00195 if (--left == 0) next_state(); 00196 y = state[inext++]; 00197 y ^= (y >> 11); 00198 y ^= (y << 7) & 0x9d2c5680L; 00199 y ^= (y << 15) & 0xefc60000L; 00200 y ^= (y >> 18); 00201 return (double)y * (1.0/4294967295.0); 00202 } 00203 00204 // generates a random number on [0,1)-real-interval 00205 public double genrand_real2() { 00206 long y; 00207 if (--left == 0) next_state(); 00208 y = state[inext++]; 00209 y ^= (y >> 11); 00210 y ^= (y << 7) & 0x9d2c5680L; 00211 y ^= (y << 15) & 0xefc60000L; 00212 y ^= (y >> 18); 00213 return (double)y * (1.0/4294967296.0); 00214 } 00215 00216 // generates a random number on (0,1)-real-interval 00217 public double genrand_real3() { 00218 long y; 00219 if (--left == 0) 00220 next_state(); 00221 y = state[inext++]; 00222 y ^= (y >> 11); 00223 y ^= (y << 7) & 0x9d2c5680L; 00224 y ^= (y << 15) & 0xefc60000L; 00225 y ^= (y >> 18); 00226 return ((double)y + 0.5) * (1.0/4294967296.0); 00227 } 00228 00229 // generates a random number on [0,1) with 53-bit resolution 00230 public double genrand_res53() { 00231 long a=genrand_int32()>>5, b=genrand_int32()>>6; 00232 return(a*67108864.0+b)*(1.0/9007199254740992.0); 00233 } 00234 00235 // generates a standardized gaussian random number 00236 public double genrand_gaussian() { 00237 int i; 00238 double a; 00239 a=0.0; 00240 for(i=0; i<6; i++) { 00241 a += genrand_real1(); 00242 a -= genrand_real1(); 00243 } 00244 return a; 00245 } 00246 00247 // returns the state of the random number generator 00248 public long[] getState() { 00249 int i; 00250 long[] savedState=new long[627]; 00251 for(i=0; i<624; i++) savedState[i] = state[i]; 00252 savedState[624] = (long) left; 00253 savedState[625] = (long) initf; 00254 savedState[626] = (long) inext; 00255 return savedState; 00256 } 00257 00258 // restores the state of the random number generator 00259 public void setState(long[] savedState) { 00260 int i; 00261 for(i=0; i<624; i++) state[i] = savedState[i]; 00262 left = (int) savedState[624]; 00263 initf = (int) savedState[625]; 00264 inext = (int) savedState[626]; 00265 } 00266 00267 /* 00268 // main program for testing as a standalone program 00269 // same variables printed as in the C program plus 00270 // 100000 gaussian draws 00271 public static void main(String[] args) { 00272 int i; 00273 long[] init={0x123, 0x234, 0x345, 0x456}; 00274 long length=4; 00275 MTwister mt0 = new MTwister(); 00276 mt0.init_by_array(init); 00277 LogManager.println("1000 outputs of genrand_int32()"); 00278 for (i=0; i<10; i++) { 00279 LogManager.println(mt0.genrand_int32()); 00280 if (i%5==4) LogManager.println("\n"); 00281 } 00282 LogManager.println("\n1000 outputs of genrand_real2()\n"); 00283 for (i=0; i<10; i++) { 00284 LogManager.println(mt0.genrand_real2()); 00285 if (i%5==4) LogManager.println("\n"); 00286 } 00287 LogManager.println("\n100000 outputs of genrand_gaussian()\n"); 00288 for (i=0; i<10; i++) { 00289 LogManager.println(mt0.genrand_gaussian()); 00290 if (i%5==4) LogManager.println("\n"); 00291 } 00292 } 00293 */ 00294 } 00295