证明任意两个自然数互质的概率为6/pi^2

证明任意两个自然数互质的概率为$\frac{6}{{\pi }^2}$

​ 设事件 A : { ′ m , n 互 质 ′ } A:\{'m,n互质'\} A:{′m,n互质′}，事件 B : { ′ a , b 互 质 ′ } B:\{'a,b互质'\} B:{′a,b互质′}，事件 C k : { ′ k 是 a , b 的 公 因 数 ′ } C_k:\{'k是a,b的公因数'\} Ck​:{′k是a,b的公因数′}， D k : { ′ k 是 a , b 的 最 大 公 因 数 ′ } D_k:\{'k是a,b的最大公因数'\} Dk​:{′k是a,b的最大公因数′}。

​ 若设任意两个自然数分别为$a,b$，那么它们互质的概率我们设为$p$,但是如果这两个自然数同时是$k$的倍数的概率为:
P ( C k ) = p { a 是 k 的 倍 数 } p { b 是 k 的 倍 数 } = p { ∀ a ∈ N , a k 余 数 是 0 } p { ∀ b ∈ N , b k 余 数 是 0 } = N / k N N / k N = 1 k 2 P(C_k)=p\{a是k的倍数\}p\{b是k的倍数\}=p\{\forall a \in N,\frac{a}{k}余数是0\}p\{\forall b \in N,\frac{b}{k}余数是0\} = \frac{N/k}{N}\frac{N/k}{N} = \frac{1}{k^2} P(Ck​)=p{a是k的倍数}p{b是k的倍数}=p{∀a∈N,ka​余数是0}p{∀b∈N,kb​余数是0}=NN/k​NN/k​=k21​
​ 那么此时我们设:
$$a=km,b=kn$$
​ 那么我们可以得到以下充要条件:
{ " m , n 互 质 " } ⇐ ⇒ { " k 是 a , b 的 最 大 公 因 数 ∣ k 是 a , b 的 公 因 数 " } \{"m,n互质"\}\Leftarrow \Rightarrow\{"k是a,b的最大公因数|k是a,b的公因数"\} {"m,n互质"}⇐⇒{"k是a,b的最大公因数∣k是a,b的公因数"}
​ 这也就是说
$$P(D_k|C_k)=P(A)$$
​ 而此时$m$和$n$同样为任意两个自然数，他们互质的概率也为$p$。
$$P(A)=P(B)=p$$
​ 而$P(C_k)=\frac{1}{k^2}$。

​ 而由贝叶斯定理我们可以得到:
$$P(D_k)=P(D_k|C_k)P(C_k)$$
​ 而$D_k|C_k{:}^{\prime }在k是a,b的公因数的情况下，k是a,b的最大公因{数}^{\prime }$。
$$P(D_k|C_k)=P(A)$$
​ 这里，我们需要明白一个事实，那就是，对于$k=1,2,...\infty$，总有一个数$k$是$a,b$的最大公因数，这件事情的概率是$1$。也就是说是确定事件。而$\forall i,j\in N,D_i\cap D_j=\varphi$。所以：
$$\sum _{k=1}^{\infty }P(D_k)=1$$
​ 带入由贝叶斯定理导入的公式可得到:
$$\sum _{k=1}^{\infty }\frac{p}{k^2}=1$$
​ 因此化简之后我们可以得到:
$$p=\frac{1}{{\sum }_{k=1}^{\infty }\frac{1}{k^2}}=\frac{6}{{\pi }^2}$$
​ 接下来我们用$Python$借助上面的原理估计一下$\pi$的大小来验证结果的正确性:

import math
import random
def isInterSushu(sx,sy):
    num = min(sx,sy)
    res = 2
    for i in range(num,1,-1):
        if (sx%i == 0)&(sy%i == 0):
            res = 0
            break
        else:
            res = 1
    return res

def calPi(numbers):
    sum = 0
    error = 0
    random.seed(10)
    for i in range(numbers):
        sx_value = int(numbers*random.random())
        sy_value = int(numbers*random.random())
        if isInterSushu(sx_value,sy_value):
            sum += 1
            error += 1**2
    Prop = sum/numbers
    sigma = math.sqrt(1/numbers*error-(sum/numbers)**2)
    epsilon = sigma*1.96/math.sqrt(numbers)
    pi = math.sqrt(6 / (Prop))
    error_min = math.sqrt(6/(Prop+epsilon))
    error_max = math.sqrt(6/(Prop-epsilon))
    print("pi在{:}的试验次数下pi的估计值为{:.3f},由0.95的概率认为真值在[{:.5f},{:.5f}]之间。".format(numbers,pi,error_min,error_max))

calPi(10000)

最后的结果是:

pi在10000的试验次数下pi的估计值为3.131,由0.95的概率认为真值在[3.10672,3.15558]之间。

Process finished with exit code 0

可以发现结论基本上是正确的。