Let's get random, extending Python's base random class

Question

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.

Source

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 and MINSTD 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 them randu and minstd, 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)

Answer 1

it is slightly cumbersome to extend the class, as the C module steals the first argument. (Somehow).

I really don't understand this. I had no trouble at all with adapting your code to use a traditional override and super()-constructor invocation. Delete your __new__, give your constructor parameters better names and typehints, and don't use the seed in your own constructor: super() will pass you back your seed to the seed() function.

Your seed function implementation is non-conformant to the base. For it to be conformant, it should be using the system time as a default, and accepting a version number.

Don't plt.show() in your class methods; just populate an existing axis object so that multiple plots can be manipulated as desired by the caller.

my concrete implementations RANDU and MINSTD are classes, but they are implemented as functions.

So.. make them classes? Basic inheritance works just fine.

Your MultipleRecursiveGenerator didn't work for me at all. At the point where you attempt a matrix multiplication @, your array dimensions are incompatible. Also your vectorize doesn't work because it's attempting to pass two arguments to random. I have not fixed this.

Suggested

import random
from time import time

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:
         modulus:    0 <  m
         multiplier: 0 <  a < m
         increment:  0 <= c < m
         seed:       0 <= X0 < m
    """

    def __init__(self, modulus: int, multiplier: int, increment: int = 0, seed: int = 0) -> None:
        super().__init__(x=seed)
        self.t = 0
        self.m = modulus
        self.a = multiplier
        self.c = increment

    def seed(self, a: int = None, version: int = 2) -> None:
        if a is None:
            a = time()
        self.X0 = a
        self.X = a

    def nextstate(self) -> int:
        self.t += 1
        self.X = int((self.a * self.X + self.c) % self.m)
        return self.X

    def random(self) -> float:
        self.nextstate()
        return self.X / (self.m + 1)

    def image(self, ax: plt.Axes, shape=(400, 400)) -> None:
        """Produces a simple grayscale image from rand()"""
        image = [[self.random() for _ in range(shape[0])] for _ in range(shape[1])]
        ax.imshow(image, cmap="gray", interpolation="none")
        ax.set_title(f'a={self.a}')
        ax.axis('off')


class RANDU(Lcg):
    def __init__(self, seed: int = 0) -> None:
        super().__init__(multiplier=65539, modulus=2**31, seed=seed)


class MINSTD(Lcg):
    """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),
    """
    def __init__(self, seed: int = 0) -> None:
        super().__init__(multiplier=7 ** 5, modulus=2**31 - 1, seed=seed)


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, modulus: int, multiplier: int, seed: int = 0, offset: int = 0) -> None:
        """Create a new PRNG with seed X0."""
        super().__init__(x=seed)

        self.a = np.array(modulus)
        self.m = multiplier
        self.c = np.array(offset)
        self.t = 0

    def seed(self, a: int = None, version: int = 2) -> None:
        if a is None:
            a = time()
        self.X0 = np.array(a)
        self.X = np.array(a)

    def nextstate(self) -> int:
        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) -> float:
        x = self.nextstate()
        return x / (self.m + 1)

    def image(self, ax: plt.Axes, shape=(400, 400)) -> None:
        """Produces a simple grayscale image from rand()"""
        image = [[self.random() for _ in range(shape[0])] for _ in range(shape[1])]
        ax.imshow(image, cmap="gray", interpolation="none")
        ax.set_title('MRG')
        ax.axis('off')


if __name__ == "__main__":
    randu = RANDU(seed=1)
    minstd = MINSTD(seed=1)
    mrg = MultipleRecursiveGenerator(multiplier=65539, modulus=2**31, offset=3)

    fig, (ax_randu, ax_minstd,
          # ax_mrg
    ) = plt.subplots(nrows=1, ncols=2)
    randu.image(ax_randu)
    minstd.image(ax_minstd)
    # mrg.image(ax_mrg)
    plt.show()

    print([randu.nextstate() for _ in range(5)])