KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > org > jahia > utils > keygenerator > MersenneTwisterFast


1 package org.jahia.utils.keygenerator;
2
3
4 import java.io.Serializable JavaDoc;
5
6 /**
7  * Mersenne Twister and MersenneTwisterFast:
8  * <P>
9  * <b>MersenneTwisterFast</b> is a drop-in subclass replacement
10  * for java.util.Random. It is properly synchronized and
11  * can be used in a multithreaded environment.
12  *
13  * <p><b>MersenneTwisterFast</b> is not a subclass of java.util.Random. It has
14  * the same public methods as Random does, however, and it is
15  * algorithmically identical to MersenneTwister. MersenneTwisterFast
16  * has hard-code inlined all of its methods directly, and made all of them
17  * final (well, the ones of consequence anyway). Further, these
18  * methods are <i>not</i> synchronized, so the same MersenneTwisterFast
19  * instance cannot be shared by multiple threads. But all this helps
20  * MersenneTwisterFast achieve over twice the speed of MersenneTwister.
21  *
22  * <p><b>About the Mersenne Twister. </b>
23  * This is a Java version of the C-program for MT19937: Integer version.
24  * next(32) generates one pseudorandom unsigned integer (32bit)
25
26  * which is uniformly distributed among 0 to 2^32-1 for each
27  * call. next(int bits) >>>'s by (32-bits) to get a value ranging
28  * between 0 and 2^bits-1 long inclusive; hope that's correct.
29  * setSeed(seed) set initial values to the working area
30  * of 624 words. For setSeed(seed), seed is any 32-bit integer
31  * <b>except for 0</b>.
32  *
33  * <p>Orignally Coded by Takuji Nishimura, considering the suggestions by
34  * Topher Cooper and Marc Rieffel in July-Aug. 1997.
35  * More information can be found
36  * <A HREF="http://www.math.keio.ac.jp/matumoto/emt.html">
37  * here. </a>
38
39  * <P>
40  * Translated to Java by Michael Lecuyer January 30, 1999
41  * Copyright (C) 1999 Michael Lecuyer
42  * <P>
43  * This library is free software; you can redistribute it and or
44
45  * modify it under the terms of the GNU Library General Public
46  * License as published by the Free Software Foundation; either
47  * version 2 of the License, or (at your option) any later
48  * version.
49  * This library is distributed in the hope that it will be useful,
50
51  * but WITHOUT ANY WARRANTY; without even the implied warranty of
52  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
53  * See the GNU Library General Public License for more details.
54  * You should have received a copy of the GNU Library General
55  * Public License along with this library; if not, write to the
56  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
57  * 02111-1307 USA
58  * <P>
59  * Makoto Matsumoto and Takuji Nishimura, the original authors
60  * ask "When you use this, send an email to: matumoto@math.keio.ac.jp
61  * with an appropriate reference to your work" You might also point
62  * out this was a translation.
63  * <P>
64  * <b>Reference. </b>
65  * M. Matsumoto and T. Nishimura,
66  * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
67  * Pseudo-Random Number Generator",
68  * <i>ACM Transactions on Modeling and Computer Simulation,</i>
69  * Vol. 8, No. 1, January 1998, pp 3--30.
70  *
71  * <p><b>About this version. </b> This is a modification of the
72  * <a HREF="http://www.theorem.com/java/index.htm#Mersenne">original
73  * code</a> made to conform to proper java.util.Random format by
74  * <a HREF="http://www.cs.umd.edu/users/seanl/">Sean Luke,</a>
75  * August 7, 1999.
76  *
77  * <p><b>Bug Fixes. </b>This implementation implements the bug fixes made
78  * in Java 1.2's version of Random, which means it can be used with
79  * earlier versions of Java. See
80  * <a HREF="http://www.javasoft.com/products/jdk/1.2/docs/api/java/util/Random.html">
81  * the JDK 1.2 java.util.Random documentation</a> for further documentation
82  * on the random-number generation contracts made. Additionally, there's
83  * an undocumented bug in the JDK java.util.Random.nextBytes() method,
84  * which this code fixes.
85  *
86  * <p><b>Important Note. </b> Just like java.util.Random, this
87  * generator accepts a long seed but doesn't use all of it. java.util.Random
88  * uses 48 bits. The Mersenne Twister instead uses 32 bits (int size).
89  * So it's best if your seed does not exceed the int range.
90  */
91
92
93
94 public class MersenneTwisterFast implements Serializable JavaDoc
95     {
96     // Period parameters
97 private static final int N = 624;
98     private static final int M = 397;
99     private static final int MATRIX_A = 0x9908b0df; // private static final * constant vector a
100 private static final int UPPER_MASK = 0x80000000; // most significant w-r bits
101 private static final int LOWER_MASK = 0x7fffffff; // least significant r bits
102
103
104     // Tempering parameters
105 private static final int TEMPERING_MASK_B = 0x9d2c5680;
106     private static final int TEMPERING_MASK_C = 0xefc60000;
107     
108     // #define TEMPERING_SHIFT_U(y) (y >>> 11)
109 // #define TEMPERING_SHIFT_S(y) (y << 7)
110 // #define TEMPERING_SHIFT_T(y) (y << 15)
111 // #define TEMPERING_SHIFT_L(y) (y >>> 18)
112
113     private int mt[]; // the array for the state vector
114 private int mti; // mti==N+1 means mt[N] is not initialized
115 private int mag01[];
116     
117     // a good initial seed (of int size, though stored in a long)
118 private static final long GOOD_SEED = 4357;
119
120     private double nextNextGaussian;
121     private boolean haveNextNextGaussian;
122
123
124     /**
125      * Constructor using the default seed.
126
127      */
128     public MersenneTwisterFast()
129         {
130         setSeed(GOOD_SEED);
131         }
132     
133     /**
134      * Constructor using a given seed. Though you pass this seed in
135      * as a long, it's best to make sure it's actually an integer.
136      *
137      * @param seed generator starting number, often the time of day.
138      */
139     public MersenneTwisterFast(long seed)
140         {
141         setSeed(seed);
142         }
143     
144     /**
145      * Initalize the pseudo random number generator.
146      * The Mersenne Twister only uses an integer for its seed;
147      * It's best that you don't pass in a long that's bigger
148      * than an int.
149      *
150      * @param seed from constructor
151      *
152      */
153     public final void setSeed(long seed)
154         {
155         haveNextNextGaussian = false;
156
157         mt = new int[N];
158         
159         // setting initial seeds to mt[N] using
160 // the generator Line 25 of Table 1 in
161 // [KNUTH 1981, The Art of Computer Programming
162 // Vol. 2 (2nd Ed.), pp102]
163
164         // the 0xffffffff is commented out because in Java
165 // ints are always 32 bits; hence i & 0xffffffff == i
166
167         mt[0]= ((int)seed); // & 0xffffffff;
168
169         for (mti = 1; mti < N; mti++)
170             mt[mti] = (69069 * mt[mti-1]); //& 0xffffffff;
171
172         // mag01[x] = x * MATRIX_A for x=0,1
173 mag01 = new int[2];
174         mag01[0] = 0x0;
175         mag01[1] = MATRIX_A;
176         }
177     
178     public final int nextInt()
179         {
180         int y;
181         
182         if (mti >= N) // generate N words at one time
183 {
184             int kk;
185             
186             for (kk = 0; kk < N - M; kk++)
187                 {
188                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
189                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
190                 }
191             for (; kk < N-1; kk++)
192                 {
193                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
194                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
195                 }
196             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
197             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
198
199             mti = 0;
200             }
201
202   
203         y = mt[mti++];
204         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
205 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
206 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
207 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
208
209         return y;
210         }
211
212
213
214     public final short nextShort()
215         {
216         int y;
217         
218         if (mti >= N) // generate N words at one time
219 {
220             int kk;
221             
222             for (kk = 0; kk < N - M; kk++)
223                 {
224                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
225                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
226                 }
227             for (; kk < N-1; kk++)
228                 {
229                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
230                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
231                 }
232             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
233             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
234
235
236             mti = 0;
237             }
238   
239         y = mt[mti++];
240         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
241 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
242 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
243 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
244
245         return (short)(y >>> 16);
246         }
247
248
249
250     public final char nextChar()
251         {
252         int y;
253         
254         if (mti >= N) // generate N words at one time
255
256             {
257             int kk;
258             
259             for (kk = 0; kk < N - M; kk++)
260                 {
261                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
262                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
263                 }
264             for (; kk < N-1; kk++)
265                 {
266                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
267                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
268                 }
269             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
270             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
271
272             mti = 0;
273             }
274   
275         y = mt[mti++];
276         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
277 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
278 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
279 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
280
281         return (char)(y >>> 16);
282         }
283
284
285     public final boolean nextBoolean()
286         {
287         int y;
288         
289         if (mti >= N) // generate N words at one time
290 {
291             int kk;
292             
293             for (kk = 0; kk < N - M; kk++)
294                 {
295                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
296                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
297                 }
298
299             for (; kk < N-1; kk++)
300                 {
301                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
302
303                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
304                 }
305             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
306             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
307
308             mti = 0;
309             }
310   
311         y = mt[mti++];
312         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
313 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
314 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
315 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
316
317         return (boolean)((y >>> 31) != 0);
318         }
319
320
321     public final byte nextByte()
322         {
323         int y;
324         
325         if (mti >= N) // generate N words at one time
326 {
327             int kk;
328             
329             for (kk = 0; kk < N - M; kk++)
330                 {
331                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
332                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
333                 }
334             for (; kk < N-1; kk++)
335                 {
336                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
337                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
338                 }
339             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
340             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
341
342             mti = 0;
343             }
344   
345         y = mt[mti++];
346         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
347 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
348 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
349 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
350
351         return (byte)(y >>> 24);
352         }
353
354
355     public final void nextBytes(byte[] bytes)
356         {
357         int y;
358         
359         for (int x=0;x<bytes.length;x++)
360             {
361             if (mti >= N) // generate N words at one time
362 {
363                 int kk;
364                 
365                 for (kk = 0; kk < N - M; kk++)
366                     {
367                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
368                     mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
369                     }
370                 for (; kk < N-1; kk++)
371                     {
372                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
373                     mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
374                     }
375                 y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
376                 mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
377                 
378                 mti = 0;
379                 }
380             
381             y = mt[mti++];
382             y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
383 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
384 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
385 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
386
387             bytes[x] = (byte)(y >>> 24);
388             }
389         }
390
391
392     public final long nextLong()
393         {
394         int y;
395         int z;
396
397         if (mti >= N) // generate N words at one time
398 {
399             int kk;
400             
401             for (kk = 0; kk < N - M; kk++)
402                 {
403                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
404                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
405                 }
406             for (; kk < N-1; kk++)
407                 {
408                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
409                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
410                 }
411             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
412             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
413
414
415             mti = 0;
416             }
417   
418         y = mt[mti++];
419         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
420 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
421 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
422 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
423
424         if (mti >= N) // generate N words at one time
425
426             {
427             int kk;
428             
429             for (kk = 0; kk < N - M; kk++)
430                 {
431                 z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
432                 mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
433                 }
434             for (; kk < N-1; kk++)
435                 {
436                 z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
437                 mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
438                 }
439             z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
440             mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
441             
442             mti = 0;
443             }
444         
445         z = mt[mti++];
446         z ^= z >>> 11; // TEMPERING_SHIFT_U(z)
447 z ^= (z << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(z)
448 z ^= (z << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(z)
449 z ^= (z >>> 18); // TEMPERING_SHIFT_L(z)
450
451         return (((long)y) << 32) + (long)z;
452         }
453
454
455     public final double nextDouble()
456         {
457         int y;
458         int z;
459
460         if (mti >= N) // generate N words at one time
461 {
462             int kk;
463             
464             for (kk = 0; kk < N - M; kk++)
465                 {
466
467                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
468                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
469                 }
470             for (; kk < N-1; kk++)
471                 {
472                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
473                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
474                 }
475             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
476             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
477
478             mti = 0;
479             }
480   
481         y = mt[mti++];
482         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
483 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
484 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
485 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
486
487         if (mti >= N) // generate N words at one time
488 {
489             int kk;
490             
491             for (kk = 0; kk < N - M; kk++)
492                 {
493                 z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
494                 mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
495                 }
496             for (; kk < N-1; kk++)
497                 {
498                 z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
499                 mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
500                 }
501             z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
502
503             mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
504             
505             mti = 0;
506             }
507         
508         z = mt[mti++];
509         z ^= z >>> 11; // TEMPERING_SHIFT_U(z)
510 z ^= (z << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(z)
511 z ^= (z << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(z)
512 z ^= (z >>> 18); // TEMPERING_SHIFT_L(z)
513
514         /* derived from nextDouble documentation in jdk 1.2 docs, see top */
515         return ((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53);
516         }
517
518
519
520
521
522     public final double nextGaussian()
523         {
524         if (haveNextNextGaussian)
525             {
526             haveNextNextGaussian = false;
527             return nextNextGaussian;
528             }
529         else
530             {
531             double v1, v2, s;
532             do
533                 {
534                 int y;
535                 int z;
536                 int a;
537                 int b;
538                     
539                     if (mti >= N) // generate N words at one time
540 {
541                         int kk;
542                         
543                         for (kk = 0; kk < N - M; kk++)
544
545                             {
546                             y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
547                             mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
548                             }
549                         for (; kk < N-1; kk++)
550                             {
551                             y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
552                             mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
553                             }
554                         y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
555
556                         mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
557                         
558                         mti = 0;
559                         }
560                 
561                 y = mt[mti++];
562                 y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
563 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
564 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
565 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
566
567                 if (mti >= N) // generate N words at one time
568 {
569                     int kk;
570                     
571                     for (kk = 0; kk < N - M; kk++)
572                         {
573                         z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
574                         mt[kk] = mt[kk+M] ^ (z >>> 1) ^ mag01[z & 0x1];
575                         }
576                     for (; kk < N-1; kk++)
577                         {
578                         z = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
579                         mt[kk] = mt[kk+(M-N)] ^ (z >>> 1) ^ mag01[z & 0x1];
580                         }
581                     z = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
582                     mt[N-1] = mt[M-1] ^ (z >>> 1) ^ mag01[z & 0x1];
583                     
584                     mti = 0;
585                     }
586                 
587                 z = mt[mti++];
588                 z ^= z >>> 11; // TEMPERING_SHIFT_U(z)
589 z ^= (z << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(z)
590 z ^= (z << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(z)
591 z ^= (z >>> 18); // TEMPERING_SHIFT_L(z)
592
593                 if (mti >= N) // generate N words at one time
594 {
595                     int kk;
596                     
597                     for (kk = 0; kk < N - M; kk++)
598                         {
599                         a = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
600                         mt[kk] = mt[kk+M] ^ (a >>> 1) ^ mag01[a & 0x1];
601                         }
602                     for (; kk < N-1; kk++)
603                         {
604
605                         a = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
606                         mt[kk] = mt[kk+(M-N)] ^ (a >>> 1) ^ mag01[a & 0x1];
607                         }
608                     a = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
609                     mt[N-1] = mt[M-1] ^ (a >>> 1) ^ mag01[a & 0x1];
610                     
611                     mti = 0;
612                     }
613                 
614                 a = mt[mti++];
615                 a ^= a >>> 11; // TEMPERING_SHIFT_U(a)
616 a ^= (a << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(a)
617 a ^= (a << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(a)
618 a ^= (a >>> 18); // TEMPERING_SHIFT_L(a)
619
620                 if (mti >= N) // generate N words at one time
621 {
622                     int kk;
623                     
624                     for (kk = 0; kk < N - M; kk++)
625                         {
626                         b = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
627                         mt[kk] = mt[kk+M] ^ (b >>> 1) ^ mag01[b & 0x1];
628                         }
629                     for (; kk < N-1; kk++)
630                         {
631                         b = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
632                         mt[kk] = mt[kk+(M-N)] ^ (b >>> 1) ^ mag01[b & 0x1];
633                         }
634                     b = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
635                     mt[N-1] = mt[M-1] ^ (b >>> 1) ^ mag01[b & 0x1];
636                     
637                     mti = 0;
638                     }
639                 
640                 b = mt[mti++];
641                 b ^= b >>> 11; // TEMPERING_SHIFT_U(b)
642 b ^= (b << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(b)
643 b ^= (b << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(b)
644 b ^= (b >>> 18); // TEMPERING_SHIFT_L(b)
645
646                 /* derived from nextDouble documentation in jdk 1.2 docs, see top */
647                 v1 = 2 *
648                     (((((long)(y >>> 6)) << 27) + (z >>> 5)) / (double)(1L << 53))
649                     - 1;
650                 v2 = 2 * (((((long)(a >>> 6)) << 27) + (b >>> 5)) / (double)(1L << 53))
651                     - 1;
652                 s = v1 * v1 + v2 * v2;
653                 } while (s >= 1);
654             double multiplier = Math.sqrt(-2 * Math.log(s)/s);
655             nextNextGaussian = v2 * multiplier;
656             haveNextNextGaussian = true;
657             return v1 * multiplier;
658             }
659         }
660     
661     
662     
663
664     
665
666
667
668
669
670
671
672
673     public final float nextFloat()
674         {
675         int y;
676         
677         if (mti >= N) // generate N words at one time
678 {
679             int kk;
680             
681             for (kk = 0; kk < N - M; kk++)
682                 {
683                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
684                 mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
685                 }
686             for (; kk < N-1; kk++)
687                 {
688                 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
689                 mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
690                 }
691             y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
692             mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
693
694             mti = 0;
695
696             }
697   
698         y = mt[mti++];
699         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
700 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
701 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
702 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
703
704         return (y >>> 8) / ((float)(1 << 24));
705         }
706
707
708
709     /** Returns an integer drawn uniformly from 0 to n-1. Suffice it to say,
710         n must be > 0, or an IllegalArgumentException is raised. */
711     public int nextInt(int n)
712         {
713         if (n<=0)
714             throw new IllegalArgumentException JavaDoc("n must be positive");
715         
716         if ((n & -n) == n) // i.e., n is a power of 2
717 {
718             int y;
719         
720             if (mti >= N) // generate N words at one time
721 {
722                 int kk;
723                 
724                 for (kk = 0; kk < N - M; kk++)
725                     {
726                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
727                     mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
728                     }
729                 for (; kk < N-1; kk++)
730                     {
731                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
732                     mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
733                     }
734                 y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
735                 mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
736                 
737                 mti = 0;
738                 }
739             
740             y = mt[mti++];
741             y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
742 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
743 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
744 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
745
746             return (int)((n * (long) (y >>> 1) ) >> 31);
747             }
748         
749         int bits, val;
750         do
751             {
752             int y;
753             
754             if (mti >= N) // generate N words at one time
755 {
756                 int kk;
757                 
758                 for (kk = 0; kk < N - M; kk++)
759                     {
760                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
761                     mt[kk] = mt[kk+M] ^ (y >>> 1) ^ mag01[y & 0x1];
762                     }
763                 for (; kk < N-1; kk++)
764                     {
765                     y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
766                     mt[kk] = mt[kk+(M-N)] ^ (y >>> 1) ^ mag01[y & 0x1];
767                     }
768                 y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
769                 mt[N-1] = mt[M-1] ^ (y >>> 1) ^ mag01[y & 0x1];
770                 
771                 mti = 0;
772                 }
773             
774             y = mt[mti++];
775             y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
776 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
777 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
778 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
779
780             bits = (y >>> 1);
781             val = bits % n;
782             } while(bits - val + (n-1) < 0);
783         return val;
784         }
785     
786
787
788
789
790
791
792
793
794     /**
795      * Tests the code.
796      */
797     public static void main(String JavaDoc args[])
798         {
799         int j;
800
801         MersenneTwisterFast r;
802
803         // UNCOMMENT THIS TO TEST FOR PROPER GAUSSIAN STATE INITIALIZATION
804
805         /*
806         System.out.println("If the gaussian state is properly initialized when setSeed() is called,\nthen #1 != #2, but #1 == #3\nIt's known that java 1.0.2 doesn't do gaussian initialization right,\nso setSeed() may result in one last gaussian drawn from the *previous* seed.");
807         r = new MersenneTwisterFast(1);
808         r.nextGaussian(); // loads the later gaussian into the state
809         System.out.println("1: " + r.nextGaussian());
810         r = new MersenneTwisterFast(1);
811         r.nextGaussian(); // loads the later gaussian into the state
812         r.setSeed(1); // should reset the gaussian state
813         System.out.println("2: " + r.nextGaussian());
814         System.out.println("3: " + r.nextGaussian());
815         */
816
817         
818         // UNCOMMENT THIS TO TEST FOR CORRECTNESS
819 // COMPARE WITH http://www.math.keio.ac.jp/~nisimura/random/int/mt19937int.out
820
821         /*
822         r = new MersenneTwisterFast(4357);
823         System.out.println("Output of MersenneTwisterFast.java");
824         for (j=0;j<1000;j++)
825             {
826             // first, convert the int from signed to "unsigned"
827             long l = (long)r.nextInt();
828
829             if (l < 0 ) l += 4294967296L; // max int value
830             String s = String.valueOf(l);
831
832             while(s.length() < 10) s = " " + s; // buffer
833             System.out.print(s + " ");
834             if (j%8==7) System.out.println();
835             }
836         */
837
838
839         // UNCOMMENT THIS TO TEST FOR SPEED
840
841         /*
842         r = new MersenneTwisterFast();
843         System.out.println("\nTime to test grabbing 10000000 ints");
844         long ms = System.currentTimeMillis();
845         int xx=0;
846         for (j = 0; j < 10000000; j++)
847             xx += r.nextInt();
848         System.out.println("Mersenne Twister: " + (System.currentTimeMillis()-ms + " Ignore this: " + xx));
849
850         Random rr = new Random(1);
851         xx = 0;
852         ms = System.currentTimeMillis();
853         for (j = 0; j < 10000000; j++)
854             xx += rr.nextInt();
855         System.out.println("java.util.Random: " + (System.currentTimeMillis()-ms + " Ignore this: " + xx));
856         */
857
858         
859         // UNCOMMENT THIS TO DO TEST DIFFERENT TYPE OUTPUTS
860 // THIS CAN BE USED TO COMPARE THE DIFFERENCE BETWEEN
861
862         // MersenneTwisterFast.java AND MersenneTwister.java
863
864         /*
865         System.out.println("\nGrab the first 1000 booleans");
866         r = new MersenneTwisterFast();
867         for (j = 0; j < 1000; j++)
868             {
869             System.out.print(r.nextBoolean() + " ");
870             if (j%8==7) System.out.println();
871             }
872         if (!(j%8==7)) System.out.println();
873
874         byte[] bytes = new byte[1000];
875         System.out.println("\nGrab the first 1000 bytes using nextBytes");
876         r = new MersenneTwisterFast();
877         r.nextBytes(bytes);
878         for (j = 0; j < 1000; j++)
879             {
880             System.out.print(bytes[j] + " ");
881             if (j%16==15) System.out.println();
882             }
883         if (!(j%16==15)) System.out.println();
884         
885         byte b;
886         System.out.println("\nGrab the first 1000 bytes -- must be same as nextBytes");
887         r = new MersenneTwisterFast();
888
889         for (j = 0; j < 1000; j++)
890             {
891             System.out.print((b = r.nextByte()) + " ");
892             if (b!=bytes[j]) System.out.print("BAD ");
893             if (j%16==15) System.out.println();
894             }
895         if (!(j%16==15)) System.out.println();
896
897         System.out.println("\nGrab the first 1000 shorts");
898         r = new MersenneTwisterFast();
899         for (j = 0; j < 1000; j++)
900             {
901             System.out.print(r.nextShort() + " ");
902             if (j%8==7) System.out.println();
903             }
904         if (!(j%8==7)) System.out.println();
905
906         System.out.println("\nGrab the first 1000 ints");
907         r = new MersenneTwisterFast();
908         for (j = 0; j < 1000; j++)
909             {
910             System.out.print(r.nextInt() + " ");
911             if (j%4==3) System.out.println();
912             }
913         if (!(j%4==3)) System.out.println();
914
915         System.out.println("\nGrab the first 1000 ints of different sizes");
916         r = new MersenneTwisterFast();
917         for (j = 0; j < 1000; j++)
918             {
919             System.out.print(r.nextInt(j+1) + " ");
920             if (j%4==3) System.out.println();
921             }
922         if (!(j%4==3)) System.out.println();
923
924         System.out.println("\nGrab the first 1000 longs");
925         r = new MersenneTwisterFast();
926         for (j = 0; j < 1000; j++)
927             {
928             System.out.print(r.nextLong() + " ");
929             if (j%3==2) System.out.println();
930             }
931         if (!(j%3==2)) System.out.println();
932
933         System.out.println("\nGrab the first 1000 floats");
934         r = new MersenneTwisterFast();
935         for (j = 0; j < 1000; j++)
936             {
937             System.out.print(r.nextFloat() + " ");
938             if (j%4==3) System.out.println();
939             }
940         if (!(j%4==3)) System.out.println();
941
942         System.out.println("\nGrab the first 1000 doubles");
943         r = new MersenneTwisterFast();
944         for (j = 0; j < 1000; j++)
945             {
946             System.out.print(r.nextDouble() + " ");
947             if (j%3==2) System.out.println();
948             }
949         if (!(j%3==2)) System.out.println();
950
951         System.out.println("\nGrab the first 1000 gaussian doubles");
952         r = new MersenneTwisterFast();
953         for (j = 0; j < 1000; j++)
954             {
955             System.out.print(r.nextGaussian() + " ");
956             if (j%3==2) System.out.println();
957             }
958         if (!(j%3==2)) System.out.println();
959         */
960         }
961
962     }
963
Popular Tags