1 | /* |
2 | A C-program for MT19937-64 (2004/9/29 version). |
3 | Coded by Takuji Nishimura and Makoto Matsumoto. |
4 | |
5 | This is a 64-bit version of Mersenne Twister pseudorandom number |
6 | generator. |
7 | |
8 | Before using, initialize the state by using init_genrand64(seed) |
9 | or init_by_array64(init_key, key_length). |
10 | |
11 | Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura, |
12 | All rights reserved. |
13 | |
14 | Redistribution and use in source and binary forms, with or without |
15 | modification, are permitted provided that the following conditions |
16 | are met: |
17 | |
18 | 1. Redistributions of source code must retain the above copyright |
19 | notice, this list of conditions and the following disclaimer. |
20 | |
21 | 2. Redistributions in binary form must reproduce the above copyright |
22 | notice, this list of conditions and the following disclaimer in the |
23 | documentation and/or other materials provided with the distribution. |
24 | |
25 | 3. The names of its contributors may not be used to endorse or promote |
26 | products derived from this software without specific prior written |
27 | permission. |
28 | |
29 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
30 | "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
31 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
32 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
33 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
34 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
35 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
36 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
37 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
38 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
39 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
40 | |
41 | References: |
42 | T. Nishimura, ``Tables of 64-bit Mersenne Twisters'' |
43 | ACM Transactions on Modeling and |
44 | Computer Simulation 10. (2000) 348--357. |
45 | M. Matsumoto and T. Nishimura, |
46 | ``Mersenne Twister: a 623-dimensionally equidistributed |
47 | uniform pseudorandom number generator'' |
48 | ACM Transactions on Modeling and |
49 | Computer Simulation 8. (Jan. 1998) 3--30. |
50 | |
51 | Any feedback is very welcome. |
52 | http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
53 | email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces) |
54 | */ |
55 | |
56 | |
57 | #include "mt19937-64.h" |
58 | #include <stdio.h> |
59 | |
60 | #define NN 312 |
61 | #define MM 156 |
62 | #define MATRIX_A 0xB5026F5AA96619E9ULL |
63 | #define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */ |
64 | #define LM 0x7FFFFFFFULL /* Least significant 31 bits */ |
65 | |
66 | |
67 | /* The array for the state vector */ |
68 | static unsigned long long mt[NN]; |
69 | /* mti==NN+1 means mt[NN] is not initialized */ |
70 | static int mti=NN+1; |
71 | |
72 | /* initializes mt[NN] with a seed */ |
73 | void init_genrand64(unsigned long long seed) |
74 | { |
75 | mt[0] = seed; |
76 | for (mti=1; mti<NN; mti++) |
77 | mt[mti] = (6364136223846793005ULL * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti); |
78 | } |
79 | |
80 | /* initialize by an array with array-length */ |
81 | /* init_key is the array for initializing keys */ |
82 | /* key_length is its length */ |
83 | void init_by_array64(unsigned long long init_key[], |
84 | unsigned long long key_length) |
85 | { |
86 | unsigned long long i, j, k; |
87 | init_genrand64(19650218ULL); |
88 | i=1; j=0; |
89 | k = (NN>key_length ? NN : key_length); |
90 | for (; k; k--) { |
91 | mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 3935559000370003845ULL)) |
92 | + init_key[j] + j; /* non linear */ |
93 | i++; j++; |
94 | if (i>=NN) { mt[0] = mt[NN-1]; i=1; } |
95 | if (j>=key_length) j=0; |
96 | } |
97 | for (k=NN-1; k; k--) { |
98 | mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 2862933555777941757ULL)) |
99 | - i; /* non linear */ |
100 | i++; |
101 | if (i>=NN) { mt[0] = mt[NN-1]; i=1; } |
102 | } |
103 | |
104 | mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */ |
105 | } |
106 | |
107 | /* generates a random number on [0, 2^64-1]-interval */ |
108 | unsigned long long genrand64_int64(void) |
109 | { |
110 | int i; |
111 | unsigned long long x; |
112 | static unsigned long long mag01[2]={0ULL, MATRIX_A}; |
113 | |
114 | if (mti >= NN) { /* generate NN words at one time */ |
115 | |
116 | /* if init_genrand64() has not been called, */ |
117 | /* a default initial seed is used */ |
118 | if (mti == NN+1) |
119 | init_genrand64(5489ULL); |
120 | |
121 | for (i=0;i<NN-MM;i++) { |
122 | x = (mt[i]&UM)|(mt[i+1]&LM); |
123 | mt[i] = mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)]; |
124 | } |
125 | for (;i<NN-1;i++) { |
126 | x = (mt[i]&UM)|(mt[i+1]&LM); |
127 | mt[i] = mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)]; |
128 | } |
129 | x = (mt[NN-1]&UM)|(mt[0]&LM); |
130 | mt[NN-1] = mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)]; |
131 | |
132 | mti = 0; |
133 | } |
134 | |
135 | x = mt[mti++]; |
136 | |
137 | x ^= (x >> 29) & 0x5555555555555555ULL; |
138 | x ^= (x << 17) & 0x71D67FFFEDA60000ULL; |
139 | x ^= (x << 37) & 0xFFF7EEE000000000ULL; |
140 | x ^= (x >> 43); |
141 | |
142 | return x; |
143 | } |
144 | |
145 | /* generates a random number on [0, 2^63-1]-interval */ |
146 | long long genrand64_int63(void) |
147 | { |
148 | return (long long)(genrand64_int64() >> 1); |
149 | } |
150 | |
151 | /* generates a random number on [0,1]-real-interval */ |
152 | double genrand64_real1(void) |
153 | { |
154 | return (genrand64_int64() >> 11) * (1.0/9007199254740991.0); |
155 | } |
156 | |
157 | /* generates a random number on [0,1)-real-interval */ |
158 | double genrand64_real2(void) |
159 | { |
160 | return (genrand64_int64() >> 11) * (1.0/9007199254740992.0); |
161 | } |
162 | |
163 | /* generates a random number on (0,1)-real-interval */ |
164 | double genrand64_real3(void) |
165 | { |
166 | return ((genrand64_int64() >> 12) + 0.5) * (1.0/4503599627370496.0); |
167 | } |
168 | |
169 | #ifdef MT19937_64_MAIN |
170 | int main(void) |
171 | { |
172 | int i; |
173 | unsigned long long init[4]={0x12345ULL, 0x23456ULL, 0x34567ULL, 0x45678ULL}, length=4; |
174 | init_by_array64(init, length); |
175 | printf("1000 outputs of genrand64_int64()\n" ); |
176 | for (i=0; i<1000; i++) { |
177 | printf("%20llu " , genrand64_int64()); |
178 | if (i%5==4) printf("\n" ); |
179 | } |
180 | printf("\n1000 outputs of genrand64_real2()\n" ); |
181 | for (i=0; i<1000; i++) { |
182 | printf("%10.8f " , genrand64_real2()); |
183 | if (i%5==4) printf("\n" ); |
184 | } |
185 | return 0; |
186 | } |
187 | #endif |
188 | |