View Javadoc
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-2024 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/" target="_blank"> https://simulation.tudelft.nl</a>. The DSOL
19   * project is distributed under a three-clause BSD-style license, which can be found at
20   * <a href="https://https://simulation.tudelft.nl/dsol/docs/latest/license.html" target="_blank">
21   * https://https://simulation.tudelft.nl/dsol/docs/latest/license.html</a>.
22   * </p>
23   * @author <a href="https://www.tudelft.nl/averbraeck" target="_blank"> Alexander Verbraeck</a>
24   * @author <a href="https://www.linkedin.com/in/peterhmjacobs">Peter Jacobs </a>
25   */
26  public class MersenneTwister extends RandomNumberGenerator
27  {
28      /** */
29      private static final long serialVersionUID = 20150426L;
30  
31      /** Period parameter N. */
32      private static final int N = 624;
33  
34      /** Period parameter M. */
35      private static final int M = 397;
36  
37      /** constant vector a. */
38      private static final int MATRIX_A = 0x9908b0df;
39  
40      /** most significant w-r bits. */
41      private static final int UPPER_MASK = 0x80000000;
42  
43      /** least significant w-r bits. */
44      private static final int LOWER_MASK = 0x7fffffff;
45  
46      /** tempering mask b. */
47      private static final int TEMPERING_MASK_B = 0x9d2c5680;
48  
49      /** tempering mask b. */
50      private static final int TEMPERING_MASK_C = 0xefc60000;
51  
52      /** unsigned mask for promoting int -> long. */
53      private static final int UMASK = (1 << 31) - 1;
54  
55      /** the array for the state vector. */
56      private int[] mt;
57  
58      /** The counter mti==N+1 means mt[N] is not initialized. */
59      private int mti;
60  
61      /** magic01. */
62      private int[] mag01;
63  
64      /**
65       * constructs a new Mersenne Twister. <code>System.currentTimeMillis()</code> is used as seed value.
66       */
67      public MersenneTwister()
68      {
69          this(System.currentTimeMillis());
70      }
71  
72      /**
73       * Constructor using a given seed.
74       * @param seed long; The initial seed.
75       */
76      public MersenneTwister(final long seed)
77      {
78          super(seed);
79      }
80  
81      /**
82       * initalizes the MersenneTwister.
83       */
84      private void initialize()
85      {
86          this.mt = new int[N];
87  
88          this.mt[0] = ((int) super.seed) & UMASK;
89          if (this.mt[0] == 0)
90          {
91              // super.seed=Integer.MAXValue --> seed & UMASK==0
92              // We set the seed again and enforce a different value.
93              this.setSeed(System.currentTimeMillis());
94          }
95          for (this.mti = 1; this.mti < N; this.mti++)
96          {
97              this.mt[this.mti] = (69069 * this.mt[this.mti - 1]);
98          }
99  
100         // mag01[x] = x * MATRIX_A for x=0,1
101         this.mag01 = new int[2];
102         this.mag01[0] = 0x0;
103         this.mag01[1] = MATRIX_A;
104     }
105 
106     /** {@inheritDoc} */
107     @Override
108     public synchronized long next(final int bits)
109     {
110         if (bits < 0 || bits > 64)
111         {
112             throw new IllegalArgumentException("bits (" + bits + ") not in range [0,64]");
113         }
114         int y;
115         if (this.mti >= N) // generate N words at one time
116         {
117             int kk;
118             for (kk = 0; kk < N - M; kk++)
119             {
120                 y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK);
121                 this.mt[kk] = this.mt[kk + M] ^ (y >>> 1) ^ this.mag01[y & 0x1];
122             }
123             for (; kk < N - 1; kk++)
124             {
125                 y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK);
126                 this.mt[kk] = this.mt[kk + (M - N)] ^ (y >>> 1) ^ this.mag01[y & 0x1];
127             }
128             y = (this.mt[N - 1] & UPPER_MASK) | (this.mt[0] & LOWER_MASK);
129             this.mt[N - 1] = this.mt[M - 1] ^ (y >>> 1) ^ this.mag01[y & 0x1];
130             this.mti = 0;
131         }
132         y = this.mt[this.mti++];
133         y ^= y >>> 11; // TEMPERING_SHIFT_U(y)
134         y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y)
135         y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y)
136         y ^= (y >>> 18); // TEMPERING_SHIFT_L(y)
137         if (bits <= 32)
138         {
139             return y >>> (32 - bits);
140         }
141         return y << 32 + this.next(bits - 32);
142     }
143 
144     /** {@inheritDoc} */
145     @Override
146     public synchronized void setSeed(final long seed)
147     {
148         super.seed = seed;
149         this.initialize();
150     }
151 }