Tweet

吉田武『素数夜曲 女王陛下のLISP』（東海大学出版会）が発売されて一部で話題になっている今日この頃ですが、この900ページ近い大著についてはまだコメントできる段階ではないので、この機会に思い出した同じ著者の『オイラーの贈物』のバグを報告しようと思います。この本を貶めるのではなく、この本のほんの少しの断片でさえもいろいろ調べて楽しむ題材になるということを伝えるのが目的です。

高校生のための数学入門書『オイラーの贈物』

『新装版 オイラーの贈物』第1版第7刷の付録A p. 314 には、モンテカルロ法を用いて円周率の近似値を求める、次のようなUBASICのコードが掲載されています（参照：UBASICの動かし方）。

10   randomize
20   input "10^n=";n
30   for i=1 to n
40   k=0:j=0
50   for j=1 to 10^i-1
60   x=rnd:y=rnd:if x*x+y*y<=1 then k=k+1
70   next j
80   print 10^i;"点中";k;"個命中";"Pi=";4*k/j
90   next i

このコードには、1つのバグといくつかの問題があります。

バグ

50行目のfor文で、jは1から(10^i-1)まで動くので、モンテカルロシミュレーションの試行回数は(10^i-1)回になります。しかし、80行目では試行回数を10^iあるいはj（ループが終了した時点でj==10^i）と、1回多く数えています。これにあわせるためには、50行目を次のように修正しなければなりません。

50   for j=0 to 10^i-1

問題

プロンプトがおかしい、とかいう細かい話はおくとして、

問題（無駄な計算）

50行目のfor文が終わると（80行目の直前）jの値は10^iになっているので、80行目で10^iを再度計算する必要はありません。80行目の末尾では、10^iの意味でjを使っているので、先頭でもそれを採用するといいでしょう。つまり、80行目は次のように書いた方がいいでしょう。

80   print j;"点中";k;"個命中";"Pi=";4*k/j

問題（無駄な代入）

40行目の「:j=0」は、その直後に「j=1」があるので無駄です。しかし、次の問題の解決に必要なので、そのままにしておいてもいいでしょう。

問題（for文の範囲）

UBASICの仕様はよくわかりませんが、手元のバージョン8.8fでは、for文の範囲（上限ではない）が2^32を超えると「数値が大き過ぎます」というエラーが出ます。ですから、せっかく大きな桁数をサポートするUBASICを使っても、このプログラムでは試行回数10^9回までしかシミュレートできません。本文のコラム（p. 226）で紹介されている実行結果もそうなっています。これではもったいないので、for文ではなくwhile文を使うようにするといいでしょう。

45   jMax=10^i
50   while j<jMax
60   x=rnd:y=rnd:if x*x+y*y<=1 then k=k+1
65   j=j+1
70   wend

問題というか要望（プログラミング言語の選択）

1993年に初版が出版されて以来、本書は売れ続け、つまり生き続けているわけですが、利用されている言語UBASICのほうは瀕死の状態、64bit版Windows 7では動きません（参考：UBASIC を Windows Vista/7 で動かす方法）。ですから、長生きしそうなプログラミング言語を採用し、改訂していただきたいものです。大きな桁の整数のサポートが必要です。C言語＋GMPなら寿命は長そうですがハードルが高いので、サポートが組み込まれているJava, Ruby, Python, Lispなどがいいでしょう。Javaに翻訳すると、下のようになります（まず紹介すべきはBigIntegerではなくintを使うバージョンですが、ここでは省略します）。

import java.util.*;
import java.math.*;

public class MonteCarlo {

  public static void main(String[] args) {
    BigInteger one = new BigInteger("1");

    System.out.print("10^n=");
    int n = new Scanner(System.in).nextInt();

    for (int i = 1; i <= n; i++) {
      BigInteger jMax = new BigInteger("10").pow(i);
      BigInteger k = new BigInteger("0");
      for (BigInteger j = k; j.compareTo(jMax) < 0; j = j.add(one)) {
        double x = Math.random();
        double y = Math.random();
        if (x * x + y * y <= 1) {
          k = k.add(one);
        }
      }
      System.out.printf("%s 点中 %s 個命中 Pi= %f\n",
          jMax, k, 4. * k.doubleValue() / jMax.doubleValue());
    }
  }
}

Javaで書くとなんかごちゃごちゃしていやですね。言語仕様が安定していない印象が強いRubyとPythonは採用しにくいので、残るのはLispですか。それなら、『素数夜曲 女王陛下のLISP』がちょうどいいですね（モンテカルロ法も扱われているようですし）。