COPASI API  4.16.103
Cmt19937.cpp
Go to the documentation of this file.
1 /* Begin CVS Header
2  $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/randomGenerator/Cmt19937.cpp,v $
3  $Revision: 1.8 $
4  $Name: $
5  $Author: shoops $
6  $Date: 2006/07/21 18:15:48 $
7  End CVS Header */
8 
9 // Copyright 2005 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc. and EML Research, gGmbH.
11 // All rights reserved.
12 
13 /**
14  * Cmt19937 class implementing the Mersenne Twister random number generator,
15  *
16  * Created for Copasi by Stefan Hoops
17  *
18  * A C-program for MT19937, with initialization improved 2002/2/10.
19  * Coded by Takuji Nishimura and Makoto Matsumoto.
20  * This is a faster version by taking Shawn Cokus's optimization,
21  * Matthe Bellew's simplification, Isaku Wada's real version.
22  *
23  * Before using, initialize the state by using init_genrand(seed)
24  * or init_by_array(init_key, key_length).
25  *
26  * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
27  * All rights reserved.
28  *
29  * Redistribution and use in source and binary forms, with or without
30  * modification, are permitted provided that the following conditions
31  * are met:
32  *
33  * 1. Redistributions of source code must retain the above copyright
34  * notice, this list of conditions and the following disclaimer.
35  *
36  * 2. Redistributions in binary form must reproduce the above copyright
37  * notice, this list of conditions and the following disclaimer in the
38  * documentation and/or other materials provided with the distribution.
39  *
40  * 3. The names of its contributors may not be used to endorse or promote
41  * products derived from this software without specific prior written
42  * permission.
43  *
44  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
45  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
46  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
47  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
48  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
49  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
50  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
51  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
52  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
53  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
54  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
55  *
56  *
57  * Any feedback is very welcome.
58  * http://www.math.keio.ac.jp/matumoto/emt.html
59  * email: matumoto@math.keio.ac.jp
60  */
61 
62 #include "copasi.h"
63 #include "CRandom.h"
64 
65 /* Period parameters */
66 #define Cmt19937_M 397
67 #define Cmt19937_MATRIX_A 0x9908b0dfUL /* constant vector a */
68 #define Cmt19937_UMASK 0x80000000UL /* most significant w-r bits */
69 #define Cmt19937_LMASK 0x7fffffffUL /* least significant r bits */
70 #define Cmt19937_MIXBITS(u,v) \
71  (((u) & Cmt19937_UMASK) | ((v) & Cmt19937_LMASK))
72 #define Cmt19937_TWIST(u,v) \
73  ((Cmt19937_MIXBITS(u,v) >> 1) ^ ((v)&1UL ? Cmt19937_MATRIX_A : 0UL))
74 
75 Cmt19937::Cmt19937(unsigned C_INT32 seed):
76  CRandom(),
77  mLeft(1),
78  mNext(NULL)
79 {
80  setModulus(0xffffffffUL);
81  initialize(seed);
82 }
83 
85 
86 /* initializes mState[Cmt19937_N] with a seed */
87 void Cmt19937::initialize(unsigned C_INT32 seed)
88 {
89  int j;
90  mState[0] = seed & 0xffffffffUL;
91  for (j = 1; j < Cmt19937_N; j++)
92  {
93  mState[j] = (1812433253UL * (mState[j - 1] ^ (mState[j - 1] >> 30)) + j);
94  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
95  /* In the previous versions, MSBs of the seed affect */
96  /* only MSBs of the array mState[]. */
97  /* 2002/01/09 modified by Makoto Matsumoto */
98  mState[j] &= 0xffffffffUL; /* for >32 bit machines */
99  }
100  mLeft = 1;
101 }
102 
103 /* initialize by an array with array-length */
104 /* init_key is the array for initializing keys */
105 /* key_length is its length */
106 void Cmt19937::init_by_array(unsigned C_INT32 init_key[],
107  C_INT32 key_length)
108 {
109  int i, j, k;
110  initialize(19650218UL);
111  i = 1;
112  j = 0;
113  k = (Cmt19937_N > key_length ? Cmt19937_N : key_length);
114  for (; k; k--)
115  {
116  mState[i] = (mState[i] ^ ((mState[i - 1] ^ (mState[i - 1] >> 30)) * 1664525UL))
117  + init_key[j] + j; /* non linear */
118  mState[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
119  i++;
120  j++;
121  if (i >= Cmt19937_N)
122  {mState[0] = mState[Cmt19937_N - 1]; i = 1;}
123  if (j >= key_length)
124  j = 0;
125  }
126  for (k = Cmt19937_N - 1; k; k--)
127  {
128  mState[i] = (mState[i] ^ ((mState[i - 1] ^ (mState[i - 1] >> 30)) * 1566083941UL))
129  - i; /* non linear */
130  mState[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
131  i++;
132  if (i >= Cmt19937_N)
133  {mState[0] = mState[Cmt19937_N - 1]; i = 1;}
134  }
135 
136  mState[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
137  mLeft = 1;
138 }
139 
141 {
142  unsigned C_INT32 *p = mState;
143  int j;
144 
145  mLeft = Cmt19937_N;
146  mNext = mState;
147 
148  for (j = Cmt19937_N - Cmt19937_M + 1; --j; p++)
149  *p = p[Cmt19937_M] ^ Cmt19937_TWIST(p[0], p[1]);
150 
151  for (j = Cmt19937_M; --j; p++)
152  *p = p[Cmt19937_M - Cmt19937_N] ^ Cmt19937_TWIST(p[0], p[1]);
153 
154  *p = p[Cmt19937_M - Cmt19937_N] ^ Cmt19937_TWIST(p[0], mState[0]);
155 }
156 
157 /* generates a random number on [0,0xffffffff]-interval */
159 {
160  if (--mLeft == 0)
161  next_state();
162  mNumberU = *mNext++;
163 
164  /* Tempering */
165  mNumberU ^= (mNumberU >> 11);
166  mNumberU ^= (mNumberU << 7) & 0x9d2c5680UL;
167  mNumberU ^= (mNumberU << 15) & 0xefc60000UL;
168  mNumberU ^= (mNumberU >> 18);
169 
170  return mNumberU;
171 }
172 
173 /* generates a random number on [0,0x7fffffff]-interval */
175 {
176  if (--mLeft == 0)
177  next_state();
178  mNumberU = *mNext++;
179 
180  /* Tempering */
181  mNumberU ^= (mNumberU >> 11);
182  mNumberU ^= (mNumberU << 7) & 0x9d2c5680UL;
183  mNumberU ^= (mNumberU << 15) & 0xefc60000UL;
184  mNumberU ^= (mNumberU >> 18);
185 
186  return mNumberU >> 1;
187 }
188 
189 /* generates a random number on [0,1]-real-interval */
191 {
192  if (--mLeft == 0)
193  next_state();
194  mNumberU = *mNext++;
195 
196  /* Tempering */
197  mNumberU ^= (mNumberU >> 11);
198  mNumberU ^= (mNumberU << 7) & 0x9d2c5680UL;
199  mNumberU ^= (mNumberU << 15) & 0xefc60000UL;
200  mNumberU ^= (mNumberU >> 18);
201 
202  return ((C_FLOAT64) mNumberU) * (1.0 / 4294967295.0);
203  /* divided by 2^32-1 */
204 }
205 
206 /* generates a random number on [0,1)-real-interval */
208 {
209  if (--mLeft == 0)
210  next_state();
211  mNumberU = *mNext++;
212 
213  /* Tempering */
214  mNumberU ^= (mNumberU >> 11);
215  mNumberU ^= (mNumberU << 7) & 0x9d2c5680UL;
216  mNumberU ^= (mNumberU << 15) & 0xefc60000UL;
217  mNumberU ^= (mNumberU >> 18);
218 
219  return ((C_FLOAT64) mNumberU) * (1.0 / 4294967296.0);
220  /* divided by 2^32 */
221 }
222 
223 /* generates a random number on (0,1)-real-interval */
225 {
226  if (--mLeft == 0)
227  next_state();
228  mNumberU = *mNext++;
229 
230  /* Tempering */
231  mNumberU ^= (mNumberU >> 11);
232  mNumberU ^= (mNumberU << 7) & 0x9d2c5680UL;
233  mNumberU ^= (mNumberU << 15) & 0xefc60000UL;
234  mNumberU ^= (mNumberU >> 18);
235 
236  return (((C_FLOAT64) mNumberU) + 0.5) * (1.0 / 4294967296.0);
237  /* divided by 2^32 */
238 }
239 
240 /* generates a random number on [0,1) with 53-bit resolution*/
242 {
243  unsigned C_INT32 a = getRandomU() >> 5, b = getRandomU() >> 6;
244  return mFloat = (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
245 }
246 /* These real versions are due to Isaku Wada, 2002/01/09 added */
247 
249  Cmt19937(seed)
250 {}
251 
252 /* generates a random number on [0,1]-real-interval */
254 {
255  unsigned C_INT32 a = getRandomU() >> 5, b = getRandomU() >> 6;
256  return mFloat = (a * 67108864.0 + b) * (1.0 / 9007199254740991.0);
257 }
258 
259 /* generates a random number on [0,1)-real-interval */
261 {
262  unsigned C_INT32 a = getRandomU() >> 5, b = getRandomU() >> 6;
263  return mFloat = (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
264 }
265 
266 /* generates a random number on (0,1)-real-interval */
268 {
269  unsigned C_INT32 a = getRandomU() >> 5, b = getRandomU() >> 6;
270  return mFloat = (a * 67108864.0 + b + 0.5) * (1.0 / 9007199254740992.0);
271 }
272 
273 /* These real versions are due to Isaku Wada, 2002/01/09 added */
274 
275 #ifdef XXXX
276 int main()
277 {
278  int i;
279  unsigned C_INT32 init[4] = {0x123, 0x234, 0x345, 0x456}, length = 4;
280  init_by_array(init, length);
281  /* This is an example of initializing by an array. */
282  /* You may use init_genrand(seed) with any 32bit integer */
283  /* as a seed for a simpler initialization */
284  printf("1000 outputs of getRandomU()\n");
285  for (i = 0; i < 1000; i++)
286  {
287  printf("%10lu ", getRandomU());
288  if (i % 5 == 4)
289  printf("\n");
290  }
291  printf("\n1000 outputs of genrand_real2()\n");
292  for (i = 0; i < 1000; i++)
293  {
294  printf("%10.8f ", getRandomCO());
295  if (i % 5 == 4)
296  printf("\n");
297  }
298 
299  return 0;
300 }
301 #endif // XXXX
Cmt19937HR(unsigned C_INT32 seed)
Definition: Cmt19937.cpp:248
C_FLOAT64 getRandomOO()
Definition: Cmt19937.cpp:267
Cmt19937(unsigned C_INT32 seed)
Definition: Cmt19937.cpp:75
C_FLOAT64 mFloat
Definition: CRandom.h:88
int main(int argc, char **argv)
C_FLOAT64 getRandomCO()
Definition: Cmt19937.cpp:260
void initialize(unsigned C_INT32 seed=CRandom::getSystemSeed())
Definition: Cmt19937.cpp:87
unsigned C_INT32 mState[624]
Definition: Cmt19937.h:74
#define C_INT32
Definition: copasi.h:90
~Cmt19937()
Definition: Cmt19937.cpp:84
void init_by_array(unsigned C_INT32 init_key[], C_INT32 key_length)
Definition: Cmt19937.cpp:106
C_FLOAT64 getRandomCO()
Definition: Cmt19937.cpp:207
C_FLOAT64 getRandomCC()
Definition: Cmt19937.cpp:253
C_FLOAT64 getRandomOO()
Definition: Cmt19937.cpp:224
#define Cmt19937_TWIST(u, v)
Definition: Cmt19937.cpp:72
unsigned C_INT32 getRandomU()
Definition: Cmt19937.cpp:158
C_FLOAT64 getRandomCC()
Definition: Cmt19937.cpp:190
void next_state()
Definition: Cmt19937.cpp:140
void setModulus(const unsigned C_INT32 &modulus)
Definition: CRandom.cpp:129
#define C_FLOAT64
Definition: copasi.h:92
unsigned C_INT32 * mNext
Definition: Cmt19937.h:78
#define Cmt19937_N
Definition: Cmt19937.h:65
C_FLOAT64 genrand_res53()
Definition: Cmt19937.cpp:241
C_INT32 getRandomS()
Definition: Cmt19937.cpp:174
C_INT mLeft
Definition: Cmt19937.h:76
unsigned C_INT32 mNumberU
Definition: CRandom.h:78
#define Cmt19937_M
Definition: Cmt19937.cpp:66