1

Define $k$ (number of faces on the fair die) and T (number of trials).

Roll it twice and record the outcomes as $x,y \in {1,2,\dots,k}$.

Let $C$ count the number of times that $GCD(x,y)=1$ after $T$ trials.

The probability of being coprime (i.e., $GCD =1)$ is

$$ \frac{C}{T} $$

However, the Probability that two random numbers are coprime is $\frac{6}{\pi^2}$

So,

$$\frac{6}{\pi^2} = \frac{C}{T}$$

Thus, an estimate of $\pi$ is

$$\pi \approx \sqrt{\frac{6T}{C}}$$

Wrote a C program to estimate $\pi$ using $10,000,000$ for the faces and trials to improve accuracy. (Approx pi value: $3.141503$)

Is this a valid method to calculate happy $\pi$ day?

vengy
  • 2,421
  • 1
    No - it will almost certainly head to the wrong limit with a $6$-sided die, no matter how many times you roll it. The statement about $\frac{6}{\pi^2}$ is the limit as the number of sides increases towards infinity (as well as the number of rolls) – Henry Mar 14 '23 at 23:25
  • 1
    There are some other monte carlo pi ideas here – Golden_Ratio Mar 14 '23 at 23:43
  • 1
    @Henry I think the OP means that the die will have 10,000,000 faces and there'll be 10,000,000 rolls. – A.J. Mar 15 '23 at 03:57
  • @A.J. Fair enough though I doubt it is accurate to more than about 3 decimal places – Henry Mar 15 '23 at 09:07
  • @Henry I think you're right about that. – A.J. Mar 15 '23 at 11:24

1 Answers1

2

This is an approximation to two distinct limits, one as the number of faces increases (with six-sided dice you might think about $\sqrt{6 \times \frac{36}{23}} \approx 3.0645$) and the other from simulation in a law of large numbers sense as the number of rolls increases.

So you will only get an approximation. C would be faster, but in R you might do something like

set.seed(2023)
faces <- 10^7
rolls <- 10^7
firstroll  <- sample(faces, rolls, replace=TRUE)
secondroll <- sample(faces, rolls, replace=TRUE)
higher <- ifelse(firstroll > secondroll, firstroll, secondroll)
lower  <- firstroll + secondroll - higher
hcf <- numeric(0)
while (length(lower) > 0){ 
   remainder <- higher %% lower
   hcf <- c(hcf, lower[remainder == 0])
   higher <- lower[remainder > 0]
   lower  <- remainder[remainder > 0]
   }
mean(hcf == 1)                   # 6/pi^2 = 0.60792710
# 0.6080553
sqrt(6 * rolls / sum(hcf == 1))  #     pi = 3.14159265
# 3.141261

which is not awful, but if you are going to use a computer then there are better and faster ways to get a good estimate of $\pi$.

Henry
  • 169,616