def square(x): return x * x
# 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 / monte_carlo(trials, dirichlet_test))

def dirichlet_test():
    return gcd(rand(), rand()) == 1

def monte_carlo(trials, experiment):
    def iter(trials_remaining, trials_passed):
        return (trials_passed / trials
                if trials_remaining == 0
                else iter(trials_remaining - 1, trials_passed + 1)
                if experiment()
                else iter(trials_remaining - 1, trials_passed))
    return iter(trials, 0)
def random_in_range(low, high):
    range = high - low
    return low + math_random() * range

def estimate_integral(pred, x1, x2, y1, y2, trials):
    area_rect = (x2 - x1) * (y2 - y1)
    return monte_carlo(trials,
                       lambda: pred(random_in_range(x1, x2),
                                    random_in_range(y1, y2))) * area_rect

def inside_unit_circle(x, y):
    return square(x) + square(y) <= 1

print(estimate_integral(inside_unit_circle, -1, 1, -1, 1, 50000))
