[Numpy-discussion] random number genration

Warren Weckesser warren.weckesser at enthought.com
Tue Mar 29 13:10:34 EDT 2011


On Tue, Mar 29, 2011 at 11:59 AM, Sturla Molden <sturla at molden.no> wrote:

> Den 29.03.2011 16:49, skrev Sturla Molden:
> >
> > This will not work. A boolean array is not compactly stored, but an
> > array of bytes. Only the first bit 0 is 1 with probability p, bits 1 to
> > 7 bits are 1 with probability 0. We thus have to do this 8 times for
> > each byte, shift left by range(8), and combine them with binary or.
> > Also the main use of random bits is crypto, which requires the use of
> > /dev/urandom or CrypGenRandom instead of Mersenne Twister
> (np.random.rand).
>
>
> Here's a cleaned-up one for those who might be interested :-)
>



How about adding it to http://www.scipy.org/Cookbook?

Warren



>
> Sturla
>
>
>
> import numpy as np
> import os
>
>
> def randombits(n, p, intention='numerical'):
>     """
>     Returns an array with packed bits drawn from n Bernoulli
>     trials with successrate p.
>     """
>     assert (intention in ('numerical','crypto'))
>     # number of bytes
>     m = (n >> 3) + 1 if n % 8 else n >> 3
>     if intention == 'numerical':
>         # Mersenne Twister
>         rflt = np.random.rand(m*8)
>     else:
>         # /dev/urandom on Linux, Apple, et al.,
>         # CryptGenRandom on Windows
>         rflt = np.frombuffer(os.urandom(m*64),dtype=np.uint64)
>         rflt = rflt / float(2**64)
>     b = (rflt.reshape((m,8))<p)
>     # pack the bits
>     b = (b<<range(8)).sum(axis=1).astype(np.uint8)
>     # zero the trailing m*8 - n bits
>     b[-1] &= (0xFF >> (m*8 - n))
>     return b
>
>
> def bitcount(a, bytewise=False):
>     """
>     Count the number of set bits in an array of np.uint8.
>     """
>     assert(a.dtype == np.uint8)
>     c = a[:,np.newaxis].repeat(8, axis=1) >> range(8) & 0x01
>     return (c.sum(axis=1) if bytewise else c.sum())
>
>
> if __name__ == '__main__':
>
>     b = randombits(int(1e6), .1, intent='numerical')
>     print bitcount(b) # should be close to 100000
>     b = randombits(int(1e6), .1, intent='crypto')
>     print bitcount(b) # should be close to 100000
>
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20110329/7726730f/attachment-0001.html>


More information about the NumPy-Discussion mailing list