Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Random/3.3/MersenneTwister.cs @ 12087

Last change on this file since 12087 was 12012, checked in by ascheibe, 10 years ago

#2212 merged r12008, r12009, r12010 back into trunk

File size: 9.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22/*
23 * C# port of the freeware implementation of the Mersenne Twister
24 * originally developed by M. Matsumoto and T. Nishimura
25 *
26 * M. Matsumoto and T. Nishimura,
27 * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
28 * Pseudo-Random Number Generator",
29 * ACM Transactions on Modeling and Computer Simulation,
30 * Vol. 8, No. 1, January 1998, pp 3-30.
31 *
32 */
33
34using System;
35using HeuristicLab.Common;
36using HeuristicLab.Core;
37using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
38
39namespace HeuristicLab.Random {
40  /// <summary>
41  /// A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator.
42  /// </summary>
43  [Item("MersenneTwister", "A high-quality pseudo random number generator which creates uniformly distributed random numbers.")]
44  [StorableClass]
45  public sealed class MersenneTwister : Item, IRandom {
46    private const int n = 624, m = 397;
47
48    private object locker = new object();
49    [Storable]
50    private uint[] state = new uint[n];
51    [Storable]
52    private int p = 0;
53    [Storable]
54    private bool init = false;
55
56    /// <summary>
57    /// Used by HeuristicLab.Persistence to initialize new instances during deserialization.
58    /// </summary>
59    /// <param name="deserializing">true, if the constructor is called during deserialization.</param>
60    [StorableConstructor]
61    private MersenneTwister(bool deserializing) : base(deserializing) { }
62    /// <summary>
63    /// Initializes a new instance from an existing one (copy constructor).
64    /// </summary>
65    /// <param name="original">The original <see cref="MersenneTwister"/> instance which is used to initialize the new instance.</param>
66    /// <param name="cloner">A <see cref="Cloner"/> which is used to track all already cloned objects in order to avoid cycles.</param>
67    private MersenneTwister(MersenneTwister original, Cloner cloner)
68      : base(original, cloner) {
69      state = (uint[])original.state.Clone();
70      p = original.p;
71      init = original.init;
72    }
73    /// <summary>
74    /// Initializes a new instance of <see cref="MersenneTwister"/>.
75    /// </summary>
76    public MersenneTwister() {
77      if (!init) Seed((uint)DateTime.Now.Ticks);
78      init = true;
79    }
80    /// <summary>
81    /// Initializes a new instance of <see cref="MersenneTwister"/>
82    /// with the given seed <paramref name="s"/>.
83    /// </summary>
84    /// <param name="s">The seed with which to initialize the random number generator.</param>
85    public MersenneTwister(uint s) {
86      Seed(s);
87      init = true;
88    }
89    /// <summary>
90    /// Initializes a new instance of <see cref="MersenneTwister"/> with the given seed array.
91    /// </summary>
92    /// <param name="array">The seed array with which to initialize the random number generator.</param>
93    public MersenneTwister(uint[] array) {
94      Seed(array);
95      init = true;
96    }
97
98    /// <summary>
99    /// Clones the current instance (deep clone).
100    /// </summary>
101    /// <param name="clonedObjects">Dictionary of all already cloned objects. (Needed to avoid cycles.)</param>
102    /// <returns>The cloned object as <see cref="MersenneTwister"/>.</returns>
103    public override IDeepCloneable Clone(Cloner cloner) {
104      return new MersenneTwister(this, cloner);
105    }
106
107    /// <summary>
108    /// Resets the current random number generator.
109    /// </summary>
110    public void Reset() {
111      lock (locker)
112        Seed((uint)DateTime.Now.Ticks);
113    }
114    /// <summary>
115    /// Resets the current random number generator with the given seed <paramref name="s"/>.
116    /// </summary>
117    /// <param name="s">The seed with which to reset the current instance.</param>
118    public void Reset(int s) {
119      lock (locker)
120        Seed((uint)s);
121    }
122
123    /// <summary>
124    /// Gets a new random number.
125    /// </summary>
126    /// <returns>A new int random number.</returns>
127    public int Next() {
128      lock (locker) {
129        return (int)(rand_int32() >> 1);
130      }
131    }
132    /// <summary>
133    /// Gets a new random number being smaller than the given <paramref name="maxVal"/>.
134    /// </summary>
135    /// <exception cref="ArgumentException">Thrown when the given maximum value is
136    /// smaller or equal to zero.</exception>
137    /// <param name="maxVal">The maximum value of the generated random number.</param>
138    /// <returns>A new int random number.</returns>
139    public int Next(int maxVal) {
140      lock (locker) {
141        if (maxVal <= 0)
142          throw new ArgumentException("The interval [0, " + maxVal + ") is empty");
143        int limit = (Int32.MaxValue / maxVal) * maxVal;
144        int value = Next();
145        while (value >= limit) value = Next();
146        return value % maxVal;
147      }
148    }
149    /// <summary>
150    /// Gets a new random number being in the given interval <paramref name="minVal"/> and
151    /// <paramref name="maxVal"/>.
152    /// </summary>
153    /// <param name="minVal">The minimum value of the generated random number.</param>
154    /// <param name="maxVal">The maximum value of the generated random number.</param>
155    /// <returns>A new int random number.</returns>
156    public int Next(int minVal, int maxVal) {
157      lock (locker) {
158        if (maxVal <= minVal)
159          throw new ArgumentException("The interval [" + minVal + ", " + maxVal + ") is empty");
160        return Next(maxVal - minVal) + minVal;
161      }
162    }
163    /// <summary>
164    /// Gets a new double random variable.
165    /// </summary>
166    /// <returns></returns>
167    public double NextDouble() {
168      lock (locker) {
169        return rand_double53();
170      }
171    }
172
173    #region Seed Methods
174    /// <summary>
175    /// Initializes current instance with random seed.
176    /// </summary>
177    /// <param name="s">A starting seed.</param>
178    public void Seed(uint s) {
179      state[0] = s & 0xFFFFFFFFU;
180      for (int i = 1; i < n; ++i) {
181        state[i] = 1812433253U * (state[i - 1] ^ (state[i - 1] >> 30)) + (uint)i;
182        state[i] &= 0xFFFFFFFFU;
183      }
184      p = n;
185    }
186    /// <summary>
187    /// Initializes current instance with random seed.
188    /// </summary>
189    /// <param name="array">A starting seed array.</param>
190    public void Seed(uint[] array) {
191      Seed(19650218U);
192      int i = 1, j = 0;
193      for (int k = ((n > array.Length) ? n : array.Length); k > 0; --k) {
194        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525U))
195          + array[j] + (uint)j;
196        state[i] &= 0xFFFFFFFFU;
197        ++j;
198        j %= array.Length;
199        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
200      }
201      for (int k = n - 1; k > 0; --k) {
202        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941U)) - (uint)i;
203        state[i] &= 0xFFFFFFFFU;
204        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
205      }
206      state[0] = 0x80000000U;
207      p = n;
208    }
209    #endregion
210
211    #region Random Number Generation Methods
212    private uint rand_int32() {
213      if (p == n) gen_state();
214      uint x = state[p++];
215      x ^= (x >> 11);
216      x ^= (x << 7) & 0x9D2C5680U;
217      x ^= (x << 15) & 0xEFC60000U;
218      return x ^ (x >> 18);
219    }
220    private double rand_double() { // interval [0, 1)
221      return ((double)rand_int32()) * (1.0 / 4294967296.0);
222    }
223    private double rand_double_closed() { // interval [0, 1]
224      return ((double)rand_int32()) * (1.0 / 4294967295.0);
225    }
226    private double rand_double_open() { // interval (0, 1)
227      return (((double)rand_int32()) + 0.5) * (1.0 / 4294967296.0);
228    }
229    private double rand_double53() { // 53 bit resolution, interval [0, 1)
230      return (((double)(rand_int32() >> 5)) * 67108864.0 +
231        ((double)(rand_int32() >> 6))) * (1.0 / 9007199254740992.0);
232    }
233    #endregion
234
235    #region Private Helper Methods
236    private uint twiddle(uint u, uint v) {
237      return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
238        ^ (((v & 1U) != 0) ? 0x9908B0DFU : 0x0U);
239    }
240    private void gen_state() {
241      for (int i = 0; i < (n - m); ++i)
242        state[i] = state[i + m] ^ twiddle(state[i], state[i + 1]);
243      for (int i = n - m; i < (n - 1); ++i)
244        state[i] = state[i + m - n] ^ twiddle(state[i], state[i + 1]);
245      state[n - 1] = state[m - 1] ^ twiddle(state[n - 1], state[0]);
246      p = 0; // reset position
247    }
248    #endregion
249  }
250}
Note: See TracBrowser for help on using the repository browser.