math/03 のバックアップ(No.14)

円周率を求めよう！

(by K, 2004.01.30)
中学生にも理解できる範囲で、円周率の計算方法を説明できないかと言われて、考えた話。
読み物としてちょっと面白いだけで、学術的価値はありません。

↑

円周率パイ

円周率パイは、中学の範囲では、おそらく次のように認識されているでしょう。
```
円周と直径の比。もしくは、直径が1の円の円周。
3.1415926535 8979323846 2643383279 5028841971 6939937510 ...
```
- パイをギリシア文字で書きたかったんですが、ブラウザによってはきれいに出ないので、カタカナにしました。ご了承ください。
ということで、普通に考えると、円周率を100桁、1万桁、・・・と計算するには、とにかく正確にまんまるの物質を探して、これの円周と直径をとにかく正確に測って、割り算すればいいと思うでしょう。
しかしこの方法では100桁を正確に求めることはできません。そこまで精度のよい測定機器（ものさし）がまだないからです。
仮に測定器の精度の問題が解決したとしても、まだ問題はあります。その丸い物体はきっと原子でできています。ということは、本当はまんまるではなく、拡大してみれば原子くらいの大きさの凹凸があるでしょう。
この凹凸が円周率の計算にどのくらい影響するのか考えてみます。直径が2万kmある、非常に大きな円盤を考えます。これは地球より大きいです（！）。円周はだいたい62831853.0718mになるでしょう。一方原子の大きさは一番小さい水素原子が0.1nmくらいです。これは0.0000000001mです。
そうなると、62831853.07179...の先はしばらく正確に続きますが、小数第10位のくらいから先は原子のでこぼこの影響のせいで、正確なまんまるとは違う結果になるでしょう。
これだと円周率はせいぜい20桁くらいしか求められないことになります。
いやいや、これからは宇宙の時代なんだから、星をいくつかつぶしてでも宇宙空間上に非常に大きい円盤を作って、それを測ればきっと何万桁でも計算できる、と思うかもしれません。なんたって宇宙は広いですからね。また、そのころには技術が発達して原子よりも小さな原子核で円盤を作れるようになるはずだ、と考えるかもしれません。なるほど、それはいいアイデアです。
では円の直径を1千億光年にしましょう。これは現在考えられている宇宙の直径よりも大きいです。1千億光年は、だいたい1兆x1兆kmです。また原子核の大きさは、原子の大きさの約1万分の1です。
そうすると、さっきの例と比べれば、小数点の下に4桁増えて、小数点の上に16桁増えることになります。結局これでも合計40桁くらいです。上に書いた50桁すら求められません。
この先もっと小さい素粒子の利用を考えるかもしれませんし、それも可能かもしれませんが、それでも100桁の達成すらどうも無理そうです。
そんなわけでこれ以上正確な円周率を計算するには、測定ではなくて、もっと数学的な方法を使わないといけないわけです。

↑

数学的な計算方法

まず1辺が1mの正方形を考えます。
この左下の頂点を原点として、下の辺をx軸に見立てて0.1ずつ目盛を付けるとします。同じく左の辺をy軸に見立てて0.1ずつ目盛を付けることにします。さて、こうすると、0.1*0.1な正方形が100個できたことになります。

次に、原点を中心として、(1.0, 0.0)から(0.0, 1.0)を円弧にするような円盤を考えます。つまり半径1の円を1/4だけ切り取ったような形です。・・・さてこの円盤の中は、一体いくつの0.1*0.1のミニ正方形が含まれているでしょうか？
「仮に」78個だとしたら、この1/4円の面積は0.78平方メートルであると近似できることになります。
そして半径1の円の面積はパイであって、この円盤はその1/4なので、
- パイ＝4*0.78=3.12
と推定できることになります。これはなかなか良い感じです。

次に、ミニ正方形がいくつ含まれているのかを判定する方法を考えます。

たとえば4つの頂点が

 A(0.4, 0.3)  B(0.5, 0.3)
 D(0.4, 0.2)  C(0.5, 0.2)

のミニ正方形ABCDは円盤の中に含まれているでしょうかいないでしょうか？
これは三平方の定理を使うと判定できるのです。 0.4*0.4+0.3*0.3=0.25 なので、A点と原点の距離は0.5ということになります。これは円の半径よりも小さいので、A点は円盤の中です。
一般に、点(x, y)が円盤の中かどうかは、x*x+y*yを計算して、それが1以下かどうかを見ればいいのです。
これで4点とも円盤の中であれば、そのミニ正方形は間違いなく円盤の中にあります。また4点とも円盤の外であれば、そのミニ正方形は間違いなく円盤の外にあります。
4点のうち、一部だけが円盤の中に入っている場合は、そのミニ正方形が円弧によって分割されていることを意味します。・・・この分割された面積を0%とみなすか100%とみなすかで、円周率の推定の下限と上限を求めることができます。

この計算方法のポイントは、円周率の推定にあたって、定規などで「測る」という作業をなくしてしまったことです。
だから0.1ずつの目盛を0.000000001ずつの目盛とかにすれば、相応に高い精度で円周率を求められます。
それどころか、目盛を十分に細かくすれば、どこまででも好きな精度で円周率を求められます。

↑

プログラムで実験

では、実際に計算したらどうなるかやってみましょう。

#include <stdio.h>

#define N 10

int main()
{
    int i, j, k, c[5];
    double x0, y0, x1, y1;
    for (k = 0; k < 5; k++)
        c[k] = 0;
    for (j = 0; j < N; j++) {
        y0 = (double) j / N;
        y1 = (double) (j + 1) / N;
        for (i = 0; i < N; i++) {
            x0 = (double) i / N;
            x1 = (double) (i + 1) / N;
            k = 0;
            if (x0 * x0 + y0 * y0 <= 1.0) k++;
            if (x1 * x1 + y0 * y0 <= 1.0) k++;
            if (x0 * x0 + y1 * y1 <= 1.0) k++;
            if (x1 * x1 + y1 * y1 <= 1.0) k++;
            c[k]++;
        }
    }
    for (k = 0; k < 5; k++)
        printf("c[%d]=%d\n", k, c[k]);
    return 0;
}

