Why doesn't my program come close to pi?

Question

Last day, Pi Day Matt Harper posted a video in which he approximated the Pi by rolling two 120-sided dice 500 times ( see video here ). Basically, for each pair of random numbers you need to check if they are mutual or not. Then the formula

pi = sqrt(6/(n_coprimes/n_cofactors)) # EDIT: Wrong premise. Misremembered the formula. 

.

His result was around 3.05, which is pretty close.

I wanted to see what happens when more rolls are made or when the range of random integers is increased. Interestingly, my program almost always gave a result of 3.05 or is close to it, regardless of how high the setting of iterations or the random range is.

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

 import random from math import gcd, sqrt def pi(cp, cf): return sqrt(6/(cf/cp)) # EDIT: Second error - switched numerator/denominator... 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)) 

Answer 1

I have not watched the video, but your formula is wrong.

The probability that two ints, randomly selected from 1 to N, are coprime, tends to 6 / pi ^ 2, when N tends to infinity. This is cp / (cf + cp), not cp / cf.

Replace pi as follows:

 def pi(cp, cf): fcp = cp / float(cp + cf) return sqrt(6/fcp) 

Gives 3.14263472915 when I run it on my machine.