Pseudo-random number generators

The classic and infamous RANDU. This is a simplified version of a version courtesy of ptomato on github.

"...its very name RANDU is enough to bring dismay into the eyes and stomachs of many computer scientists!"

Donald E. Knuth, The Art of Computer Programming

In [1]:
import numpy as np

class RANDU:
    "Pseudorandom number generator RANDU, a (flawed) linear congruential PRNG."
    def __init__(self, seed=None):
        self.seed(seed)

    def seed(self, s):
        self._state = hash(s) % 2147483648  # == 2^31

    def _random(self):
        self._state = (65539 * self._state) % 2147483648
        return self._state / 2147483648
    
    def random(self, size=None):
        "Return one random number or an array of them, with `size` elements."
        if size is None:
            return self._random()
        else:
            r = self._random
            return np.array([r() for i in range(size)])
In [2]:
RANDU().random(10)
Out[2]:
array([ 0.68546926,  0.9696463 ,  0.64865447,  0.16511013,  0.15277057,
        0.43063228,  0.2088585 ,  0.37746052,  0.38503659,  0.91307486])
In [3]:
%matplotlib notebook
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
In [4]:
plt.rcParams['figure.figsize'] = (10, 6)
plt.style.use('seaborn')
sns.set_context('talk', font_scale=1.4)
In [5]:
npts = 1_000

pts = RANDU().random(npts*3).reshape(npts, 3)
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=5);
plt.title("RANDU");

This is what similar output looks like for Numpy, which uses the Mersenne Twister PRNG, whose period is $2^{19937}-1 (\approx \cal{O}(10^{6001}))$ and which does not suffer from known serial correlations:

In [6]:
pts = np.random.uniform(size=(npts, 3))
fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], s=5)
plt.title("Numpy Random");

A terrible PRNG

Let's implement a simple linear congruential PRNG of the type

$$ X_{k+1} = a X_k \mod M $$

with particularly bad constants $a$, $M$ and seed $X_0$:

In [7]:
a = 3
M = 64
seed = 17
In [8]:
state = seed
numbers = []
for i in range(100):
    state = (a * state) % M
    numbers.append(state)
numbers = np.array(numbers)
print(numbers)

[51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19
 57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25
 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43
  1  3  9 27 17 51 25 11 33 35 41 59 49 19 57 43  1  3  9 27 17 51 25 11 33]
In [9]:
sorted(set(numbers))
Out[9]:
[1, 3, 9, 11, 17, 19, 25, 27, 33, 35, 41, 43, 49, 51, 57, 59]
In [10]:
len(_)
Out[10]:
16
In [11]:
np.nonzero(numbers==1)
Out[11]:
(array([11, 27, 43, 59, 75, 91]),)
In [12]:
np.diff(_)
Out[12]:
array([[16, 16, 16, 16, 16]])
In [14]:
%matplotlib inline
In [15]:
randx = numbers[:-1]
randy = numbers[1:]
plt.scatter(randx, randy, facecolors='none', edgecolors='b', linewidth=2);