1 package nl.tudelft.simulation.jstats.streams;
2
3 /**
4 * A java implementation of the Mersenne Twister pseudo random number generator.
5 * <p>
6 * This is a Java version of the C-program for MT19937: Integer version. genrand() generates one pseudorandom unsigned integer
7 * (32bit) which is uniformly distributed among 0 to 2^32-1 for each call. sgenrand(seed) set initial values to the working area
8 * of 624 words. (seed is any 32-bit integer except for 0).
9 * </p>
10 * <p>
11 * Orignally Coded by Takuji Nishimura, considering the suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997. More
12 * information can be found at <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">
13 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf</a> and
14 * <a href="https://en.wikipedia.org/wiki/Mersenne_Twister"> https://en.wikipedia.org/wiki/Mersenne_Twister</a>.
15 * </p>
16 * <p>
17 * Copyright (c) 2002-2025 Delft University of Technology, Jaffalaan 5, 2628 BX Delft, the Netherlands. All rights reserved. See
18 * for project information <a href="https://simulation.tudelft.nl/dsol/manual/" target="_blank">DSOL Manual</a>. The DSOL
19 * project is distributed under a three-clause BSD-style license, which can be found at
20 * <a href="https://simulation.tudelft.nl/dsol/docs/latest/license.html" target="_blank">DSOL License</a>.
21 * </p>
22 * @author <a href="https://github.com/averbraeck" target="_blank"> Alexander Verbraeck</a>
23 * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
24 */
25 public class MersenneTwister extends RandomNumberGenerator
26 {
27 /** Period parameter N. */
28 private static final int N = 624;
29
30 /** Period parameter M. */
31 private static final int M = 397;
32
33 /** constant vector a. */
34 private static final int MATRIX_A = 0x9908b0df;
35
36 /** most significant w-r bits. */
37 private static final int UPPER_MASK = 0x80000000;
38
39 /** least significant w-r bits. */
40 private static final int LOWER_MASK = 0x7fffffff;
41
42 /** tempering mask b. */
43 private static final int TEMPERING_MASK_B = 0x9d2c5680;
44
45 /** tempering mask b. */
46 private static final int TEMPERING_MASK_C = 0xefc60000;
47
48 /** unsigned mask for promoting int -> long. */
49 private static final int UMASK = (1 << 31) - 1;
50
51 /** the array for the state vector. */
52 private int[] mt;
53
54 /** The counter mti==N+1 means mt[N] is not initialized. */
55 private int mti;
56
57 /** magic01. */
58 private int[] mag01;
59
60 /**
61 * constructs a new Mersenne Twister. <code>System.currentTimeMillis()</code> is used as seed value.
62 */
63 public MersenneTwister()
64 {
65 this(System.currentTimeMillis());
66 }
67
68 /**
69 * Constructor using a given seed.
70 * @param seed The initial seed.
71 */
72 public MersenneTwister(final long seed)
73 {
74 super(seed);
75 }
76
77 /**
78 * initalizes the MersenneTwister.
79 */
80 private void initialize()
81 {
82 this.mt = new int[N];
83
84 this.mt[0] = ((int) super.seed) & UMASK;
85 if (this.mt[0] == 0)
86 {
87 // super.seed=Integer.MAXValue --> seed & UMASK==0
88 // We set the seed again and enforce a different value.
89 this.setSeed(System.currentTimeMillis());
90 }
91 for (this.mti = 1; this.mti < N; this.mti++)
92 {
93 this.mt[this.mti] = (69069 * this.mt[this.mti - 1]);
94 }
95
96 // mag01[x] = x * MATRIX_A for x=0,1
97 this.mag01 = new int[2];
98 this.mag01[0] = 0x0;
99 this.mag01[1] = MATRIX_A;
100 }
101
102 @Override
103 public synchronized long next(final int bits)
104 {
105 if (bits < 0 || bits > 64)
106 {
107 throw new IllegalArgumentException("bits (" + bits + ") not in range [0,64]");
108 }
109 int y;
110 if (this.mti >= N) // generate N words at one time
111 {
112 int kk;
113 for (kk = 0; kk < N - M; kk++)
114 {
115 y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK);
116 this.mt[kk] = this.mt[kk + M] ^ (y >>> 1) ^ this.mag01[y & 0x1];
117 }
118 for (; kk < N - 1; kk++)
119 {
120 y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK);
121 this.mt[kk] = this.mt[kk + (M - N)] ^ (y >>> 1) ^ this.mag01[y & 0x1];
122 }
123 y = (this.mt[N - 1] & UPPER_MASK) | (this.mt[0] & LOWER_MASK);
124 this.mt[N - 1] = this.mt[M - 1] ^ (y >>> 1) ^ this.mag01[y & 0x1];
125 this.mti = 0;
126 }
127 y = this.mt[this.mti++];
128 y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
129 y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
130 y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
131 y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
132 if (bits <= 32)
133 {
134 return y >>> (32 - bits);
135 }
136 return y << 32 + this.next(bits - 32);
137 }
138
139 @Override
140 public synchronized void setSeed(final long seed)
141 {
142 super.seed = seed;
143 this.initialize();
144 }
145 }