#!/usr/bin/env python3
from mpi4py import MPI
import numpy as np
import math

comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()

total = int(1e7)
n = total // size

f = lambda x,y: math.exp(-0.5*(x**2 + y**2))
xy = np.random.uniform(size=n*2, low=-1, high=1)
z = np.random.uniform(size=n, low=0, high=1)

lower = 0
for i in range(n):
        if z[i] < f(xy[i], xy[n+i]):
                lower += 1

l = np.array(lower, dtype=np.int)
t = np.zeros(size, dtype=np.int)
#comm.Gather(l, t, root=0)
comm.Reduce(l, t, op=MPI.SUM, root=0)

if rank == 0:
#       print(2*2*(sum(t)/total))
        print(2*2*(t[0]/total))