In this article, we will examine some interesting relationships between $\pi$ and prime numbers. This is a special article that I prepared to celebrate the Pi day, which is March 14th.
We start from a question that is simple to formulate, but difficult to answer:
What is the probability that two integers are relatively prime?
Similar to other articles that I wrote, we will start by using the violent approach, to calculate the probability brute force with the code below:
# we first define a coprime function to check whether two numbers are relatively prime
def coprime (a, b):
noncoprime = 1
smaller = a
if (a>b):
smaller = b
for i in range (2, smaller+1): #range from 2 to the smaller of a and b
if (a % i == 0 and b % i == 0):
noncoprime = 0
break #we break here because we do not really need the count, just one example would work
return noncoprime
import random
random.seed(123)
count = 0
max = 20000
for i in range (2, max+2): #test this max times (it can be any other number, I just used max for convenience)
a = random.randint(2, max)
b = random.randint(2, max)
if (coprime(a, b)):
count += 1
if (i % 1000 == 0):
print ("current count is", count, " percentage is ", count/(i-1))
print ("Final count is ", count, " percentage is ", count/max)
current count is 618 percentage is 0.6186186186186187 current count is 1216 percentage is 0.608304152076038 current count is 1826 percentage is 0.608869623207736 current count is 2450 percentage is 0.6126531632908228 current count is 3058 percentage is 0.6117223444688937 current count is 3665 percentage is 0.6109351558593099 current count is 4278 percentage is 0.6112301757393913 current count is 4858 percentage is 0.6073259157394675 current count is 5460 percentage is 0.6067340815646183 current count is 6076 percentage is 0.6076607660766077 current count is 6688 percentage is 0.608055277752523 current count is 7322 percentage is 0.6102175181265106 current count is 7920 percentage is 0.6092776367412878 current count is 8524 percentage is 0.6089006357596971 current count is 9100 percentage is 0.6067071138075871 current count is 9721 percentage is 0.6076004750296894 current count is 10326 percentage is 0.6074474969115831 current count is 10907 percentage is 0.6059781098949941 current count is 11543 percentage is 0.6075582925417127 current count is 12142 percentage is 0.6071303565178259 Final count is 12143 percentage is 0.60715
So after running the program above, it seems that the probability that any two integers are relatively prime is 0.607. We only sampled integers up to 20000 for speed reasons, but you can test other values, and it should yield similar answers.
One concern that some readers may have is that the density of prime numbers will decrease when the maximum range of random sampling becomes larger (for example, there are 4 prime numbers within 10, but only 25 prime numbers within 100). In fact, the density of prime numbers around an integer $x$ is shown to be $1/log(x)$. This does not really impact our findings above since we are only considering relatively prime, not prime by itself. Again you can try different numbers for the minimum and maximum range of integer sampling in the code above and test the program yourself.
To calculate the exact probability analytically, we can consider the following scenarios one by one, and then multiply them:
1, given two integer $a$ and $b$, the probability that $a$ can be divided by 2 is 1/2, and the probability that $b$ can be divided by 2 is also 1/2. Therefore, the probability that $a$ and $b$ can both be divided by 2 is (1/2)*(1/2)=1/4. Then the probability that $a$ and $b$ do not have 2 as common divider (or, relatively prime with respect to 2) is 1-1/4.
2, similarly, the probability that $a$ and $b$ do not have 3 as common divider is 1-1/3*1/3=1=1/9.
3, let $p_k$ be the k-th prime number ($p_1$=2, $p_2$=3, $p_3$=5, ...). Then the probability that $a$ and $b$ do not have $p_k$ as common divider is $1-1/p_k*1/p_k=1-1/p_k^2$.
4, then the probability that any two integers are relatively prime is the multiplication of all these probabilities as
$P=\displaystyle\prod_{k=1}^{\infty} (1-1/p_k^2)$
The only difference with the Python code earlier is that we put a maximum boundary of sampling the two numbers ($a$ and $b$), so that the computation can be completed in short period of time. Here there is no upper boundary how much $a$ and $b$ can be.
5, Because Euler Product formula says that
$\displaystyle\prod_{k=1}^{\infty} \frac{1}{1-1/p_k^s} = \displaystyle\sum_{k=1}^{\infty} \frac{1}{n^s} = \zeta(s)$
We have $P = \frac{1}{\zeta(2)} = \frac{6}{\pi^2}$
The Zeta function $\zeta$ is defined in Riemann Hypothesis, which is arguably the most important problem in mathematics and human history. If it is wrong, then the entire infrastructure of modern mathematics will crash. Every week, somebody will publish a paper claiming that they proved the hypothesis, yet none has gained widespread recognition or approval yet. So the hypothesis is better not to be wrong.
import math
x = 6/math.pi**2
x
0.6079271018540267
So the theoretical answer (0.608) is quite close to our data simulation above.
Since the zeta function is also available in SciPy, we can calculate it and validate its values:
from scipy.special import zeta
1/zeta(2)
0.6079271018540267
In the situation above we considered two integers that are relatively prime. Now we extend the problem to "what is the probability that m integers that are relatively prime". The question is a little tricky, because we only require that the k numbers do not share a common divider, yet it is totally fine that k-1 numbers share a common divider.
Intuitively, since the answer to the qestion above is $1/\zeta(2)$, the answer to this question is likely $1/\zeta(m)$. This is correct.
The same reasoning can be used here. Say for example if we have 3 numbers, and we want to calculate the probability that they are relatively prime with respect to 2 (that is, they do not share a common divider of 2). It is possible that 2 of the 3 numbers are even and the third one is odd. Since the probability that all 3 numbers are even is 1/8, then the probability of the reverse is 1-1/8.
So the final $P=\displaystyle\prod_{k=1}^{\infty} (1-1/p_k^m) = 1/\zeta(m)$
Based on the webpage, the numbers can be reprented as a function of $\pi$ for $m=2$ and $m=4$, but not $m=3$ and other values.
We can use SciPy to calculate the values when m=2, 3, 4. Note that m=2 is a special case of the previous question.
from scipy.special import zeta
1/zeta(2), 1/zeta(3), 1/zeta(4)
(0.6079271018540267, 0.8319073725807075, 0.9239384029215902)
This is a slightly different problem than the above one, because we require "mutually relatively prime". This means no two numbers in the $n$ numbers can share a common divider, yet in the previous question it is okay for $n-1$ numbers to share a common divider.
To solve this problem, we still start our thinking from 2 as a divider. Say if we have 3 numbers, to ensure that they are mutually relatively prime, then at most one number can be even, and the other two numbers must be odd. Or all the 3 numbers must be odd. Similarly, if we have 4 numbers, at most one can be even and all other three numbers must be odd.
The general formula therefore seems to be $(1-1/2)^{n}+(1-1/2)^{n-1} n/2$. The first part calculates the probability that all numbers are odd, and the second one calculate the probability that one number is even, and the other ones are odd.
If we then consider the prime number 3 as a divider, we can get similar conclusion that the probability is $(1-1/3)^{n}+(1-1/3)^{n-1} n/3$.
Now consider all prime numbers $p \in \mathbb{p}$. The probability is
$P=\displaystyle\prod_{k=1}^{\infty} ((1-1/p_k)^n+(1-1/p_k)^{n-1}n/p_k)$
It cannot be represented by a zeta function.
The code below tests n=2, n=3 and so on.
prob = 1
n = 2
for i in range(2, 1000): #cycle up to 1000
for j in range(2, i):
if i % j == 0:
break
else:
sum = (1-1/i)**n + (1-1/i)**(n-1)*n/i
prob *= sum
print ("total P = ", prob)
total P = 0.6080043073061262
prob = 1
n = 3
for i in range(2, 1000): #cycle up to 1000
for j in range(2, i):
if i % j == 0:
break
else:
sum = (1-1/i)**n + (1-1/i)**(n-1)*n/i
prob *= sum
print ("total P = ", prob)
total P = 0.28685665267717797
prob = 1
n = 4
for i in range(2, 1000): #cycle up to 1000
for j in range(2, i):
if i % j == 0:
break
else:
sum = (1-1/i)**n + (1-1/i)**(n-1)*n/i
prob *= sum
print ("total P = ", prob)
total P = 0.11497155024232823