1 | // this program is a JavaScript version of Mersenne Twister, with concealment and encapsulation in class,
|
2 | // an almost straight conversion from the original program, mt19937ar.c,
|
3 | // translated by y. okada on July 17, 2006.
|
4 | // and modified a little at july 20, 2006, but there are not any substantial differences.
|
5 | // in this program, procedure descriptions and comments of original source code were not removed.
|
6 | // lines commented with //c// were originally descriptions of c procedure. and a few following lines are appropriate JavaScript descriptions.
|
7 | // lines commented with /* and */ are original comments.
|
8 | // lines commented with // are additional comments in this JavaScript version.
|
9 | // before using this version, create at least one instance of MersenneTwister19937 class, and initialize the each state, given below in c comments, of all the instances.
|
10 | /*
|
11 | A C-program for MT19937, with initialization improved 2002/1/26.
|
12 | Coded by Takuji Nishimura and Makoto Matsumoto.
|
13 |
|
14 | Before using, initialize the state by using init_genrand(seed)
|
15 | or init_by_array(init_key, key_length).
|
16 |
|
17 | Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
|
18 | All rights reserved.
|
19 |
|
20 | Redistribution and use in source and binary forms, with or without
|
21 | modification, are permitted provided that the following conditions
|
22 | are met:
|
23 |
|
24 | 1. Redistributions of source code must retain the above copyright
|
25 | notice, this list of conditions and the following disclaimer.
|
26 |
|
27 | 2. Redistributions in binary form must reproduce the above copyright
|
28 | notice, this list of conditions and the following disclaimer in the
|
29 | documentation and/or other materials provided with the distribution.
|
30 |
|
31 | 3. The names of its contributors may not be used to endorse or promote
|
32 | products derived from this software without specific prior written
|
33 | permission.
|
34 |
|
35 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
36 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
37 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
38 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
|
39 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
40 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
41 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
42 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
43 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
44 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
45 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
46 |
|
47 |
|
48 | Any feedback is very welcome.
|
49 | http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
|
50 | email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
|
51 | */
|
52 |
|
53 | function MersenneTwister19937()
|
54 | {
|
55 | /* constants should be scoped inside the class */
|
56 | var N, M, MATRIX_A, UPPER_MASK, LOWER_MASK;
|
57 | /* Period parameters */
|
58 | //c//#define N 624
|
59 | //c//#define M 397
|
60 | //c//#define MATRIX_A 0x9908b0dfUL /* constant vector a */
|
61 | //c//#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
|
62 | //c//#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
|
63 | N = 624;
|
64 | M = 397;
|
65 | MATRIX_A = 0x9908b0df; /* constant vector a */
|
66 | UPPER_MASK = 0x80000000; /* most significant w-r bits */
|
67 | LOWER_MASK = 0x7fffffff; /* least significant r bits */
|
68 | //c//static unsigned long mt[N]; /* the array for the state vector */
|
69 | //c//static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
|
70 | var mt = new Array(N); /* the array for the state vector */
|
71 | var mti = N+1; /* mti==N+1 means mt[N] is not initialized */
|
72 |
|
73 | function unsigned32 (n1) // returns a 32-bits unsiged integer from an operand to which applied a bit operator.
|
74 | {
|
75 | return n1 < 0 ? (n1 ^ UPPER_MASK) + UPPER_MASK : n1;
|
76 | }
|
77 |
|
78 | function subtraction32 (n1, n2) // emulates lowerflow of a c 32-bits unsiged integer variable, instead of the operator -. these both arguments must be non-negative integers expressible using unsigned 32 bits.
|
79 | {
|
80 | return n1 < n2 ? unsigned32((0x100000000 - (n2 - n1)) & 0xffffffff) : n1 - n2;
|
81 | }
|
82 |
|
83 | function addition32 (n1, n2) // emulates overflow of a c 32-bits unsiged integer variable, instead of the operator +. these both arguments must be non-negative integers expressible using unsigned 32 bits.
|
84 | {
|
85 | return unsigned32((n1 + n2) & 0xffffffff)
|
86 | }
|
87 |
|
88 | function multiplication32 (n1, n2) // emulates overflow of a c 32-bits unsiged integer variable, instead of the operator *. these both arguments must be non-negative integers expressible using unsigned 32 bits.
|
89 | {
|
90 | var sum = 0;
|
91 | for (var i = 0; i < 32; ++i){
|
92 | if ((n1 >>> i) & 0x1){
|
93 | sum = addition32(sum, unsigned32(n2 << i));
|
94 | }
|
95 | }
|
96 | return sum;
|
97 | }
|
98 |
|
99 | /* initializes mt[N] with a seed */
|
100 | //c//void init_genrand(unsigned long s)
|
101 | this.init_genrand = function (s)
|
102 | {
|
103 | //c//mt[0]= s & 0xffffffff;
|
104 | mt[0]= unsigned32(s & 0xffffffff);
|
105 | for (mti=1; mti<N; mti++) {
|
106 | mt[mti] =
|
107 | //c//(1812433253 * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
|
108 | addition32(multiplication32(1812433253, unsigned32(mt[mti-1] ^ (mt[mti-1] >>> 30))), mti);
|
109 | /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
|
110 | /* In the previous versions, MSBs of the seed affect */
|
111 | /* only MSBs of the array mt[]. */
|
112 | /* 2002/01/09 modified by Makoto Matsumoto */
|
113 | //c//mt[mti] &= 0xffffffff;
|
114 | mt[mti] = unsigned32(mt[mti] & 0xffffffff);
|
115 | /* for >32 bit machines */
|
116 | }
|
117 | }
|
118 |
|
119 | /* initialize by an array with array-length */
|
120 | /* init_key is the array for initializing keys */
|
121 | /* key_length is its length */
|
122 | /* slight change for C++, 2004/2/26 */
|
123 | //c//void init_by_array(unsigned long init_key[], int key_length)
|
124 | this.init_by_array = function (init_key, key_length)
|
125 | {
|
126 | //c//int i, j, k;
|
127 | var i, j, k;
|
128 | //c//init_genrand(19650218);
|
129 | this.init_genrand(19650218);
|
130 | i=1; j=0;
|
131 | k = (N>key_length ? N : key_length);
|
132 | for (; k; k--) {
|
133 | //c//mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525))
|
134 | //c// + init_key[j] + j; /* non linear */
|
135 | mt[i] = addition32(addition32(unsigned32(mt[i] ^ multiplication32(unsigned32(mt[i-1] ^ (mt[i-1] >>> 30)), 1664525)), init_key[j]), j);
|
136 | mt[i] =
|
137 | //c//mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
|
138 | unsigned32(mt[i] & 0xffffffff);
|
139 | i++; j++;
|
140 | if (i>=N) { mt[0] = mt[N-1]; i=1; }
|
141 | if (j>=key_length) {j=0;}
|
142 | }
|
143 | for (k=N-1; k; k--) {
|
144 | //c//mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941))
|
145 | //c//- i; /* non linear */
|
146 | mt[i] = subtraction32(unsigned32((dbg=mt[i]) ^ multiplication32(unsigned32(mt[i-1] ^ (mt[i-1] >>> 30)), 1566083941)), i);
|
147 | //c//mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
|
148 | mt[i] = unsigned32(mt[i] & 0xffffffff);
|
149 | i++;
|
150 | if (i>=N) { mt[0] = mt[N-1]; i=1; }
|
151 | }
|
152 | mt[0] = 0x80000000; /* MSB is 1; assuring non-zero initial array */
|
153 | }
|
154 |
|
155 | /* moved outside of genrand_int32() by jwatte 2010-11-17; generate less garbage */
|
156 | var mag01 = [0x0, MATRIX_A];
|
157 |
|
158 | /* generates a random number on [0,0xffffffff]-interval */
|
159 | //c//unsigned long genrand_int32(void)
|
160 | this.genrand_int32 = function ()
|
161 | {
|
162 | //c//unsigned long y;
|
163 | //c//static unsigned long mag01[2]={0x0UL, MATRIX_A};
|
164 | var y;
|
165 | /* mag01[x] = x * MATRIX_A for x=0,1 */
|
166 |
|
167 | if (mti >= N) { /* generate N words at one time */
|
168 | //c//int kk;
|
169 | var kk;
|
170 |
|
171 | if (mti == N+1) /* if init_genrand() has not been called, */
|
172 | //c//init_genrand(5489); /* a default initial seed is used */
|
173 | {this.init_genrand(5489);} /* a default initial seed is used */
|
174 |
|
175 | for (kk=0;kk<N-M;kk++) {
|
176 | //c//y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
|
177 | //c//mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
|
178 | y = unsigned32((mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK));
|
179 | mt[kk] = unsigned32(mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1]);
|
180 | }
|
181 | for (;kk<N-1;kk++) {
|
182 | //c//y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
|
183 | //c//mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
|
184 | y = unsigned32((mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK));
|
185 | mt[kk] = unsigned32(mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1]);
|
186 | }
|
187 | //c//y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
|
188 | //c//mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
|
189 | y = unsigned32((mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK));
|
190 | mt[N-1] = unsigned32(mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1]);
|
191 | mti = 0;
|
192 | }
|
193 |
|
194 | y = mt[mti++];
|
195 |
|
196 | /* Tempering */
|
197 | //c//y ^= (y >> 11);
|
198 | //c//y ^= (y << 7) & 0x9d2c5680;
|
199 | //c//y ^= (y << 15) & 0xefc60000;
|
200 | //c//y ^= (y >> 18);
|
201 | y = unsigned32(y ^ (y >>> 11));
|
202 | y = unsigned32(y ^ ((y << 7) & 0x9d2c5680));
|
203 | y = unsigned32(y ^ ((y << 15) & 0xefc60000));
|
204 | y = unsigned32(y ^ (y >>> 18));
|
205 |
|
206 | return y;
|
207 | }
|
208 |
|
209 | /* generates a random number on [0,0x7fffffff]-interval */
|
210 | //c//long genrand_int31(void)
|
211 | this.genrand_int31 = function ()
|
212 | {
|
213 | //c//return (genrand_int32()>>1);
|
214 | return (this.genrand_int32()>>>1);
|
215 | }
|
216 |
|
217 | /* generates a random number on [0,1]-real-interval */
|
218 | //c//double genrand_real1(void)
|
219 | this.genrand_real1 = function ()
|
220 | {
|
221 | //c//return genrand_int32()*(1.0/4294967295.0);
|
222 | return this.genrand_int32()*(1.0/4294967295.0);
|
223 | /* divided by 2^32-1 */
|
224 | }
|
225 |
|
226 | /* generates a random number on [0,1)-real-interval */
|
227 | //c//double genrand_real2(void)
|
228 | this.genrand_real2 = function ()
|
229 | {
|
230 | //c//return genrand_int32()*(1.0/4294967296.0);
|
231 | return this.genrand_int32()*(1.0/4294967296.0);
|
232 | /* divided by 2^32 */
|
233 | }
|
234 |
|
235 | /* generates a random number on (0,1)-real-interval */
|
236 | //c//double genrand_real3(void)
|
237 | this.genrand_real3 = function ()
|
238 | {
|
239 | //c//return ((genrand_int32()) + 0.5)*(1.0/4294967296.0);
|
240 | return ((this.genrand_int32()) + 0.5)*(1.0/4294967296.0);
|
241 | /* divided by 2^32 */
|
242 | }
|
243 |
|
244 | /* generates a random number on [0,1) with 53-bit resolution*/
|
245 | //c//double genrand_res53(void)
|
246 | this.genrand_res53 = function ()
|
247 | {
|
248 | //c//unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
|
249 | var a=this.genrand_int32()>>>5, b=this.genrand_int32()>>>6;
|
250 | return(a*67108864.0+b)*(1.0/9007199254740992.0);
|
251 | }
|
252 | /* These real versions are due to Isaku Wada, 2002/01/09 added */
|
253 | }
|
254 |
|
255 | // Exports: Public API
|
256 |
|
257 | // Export the twister class
|
258 | exports.MersenneTwister19937 = MersenneTwister19937;
|