metapostRNG.py

#!/bin/env/py
"""This is a quick Pythonification of the Pascal code in pdftex.web

I know I should use a Python generator rather than variables with
global scope...

Anyway, brief testing showed it does match exactly with actual
\pdfuniformedeviate/\pdfsetrandomseed, and that's what I wanted
in case I attempt some statistics on this PRNG.

May 2, 2018. JF B.
"""

fraction_one = 2**28
fraction_half = 2**27
randoms = [0] * 55
j_random = 0

def next_random():
    global j_random
    if j_random==0:
        new_randoms()
    else:
        j_random -=1


def new_randoms():
    global j_random
    global randoms
    for k in range(24):
        randoms[k] = (randoms[k] - randoms[k+31]) % fraction_one
    for k in range(24, 55):
        randoms[k] = (randoms[k] - randoms[k-24]) % fraction_one
    j_random = 54


def init_randoms(seed):
    global randoms
    j = abs(seed)
    while j >= fraction_one:
        j //= 2
    k = 1
    for i in range(55):
        jj = k
        k = (j - k) % fraction_one
        j = jj
        randoms[(i * 21) % 55] = j
    new_randoms()
    new_randoms()
    new_randoms()


def unif_rand(x):
    global j_random
    global randoms
    next_random()
    # here we must round, but Python // truncates, so we are in exact
    # opposite as when we want to truncate in a \numexpr which rounds ;-)
    y = ( abs(x) * randoms[j_random] + fraction_half ) // fraction_one
    if y == abs(x):
        return 0
    else:
        if x > 0:
            return y
        else:
            return -y

# all the above is quite direct Pythonification of code from pdftex.web
# it definitely should be using a generator ("yield") to be more Pythonesque...

# now let's start using it

def test(seed, N, reps):
    """Create a stream of random integers which we can compare with similar
    one from tex code using \pdfuniformdeviate N
    
    I did it with 1000 draws and seeds 12345 and 12346 and match was perfect.
    """
    init_randoms(seed)
    with open("randomints_py_" + str(seed) + ".txt", "w") as f:
        for k in range(reps):
            f.write(str(unif_rand(N)) + '\n')


# compute a chi2 and using scipy.stats a p-value

from scipy.stats import chisquare
        
def test_uniform(seed, N, reps):
    init_randoms(seed)
    occ = { i:0 for i in range(N) }
    for k in range(reps):
        n = unif_rand(N)
        occ[n] += 1

    E = reps/N  # expected number of occurences for each value

    # we compute chi2 also ourselve for a double-check
    # chi2 = sum((occ[i] - E)**2/E for i in range(N))

    # it is more precise using (\sum O_i^2)/E  - reps which we can do
    # as Python integers are not bounded. (add this to xint?)

    chi2 = N * sum(occ[i]**2 for i in range(N))/reps - reps

    print(chi2)  # for checking

    # let's hand over to more powerful tools (which I need to learn next)
    # so we get the "p-value". This call retuns chi2 and "p-value"

    return chisquare(list(occ.values()))