# A very simple rand_update function computes a number
# from 0 (inclusive) to 200560490131 (a large prime)
# from a value x by multiplying it with a constant a,
# adding a constant b, and then taking the remainder
# of dividing it by the large prime. We used it here
# for illustration only, and do not claim any
# statistical properties.
m = 200560490131
a = 1103515245
b = 12345

def rand_update(x):
    return (a * x + b) % m
random_init = 123456789
def make_rand():
    x = random_init
    def rand():
        nonlocal x
        x = rand_update(x)
        return x
    return rand

rand = make_rand()
def gcd(a, b):
    return a if b == 0 else gcd(b, a % b)
def estimate_pi(trials):
    return math_sqrt(6 / random_gcd_test(trials, random_init))

def random_gcd_test(trials, initial_x):
    def iter(trials_remaining, trials_passed, x):
        x1 = rand_update(x)
        x2 = rand_update(x1)
        return (trials_passed / trials
                if trials_remaining == 0
                else iter(trials_remaining - 1, trials_passed + 1, x2)
                if gcd(x1, x2) == 1
                else iter(trials_remaining - 1, trials_passed, x2))
    return iter(trials, 0, initial_x)

print(estimate_pi(10000))

# expected: 3.1699476461053324
