mima2 site


円周率



弟の大学の講義で、円周率を点以下8桁までプログラミングで求めなさい、という課題が出ていました。
そんな面白そうな問題を兄である私が放っておくわけもなく・・・




円周率π

円周率は、円の周を直径で割った値で、

3. 14159 26535 89793 23846 26433 73279 50288 ...

といった具合に無限に続く無理数です。
通常、π(パイ)を使って表されます。




プログラミング

さて、プログラミングでπを求めるにはどうすればいいでしょう。

#include <stdio.h>
#include <math.h>

int main(void)
{
  printf("PI = %f\n", 2.0 * acos(0.0));
  return 0;
}

・・・いや、冗談です。

ライブラリ関数は、使わないことにしましょう。




遅い方法

πをプログラミングで求める方法ですぐに思いつく(?)のは、
arctan(x)のマクローリン展開にx = 1を代入して、

π/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...

というのがありました。(もちろん無限個は足せないので、途中で打ち切ります。)

しかしこの計算、集まりが悪すぎて正直使えませんw

というのも、精度が線形的にしか上がらないからです。例えば、500項くらいまで計算して、やっと3.14が求まります。50000項で3.1416まで、5百万項で3.1415927まで・・・と、非常に収束が遅いのがわかります。

いちおう、プログラムを載せておきます。
calculate_pi_1.c




速い方法

プログラミング向きの方法として、arcsin(x)のマクローリン展開にx=1/2を代入する方法が知られています。



ここでは、この式は正しいという前提で話を進めます。詳しいことは数学科の人にでもきいてください。




実際にコードを書いてみる。

8桁まで求めるのであれば、小数の表現はdouble型で十分でしょう。(double型の精度は10進で約15桁あります。)とはいっても、上のシグマの中の項を素直に計算するのはちょっと精度や計算時間の観点からきびしいかもしれないので、少し工夫してみます。

第(m-1)項から第m項にうつるときにどう変化するかを見てみます。(ただしmは1以上の整数)
     (2n)!・・・2m(2m-1)倍になる。
     2^(4n+1)・・・16倍になる。
     (n!)^2・・・m^2倍になる。
     (2n+1)・・・(2m+1)/(2m-1)倍になる。
つまり、まとめると、
     (2m-1)^2 / 8 m (2m+1)倍になる。
となります。

第0項は直接計算すると1/2であることがわかります。ここから第1項、第2項・・・を求めて次々に足していきます。8桁の精度で求めることになっているので、項の値がその100分の1、10-10より小さくなったら終了ということにしましょう。(最後に6倍するので、気持ち多めに計算しています。)

このロジックを元に作ったプログラムを載せておきます。
calculate_pi_2.c




発展:点以下1万桁まで求めてみよう!

つづく。