モンテカルロ法

あんまりよくわかっていないけど、ランダムな 2つの自然数の
最大公約数が 1になる確率は 6 / π^2 らしい。そんな話を
昔聞いたような聞かなかったような。忘れた。

schemeで実装してみた。

#!/usr/bin/env gosh
(use srfi-27)

(define (get-pi experiments)
  (let ((random-size 10000))
    (define (monte test-count count)
      (let ((a (+ (random-integer random-size) 1))
            (b (+ (random-integer random-size) 1)))
        (if (= test-count 0)
            (calc-pi (/ count experiments))
            (if (= (gcd a b) 1)
                (monte (- test-count 1) (+ count 1))
                (monte (- test-count 1) count)))))
    (monte experiments 0)))

(define (calc-pi probability)
  (sqrt (/ 6 probability)))

(define (gcd x y)
  (if (= (modulo x y) 0)
      y
      (gcd y (modulo x y))))

(get-pi 100000) ;; 3.1440685912455315

perlでも実装してみた。

#!/usr/bin/env perl
use strict;
use warnings;

print "pi is ", get_pi(500000), "\n"; #  3.14165751518437

sub get_pi {
    my $count = shift;

    my $num = 0;
    foreach (1..$count) {
        my $a = int(rand(10000)) + 1;
        my $b = int(rand(10000)) + 1;

        if (gcd($a, $b) == 1) {
            $num++;
        }
    }

    my $probability = $num / $count;

    my $pi = sqrt(6 / $probability);

    return $pi;
}

sub gcd {
    my ($x, $y) = @_;

    if ($x % $y == 0) {
        return $y;
    } else {
        gcd($y, $x % $y);
    }
}

回数を増やすと piの値に近づくはず。
こういうのって Schemeのような関数型言語で考えるといいんだろうけど、
なかなかそういう頭になれません。勉強しないとな。