1 /* -*- Mode: c; c-basic-offset: 2 -*-
2 *
3 * mt.c - Mersenne Twister functions
4 *
5 * This is free and unencumbered software released into the public domain.
6 *
7 * Anyone is free to copy, modify, publish, use, compile, sell, or
8 * distribute this software, either in source code form or as a compiled
9 * binary, for any purpose, commercial or non-commercial, and by any
10 * means.
11 *
12 * In jurisdictions that recognize copyright laws, the author or authors
13 * of this software dedicate any and all copyright interest in the
14 * software to the public domain. We make this dedication for the benefit
15 * of the public at large and to the detriment of our heirs and
16 * successors. We intend this dedication to be an overt act of
17 * relinquishment in perpetuity of all present and future rights to this
18 * software under copyright law.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
22 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
23 * IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
24 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
25 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
26 * OTHER DEALINGS IN THE SOFTWARE.
27 *
28 * For more information, please refer to <http://unlicense.org/>
29 *
30 */
31
32
33 #ifdef MTWIST_CONFIG
34 #include <mtwist_config.h>
35 #endif
36
37 #include <stdio.h>
38
39 #ifdef HAVE_STDLIB_H
40 #include <stdlib.h>
41 #endif
42
43 #ifdef HAVE_STDINT_H
44 /* To get UINT32_C */
45 #define __STDC_CONSTANT_MACROS 1
46 #include <stdint.h>
47 #else
48 /* pessimistic - use an unsigned long */
49 typedef unsigned long uint32_t;
50 #define UINT32_C(v) (v ## UL)
51 #endif
52
53
54 #include <mtwist.h>
55
56 #include <mtwist_internal.h>
57
58 /*
59 * Mersenne Twister (MT19937) algorithm
60 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
61 *
62 * http://en.wikipedia.org/wiki/Mersenne_twister
63 *
64 */
65
66 #define MTWIST_UPPER_MASK UINT32_C(0x80000000)
67 #define MTWIST_LOWER_MASK UINT32_C(0x7FFFFFFF)
68 #define MTWIST_FULL_MASK UINT32_C(0xFFFFFFFF)
69
70 #define MTWIST_MATRIX_A UINT32_C(0x9908B0DF)
71
72 #define MTWIST_MIXBITS(u, v) ( ( (u) & MTWIST_UPPER_MASK) | ( (v) & MTWIST_LOWER_MASK) )
73 #define MTWIST_TWIST(u, v) ( (MTWIST_MIXBITS(u, v) >> 1) ^ ( (v) & UINT32_C(1) ? MTWIST_MATRIX_A : UINT32_C(0)) )
74
75
76 /**
77 * mtwist_new:
78 *
79 * Construct a Mersenne Twister object
80 *
81 * Return value: new MT object or NULL on failure
82 */
83 mtwist*
mtwist_new(void)84 mtwist_new(void)
85 {
86 mtwist* mt;
87
88 mt = (mtwist*)calloc(1, sizeof(*mt));
89 if(!mt)
90 return NULL;
91
92 mt->remaining = 0;
93 mt->next = NULL;
94 mt->seeded = 0;
95
96 return mt;
97 }
98
99
100 /**
101 * mtwist_free:
102 * @mt: mt object
103 *
104 * Destroy a Mersenne Twister object
105 */
106 void
mtwist_free(mtwist * mt)107 mtwist_free(mtwist* mt)
108 {
109 if(mt)
110 free(mt);
111 }
112
113
114 /**
115 * mtwist_init:
116 * @mt: mt object
117 * @seed: seed (lower 32 bits used)
118 *
119 * Initialise a Mersenne Twister with an unsigned 32 bit int seed
120 */
121 void
mtwist_init(mtwist * mt,unsigned long seed)122 mtwist_init(mtwist* mt, unsigned long seed)
123 {
124 int i;
125
126 if(!mt)
127 return;
128
129 mt->state[0] = (uint32_t)(seed & MTWIST_FULL_MASK);
130 for(i = 1; i < MTWIST_N; i++) {
131 mt->state[i] = (UINT32_C(1812433253) * (mt->state[i - 1] ^ (mt->state[i - 1] >> 30)) + i);
132 mt->state[i] &= MTWIST_FULL_MASK;
133 }
134
135 mt->remaining = 0;
136 mt->next = NULL;
137
138 mt->seeded = 1;
139 }
140
141
142 static void
mtwist_update_state(mtwist * mt)143 mtwist_update_state(mtwist* mt)
144 {
145 int count;
146 uint32_t *p = mt->state;
147
148 for(count = (MTWIST_N - MTWIST_M + 1); --count; p++)
149 *p = p[MTWIST_M] ^ MTWIST_TWIST(p[0], p[1]);
150
151 for(count = MTWIST_M; --count; p++)
152 *p = p[MTWIST_M - MTWIST_N] ^ MTWIST_TWIST(p[0], p[1]);
153
154 *p = p[MTWIST_M - MTWIST_N] ^ MTWIST_TWIST(p[0], mt->state[0]);
155
156 mt->remaining = MTWIST_N;
157 mt->next = mt->state;
158 }
159
160
161 /**
162 * mtwist_u32rand:
163 * @mt: mt object
164 *
165 * Get a random unsigned 32 bit integer from the random number generator
166 *
167 * Return value: unsigned long with 32 valid bits
168 */
169 unsigned long
mtwist_u32rand(mtwist * mt)170 mtwist_u32rand(mtwist* mt)
171 {
172 uint32_t r;
173
174 if(!mt)
175 return 0UL;
176
177 if(!mt->seeded)
178 mtwist_init(mt, mtwist_seed_from_system(mt));
179
180 if(!mt->remaining)
181 mtwist_update_state(mt);
182
183 r = *mt->next++;
184 mt->remaining--;
185
186 /* Tempering */
187 r ^= (r >> 11);
188 r ^= (r << 7) & UINT32_C(0x9D2C5680);
189 r ^= (r << 15) & UINT32_C(0xEFC60000);
190 r ^= (r >> 18);
191
192 r &= MTWIST_FULL_MASK;
193
194 return (unsigned long)r;
195 }
196
197
198 /**
199 * mtwist_drand:
200 * @mt: mt object
201 *
202 * Get a random double from the random number generator
203 *
204 * Return value: random double in the range 0.0 inclusive to 1.0 exclusive; [0.0, 1.0) */
205 double
mtwist_drand(mtwist * mt)206 mtwist_drand(mtwist* mt)
207 {
208 unsigned long r;
209 double d;
210
211 if(!mt)
212 return 0.0;
213
214 r = mtwist_u32rand(mt);
215
216 d = r / 4294967296.0; /* 2^32 */
217
218 return d;
219 }
220