Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Persistence Test/HeuristicLab.Random/3.3/MersenneTwister.cs @ 4062

Last change on this file since 4062 was 1823, checked in by epitzer, 15 years ago

Namespace refactoring: rename formatters & decomposers -> primitive and composite serializers. (#603)

File size: 8.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2008 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 System.Collections.Generic;
36using System.Text;
37using System.Xml;
38using HeuristicLab.Core;
39using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
40
41namespace HeuristicLab.Random {
42  /// <summary>
43  /// A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator.
44  /// </summary>
45  public class MersenneTwister : ItemBase, 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    /// Initializes a new instance of <see cref="MersenneTwister"/>.
58    /// </summary>
59    public MersenneTwister() {
60      if (!init) seed((uint)DateTime.Now.Ticks);
61      init = true;
62    }
63    /// <summary>
64    /// Initializes a new instance of <see cref="MersenneTwister"/>
65    /// with the given seed <paramref name="s"/>.
66    /// </summary>
67    /// <param name="s">The seed with which to initialize the random number generator.</param>
68    public MersenneTwister(uint s) {
69      seed(s);
70      init = true;
71    }
72    /// <summary>
73    /// Initializes a new instance of <see cref="MersenneTwister"/> with the given seed array.
74    /// </summary>
75    /// <param name="array">The seed array with which to initialize the random number generator.</param>
76    public MersenneTwister(uint[] array) {
77      seed(array);
78      init = true;
79    }
80
81    /// <summary>
82    /// Clones the current instance (deep clone).
83    /// </summary>
84    /// <param name="clonedObjects">Dictionary of all already cloned objects. (Needed to avoid cycles.)</param>
85    /// <returns>The cloned object as <see cref="MersenneTwister"/>.</returns>
86    public override object Clone(IDictionary<Guid, object> clonedObjects) {
87      MersenneTwister clone = new MersenneTwister();
88      clonedObjects.Add(Guid, clone);
89      clone.state = (uint[])state.Clone();
90      clone.p = p;
91      clone.init = init;
92      return clone;
93    }
94
95    /// <summary>
96    /// Resets the current random number generator.
97    /// </summary>
98    public void Reset() {
99      lock (locker)
100        seed((uint)DateTime.Now.Ticks);
101    }
102    /// <summary>
103    /// Resets the current random number generator with the given seed <paramref name="s"/>.
104    /// </summary>
105    /// <param name="s">The seed with which to reset the current instance.</param>
106    public void Reset(int s) {
107      lock (locker)
108        seed((uint)s);
109    }
110
111    /// <summary>
112    /// Gets a new random number.
113    /// </summary>
114    /// <returns>A new int random number.</returns>
115    public int Next() {
116      lock (locker) {
117        return (int)(rand_int32() >> 1);
118      }
119    }
120    /// <summary>
121    /// Gets a new random number being smaller than the given <paramref name="maxVal"/>.
122    /// </summary>
123    /// <exception cref="ArgumentException">Thrown when the given maximum value is
124    /// smaller or equal to zero.</exception>
125    /// <param name="maxVal">The maximum value of the generated random number.</param>
126    /// <returns>A new int random number.</returns>
127    public int Next(int maxVal) {
128      lock (locker) {
129        if (maxVal <= 0)
130          throw new ArgumentException("The interval [0, " + maxVal + ") is empty");
131        int limit = (Int32.MaxValue / maxVal) * maxVal;
132        int value = Next();
133        while (value >= limit) value = Next();
134        return value % maxVal;
135      }
136    }
137    /// <summary>
138    /// Gets a new random number being in the given interval <paramref name="minVal"/> and
139    /// <paramref name="maxVal"/>.
140    /// </summary>
141    /// <param name="minVal">The minimum value of the generated random number.</param>
142    /// <param name="maxVal">The maximum value of the generated random number.</param>
143    /// <returns>A new int random number.</returns>
144    public int Next(int minVal, int maxVal) {
145      lock (locker) {
146        if (maxVal <= minVal)
147          throw new ArgumentException("The interval [" + minVal + ", " + maxVal + ") is empty");
148        return Next(maxVal - minVal) + minVal;
149      }
150    }
151    /// <summary>
152    /// Gets a new double random variable.
153    /// </summary>
154    /// <returns></returns>
155    public double NextDouble() {
156      lock (locker) {
157        return rand_double53();
158      }
159    }
160
161    #region Seed Methods
162    /// <summary>
163    /// Initializes current instance with random seed.
164    /// </summary>
165    /// <param name="s">A starting seed.</param>
166    public void seed(uint s) {
167      state[0] = s & 0xFFFFFFFFU;
168      for (int i = 1; i < n; ++i) {
169        state[i] = 1812433253U * (state[i - 1] ^ (state[i - 1] >> 30)) + (uint)i;
170        state[i] &= 0xFFFFFFFFU;
171      }
172      p = n;
173    }
174    /// <summary>
175    /// Initializes current instance with random seed.
176    /// </summary>
177    /// <param name="array">A starting seed array.</param>
178    public void seed(uint[] array) {
179      seed(19650218U);
180      int i = 1, j = 0;
181      for (int k = ((n > array.Length) ? n : array.Length); k > 0; --k) {
182        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525U))
183          + array[j] + (uint)j;
184        state[i] &= 0xFFFFFFFFU;
185        ++j;
186        j %= array.Length;
187        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
188      }
189      for (int k = n - 1; k > 0; --k) {
190        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941U)) - (uint)i;
191        state[i] &= 0xFFFFFFFFU;
192        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
193      }
194      state[0] = 0x80000000U;
195      p = n;
196    }
197    #endregion
198
199    #region Random Number Generation Methods
200    private uint rand_int32() {
201      if (p == n) gen_state();
202      uint x = state[p++];
203      x ^= (x >> 11);
204      x ^= (x << 7) & 0x9D2C5680U;
205      x ^= (x << 15) & 0xEFC60000U;
206      return x ^ (x >> 18);
207    }
208    private double rand_double() { // interval [0, 1)
209      return ((double)rand_int32()) * (1.0 / 4294967296.0);
210    }
211    private double rand_double_closed() { // interval [0, 1]
212      return ((double)rand_int32()) * (1.0 / 4294967295.0);
213    }
214    private double rand_double_open() { // interval (0, 1)
215      return (((double)rand_int32()) + 0.5) * (1.0 / 4294967296.0);
216    }
217    private double rand_double53() { // 53 bit resolution, interval [0, 1)
218      return (((double)(rand_int32() >> 5)) * 67108864.0 +
219        ((double)(rand_int32() >> 6))) * (1.0 / 9007199254740992.0);
220    }
221    #endregion
222
223    #region Private Helper Methods
224    private uint twiddle(uint u, uint v) {
225      return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
226        ^ (((v & 1U) != 0) ? 0x9908B0DFU : 0x0U);
227    }
228    private void gen_state() {
229      for (int i = 0; i < (n - m); ++i)
230        state[i] = state[i + m] ^ twiddle(state[i], state[i + 1]);
231      for (int i = n - m; i < (n - 1); ++i)
232        state[i] = state[i + m - n] ^ twiddle(state[i], state[i + 1]);
233      state[n - 1] = state[m - 1] ^ twiddle(state[n - 1], state[0]);
234      p = 0; // reset position
235    }
236    #endregion
237  }
238}
Note: See TracBrowser for help on using the repository browser.