Čistě multiplikativní lineární kongruentní generátor

import numpy as np
import matplotlib.pyplot as plt

$I_{j+1}=aI_j\; (\text{mod } m)$
$a = 16807$, $m = 2^{31}-1$

def ran0(i):
    a=16807
    m=2147483647
    i=(a*i) % m
    return i,i/m

i=1243546789 #seminko
N=1000000 #počet simulovaných hodnot
r=np.empty(N)
for j in range(N):
    i,r[j]=ran0(i)

plt.hist(r,bins=1000,density=True)
plt.ylim(0.85,1.15)
plt.xlabel("x")
plt.ylabel("y")
plt.show()