Aluode/PerceptionLabPortable
0
1/*
2 * PRNG support.
3 */
4
5#ifdef _MSC_VER
6#define HAVE_PTHREAD_ATFORK 0
7#else
8#define HAVE_PTHREAD_ATFORK 1
9#include <pthread.h>
10#endif
11
12
13/* Magic Mersenne Twister constants */
14#define MT_N 624
15#define MT_M 397
16#define MT_MATRIX_A 0x9908b0dfU
17#define MT_UPPER_MASK 0x80000000U
18#define MT_LOWER_MASK 0x7fffffffU
19
20/*
21 * Note this structure is accessed in numba.targets.randomimpl,
22 * any changes here should be reflected there too.
23 */
24typedef struct {
25 int index;
26 /* unsigned int is sufficient on modern machines as we only need 32 bits */
27 unsigned int mt[MT_N];
28 int has_gauss;
29 double gauss;
30 int is_initialized;
31} rnd_state_t;
32
33/* Some code portions below from CPython's _randommodule.c, some others
34 from Numpy's and Jean-Sebastien Roy's randomkit.c. */
35
36NUMBA_EXPORT_FUNC(void)
37numba_rnd_shuffle(rnd_state_t *state)
38{
39 int i;
40 unsigned int y;
41
42 for (i = 0; i < MT_N - MT_M; i++) {
43 y = (state->mt[i] & MT_UPPER_MASK) | (state->mt[i+1] & MT_LOWER_MASK);
44 state->mt[i] = state->mt[i+MT_M] ^ (y >> 1) ^
45 (-(int) (y & 1) & MT_MATRIX_A);
46 }
47 for (; i < MT_N - 1; i++) {
48 y = (state->mt[i] & MT_UPPER_MASK) | (state->mt[i+1] & MT_LOWER_MASK);
49 state->mt[i] = state->mt[i+(MT_M-MT_N)] ^ (y >> 1) ^
50 (-(int) (y & 1) & MT_MATRIX_A);
51 }
52 y = (state->mt[MT_N - 1] & MT_UPPER_MASK) | (state->mt[0] & MT_LOWER_MASK);
53 state->mt[MT_N - 1] = state->mt[MT_M - 1] ^ (y >> 1) ^
54 (-(int) (y & 1) & MT_MATRIX_A);
55}
56
57/* Initialize mt[] with an integer seed */
58NUMBA_EXPORT_FUNC(void)
59numba_rnd_init(rnd_state_t *state, unsigned int seed)
60{
61 unsigned int pos;
62 seed &= 0xffffffffU;
63
64 /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
65 for (pos = 0; pos < MT_N; pos++) {
66 state->mt[pos] = seed;
67 seed = (1812433253U * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffU;
68 }
69 state->index = MT_N;
70 state->has_gauss = 0;
71 state->gauss = 0.0;
72 state->is_initialized = 1;
73}
74
75/* Perturb mt[] with a key array */
76static void
77rnd_init_by_array(rnd_state_t *state, unsigned int init_key[], size_t key_length)
78{
79 size_t i, j, k;
80 unsigned int *mt = state->mt;
81
82 numba_rnd_init(state, 19650218U);
83 i = 1; j = 0;
84 k = (MT_N > key_length ? MT_N : key_length);
85 for (; k; k--) {
86 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
87 + init_key[j] + (unsigned int) j; /* non linear */
88 mt[i] &= 0xffffffffU;
89 i++; j++;
90 if (i >= MT_N) { mt[0] = mt[MT_N - 1]; i = 1; }
91 if (j >= key_length) j = 0;
92 }
93 for (k = MT_N - 1; k; k--) {
94 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))
95 - (unsigned int) i; /* non linear */
96 mt[i] &= 0xffffffffU;
97 i++;
98 if (i >= MT_N) { mt[0] = mt[MT_N - 1]; i=1; }
99 }
100
101 mt[0] = 0x80000000U; /* MSB is 1; ensuring non-zero initial array */
102 state->index = MT_N;
103 state->has_gauss = 0;
104 state->gauss = 0.0;
105 state->is_initialized = 1;
106}
107
108/*
109 * Management of thread-local random state.
110 */
111
112static int rnd_globally_initialized;
113
114#ifdef _MSC_VER
115#define THREAD_LOCAL(ty) __declspec(thread) ty
116#else
117/* Non-standard C99 extension that's understood by gcc and clang */
118#define THREAD_LOCAL(ty) __thread ty
119#endif
120
121static THREAD_LOCAL(rnd_state_t) numba_py_random_state;
122static THREAD_LOCAL(rnd_state_t) numba_np_random_state;
123static THREAD_LOCAL(rnd_state_t) numba_internal_random_state;
124
125/* Seed the state with random bytes */
126static int
127rnd_seed_with_bytes(rnd_state_t *state, Py_buffer *buf)
128{
129 unsigned int *keys;
130 unsigned char *bytes;
131 size_t i, nkeys;
132
133 nkeys = buf->len / sizeof(unsigned int);
134 keys = (unsigned int *) PyMem_Malloc(nkeys * sizeof(unsigned int));
135 if (keys == NULL) {
136 PyBuffer_Release(buf);
137 return -1;
138 }
139 bytes = (unsigned char *) buf->buf;
140 /* Convert input bytes to int32 keys, without violating alignment
141 * constraints.
142 */
143 for (i = 0; i < nkeys; i++, bytes += 4) {
144 keys[i] =
145 ((unsigned int)bytes[3] << 24) +
146 ((unsigned int)bytes[2] << 16) +
147 ((unsigned int)bytes[1] << 8) +
148 ((unsigned int)bytes[0] << 0);
149 }
150 PyBuffer_Release(buf);
151 rnd_init_by_array(state, keys, nkeys);
152 PyMem_Free(keys);
153 return 0;
154}
155
156#if HAVE_PTHREAD_ATFORK
157/* After a fork(), the child should reseed its random states.
158 * Since only the main thread survives in the child, it's enough to mark
159 * the current thread-local states as uninitialized.
160 */
161static void
162rnd_atfork_child(void)
163{
164 numba_py_random_state.is_initialized = 0;
165 numba_np_random_state.is_initialized = 0;
166 numba_internal_random_state.is_initialized = 0;
167}
168#endif
169
170/* Global initialization routine. It must be called as early as possible.
171 */
172NUMBA_EXPORT_FUNC(void)
173numba_rnd_ensure_global_init(void)
174{
175 if (!rnd_globally_initialized) {
176#if HAVE_PTHREAD_ATFORK
177 pthread_atfork(NULL, NULL, rnd_atfork_child);
178#endif
179 numba_py_random_state.is_initialized = 0;
180 numba_np_random_state.is_initialized = 0;
181 numba_internal_random_state.is_initialized = 0;
182 rnd_globally_initialized = 1;
183 }
184}
185
186/* First-time init a random state */
187static void
188rnd_implicit_init(rnd_state_t *state)
189{
190 /* Initialize with random bytes. The easiest way to get good-quality
191 * cross-platform random bytes is still to call os.urandom()
192 * using the Python interpreter...
193 */
194 PyObject *module, *bufobj;
195 Py_buffer buf;
196 PyGILState_STATE gilstate = PyGILState_Ensure();
197
198 module = PyImport_ImportModule("os");
199 if (module == NULL)
200 goto error;
201 /* Read as many bytes as necessary to get the full entropy
202 * exploitable by the MT generator.
203 */
204 bufobj = PyObject_CallMethod(module, "urandom", "i",
205 (int) (MT_N * sizeof(unsigned int)));
206 Py_DECREF(module);
207 if (bufobj == NULL)
208 goto error;
209 if (PyObject_GetBuffer(bufobj, &buf, PyBUF_SIMPLE))
210 goto error;
211 Py_DECREF(bufobj);
212 if (rnd_seed_with_bytes(state, &buf))
213 goto error;
214 /* state->is_initialized is set now */
215
216 PyGILState_Release(gilstate);
217 return;
218
219error:
220 /* In normal conditions, os.urandom() and PyMem_Malloc() shouldn't fail,
221 * and we don't want the caller to deal with errors, so just bail out.
222 */
223 if (PyErr_Occurred())
224 PyErr_Print();
225 Py_FatalError(NULL);
226}
227
228/* Functions returning the thread-local random state pointer.
229 * The LLVM JIT doesn't support thread-local variables so we rely
230 * on the C compiler instead.
231 */
232
233NUMBA_EXPORT_FUNC(rnd_state_t *)
234numba_get_py_random_state(void)
235{
236 rnd_state_t *state = &numba_py_random_state;
237 if (!state->is_initialized)
238 rnd_implicit_init(state);
239 return state;
240}
241
242NUMBA_EXPORT_FUNC(rnd_state_t *)
243numba_get_np_random_state(void)
244{
245 rnd_state_t *state = &numba_np_random_state;
246 if (!state->is_initialized)
247 rnd_implicit_init(state);
248 return state;
249}
250
251NUMBA_EXPORT_FUNC(rnd_state_t *)
252numba_get_internal_random_state(void)
253{
254 rnd_state_t *state = &numba_internal_random_state;
255 if (!state->is_initialized)
256 rnd_implicit_init(state);
257 return state;
258}
259
260/*
261 * Python-exposed helpers for state management and testing.
262 */
263static int
264rnd_state_converter(PyObject *obj, rnd_state_t **state)
265{
266 *state = (rnd_state_t *) PyLong_AsVoidPtr(obj);
267 return (*state != NULL || !PyErr_Occurred());
268}
269
270NUMBA_EXPORT_FUNC(PyObject *)
271_numba_rnd_get_py_state_ptr(PyObject *self)
272{
273 return PyLong_FromVoidPtr(numba_get_py_random_state());
274}
275
276NUMBA_EXPORT_FUNC(PyObject *)
277_numba_rnd_get_np_state_ptr(PyObject *self)
278{
279 return PyLong_FromVoidPtr(numba_get_np_random_state());
280}
281
282NUMBA_EXPORT_FUNC(PyObject *)
283_numba_rnd_shuffle(PyObject *self, PyObject *arg)
284{
285 rnd_state_t *state;
286 if (!rnd_state_converter(arg, &state))
287 return NULL;
288 numba_rnd_shuffle(state);
289 Py_RETURN_NONE;
290}
291
292NUMBA_EXPORT_FUNC(PyObject *)
293_numba_rnd_set_state(PyObject *self, PyObject *args)
294{
295 int i, index;
296 rnd_state_t *state;
297 PyObject *tuplearg, *intlist;
298
299 if (!PyArg_ParseTuple(args, "O&O!:rnd_set_state",
300 rnd_state_converter, &state,
301 &PyTuple_Type, &tuplearg))
302 return NULL;
303 if (!PyArg_ParseTuple(tuplearg, "iO!", &index, &PyList_Type, &intlist))
304 return NULL;
305 if (PyList_GET_SIZE(intlist) != MT_N) {
306 PyErr_SetString(PyExc_ValueError, "list object has wrong size");
307 return NULL;
308 }
309 state->index = index;
310 for (i = 0; i < MT_N; i++) {
311 PyObject *v = PyList_GET_ITEM(intlist, i);
312 unsigned long x = PyLong_AsUnsignedLong(v);
313 if (x == (unsigned long) -1 && PyErr_Occurred())
314 return NULL;
315 state->mt[i] = (unsigned int) x;
316 }
317 state->has_gauss = 0;
318 state->gauss = 0.0;
319 state->is_initialized = 1;
320 Py_RETURN_NONE;
321}
322
323NUMBA_EXPORT_FUNC(PyObject *)
324_numba_rnd_get_state(PyObject *self, PyObject *arg)
325{
326 PyObject *intlist;
327 int i;
328 rnd_state_t *state;
329 if (!rnd_state_converter(arg, &state))
330 return NULL;
331
332 intlist = PyList_New(MT_N);
333 if (intlist == NULL)
334 return NULL;
335 for (i = 0; i < MT_N; i++) {
336 PyObject *v = PyLong_FromUnsignedLong(state->mt[i]);
337 if (v == NULL) {
338 Py_DECREF(intlist);
339 return NULL;
340 }
341 PyList_SET_ITEM(intlist, i, v);
342 }
343 return Py_BuildValue("iN", state->index, intlist);
344}
345
346NUMBA_EXPORT_FUNC(PyObject *)
347_numba_rnd_seed(PyObject *self, PyObject *args)
348{
349 unsigned int seed;
350 rnd_state_t *state;
351
352 if (!PyArg_ParseTuple(args, "O&I:rnd_seed",
353 rnd_state_converter, &state, &seed)) {
354 /* rnd_seed_*(bytes-like object) */
355 Py_buffer buf;
356
357 PyErr_Clear();
358 if (!PyArg_ParseTuple(args, "O&s*:rnd_seed",
359 rnd_state_converter, &state, &buf))
360 return NULL;
361
362 if (rnd_seed_with_bytes(state, &buf))
363 return NULL;
364 else
365 Py_RETURN_NONE;
366 }
367 else {
368 /* rnd_seed_*(int32) */
369 numba_rnd_init(state, seed);
370 Py_RETURN_NONE;
371 }
372}
373
374/*
375 * Random distribution helpers.
376 * Most code straight from Numpy's distributions.c.
377 */
378
379#ifndef M_PI
380#define M_PI 3.14159265358979323846264338328
381#endif
382
383NUMBA_EXPORT_FUNC(unsigned int)
384get_next_int32(rnd_state_t *state)
385{
386 unsigned int y;
387
388 if (state->index == MT_N) {
389 numba_rnd_shuffle(state);
390 state->index = 0;
391 }
392 y = state->mt[state->index++];
393 /* Tempering */
394 y ^= (y >> 11);
395 y ^= (y << 7) & 0x9d2c5680U;
396 y ^= (y << 15) & 0xefc60000U;
397 y ^= (y >> 18);
398 return y;
399}
400
401NUMBA_EXPORT_FUNC(double)
402get_next_double(rnd_state_t *state)
403{
404 double a = get_next_int32(state) >> 5;
405 double b = get_next_int32(state) >> 6;
406 return (a * 67108864.0 + b) / 9007199254740992.0;
407}
408
409NUMBA_EXPORT_FUNC(double)
410loggam(double x)
411{
412 double x0, x2, xp, gl, gl0;
413 long k, n;
414
415 static double a[10] = {8.333333333333333e-02,-2.777777777777778e-03,
416 7.936507936507937e-04,-5.952380952380952e-04,
417 8.417508417508418e-04,-1.917526917526918e-03,
418 6.410256410256410e-03,-2.955065359477124e-02,
419 1.796443723688307e-01,-1.39243221690590e+00};
420 x0 = x;
421 n = 0;
422 if ((x == 1.0) || (x == 2.0))
423 {
424 return 0.0;
425 }
426 else if (x <= 7.0)
427 {
428 n = (long)(7 - x);
429 x0 = x + n;
430 }
431 x2 = 1.0/(x0*x0);
432 xp = 2*M_PI;
433 gl0 = a[9];
434 for (k=8; k>=0; k--)
435 {
436 gl0 *= x2;
437 gl0 += a[k];
438 }
439 gl = gl0/x0 + 0.5*log(xp) + (x0-0.5)*log(x0) - x0;
440 if (x <= 7.0)
441 {
442 for (k=1; k<=n; k++)
443 {
444 gl -= log(x0-1.0);
445 x0 -= 1.0;
446 }
447 }
448 return gl;
449}
450
451
452NUMBA_EXPORT_FUNC(int64_t)
453numba_poisson_ptrs(rnd_state_t *state, double lam)
454{
455 /* This method is invoked only if the parameter lambda of this
456 * distribution is big enough ( >= 10 ). The algorithm used is
457 * described in "Hörmann, W. 1992. 'The Transformed Rejection
458 * Method for Generating Poisson Random Variables'.
459 * The implementation comes straight from Numpy.
460 */
461 int64_t k;
462 double U, V, slam, loglam, a, b, invalpha, vr, us;
463
464 slam = sqrt(lam);
465 loglam = log(lam);
466 b = 0.931 + 2.53*slam;
467 a = -0.059 + 0.02483*b;
468 invalpha = 1.1239 + 1.1328/(b-3.4);
469 vr = 0.9277 - 3.6224/(b-2);
470
471 while (1)
472 {
473 U = get_next_double(state) - 0.5;
474 V = get_next_double(state);
475 us = 0.5 - fabs(U);
476 k = (int64_t) floor((2*a/us + b)*U + lam + 0.43);
477 if ((us >= 0.07) && (V <= vr))
478 {
479 return k;
480 }
481 if ((k < 0) ||
482 ((us < 0.013) && (V > us)))
483 {
484 continue;
485 }
486 if ((log(V) + log(invalpha) - log(a/(us*us)+b)) <=
487 (-lam + (double) k*loglam - loggam((double) k+1)))
488 {
489 return k;
490 }
491 }
492}
493