So Python's base random class uses the Mersenne Twister as its base for seeded random number generation. If available, it can use /urandom/
, but that is beside the point.
There have been several claims that this generator is not ideal, and that other generators should be used instead.
Extending the base class of Random
's library, should, in theory be simple. As stated in the source code
class Random(_random.Random):
"""Random number generator base class used by bound module functions.
Used to instantiate instances of Random to get generators that don't
share state.
Class Random can also be subclassed if you want to use a different basic
generator of your own devising: in that case, override the following
methods: random(), seed(), getstate(), and setstate().
Optionally, implement a getrandbits() method so that randrange()
can cover arbitrarily large ranges.
"""
However, things are not as simple as they seem. As can be seen from this SO post.
So it is slightly cumbersome to extend the class, as the C module steals the first argument. (Somehow).
I implemented the most basic LCG (Linear congruential PRNG) as follows
lcg.py
import random
import numpy as np
import matplotlib.pyplot as plt
class Lcg(random.Random):
"""A simple linear congruential Pseudo-Random Number Generator.
Based on the following recurrence relation
X[n + 1] = ( a X[n] + c ) mod m
Args:
m: 0 < m - the modulus
a: 0 < a < m - the multiplier
c: 0 <= c < m - the increment
X0: 0 <= X0 < m - the seed or start value
"""
def __new__(cls, m, a, c=0, X0=0):
cls.m = m
cls.a = a
cls.c = c
cls.X0 = X0
return super().__new__(cls)
def __init__(self, m, a, c=0, X0=0):
"""Create a new PRNG with seed X0."""
super().__init__(x=None)
self.t = 0
self.m = m
self.a = a
self.c = c
self.X0 = X0
self.X = X0
def seed(self, seed=None):
if seed is not None:
self.X0 = seed
self.X = self.X0
def nextstate(self) -> int:
self.t += 1
self.X = int((self.a * self.X + self.c) % self.m)
return self.X
def random(self):
self.nextstate()
return self.X / (self.m + 1)
def image(self, shape=(400, 400)):
"""Produces a simple grayscale image from rand()"""
image = [[self.random() for _ in range(shape[0])] for _ in range(shape[1])]
plt.figure(figsize=(8, 8))
plt.imshow(image, cmap="gray", interpolation="none")
plt.axis("off")
plt.show()
def RANDU(seed=0):
return Lcg(a=65539, m=2 ** 31, X0=seed)
def MINSTD(seed=0):
"""A particular version of the lehmer random number generator
X[k + 1] = a X[x] mod m
In 1988, Park and Miller[3] suggested a Lehmer RNG with particular
parameters m = 2^31 − 1 = 2,147,483,647 (a Mersenne prime M31) and
a = 75 = 16,807 (a primitive root modulo M31),
"""
return Lcg(a=7 ** 5, m=2 ** 31 - 1, X0=seed)
if __name__ == "__main__":
randu = RANDU(1)
RANDU(1).image()
# MINSTD(1).image()
# print([randu.nextstate() for _ in range(5)])
Similarly MRG (Multiple LCG) was implemented as
mrg.py
import math
import random
import matplotlib.pyplot as plt
import numpy as np
from lcg import Lcg
class MultipleRecursiveGenerator(random.Random):
"""A multiple recursive (Pseudo-Random Number) Generator.
A multiple recursive generator (MRG) of order k is defined by the
linear recurrence:
X[n] = ( a1 X[n-1] + ... + ak X[n-k] + c ) mod m
"""
def __new__(cls, a, m, X0, c=0, x=None):
cls.a = a
cls.m = m
cls.X0 = X0
cls.c = c
cls.x = x
return super().__new__(cls)
def __init__(self, a, m, X0, c=0, x=None):
"""Create a new PRNG with seed X0."""
self.a = np.array(a)
self.m = m
self.c = np.array(c)
self.t = 0
self.X0 = np.array(X0)
self.X = np.array(X0)
super().__init__(x)
def nextstate(self):
self.t += 1
x = self.X
next_x = (self.a @ self.X + self.c) % self.m # [a b] @ [c d] = ab+bd
self.X = np.roll(self.X, 1) # roll([1 2 3], 1) -> [3 1 2]
self.X[0] = next_x
return x
def random(self):
x = self.nextstate()
return x / (self.m + 1)
def image(self, shape=(400, 400)):
"""Produces a simple grayscale image from rand()"""
samples = np.fromfunction(np.vectorize(self.random), shape=shape, dtype=int)
image = samples(shape)
plt.figure(figsize=(8, 8))
plt.imshow(image, cmap="gray", interpolation="none")
plt.axis("off")
plt.show()
The code is working as intended and is an extension of the random class. E.g. the extension above uses either the LCG or MRG to manage the state, instead of the default Mersenne Twister (Note that LCG / MRG are much, much weaker and only implemented as a proof of concept).
I could have (and have) implemented much more sophisticated methods, but I want to get feedback on my implementation before I proceed.
- That random steals my arguments to use as a seed, hampers my ability to write clean code. I feel like I have to initiate every class with a
__new__
- I tried to create a buffer class to store common methods and override
random.Random
, but was unsuccessful (same reasons as above) - Similarly my concrete implementations
RANDU
andMINSTD
are classes, but they are implemented as functions. I wanted to implement them as functions, but as before I rant into issues with random stealing my arguments, and non hashable entries. Is it fine having functions with capitalized names as this is their official names? I guess I could have named themrandu
andminstd
, but then it would be confusing why they are class generators.
Anyway that is just my rambling, any improvements on how to rip out the standard Mersenne Prime generator from random and implement our own, is more than welcome =)
EDIT: For reference here is a non working result, due to values being passed directly to _random
. The problem goes away if the lists are changed into tupples, but for the more sophisticated generators this is not a solution
import numpy as np
import random
class MultipleRecursiveGenerator(random.Random):
"""A multiple recursive (Pseudo-Random Number) Generator.
A multiple recursive generator (MRG) of order k is defined by the
linear recurrence:
X[n] = ( a1 X[n-1] + ... + ak X[n-k] c ) mod m
"""
def __init__(self, a, m, X0, c=0):
"""Create a new PRNG with seed X0."""
self.a = np.array(a)
self.m = m
self.c = np.array(c)
self.t = 0
self.X0 = np.array(X0)
self.X = np.array(X0)
super(random.Random, self).__init__()
def nextstate(self):
self.t += 1
x = self.X
next_x = (self.a @ self.X + self.c) % self.m
# roll([1,2,3], 1) -> [3, 1, 2]
self.X = np.roll(self.X, 1)
self.X[0] = next_x
return x
def random(self):
x = self.nextstate()
return x / (self.m + 1)
if __name__ == "__main__":
seed = 123456
m1 = (1 << 32) - 209
a = np.array([0, 1403580, -810728])
X0 = [seed, 0, 0]
mrg1 = MultipleRecursiveGenerator(a, m1, X0)
m2 = (1 << 32) - 22853
b = np.array([527612, 0, -1370589])
Y0 = [seed, 0, 0]
mrg2 = MultipleRecursiveGenerator(b, m2, Y0)