Welcome, guest | Sign In | My Account | Store | Cart

Long period random number generator (Python recipe) by Raymond Hettinger
ActiveState Code (http://code.activestate.com/recipes/576707/)

Implements a complementary-multiply-with-carry psuedo-random-number-generator. Period is 3636507990 * 2 ** 43487 (approximately 10 ** 13101).

Python, 49 lines
Download
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import random

class CMWC(random.Random):
    'Long period random number generator: Complementary Multiply with Carry'
    # http://en.wikipedia.org/wiki/Multiply-with-carry

    a = 3636507990
    logb = 32
    b = 2 ** logb
    r = 1359

    def _gen_word(self):
        i = self.i
        xc, self.c = divmod(self.a * self.Q[i] + self.c, self.b)
        x = self.Q[i] = self.b - 1 - xc
        self.i = 0 if i + 1 == self.r else i + 1
        return x

    def getrandbits(self, k):
        while self.bits < k:
            self.f = (self.f << self.logb) | self._gen_word()
            self.bits += self.logb
        x = self.f & ((1 << k) - 1)
        self.f >>= k;  self.bits -= k
        return x

    def random(self, RECIP_BPF=random.RECIP_BPF, BPF=random.BPF):
        return self.getrandbits(BPF) * RECIP_BPF

    def seed(self, seed=None):
        seeder = random.Random(seed)
        Q = [seeder.randrange(0x100000000) for i in range(self.r)]
        c = seeder.randrange(0x100000000)
        self.setstate((0, 0, 0, c, Q))

    def getstate(self):
        return self.f, self.bits, self.i, self.c, tuple(self.Q)

    def setstate(self, (f, bits, i, c, Q)):
        self.f, self.bits, self.i, self.c, self.Q = f, bits, i, c, list(Q)


if __name__ == '__main__':
    prng = CMWC(134123413541344)
    for i in range(20):
        print prng.random()
    print
    for i in range(20):
        print normalvariate(mu=5.0, sigma=2.2)

Here's a fast generator version of random() optimized for parameter sets where b is an exact power-of-two. It uses the Mersenne Twister generator for seeding the Q array and the initial carry value.

def cmwc_random(seed=None, a=3636507990, b=2**32, logb=32, r=1359):
    seeder = random.Random(seed)
    Q = [seeder.randrange(b) for i in range(r)]
    c = seeder.randrange(b)
    f = bits = 0
    for i in itertools.cycle(range(r)):
        t = a * Q[i] + c
        c = t & (b - 1)
        x = Q[i] = b - 1 - (t >> logb)
        f = (f << logb) | x;  bits += logb
        if bits >= 53:            
            yield (f & (2 ** 53 - 1)) * (2 ** -53)
            f >>= 53;  bits -= 53
Tags: generator, number, random

5 comments

a 13 years, 6 months ago  # | flag

echo "1/(32768/($RANDOM+1))" | bc -l

Craig McQueen 13 years, 5 months ago  # | flag

getrandbits() appears to insert new bits at the least-significant-bits end, and also take bits from the least-significant-bits end. That seems a little unusual. Wouldn't it be better to shift bits through like a queue, i.e. insert at one end, remove from the other end? e.g.

    def getrandbits(self, k):
        while self.bits < k:
            self.f = (self.f << self.logb) | self._gen_word()
            self.bits += self.logb
        self.bits -= k
        x = self.f >> self.bits
        self.f &= ((1 << self.bits) - 1)
        return x
Craig McQueen 13 years, 5 months ago  # | flag

getstate() and setstate() seem to have parameters f and bits in the opposite order. Is this right?

Stamatis Karlos 9 years, 6 months ago  # | flag

Could you provide me the code for stoping the cmwc_random generator faster?

Ron Charlton 8 years ago  # | flag

Thanks for creating and sharing a clear and concise implementation of a CMWC PRNG.

I noticed a problem with integer return values from getrandbits while I was testing output with John Walker's ent.exe program: getrandbits produced a 13,107,200 byte file that had about 2.6 bits of entropy per byte. That is obviously too low. It is repeatable. Most PRNGs will have close to 8 bits of entropy per byte. function cmwc_random has the same problem (I modified it slightly by not multiplying by 2 ** -53).

I am very glad you gave the Wikipedia source for CMWC's algorithm. I found that in _gen_word the assigned-to variables in the divmod() statement are swapped. Reversing their order fixed the problem. I now get a 13,107,200 byte file with 7.999987 bits per byte.

The corrected _gen_word follows:

def _gen_word(self):
    i = self.i
    self.c, xc = divmod(self.a * self.Q[i] + self.c, self.b)
    x = self.Q[i] = (self.b - 1 - xc)
    self.i = (0 if i + 1 == self.r else i + 1)
    return x

The corrected section of cmwc_random code:

for i in itertools.cycle(list(range(r))):
    t = a * Q[i] + c
    c = b - 1 - (t >> logb)
    x = Q[i] = t & (b - 1)
    f = (f << logb) | x;  bits += logb

Thanks a billion for your code, Raymond. The corrected code is quite fast.