Random number cycles

Random number generation is often done by applying a function recursively to a seed number. The function needs to be bounded to a certain domain. Obtaining a function that is bounded is easily done by applying the modulo operation on the base function. The function must be deterministic. For example the following function

\[ f(x) = (x + 1) % 3 \]

will supply us with the sequence \(0,1,2\) looping indefinitely if we take seed value \(0\). If we want to detect the cycle in a sequence we can use the equation \(X_n = X_{2n}\) for some \(n\). Proofing this relation is left as an exercise. Now we can create some code to detect cycles for a function:

def period_fun(f, x0):
    """return length and start of recursively
    applied function 'f' with start value 'x'"""
    x1 = f(x0)
    x2 = f(f(x0))
    n = 1
    while x1 != x2:
        x1 = f(x1)
        x2 = f(f(x2))
        n += 1

    xl = f(x1)
    l = 1
    while xl != x2:
        xl = f(xl)
        l += 1

    xm = x0
    m = 0
    while xm != x1:
        xm = f(xm)
        x1 = f(x1)
        m += 1

    return l, m, xm

In this code the l indicates the length of a period, and m is the index of the first digit of the cycle. To test if this works, let’s define a recursive function. A famous random number generating function by Von Neumann is the Middle-square method. It is here implemented for 2-digit numbers.

def middle_square2(d):
    return int(str(d * d).zfill(4)[1:3])

Here is a image from the wikipedia, showing the possible cycles:

Let’s see if our results compare with a simple test.

from numpy import zeros
n = 100
result = zeros((n, 4))
for i in range(100):
    l, m, e = period_fun(middle_square2, i)
    result[i] = [i, l, m, e]

Now visualize our results with matplotlib, the default Python workflow.

import maptlotlib.pyplot as plt
fig, ax = plt.subplots(1, 3)
ax[0].hist(result[:,1])
ax[0].set_title("cycle length")
ax[1].hist(result[:,2])
ax[1].set_title("init sequence length")
ax[2].hist(result[:,3])
ax[2].set_title("end number")
plt.show()

It looks like it works, you can also rewrite the period_fun to take sequences instead of a function. Try to implement it yourself. Also test some other functions to see what happens.