PageRenderTime 716ms CodeModel.GetById 124ms app.highlight 249ms RepoModel.GetById 137ms app.codeStats 0ms

/Modules/_randommodule.c

http://unladen-swallow.googlecode.com/
C | 580 lines | 391 code | 60 blank | 129 comment | 89 complexity | f79f21b2c33a3f4dda811a0cc29d6f17 MD5 | raw file
  1/* Random objects */
  2
  3/* ------------------------------------------------------------------
  4   The code in this module was based on a download from:
  5	  http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html
  6
  7   It was modified in 2002 by Raymond Hettinger as follows:
  8
  9	* the principal computational lines untouched except for tabbing.
 10
 11	* renamed genrand_res53() to random_random() and wrapped
 12	  in python calling/return code.
 13
 14	* genrand_int32() and the helper functions, init_genrand()
 15	  and init_by_array(), were declared static, wrapped in
 16	  Python calling/return code.  also, their global data
 17	  references were replaced with structure references.
 18
 19	* unused functions from the original were deleted.
 20	  new, original C python code was added to implement the
 21	  Random() interface.
 22
 23   The following are the verbatim comments from the original code:
 24
 25   A C-program for MT19937, with initialization improved 2002/1/26.
 26   Coded by Takuji Nishimura and Makoto Matsumoto.
 27
 28   Before using, initialize the state by using init_genrand(seed)
 29   or init_by_array(init_key, key_length).
 30
 31   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
 32   All rights reserved.
 33
 34   Redistribution and use in source and binary forms, with or without
 35   modification, are permitted provided that the following conditions
 36   are met:
 37
 38     1. Redistributions of source code must retain the above copyright
 39	notice, this list of conditions and the following disclaimer.
 40
 41     2. Redistributions in binary form must reproduce the above copyright
 42	notice, this list of conditions and the following disclaimer in the
 43	documentation and/or other materials provided with the distribution.
 44
 45     3. The names of its contributors may not be used to endorse or promote
 46	products derived from this software without specific prior written
 47	permission.
 48
 49   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 50   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 51   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 52   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 53   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 54   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 55   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 56   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 57   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 58   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 59   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 60
 61
 62   Any feedback is very welcome.
 63   http://www.math.keio.ac.jp/matumoto/emt.html
 64   email: matumoto@math.keio.ac.jp
 65*/
 66
 67/* ---------------------------------------------------------------*/
 68
 69#include "Python.h"
 70#include <time.h>		/* for seeding to current time */
 71
 72/* Period parameters -- These are all magic.  Don't change. */
 73#define N 624
 74#define M 397
 75#define MATRIX_A 0x9908b0dfUL	/* constant vector a */
 76#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
 77#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
 78
 79typedef struct {
 80	PyObject_HEAD
 81	unsigned long state[N];
 82	int index;
 83} RandomObject;
 84
 85static PyTypeObject Random_Type;
 86
 87#define RandomObject_Check(v)	   (Py_TYPE(v) == &Random_Type)
 88
 89
 90/* Random methods */
 91
 92
 93/* generates a random number on [0,0xffffffff]-interval */
 94static unsigned long
 95genrand_int32(RandomObject *self)
 96{
 97	unsigned long y;
 98	static unsigned long mag01[2]={0x0UL, MATRIX_A};
 99	/* mag01[x] = x * MATRIX_A  for x=0,1 */
100	unsigned long *mt;
101
102	mt = self->state;
103	if (self->index >= N) { /* generate N words at one time */
104		int kk;
105
106		for (kk=0;kk<N-M;kk++) {
107			y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
108			mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
109		}
110		for (;kk<N-1;kk++) {
111			y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
112			mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
113		}
114		y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
115		mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
116
117		self->index = 0;
118	}
119
120    y = mt[self->index++];
121    y ^= (y >> 11);
122    y ^= (y << 7) & 0x9d2c5680UL;
123    y ^= (y << 15) & 0xefc60000UL;
124    y ^= (y >> 18);
125    return y;
126}
127
128/* random_random is the function named genrand_res53 in the original code;
129 * generates a random number on [0,1) with 53-bit resolution; note that
130 * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
131 * multiply-by-reciprocal in the (likely vain) hope that the compiler will
132 * optimize the division away at compile-time.  67108864 is 2**26.  In
133 * effect, a contains 27 random bits shifted left 26, and b fills in the
134 * lower 26 bits of the 53-bit numerator.
135 * The orginal code credited Isaku Wada for this algorithm, 2002/01/09.
136 */
137static PyObject *
138random_random(RandomObject *self)
139{
140	unsigned long a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
141    	return PyFloat_FromDouble((a*67108864.0+b)*(1.0/9007199254740992.0));
142}
143
144/* initializes mt[N] with a seed */
145static void
146init_genrand(RandomObject *self, unsigned long s)
147{
148	int mti;
149	unsigned long *mt;
150
151	mt = self->state;
152	mt[0]= s & 0xffffffffUL;
153	for (mti=1; mti<N; mti++) {
154		mt[mti] =
155		(1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
156		/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
157		/* In the previous versions, MSBs of the seed affect   */
158		/* only MSBs of the array mt[]. 		       */
159		/* 2002/01/09 modified by Makoto Matsumoto	       */
160		mt[mti] &= 0xffffffffUL;
161		/* for >32 bit machines */
162	}
163	self->index = mti;
164	return;
165}
166
167/* initialize by an array with array-length */
168/* init_key is the array for initializing keys */
169/* key_length is its length */
170static PyObject *
171init_by_array(RandomObject *self, unsigned long init_key[], unsigned long key_length)
172{
173	unsigned int i, j, k;	/* was signed in the original code. RDH 12/16/2002 */
174	unsigned long *mt;
175
176	mt = self->state;
177	init_genrand(self, 19650218UL);
178	i=1; j=0;
179	k = (N>key_length ? N : key_length);
180	for (; k; k--) {
181		mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
182			 + init_key[j] + j; /* non linear */
183		mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
184		i++; j++;
185		if (i>=N) { mt[0] = mt[N-1]; i=1; }
186		if (j>=key_length) j=0;
187	}
188	for (k=N-1; k; k--) {
189		mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
190			 - i; /* non linear */
191		mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
192		i++;
193		if (i>=N) { mt[0] = mt[N-1]; i=1; }
194	}
195
196    mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
197    Py_INCREF(Py_None);
198    return Py_None;
199}
200
201/*
202 * The rest is Python-specific code, neither part of, nor derived from, the
203 * Twister download.
204 */
205
206static PyObject *
207random_seed(RandomObject *self, PyObject *args)
208{
209	PyObject *result = NULL;	/* guilty until proved innocent */
210	PyObject *masklower = NULL;
211	PyObject *thirtytwo = NULL;
212	PyObject *n = NULL;
213	unsigned long *key = NULL;
214	unsigned long keymax;		/* # of allocated slots in key */
215	unsigned long keyused;		/* # of used slots in key */
216	int err;
217
218	PyObject *arg = NULL;
219
220	if (!PyArg_UnpackTuple(args, "seed", 0, 1, &arg))
221		return NULL;
222
223	if (arg == NULL || arg == Py_None) {
224		time_t now;
225
226		time(&now);
227		init_genrand(self, (unsigned long)now);
228		Py_INCREF(Py_None);
229		return Py_None;
230	}
231	/* If the arg is an int or long, use its absolute value; else use
232	 * the absolute value of its hash code.
233	 */
234	if (PyInt_Check(arg) || PyLong_Check(arg))
235		n = PyNumber_Absolute(arg);
236	else {
237		long hash = PyObject_Hash(arg);
238		if (hash == -1)
239			goto Done;
240		n = PyLong_FromUnsignedLong((unsigned long)hash);
241	}
242	if (n == NULL)
243		goto Done;
244
245	/* Now split n into 32-bit chunks, from the right.  Each piece is
246	 * stored into key, which has a capacity of keymax chunks, of which
247	 * keyused are filled.  Alas, the repeated shifting makes this a
248	 * quadratic-time algorithm; we'd really like to use
249	 * _PyLong_AsByteArray here, but then we'd have to break into the
250	 * long representation to figure out how big an array was needed
251	 * in advance.
252	 */
253	keymax = 8; 	/* arbitrary; grows later if needed */
254	keyused = 0;
255	key = (unsigned long *)PyMem_Malloc(keymax * sizeof(*key));
256	if (key == NULL)
257		goto Done;
258
259	masklower = PyLong_FromUnsignedLong(0xffffffffU);
260	if (masklower == NULL)
261		goto Done;
262	thirtytwo = PyInt_FromLong(32L);
263	if (thirtytwo == NULL)
264		goto Done;
265	while ((err=PyObject_IsTrue(n))) {
266		PyObject *newn;
267		PyObject *pychunk;
268		unsigned long chunk;
269
270		if (err == -1)
271			goto Done;
272		pychunk = PyNumber_And(n, masklower);
273		if (pychunk == NULL)
274			goto Done;
275		chunk = PyLong_AsUnsignedLong(pychunk);
276		Py_DECREF(pychunk);
277		if (chunk == (unsigned long)-1 && PyErr_Occurred())
278			goto Done;
279		newn = PyNumber_Rshift(n, thirtytwo);
280		if (newn == NULL)
281			goto Done;
282		Py_DECREF(n);
283		n = newn;
284		if (keyused >= keymax) {
285			unsigned long bigger = keymax << 1;
286			if ((bigger >> 1) != keymax) {
287				PyErr_NoMemory();
288				goto Done;
289			}
290			key = (unsigned long *)PyMem_Realloc(key,
291						bigger * sizeof(*key));
292			if (key == NULL)
293				goto Done;
294			keymax = bigger;
295		}
296		assert(keyused < keymax);
297		key[keyused++] = chunk;
298	}
299
300	if (keyused == 0)
301		key[keyused++] = 0UL;
302	result = init_by_array(self, key, keyused);
303Done:
304	Py_XDECREF(masklower);
305	Py_XDECREF(thirtytwo);
306	Py_XDECREF(n);
307	PyMem_Free(key);
308	return result;
309}
310
311static PyObject *
312random_getstate(RandomObject *self)
313{
314	PyObject *state;
315	PyObject *element;
316	int i;
317
318	state = PyTuple_New(N+1);
319	if (state == NULL)
320		return NULL;
321	for (i=0; i<N ; i++) {
322		element = PyLong_FromUnsignedLong(self->state[i]);
323		if (element == NULL)
324			goto Fail;
325		PyTuple_SET_ITEM(state, i, element);
326	}
327	element = PyLong_FromLong((long)(self->index));
328	if (element == NULL)
329		goto Fail;
330	PyTuple_SET_ITEM(state, i, element);
331	return state;
332
333Fail:
334	Py_DECREF(state);
335	return NULL;
336}
337
338static PyObject *
339random_setstate(RandomObject *self, PyObject *state)
340{
341	int i;
342	unsigned long element;
343	long index;
344
345	if (!PyTuple_Check(state)) {
346		PyErr_SetString(PyExc_TypeError,
347			"state vector must be a tuple");
348		return NULL;
349	}
350	if (PyTuple_Size(state) != N+1) {
351		PyErr_SetString(PyExc_ValueError,
352			"state vector is the wrong size");
353		return NULL;
354	}
355
356	for (i=0; i<N ; i++) {
357		element = PyLong_AsUnsignedLong(PyTuple_GET_ITEM(state, i));
358		if (element == (unsigned long)-1 && PyErr_Occurred())
359			return NULL;
360		self->state[i] = element & 0xffffffffUL; /* Make sure we get sane state */
361	}
362
363	index = PyLong_AsLong(PyTuple_GET_ITEM(state, i));
364	if (index == -1 && PyErr_Occurred())
365		return NULL;
366	self->index = (int)index;
367
368	Py_INCREF(Py_None);
369	return Py_None;
370}
371
372/*
373Jumpahead should be a fast way advance the generator n-steps ahead, but
374lacking a formula for that, the next best is to use n and the existing
375state to create a new state far away from the original.
376
377The generator uses constant spaced additive feedback, so shuffling the
378state elements ought to produce a state which would not be encountered
379(in the near term) by calls to random().  Shuffling is normally
380implemented by swapping the ith element with another element ranging
381from 0 to i inclusive.  That allows the element to have the possibility
382of not being moved.  Since the goal is to produce a new, different
383state, the swap element is ranged from 0 to i-1 inclusive.  This assures
384that each element gets moved at least once.
385
386To make sure that consecutive calls to jumpahead(n) produce different
387states (even in the rare case of involutory shuffles), i+1 is added to
388each element at position i.  Successive calls are then guaranteed to
389have changing (growing) values as well as shuffled positions.
390
391Finally, the self->index value is set to N so that the generator itself
392kicks in on the next call to random().	This assures that all results
393have been through the generator and do not just reflect alterations to
394the underlying state.
395*/
396
397static PyObject *
398random_jumpahead(RandomObject *self, PyObject *n)
399{
400	long i, j;
401	PyObject *iobj;
402	PyObject *remobj;
403	unsigned long *mt, tmp;
404
405	if (!PyInt_Check(n) && !PyLong_Check(n)) {
406		PyErr_Format(PyExc_TypeError, "jumpahead requires an "
407			     "integer, not '%s'",
408			     Py_TYPE(n)->tp_name);
409		return NULL;
410	}
411
412	mt = self->state;
413	for (i = N-1; i > 1; i--) {
414		iobj = PyInt_FromLong(i);
415		if (iobj == NULL)
416			return NULL;
417		remobj = PyNumber_Remainder(n, iobj);
418		Py_DECREF(iobj);
419		if (remobj == NULL)
420			return NULL;
421		j = PyInt_AsLong(remobj);
422		Py_DECREF(remobj);
423		if (j == -1L && PyErr_Occurred())
424			return NULL;
425		tmp = mt[i];
426		mt[i] = mt[j];
427		mt[j] = tmp;
428	}
429
430	for (i = 0; i < N; i++)
431		mt[i] += i+1;
432
433	self->index = N;
434	Py_INCREF(Py_None);
435	return Py_None;
436}
437
438static PyObject *
439random_getrandbits(RandomObject *self, PyObject *args)
440{
441	int k, i, bytes;
442	unsigned long r;
443	unsigned char *bytearray;
444	PyObject *result;
445
446	if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
447		return NULL;
448
449	if (k <= 0) {
450		PyErr_SetString(PyExc_ValueError,
451				"number of bits must be greater than zero");
452		return NULL;
453	}
454
455	bytes = ((k - 1) / 32 + 1) * 4;
456	bytearray = (unsigned char *)PyMem_Malloc(bytes);
457	if (bytearray == NULL) {
458		PyErr_NoMemory();
459		return NULL;
460	}
461
462	/* Fill-out whole words, byte-by-byte to avoid endianness issues */
463	for (i=0 ; i<bytes ; i+=4, k-=32) {
464		r = genrand_int32(self);
465		if (k < 32)
466			r >>= (32 - k);
467		bytearray[i+0] = (unsigned char)r;
468		bytearray[i+1] = (unsigned char)(r >> 8);
469		bytearray[i+2] = (unsigned char)(r >> 16);
470		bytearray[i+3] = (unsigned char)(r >> 24);
471	}
472
473	/* little endian order to match bytearray assignment order */
474	result = _PyLong_FromByteArray(bytearray, bytes, 1, 0);
475	PyMem_Free(bytearray);
476	return result;
477}
478
479static PyObject *
480random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
481{
482	RandomObject *self;
483	PyObject *tmp;
484
485	if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds))
486		return NULL;
487
488	self = (RandomObject *)type->tp_alloc(type, 0);
489	if (self == NULL)
490		return NULL;
491	tmp = random_seed(self, args);
492	if (tmp == NULL) {
493		Py_DECREF(self);
494		return NULL;
495	}
496	Py_DECREF(tmp);
497	return (PyObject *)self;
498}
499
500static PyMethodDef random_methods[] = {
501	{"random",	(PyCFunction)random_random,  METH_NOARGS,
502		PyDoc_STR("random() -> x in the interval [0, 1).")},
503	{"seed",	(PyCFunction)random_seed,  METH_VARARGS,
504		PyDoc_STR("seed([n]) -> None.  Defaults to current time.")},
505	{"getstate",	(PyCFunction)random_getstate,  METH_NOARGS,
506		PyDoc_STR("getstate() -> tuple containing the current state.")},
507	{"setstate",	  (PyCFunction)random_setstate,  METH_O,
508		PyDoc_STR("setstate(state) -> None.  Restores generator state.")},
509	{"jumpahead",	(PyCFunction)random_jumpahead,	METH_O,
510		PyDoc_STR("jumpahead(int) -> None.  Create new state from "
511			  "existing state and integer.")},
512	{"getrandbits",	(PyCFunction)random_getrandbits,  METH_VARARGS,
513		PyDoc_STR("getrandbits(k) -> x.  Generates a long int with "
514			  "k random bits.")},
515	{NULL,		NULL}		/* sentinel */
516};
517
518PyDoc_STRVAR(random_doc,
519"Random() -> create a random number generator with its own internal state.");
520
521static PyTypeObject Random_Type = {
522	PyVarObject_HEAD_INIT(NULL, 0)
523	"_random.Random",		/*tp_name*/
524	sizeof(RandomObject),		/*tp_basicsize*/
525	0,				/*tp_itemsize*/
526	/* methods */
527	0,				/*tp_dealloc*/
528	0,				/*tp_print*/
529	0,				/*tp_getattr*/
530	0,				/*tp_setattr*/
531	0,				/*tp_compare*/
532	0,				/*tp_repr*/
533	0,				/*tp_as_number*/
534	0,				/*tp_as_sequence*/
535	0,				/*tp_as_mapping*/
536	0,				/*tp_hash*/
537	0,				/*tp_call*/
538	0,				/*tp_str*/
539	PyObject_GenericGetAttr,	/*tp_getattro*/
540	0,				/*tp_setattro*/
541	0,				/*tp_as_buffer*/
542	Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,	/*tp_flags*/
543	random_doc,			/*tp_doc*/
544	0,				/*tp_traverse*/
545	0,				/*tp_clear*/
546	0,				/*tp_richcompare*/
547	0,				/*tp_weaklistoffset*/
548	0,				/*tp_iter*/
549	0,				/*tp_iternext*/
550	random_methods, 		/*tp_methods*/
551	0,				/*tp_members*/
552	0,				/*tp_getset*/
553	0,				/*tp_base*/
554	0,				/*tp_dict*/
555	0,				/*tp_descr_get*/
556	0,				/*tp_descr_set*/
557	0,				/*tp_dictoffset*/
558	0,				/*tp_init*/
559	0,				/*tp_alloc*/
560	random_new,			/*tp_new*/
561	_PyObject_Del,			/*tp_free*/
562	0,				/*tp_is_gc*/
563};
564
565PyDoc_STRVAR(module_doc,
566"Module implements the Mersenne Twister random number generator.");
567
568PyMODINIT_FUNC
569init_random(void)
570{
571	PyObject *m;
572
573	if (PyType_Ready(&Random_Type) < 0)
574		return;
575	m = Py_InitModule3("_random", NULL, module_doc);
576	if (m == NULL)
577		return;
578	Py_INCREF(&Random_Type);
579	PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
580}