円周率 - peanutsjamjam's diary

id:Jagdpanther さんが Gauche で円周率を計算しているのをみて、私もやってみる。

ガウチェ・モンテール円周率 - Jagdpantherの日記: http://d.hatena.ne.jp/Jagdpanther/20090126/1232921485

id:Jagdpanther さんは Gauche に標準で付いてくる Mersenne Twister を使ったモンテカルロ法で円周率を計算していたけど、こちらは別な方法で。

Gauche には Mersenne Twister も標準で付いているけれど、gcd という手続きも標準で付いてくるのだ。最大公約数。これは使える。なぜかというと、Wikipedia で、こんな記述を見つけてしまったから。

自然数から無作為に2つを取り出した時、その2つが互いに素である確率は 6/π² となる。

円周率 - Wikipedia

http://ja.wikipedia.org/wiki/%E5%86%86%E5%91%A8%E7%8E%87

ま、まじで？

と思ってしまったので、早速 Scheme で確認してみましょう。まずは手始めに、こんなのを。（まだ gcd は出てこない）：

(define (rc lim)
  (let loop-n ((m 1))
    (let loop-m ((n 1))

      (format #t "m=~d\tn=~d\n" m n)

      (when (< n m) (loop-m (+ n 1))))
    (when (< m lim) (loop-n (+ m 1)))))

これを、たとえば (rc 4) と呼び出してやる。実行するとこうなる。

m=1     n=1
m=2     n=1
m=2     n=2
m=3     n=1
m=3     n=2
m=3     n=3
m=4     n=1
m=4     n=2
m=4     n=3
m=4     n=4

やっていることは、集合 {1, 2, 3, 4} から、2 つを選び出しているだけ。「4 つのものから重複を許して 2 つを選ぶ」という重複組合せ（repeated combination）だ。出力の中に {1, 1} とか {3, 3} っていう行があるので、ちゃんと重複を許しているのがわかる。実際に動かすときには 4 ではなくて、もっと大きな数を与えてやる。

これを使って、「数のペアをひとつ決めるたびに、その 2 つの数が互いに素であるかどうか」を調べたい。そして、「互いに素だった場合」と「互いに素でなかった場合」の数を記録する。互いに素かどうか、を調べるには gcd（最大公約数）を調べればよい。ふたつの数が互いに素であるということは、最大公約数が 1 だ、ということだ。

でも、まず、gcd を使う前にちょっとリファクタリング。

さっき書いたコードを修正して、「数のペアをひとつ決めるたびに、何かする」という形に書き換える。その「何か」の部分に、さっきと同じ m と n を表示する部分を入れてやる。

#! /usr/local/bin/gosh

(define (rc proc lim)  ; lim だけじゃなく、手続きもひとつもらうようにした。
  (let loop-n ((m 1))
    (let loop-m ((n 1))

      (proc m n)       ; 手続きを m n のペアに適用する。

      (when (< n m) (loop-m (+ n 1))))
    (when (< m lim) (loop-n (+ m 1)))))

(define (main args)

  (define (print-m-n m n)
    (format #t "m=~d\tn=~d\n" m n))

  (rc print-m-n 4)    ; print-m-n という手続きと、4 を引数にして rc を呼ぶ。

  0)

これが、さっきのとまったく同じ出力になることを確認してから、gcd を使って「互いに素かどうか」を調べる部分を書く。ついでに、πも計算。

#! /usr/local/bin/gosh

(define (rc proc lim)
  (let loop-n ((m 1))
    (let loop-m ((n 1))

      (proc m n)

      (when (< n m) (loop-m (+ n 1))))
    (when (< m lim) (loop-n (+ m 1)))))

(define (main args)
  (define coprime 0)     ; 互いに素だった時に 1 増やす。
  (define non-coprime 0) ; 互いに素ではなかった時に 1 増やす。

  (define (calc-pi coprime non-coprime)
    (sqrt (/ 6 (/ coprime (+ coprime non-coprime))))) ; √(6/(互いに素が占める割合))

  (define (count m n)
    (if (= 1 (gcd m n))
      (set! coprime (+ coprime 1))
      (set! non-coprime (+ non-coprime 1)))
    (format #t "m=~d n=~d  gcd(~d, ~d)=~d" m n m n (gcd m n))
    (format #t "  coprime=~d  non-coprime=~d" coprime non-coprime)
    (format #t "  pi=~d\n" (calc-pi coprime non-coprime)))

  (rc count 4)

  0)

実行する。

m=1 n=1  gcd(1, 1)=1  coprime=1  non-coprime=0  pi=2.449489742783178
m=2 n=1  gcd(2, 1)=1  coprime=2  non-coprime=0  pi=2.449489742783178
m=2 n=2  gcd(2, 2)=2  coprime=2  non-coprime=1  pi=3.0
m=3 n=1  gcd(3, 1)=1  coprime=3  non-coprime=1  pi=2.8284271247461903
m=3 n=2  gcd(3, 2)=1  coprime=4  non-coprime=1  pi=2.7386127875258306
m=3 n=3  gcd(3, 3)=3  coprime=4  non-coprime=2  pi=3.0
m=4 n=1  gcd(4, 1)=1  coprime=5  non-coprime=2  pi=2.898275349237888
m=4 n=2  gcd(4, 2)=2  coprime=5  non-coprime=3  pi=3.0983866769659336
m=4 n=3  gcd(4, 3)=1  coprime=6  non-coprime=3  pi=3.0
m=4 n=4  gcd(4, 4)=4  coprime=6  non-coprime=4  pi=3.1622776601683795

たった 10 個のペアを与えただけで、もう 3.1 になってる。(rc count 4) の 4 を大きな数にしてあげると、じわじわとπに近づいていくのがわかる。
そのへんにいくらでも転がっている自然数を 2 つもってきて、「互いに素」かどうかを調べていくだけで円周率が現れてくれるってのが驚き。