Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Random/MersenneTwister.cs @ 143

Last change on this file since 143 was 2, checked in by swagner, 17 years ago

Added HeuristicLab 3.0 sources from former SVN repository at revision 52

File size: 7.0 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;
39
40namespace HeuristicLab.Random {
41  public class MersenneTwister : ItemBase, IRandom {
42    private const int n = 624, m = 397;
43
44    private object locker = new object();
45    private uint[] state = new uint[n];
46    private int p = 0;
47    private bool init = false;
48
49    public MersenneTwister() {
50      if (!init) seed((uint)DateTime.Now.Ticks);
51      init = true;
52    }
53    public MersenneTwister(uint s) {
54      seed(s);
55      init = true;
56    }
57    public MersenneTwister(uint[] array) {
58      seed(array);
59      init = true;
60    }
61
62    public override object Clone(IDictionary<Guid, object> clonedObjects) {
63      MersenneTwister clone = new MersenneTwister();
64      clonedObjects.Add(Guid, clone);
65      clone.state = (uint[])state.Clone();
66      clone.p = p;
67      clone.init = init;
68      return clone;
69    }
70
71    public void Reset() {
72      lock (locker)
73        seed((uint)DateTime.Now.Ticks);
74    }
75    public void Reset(int s) {
76      lock (locker)
77        seed((uint)s);
78    }
79
80    public int Next() {
81      lock (locker) {
82        return (int)(rand_int32() >> 1);
83      }
84    }
85    public int Next(int maxVal) {
86      lock (locker) {
87        if (maxVal <= 0)
88          throw new ArgumentException("The interval [0, " + maxVal + ") is empty");
89        int limit = (Int32.MaxValue / maxVal) * maxVal;
90        int value = Next();
91        while (value >= limit) value = Next();
92        return value % maxVal;
93      }
94    }
95    public int Next(int minVal, int maxVal) {
96      lock (locker) {
97        if (maxVal <= minVal)
98          throw new ArgumentException("The interval [" + minVal + ", " + maxVal + ") is empty");
99        return Next(maxVal - minVal) + minVal;
100      }
101    }
102    public double NextDouble() {
103      lock (locker) {
104        return rand_double53();
105      }
106    }
107
108    #region Persistence Methods
109    public override XmlNode GetXmlNode(string name, XmlDocument document, IDictionary<Guid,IStorable> persistedObjects) {
110      XmlNode node = base.GetXmlNode(name, document, persistedObjects);
111
112      StringBuilder builder = new StringBuilder();
113      builder.Append(state[0]);
114      for (int i = 1; i < state.Length; i++) {
115        builder.Append(';');
116        builder.Append(state[i]);
117      }
118      XmlNode stateNode = document.CreateNode(XmlNodeType.Element, "State", null);
119      stateNode.InnerText = builder.ToString();
120      node.AppendChild(stateNode);
121
122      XmlNode pNode = document.CreateNode(XmlNodeType.Element, "P", null);
123      pNode.InnerText = p.ToString();
124      node.AppendChild(pNode);
125
126      XmlNode initNode = document.CreateNode(XmlNodeType.Element, "Init", null);
127      initNode.InnerText = init.ToString();
128      node.AppendChild(initNode);
129
130      return node;
131    }
132    public override void Populate(XmlNode node, IDictionary<Guid,IStorable> restoredObjects) {
133      base.Populate(node, restoredObjects);
134
135      string stateString = node.SelectSingleNode("State").InnerText;
136      string[] tokens = stateString.Split(';');
137      for (int i = 0; i < tokens.Length; i++)
138        state[i] = uint.Parse(tokens[i]);
139      p = int.Parse(node.SelectSingleNode("P").InnerText);
140      init = bool.Parse(node.SelectSingleNode("Init").InnerText);
141    }
142    #endregion
143
144    #region Seed Methods
145    public void seed(uint s) {
146      state[0] = s & 0xFFFFFFFFU;
147      for (int i = 1; i < n; ++i) {
148        state[i] = 1812433253U * (state[i - 1] ^ (state[i - 1] >> 30)) + (uint)i;
149        state[i] &= 0xFFFFFFFFU;
150      }
151      p = n;
152    }
153    public void seed(uint[] array) {
154      seed(19650218U);
155      int i = 1, j = 0;
156      for (int k = ((n > array.Length) ? n : array.Length); k > 0; --k) {
157        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525U))
158          + array[j] + (uint)j;
159        state[i] &= 0xFFFFFFFFU;
160        ++j;
161        j %= array.Length;
162        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
163      }
164      for (int k = n - 1; k > 0; --k) {
165        state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941U)) - (uint)i;
166        state[i] &= 0xFFFFFFFFU;
167        if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
168      }
169      state[0] = 0x80000000U;
170      p = n;
171    }
172    #endregion
173
174    #region Random Number Generation Methods
175    private uint rand_int32() {
176      if (p == n) gen_state();
177      uint x = state[p++];
178      x ^= (x >> 11);
179      x ^= (x << 7) & 0x9D2C5680U;
180      x ^= (x << 15) & 0xEFC60000U;
181      return x ^ (x >> 18);
182    }
183    private double rand_double() { // interval [0, 1)
184      return ((double)rand_int32()) * (1.0 / 4294967296.0);
185    }
186    private double rand_double_closed() { // interval [0, 1]
187      return ((double)rand_int32()) * (1.0 / 4294967295.0);
188    }
189    private double rand_double_open() { // interval (0, 1)
190      return (((double)rand_int32()) + 0.5) * (1.0 / 4294967296.0);
191    }
192    private double rand_double53() { // 53 bit resolution, interval [0, 1)
193      return (((double)(rand_int32() >> 5)) * 67108864.0 +
194        ((double)(rand_int32() >> 6))) * (1.0 / 9007199254740992.0);
195    }
196    #endregion
197
198    #region Private Helper Methods
199    private uint twiddle(uint u, uint v) {
200      return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
201        ^ (((v & 1U) != 0) ? 0x9908B0DFU : 0x0U);
202    }
203    private void gen_state() {
204      for (int i = 0; i < (n - m); ++i)
205        state[i] = state[i + m] ^ twiddle(state[i], state[i + 1]);
206      for (int i = n - m; i < (n - 1); ++i)
207        state[i] = state[i + m - n] ^ twiddle(state[i], state[i + 1]);
208      state[n - 1] = state[m - 1] ^ twiddle(state[n - 1], state[0]);
209      p = 0; // reset position
210    }
211    #endregion
212  }
213}
Note: See TracBrowser for help on using the repository browser.