为什么我的程序无法逼近圆周率？

Question

针对昨天的“圆周率日”，Matt Harper发布了一段视频，其中他通过掷500次两个120面骰子来近似计算圆周率(点击此处查看视频)。基本上，对于每组随机数，您需要检查它们是否互质。然后，使用以下公式进行计算：

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

被计算出来。

他的结果约为3.05，非常接近。

我想看看当更多的掷骰子或随机整数范围增加时会发生什么。有趣的是，无论我设置迭代次数或随机范���有多高，我的程序几乎总是给出一个接近3.05的结果。

这是我的程序。我在Python 3.6（Win64）上运行它。Python使用的随机数生成器应该是非常好的，所以也许我在程序中犯了错误？

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))

- Tim Pietzcker

---

昨天看了这个视频，甚至没有想过自己尝试。谢谢你的创意！至于为什么它不能提高准确性，我不知道。我最好的猜测是生成器的限制，但我对Python的生成器知之甚少。 - Carcigenicate

Answer 1

我还没有看过这个视频，但是你的公式是错误的。

从1到N中随机选取两个整数互质的概率趋近于6/pi^2，当N趋向于无穷大时。这是cp/(cf+cp)而不是cp/cf。

将你的pi替换为这个:

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

在我的机器上运行时，它给出3.14263472915。

- Paul Hankin

---

@Random832：那是因为还有另一个错误。应该是6/(cp/total)，而不是6/(cf/cp)或6/(cf/total)。 - user2357112

你能引用一下你最终公式的参考文献吗？我在查找正式推导，但它们都使用了测度论，我无法将其与互质数和余子式联系起来。 - Akshat Mahajan

https://en.wikipedia.org/wiki/Coprime_integers#Probabilities展示了一直到无限和的推导过程。展示无限和为pi^2/6则稍微困难一些，如果我没记错的话，你需要使用atan或类似函数的幂级数，在某个值处进行求值。这有点模糊，但希望能指引你朝着正确的方向前进。 - Paul Hankin

你程序中的一个小错误是，你的随机范围应该从1开始，而不是0。这会添加一个渐近趋近于零但不为零的误差。 - Paul Hankin

@Aaron 当然，我包含它是因为那是我在自己的电脑上运行的版本，而我的电脑没有py3.6。 - Paul Hankin