import numpy.random as rd

def f(n,i) :
    return (n//(2**i))%2
    
def simulation(N,p,m) :
    L = [0]*m
    for _ in range(N) :
        X = rd.geometric(p) - 1
        for i in range(m) :
            L[i] += f(X,i)
    return [k/N for k in L]

for p in [1/2,1/10,1/100,1/1000] :
    print("pour p = ",p,simulation(10000,p,15))