Implements a complementary-multiply-with-carry psuedo-random-number-generator. Period is 3636507990 * 2 ** 43487 (approximately 10 ** 13101).
Python, 49 lines
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
echo "1/(32768/($RANDOM+1))" | bc -l
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.
setstate()seem to have parameters
bitsin the opposite order. Is this right?
Could you provide me the code for stoping the cmwc_random generator faster?
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:
The corrected section of cmwc_random code:
Thanks a billion for your code, Raymond. The corrected code is quite fast.