とても単純なプログラムです。<math.h>すら使っていません。ミニ正方形の4つの頂点が円盤の中にあるかどうかを判定して、分類して集計しているだけです。
実行するとこんな結果になります。
```
c[0]=12, c[1]=5, c[2]=8, c[3]=6, c[4]=69
```
これはつまり、完全に円盤の中に入っていたものが69個で、頂点が3つだけ円盤の中に入ってたのが6個で、2つ入っていたのが8個で、1つ入っていたのが5つで、完全に円盤の外にあったものが12個あった、ということです。
c[1～3]の面積を0とみなしてパイを推定すると、2.76にしかなりません。これはかなり3.14から遠いです。がっかりですね。・・・逆にc[1～3]を100%のミニ正方形だと仮定してパイを推定すると、3.52となり、これも3.14から遠いです。
でもとりあえず上限と下限は分かりました。Nを大きくすれば、より誤差を狭めることが出来そうです。
ここでN=100へ進んでもいいのですが、c[1]を25%、c[2]を50%、c[3]を75%として計算すると、
```
69+5*0.25+8*0.50+6*0.75=78.75
```
となって、円盤の中にはミニ正方形が78.75個相当あると推定することもできます。この推定値からパイを計算すればなんと3.15となり、N=10の計算だけでも3.14にかなり近づくことができます！

ではN=100で計算してみます。

c[0]=2047, c[1]=58, c[2]=82, c[3]=59, c[4]=7754

これを使うとパイの値は、3.1016～3.1812だということが分かります。
さらに先の25%～75%の方法を使えば、パイの値の近似値は3.1415になります。

もっと先に進んでN=1000にしてみます。

c[0]=213614, c[1]=585, c[2]=828, c[3]=586, c[4]=784387
上限下限:3.137548～3.145544, (25%～75%で近似)=3.141547

さらにN=1万をやってみます。

c[0]=21450238, c[1]=5857, c[2]=8284, c[3]=5858, c[4]=78529763
上限下限:3.14119052～3.14199048, (25%～75%で近似)=3.14159051

↑

考察

上限下限方式は数学的な意味が一番高いと思います。いくつより大きくていくつより小さいかはっきりするからです。たとえばN=1万の結果を見れば 3.14119052～3.14199048 なので、3.141までは確定させることができます。しかしミニ正方形を1億個も作って比較して分類したのに、たったの4桁しか分からないというのは、なかなかさみしい気もします。
これを改善する方法としては、まず上限値と下限値を足して2で割ってみる、という単純な方法があります。もう一つは、上でやったように、25%～75%で近似するという方法もあります。僕の直感では、単純に足して2で割るよりは、25%～75%で近似するほうが正解に近いかなと思っているのですが、さてどうでしょうか？
- 推定したパイが3.14159265からどのくらいずれているのか？
  N=10 N=100 N=1000 N=1万
  単純に足して2で割る -0.00159265 -0.00019265 -0.00004665 -0.00000215
  25%～75%で近似 +0.00840735 -0.00009265 -0.00004565 -0.00000214
こうして見ると、25%～75%の近似が特に優れているとかそういうことはなくて、適当にやった上限と下限の単純平均でもそんなに悪くなさそうです。

（prog/01に、このプログラムの高速化の話題があります）

↑

リンク

http://www1.coralnet.or.jp/kusuto/PI/
- 円周率計算記録の歴史などが紹介されています。
- 16000桁、13万桁の計算結果も見られます。

↑

雑談

小学校のゆとり教育が始まったときに、円周率を3に近似するのは円を円に内接する正六角形に近似するのに相当するという批判があった。実は愚かにも僕もそれをすっかり信じた。しかしこれは少々おかしい気がする。正六角形の面積は、正三角形x6に等しいが、斜辺が1の正三角形の面積はsqrt(3)/4であって、これを6倍しても2.6にすら届かない（半径1の円の面積は円周率に等しい）。
いやおかしくなかった。円周で考えると、正六角形の周の長さは6だから、直径の2で割れば3なのか。なるほど。

となると日能研で紹介されていた東京大学の入試問題の、円周率が3.05以上であることを証明せよというのは、正12角形の周の長さを出せばいいのかな。それでも足りないのかな。足りなければ24角形だなあ。・・・いつか暇になったら計算しよう。

↑

こめんと欄

ルール違反でページを編集する荒らしが何度もきたので、当分の間このページは凍結します。コメントを書きたい場合はimpressionsへ。 -- K 2008-04-18 (金) 17:39:07
ちなみに凍結中でもこのようにコメントは書けますが、書き損じた場合に修正ができません。書き損じない自信があればここに書いてもいいです。 -- K 2008-04-18 (金) 17:42:59
書き始めから数えてもう10年以上になるのですが、ずっと書き途中になっていました。さすがに恥ずかしくなったので、予定していた内容を最後まで書きました。 -- K 2015-11-10 (火) 18:41:34

コメント	お名前	NameLink

	N=10	N=100	N=1000	N=1万
単純に足して2で割る	-0.00159265	-0.00019265	-0.00004665	-0.00000215
25%～75%で近似	+0.00840735	-0.00009265	-0.00004565	-0.00000214

math/03 のバックアップ(No.14)

円周率を求めよう！

円周率 パイ

数学的な計算方法

プログラムで実験

考察

リンク

雑談

こめんと欄

円周率パイ