rng.mlcg = function(n, a1=20, a0=0, m=53, s0=21)
{
    s = numeric(n)
    s[1] = s0
    for (i in 1:(n-1))
        s[i+1] = (a1*s[i] + a0) %% m
        
    return (s)
}