mercredi 15 mars 2017

Why doesn't my program approximate pi?

For yesterday's Pi Day, Matt Harper published a video in which he approximated Pi by rolling two 120-sided dice 500 times (see the video here). Basically, for each pair of random numbers, you have to check whether they are coprime or not. Then, the formula

pi = sqrt(6/(n_coprimes/n_cofactors))

is calculated.

His result was about 3.05 which is rather close.

I wanted to see what happens when more rolls are done or when the range of random integers is increased. Interestingly, my program nearly always gave a result of 3.05 or close to it, no matter how high I set the iterations or the random range.

Here is my program. I ran it on Python 3.6 (Win64). The random number generator Python uses is supposed to be very good, so maybe I've made a mistake in my program?

import random
from math import gcd, sqrt

def pi(cp, cf):
    return sqrt(6/(cf/cp))

coprime = 0
cofactor = 0

iterations = 1000000

for i in range(iterations):
    x = random.randint(0,1000000)
    y = random.randint(0,1000000)
    if gcd(x,y) > 1: 
        cofactor += 1
    else:
        coprime += 1

print(pi(coprime, cofactor))




Aucun commentaire:

Enregistrer un commentaire