← Back to branch summary

~ubuntu-branches/ubuntu/maverick/python3.1/maverick

« back to all changes in this revision

Viewing changes to Modules/_randommodule.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2009-03-23 00:01:27 UTC
  • Revision ID: james.westby@ubuntu.com-20090323000127-5fstfxju4ufrhthq
Tags: upstream-3.1~a1+20090322
Import upstream version 3.1~a1+20090322

Show diffs side-by-side

added added

removed removed

Lines of Context:
Modules/_randommodule.c
 
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
 
 
79
typedef struct {
 
80
        PyObject_HEAD
 
81
        unsigned long state[N];
 
82
        int index;
 
83
} RandomObject;
 
84
 
 
85
static 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 */
 
94
static unsigned long
 
95
genrand_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
 */
 
137
static PyObject *
 
138
random_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 */
 
145
static void
 
146
init_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 */
 
170
static PyObject *
 
171
init_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
 
 
206
static PyObject *
 
207
random_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 (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 = PyLong_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);
 
303
Done:
 
304
        Py_XDECREF(masklower);
 
305
        Py_XDECREF(thirtytwo);
 
306
        Py_XDECREF(n);
 
307
        PyMem_Free(key);
 
308
        return result;
 
309
}
 
310
 
 
311
static PyObject *
 
312
random_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
 
 
333
Fail:
 
334
        Py_DECREF(state);
 
335
        return NULL;
 
336
}
 
337
 
 
338
static PyObject *
 
339
random_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 == -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
static PyObject *
 
373
random_getrandbits(RandomObject *self, PyObject *args)
 
374
{
 
375
        int k, i, bytes;
 
376
        unsigned long r;
 
377
        unsigned char *bytearray;
 
378
        PyObject *result;
 
379
 
 
380
        if (!PyArg_ParseTuple(args, "i:getrandbits", &k))
 
381
                return NULL;
 
382
 
 
383
        if (k <= 0) {
 
384
                PyErr_SetString(PyExc_ValueError,
 
385
                                "number of bits must be greater than zero");
 
386
                return NULL;
 
387
        }
 
388
 
 
389
        bytes = ((k - 1) / 32 + 1) * 4;
 
390
        bytearray = (unsigned char *)PyMem_Malloc(bytes);
 
391
        if (bytearray == NULL) {
 
392
                PyErr_NoMemory();
 
393
                return NULL;
 
394
        }
 
395
 
 
396
        /* Fill-out whole words, byte-by-byte to avoid endianness issues */
 
397
        for (i=0 ; i<bytes ; i+=4, k-=32) {
 
398
                r = genrand_int32(self);
 
399
                if (k < 32)
 
400
                        r >>= (32 - k);
 
401
                bytearray[i+0] = (unsigned char)r;
 
402
                bytearray[i+1] = (unsigned char)(r >> 8);
 
403
                bytearray[i+2] = (unsigned char)(r >> 16);
 
404
                bytearray[i+3] = (unsigned char)(r >> 24);
 
405
        }
 
406
 
 
407
        /* little endian order to match bytearray assignment order */
 
408
        result = _PyLong_FromByteArray(bytearray, bytes, 1, 0);
 
409
        PyMem_Free(bytearray);
 
410
        return result;
 
411
}
 
412
 
 
413
static PyObject *
 
414
random_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
 
415
{
 
416
        RandomObject *self;
 
417
        PyObject *tmp;
 
418
 
 
419
        if (type == &Random_Type && !_PyArg_NoKeywords("Random()", kwds))
 
420
                return NULL;
 
421
 
 
422
        self = (RandomObject *)type->tp_alloc(type, 0);
 
423
        if (self == NULL)
 
424
                return NULL;
 
425
        tmp = random_seed(self, args);
 
426
        if (tmp == NULL) {
 
427
                Py_DECREF(self);
 
428
                return NULL;
 
429
        }
 
430
        Py_DECREF(tmp);
 
431
        return (PyObject *)self;
 
432
}
 
