random number generation from cython

[Flavio Coelho]

Hi,

anyone know how could I access the Python's builtin mersenne twister
random number generator from Cython?

I'd like to generate random numbers in the interval [0,1[ faster than
by calling random.random(). I know the implementation of random.random
is in C, but there is the cost of a python call. I want something I
can define as a "cdef extern" in Cython.

thanks,

Flávio

[Jason Grout]

You might read the (very good) documentation and code in
sage/misc/randstate.pyx.

I just noticed that this file has not been converted to rest, so it is
not in the reference manual. Ouch; it has such nice documentation!

Thanks,

Jason

--
Jason Grout

[Flavio Coelho]

Thanks for the pointer,

but randstate.pyx, which allows one to choose between differents RNGs,
offers the built-in python RNG as a python object.

on line 561 it does a
import random
rand = random.Random()
return rand

What I am looking for is a way to call the C function behind
random.random directly from Cython, thus saving the overhead of a
python call. check the performance difference of the two calls from
the cython snippets below:

%cython
%timeit
from random import random
cdef int i
for i in xrange(100000000):
random()

CPU time: 11.78 s, Wall time: 13.13 s

%cython
%timeit
cdef extern from "stdlib.h":
long c_libc_random "random"()
void c_libc_srandom "srandom"(unsigned int seed)
cdef int i
for i in xrange(100000000):
c_libc_random()

CPU time: 1.85 s, Wall time: 3.11 s

Now, the above code would suffice me if the documentation of
randstate.pyx didn't mention explicitly that the RNG from libc is of
poor quality and to be avoided. Moreover, by using the RNG built into
Python, I guarantee reproducibility across platforms (I think) and
avoid external dependencies...

[Dag Sverre Seljebotn]

If you look at the sources for the random module in Python you'll see
that this essentially can't be done -- it's all declared static in the C
file creating the "_random" module.

Of course, if it is for some reason *very* important to have the Python
RNG, you could copy&paste the source code from Python. Then you have
full control. (But I think you could probably just as well use another
of Sage's PRNG).

Then there's NumPy's PRNG (yet another one) where you can draw as many
numbers at the time as you wish in an array, then access them afterwards
using Cython's fast array access. Then the Python call is only done once
regardless of N (but you need enough memory...)

randnums = numpy.random.uniform(size=100000000)

--
Dag Sverre