使用蒙特卡洛方法求PI的位数

Question

我尝试了很多用蒙特卡罗方法计算π的算法。这里有一个用Python写的解决方案：

def calc_PI():
    n_points = 1000000
    hits = 0

    for i in range(1, n_points):
        x, y = uniform(0.0, 1.0), uniform(0.0, 1.0)

        if (x**2 + y**2) <= 1.0:
            hits += 1

    print "Calc2: PI result", 4.0 * float(hits) / n_points

可悲的是，即使用到10亿次计算，结果的精度还是非常差（只有3.141...）。

这个方法能提供的最大精度就是这样吗？我选择蒙特卡罗方法是因为它很容易分成多个部分来并行计算。有没有其他的计算π的算法，也容易分成小块来计算呢？

Answer 1

使用一种叫做准随机数生成器的东西，代替普通的伪随机数生成器。准随机数生成器生成的数字在你需要计算的区域里分布得更均匀，这样可以让你在做蒙特卡罗积分时，结果更准确，收敛得更好。

Answer 2

你的误差是通过 sqrt(N)/N = 1/sqrt(N) 来计算的，所以这种方法获取准确估计的效率很低。这是由测量的统计特性决定的，无法改善。

对于 N 次投掷，你应该能得到大约 floor(log_10(N))/2-1 位的良好精度。为了安全起见，可能还要减去 -2。

即使这样，也假设你使用的是一个真正的随机数生成器（RNG）或者一个足够好的伪随机数生成器（PRNG）。

Answer 3

这是一个经典的蒙特卡洛方法的例子。不过，如果你想把计算圆周率的过程分成多个部分来并行处理，为什么不直接使用一个无限级数呢？你可以让每个处理器负责一个范围，然后在计算的过程中把结果加起来。

http://mathworld.wolfram.com/PiFormulas.html