Nim编程早茶
我们使用 Nim 语言根据蒙特卡罗算法估算圆周率。
蒙特卡罗算法简介
蒙特卡罗是一种随机模拟的算法，根据大数定律，蒙特卡洛算法的采样越多，越接近问题的最优解。但是蒙特卡洛算法只能保证尽量找到好的，不保证找到最优解。
用蒙特卡罗算法求圆周率
首先，我们考虑单位圆，单位圆的面积为 Pi / 4，蒙特卡罗算法的思路就是，以圆心为中心构建边长为 1 的正方形。随机生成大量(-1.0 .. 1.0, -1.0 .. 1.0) 的点，这些点落在正方形内的概率为 1，而落在单位圆内的概率为 Pi / 4。所以说，我们只要生成足够多的随机点，根据大数定律，落在单位圆内这一事件发生的频数除以随机点的总数，就近似于它发生的概率 Pi / 4。从而，我们可以计算出 Pi 的近似值。
首先，我们写一个函数，用来判断随机点是否落在单位圆内。
import math, random
proc inCircle(x: float, y: float, radius: float = 1.0): bool =
if x ^ 2 + y ^ 2 < radius ^ 2:
return true
else:
return false
我们可以自己写一个简易的随机数发生器，其中 n 为事件发生的次数。当随机点落到圆内，次数加一。最后乘 4 除以总数，就可以得到圆周率的近似解。
proc monteCarlo*(a: int = 1103515245,
c: int = 12345,
m: int = 1429496829,
z: var int,
n: int = 1000_0000): float =
var
x: float
y: float
hits = 0
for i in 1 .. n:
z = (a * z + c) mod m
x = z / m
z = (a * z + c) mod m
y = z / m
if inCircle(x, y):
hits += 1
return 4 * hits / n
我们还可以借助于 nim 的随机数模块
proc monteCarloRand*(n: int = 1000_0000): float =
randomize()
var
x: float
y: float
hits = 0
for i in 1 .. n:
x = rand(-1.0 .. 1.0)
y = rand(-1.0 .. 1.0)
if inCircle(x, y):
hits += 1
return 4 * hits / n
最后让我们来测试程序，可以使用命令 nimble install timeit 。
when isMainModule:
import timeit
block:
var z = 1
timeOnce("z = 1"):
echo monteCarlo(z = z)
timeOnce("rand"):
echo monteCarloRand()
block:
var z = 10
timeOnce("z = 10"):
echo monteCarlo(z = z)
block:
var z = 100
timeonce("z = 100"):
echo monteCarlo(z = z)
我们使用 -d:release 编译，结果还是有很大误差，主要是 n 的数值不足够大
3.1411848
z = 1 -> [433ms 839μs 100.00ns]
3.141876
rand -> [220ms 666μs 300.00ns]
3.1413456
z = 10 -> [480ms 766μs 500.00ns]
3.1413044
z = 100 -> [530ms 436μs 600.00ns]
我们可以增大 n 的数值，近似解变成了 3.141512232(由于随机的缘故，每次生成的结果都不同，但会稳定在 3.1415 左右)，但是相应的，执行时间也增加到了 20 s。
when isMainModule:
import timeit
block:
timeOnce("rand"):
echo monteCarloRand(10_0000_0000)
