蒙地卡羅法求 PI

蒙地卡羅為摩洛哥王國之首都，該國位於法國與義大利國境，以賭博聞名。蒙地卡羅法是一種隨機演算，是指使用隨機亂數來解決計算問題的方法。

例如，以亂數散佈點配合面積公式，可以求得近似的 PI，是簡單的蒙地卡羅法運用。隨機演算雖然會有精確度等方面的疑慮，然而有些需求可能無最佳解，或者需要耗費大量運算才能得到理想解答，若隨機演算可以快速求得某個可接受的結果，就可能採取這類演算。

在圍棋上贏得人機之戰的 AlphaGo，採用的決策過程之一，就有蒙特卡羅樹搜尋（Monte Carlo Tree Search，MCTS），而 MCTS 存在目的是為了解決搜尋空間過於巨大，難以窮舉全部子樹的問題，MCTS 子節點展開後隨機進行遊戲，根據勝負結果來更新沿路至根節點的資訊，作為下次搜尋選取節點的評估之用。

解法思路

假設圓半徑為 1，四分之一圓面積就是 PI，能包含此四分之一圓的正方形面積就為 1，如果隨意地在正方形中散佈點，這些點有些會落於四分之一圓內，假設散佈了 n 個點，在圓內的點有 c 個，依比例來算，就會得到以下的公式：

蒙地卡羅法求 PI

要判斷產生的點是否落於圓內，可以令亂數產生 X 與 Y 兩個數值，如果 X² + Y² 小於 1 就是落在圓內。

程式實作

#include <stdio.h> 
#include <stdlib.h> 
#include <time.h> 

#define N 50001

int main(void) { 
    srand(time(NULL)); 
    
    int sum = 0;
    
    int i;
    for(i = 1; i < N; i++) { 
        double x = (double) rand() / RAND_MAX; 
        double y = (double) rand() / RAND_MAX; 
        if((x * x + y * y) < 1) {
            sum++; 
        }
    } 

    printf("PI = %f\n", (double) 4 * sum / (N - 1));
    
    return 0; 
}

import static java.lang.Math.*;
public class MonteCarlo {
    public static void main(String[] args) {
        final int N = 50001;
        int sum = 0; 
        for(int i = 1; i < N; i++) if(pow(random(), 2) + pow(random(), 2) < 1) {
            sum++;
        }
        System.out.printf("PI = %f%n", 4.0 * sum / (N - 1));  
    }
}

from random import random
N = 50001
print("PI =", 4 * len([1 for i in range(1, N) 
        if random() ** 2 + random() ** 2 < 1]) / (N - 1))

import java.lang.Math._
val N = 50000
printf("PI = %f%n", 4.0 * (for(i <- 1 to N
   if(pow(random(), 2) + pow(random(), 2) < 1)) yield 1).size / N)

N = 50000
print "PI = ", 4.0 * (1...N).map { 
    rand ** 2 + rand ** 2 < 1 ? 1 : 0 }.reduce(:+) / N

Array.prototype.reduce = function(init, f) {
    var value = init;
    for(var i = 0; i < this.length; i++) {
        value = f(value, this[i]);
    }
    return value;
};

function range(n) {
    var r = [];
    for(var i = 0; i < n; i++) {
        r[i] = i;
    }
    return r;
}

var n = 50000;
print(4 * range(n).map(function() {
    var x = Math.random();
    var y = Math.random();
    return x * x + y * y < 1 ? 1 : 0;
}).reduce(0, function(ac, elem) {
    return ac + elem;
}) / n);

import System.Random

rand gen n= take n $ randomRs (0.0, 1.0) gen::[Float]

main = do
    gen1 <- getStdGen
    gen2 <- newStdGen
    let n = 50000
        ic = length [1 | (x, y) <- 
                 zip (rand gen1 n) (rand gen2 n), (x ** 2 + y ** 2) < 1]
    print (4 * fromIntegral ic / fromIntegral n)

from '/lib/math' import random
from '/lib/math' import pow

N = 50001
(sum = iterate(0, N)
         .select(_ -> pow(random(), 2) + pow(random(), 2) < 1)
         .length())

println('PI = {0}'.format(4 * sum / (N - 1)))