Mercurial > repos > rhope
comparison mt19937ar.c @ 0:76568becd6d6
Rhope Alpha 2a source import
author | Mike Pavone <pavone@retrodev.com> |
---|---|
date | Tue, 28 Apr 2009 23:06:07 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:76568becd6d6 |
---|---|
1 /* | |
2 A C-program for MT19937, with initialization improved 2002/1/26. | |
3 Coded by Takuji Nishimura and Makoto Matsumoto. | |
4 | |
5 Before using, initialize the state by using init_genrand(seed) | |
6 or init_by_array(init_key, key_length). | |
7 | |
8 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | |
9 All rights reserved. | |
10 Copyright (C) 2005, Mutsuo Saito, | |
11 All rights reserved. | |
12 | |
13 Redistribution and use in source and binary forms, with or without | |
14 modification, are permitted provided that the following conditions | |
15 are met: | |
16 | |
17 1. Redistributions of source code must retain the above copyright | |
18 notice, this list of conditions and the following disclaimer. | |
19 | |
20 2. Redistributions in binary form must reproduce the above copyright | |
21 notice, this list of conditions and the following disclaimer in the | |
22 documentation and/or other materials provided with the distribution. | |
23 | |
24 3. The names of its contributors may not be used to endorse or promote | |
25 products derived from this software without specific prior written | |
26 permission. | |
27 | |
28 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | |
29 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | |
30 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | |
31 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR | |
32 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, | |
33 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | |
34 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | |
35 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | |
36 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | |
37 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | |
38 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
39 | |
40 | |
41 Any feedback is very welcome. | |
42 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html | |
43 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) | |
44 */ | |
45 | |
46 #include <stdio.h> | |
47 #include "mt19937ar.h" | |
48 | |
49 /* Period parameters */ | |
50 #define N 624 | |
51 #define M 397 | |
52 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
53 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */ | |
54 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */ | |
55 | |
56 static unsigned long mt[N]; /* the array for the state vector */ | |
57 static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ | |
58 | |
59 /* initializes mt[N] with a seed */ | |
60 void init_genrand(unsigned long s) | |
61 { | |
62 mt[0]= s & 0xffffffffUL; | |
63 for (mti=1; mti<N; mti++) { | |
64 mt[mti] = | |
65 (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); | |
66 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ | |
67 /* In the previous versions, MSBs of the seed affect */ | |
68 /* only MSBs of the array mt[]. */ | |
69 /* 2002/01/09 modified by Makoto Matsumoto */ | |
70 mt[mti] &= 0xffffffffUL; | |
71 /* for >32 bit machines */ | |
72 } | |
73 } | |
74 | |
75 /* initialize by an array with array-length */ | |
76 /* init_key is the array for initializing keys */ | |
77 /* key_length is its length */ | |
78 /* slight change for C++, 2004/2/26 */ | |
79 void init_by_array(unsigned long init_key[], int key_length) | |
80 { | |
81 int i, j, k; | |
82 init_genrand(19650218UL); | |
83 i=1; j=0; | |
84 k = (N>key_length ? N : key_length); | |
85 for (; k; k--) { | |
86 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) | |
87 + init_key[j] + j; /* non linear */ | |
88 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ | |
89 i++; j++; | |
90 if (i>=N) { mt[0] = mt[N-1]; i=1; } | |
91 if (j>=key_length) j=0; | |
92 } | |
93 for (k=N-1; k; k--) { | |
94 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) | |
95 - i; /* non linear */ | |
96 mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ | |
97 i++; | |
98 if (i>=N) { mt[0] = mt[N-1]; i=1; } | |
99 } | |
100 | |
101 mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ | |
102 } | |
103 | |
104 /* generates a random number on [0,0xffffffff]-interval */ | |
105 unsigned long genrand_int32(void) | |
106 { | |
107 unsigned long y; | |
108 static unsigned long mag01[2]={0x0UL, MATRIX_A}; | |
109 /* mag01[x] = x * MATRIX_A for x=0,1 */ | |
110 | |
111 if (mti >= N) { /* generate N words at one time */ | |
112 int kk; | |
113 | |
114 if (mti == N+1) /* if init_genrand() has not been called, */ | |
115 init_genrand(5489UL); /* a default initial seed is used */ | |
116 | |
117 for (kk=0;kk<N-M;kk++) { | |
118 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); | |
119 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
120 } | |
121 for (;kk<N-1;kk++) { | |
122 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); | |
123 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
124 } | |
125 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); | |
126 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; | |
127 | |
128 mti = 0; | |
129 } | |
130 | |
131 y = mt[mti++]; | |
132 | |
133 /* Tempering */ | |
134 y ^= (y >> 11); | |
135 y ^= (y << 7) & 0x9d2c5680UL; | |
136 y ^= (y << 15) & 0xefc60000UL; | |
137 y ^= (y >> 18); | |
138 | |
139 return y; | |
140 } | |
141 | |
142 /* generates a random number on [0,0x7fffffff]-interval */ | |
143 long genrand_int31(void) | |
144 { | |
145 return (long)(genrand_int32()>>1); | |
146 } | |
147 | |
148 /* generates a random number on [0,1]-real-interval */ | |
149 double genrand_real1(void) | |
150 { | |
151 return genrand_int32()*(1.0/4294967295.0); | |
152 /* divided by 2^32-1 */ | |
153 } | |
154 | |
155 /* generates a random number on [0,1)-real-interval */ | |
156 double genrand_real2(void) | |
157 { | |
158 return genrand_int32()*(1.0/4294967296.0); | |
159 /* divided by 2^32 */ | |
160 } | |
161 | |
162 /* generates a random number on (0,1)-real-interval */ | |
163 double genrand_real3(void) | |
164 { | |
165 return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0); | |
166 /* divided by 2^32 */ | |
167 } | |
168 | |
169 /* generates a random number on [0,1) with 53-bit resolution*/ | |
170 double genrand_res53(void) | |
171 { | |
172 unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; | |
173 return(a*67108864.0+b)*(1.0/9007199254740992.0); | |
174 } | |
175 /* These real versions are due to Isaku Wada, 2002/01/09 added */ |