433
 
 
434
static PyMethodDef random_methods[] = {
 
435
        {"random",      (PyCFunction)random_random,  METH_NOARGS,
 
436
                PyDoc_STR("random() -> x in the interval [0, 1).")},
 
437
        {"seed",        (PyCFunction)random_seed,  METH_VARARGS,
 
438
                PyDoc_STR("seed([n]) -> None.  Defaults to current time.")},
 
439
        {"getstate",    (PyCFunction)random_getstate,  METH_NOARGS,
 
440
                PyDoc_STR("getstate() -> tuple containing the current state.")},
 
441
        {"setstate",      (PyCFunction)random_setstate,  METH_O,
 
442
                PyDoc_STR("setstate(state) -> None.  Restores generator state.")},
 
443
        {"getrandbits", (PyCFunction)random_getrandbits,  METH_VARARGS,
 
444
                PyDoc_STR("getrandbits(k) -> x.  Generates a long int with "
 
445
                          "k random bits.")},
 
446
        {NULL,          NULL}           /* sentinel */
 
447
};
 
448
 
 
449
PyDoc_STRVAR(random_doc,
 
450
"Random() -> create a random number generator with its own internal state.");
 
451
 
 
452
static PyTypeObject Random_Type = {
 
453
        PyVarObject_HEAD_INIT(NULL, 0)
 
454
        "_random.Random",               /*tp_name*/
 
455
        sizeof(RandomObject),           /*tp_basicsize*/
 
456
        0,                              /*tp_itemsize*/
 
457
        /* methods */
 
458
        0,                              /*tp_dealloc*/
 
459
        0,                              /*tp_print*/
 
460
        0,                              /*tp_getattr*/
 
461
        0,                              /*tp_setattr*/
 
462
        0,                              /*tp_reserved*/
 
463
        0,                              /*tp_repr*/
 
464
        0,                              /*tp_as_number*/
 
465
        0,                              /*tp_as_sequence*/
 
466
        0,                              /*tp_as_mapping*/
 
467
        0,                              /*tp_hash*/
 
468
        0,                              /*tp_call*/
 
469
        0,                              /*tp_str*/
 
470
        PyObject_GenericGetAttr,        /*tp_getattro*/
 
471
        0,                              /*tp_setattro*/
 
472
        0,                              /*tp_as_buffer*/
 
473
        Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,       /*tp_flags*/
 
474
        random_doc,                     /*tp_doc*/
 
475
        0,                              /*tp_traverse*/
 
476
        0,                              /*tp_clear*/
 
477
        0,                              /*tp_richcompare*/
 
478
        0,                              /*tp_weaklistoffset*/
 
479
        0,                              /*tp_iter*/
 
480
        0,                              /*tp_iternext*/
 
481
        random_methods,                 /*tp_methods*/
 
482
        0,                              /*tp_members*/
 
483
        0,                              /*tp_getset*/
 
484
        0,                              /*tp_base*/
 
485
        0,                              /*tp_dict*/
 
486
        0,                              /*tp_descr_get*/
 
487
        0,                              /*tp_descr_set*/
 
488
        0,                              /*tp_dictoffset*/
 
489
        0,                              /*tp_init*/
 
490
        0,                              /*tp_alloc*/
 
491
        random_new,                     /*tp_new*/
 
492
        PyObject_Free,                  /*tp_free*/
 
493
        0,                              /*tp_is_gc*/
 
494
};
 
495
 
 
496
PyDoc_STRVAR(module_doc,
 
497
"Module implements the Mersenne Twister random number generator.");
 
498
 
 
499
 
 
500
static struct PyModuleDef _randommodule = {
 
501
        PyModuleDef_HEAD_INIT,
 
502
        "_random",
 
503
        module_doc,
 
504
        -1,
 
505
        NULL,
 
506
        NULL,
 
507
        NULL,
 
508
        NULL,
 
509
        NULL
 
510
};
 
511
 
 
512
PyMODINIT_FUNC
 
513
PyInit__random(void)
 
514
{
 
515
        PyObject *m;
 
516
 
 
517
        if (PyType_Ready(&Random_Type) < 0)
 
518
                return NULL;
 
519
        m = PyModule_Create(&_randommodule);
 
520
        if (m == NULL)
 
521
                return NULL;
 
522
        Py_INCREF(&Random_Type);
 
523
        PyModule_AddObject(m, "Random", (PyObject *)&Random_Type);
 
524
        return m;
 
525
}