3月14日は円周率の日！ということで円周率を求めるプログラムを書こう！

3月14日

3月14日は何の日ですか？
え、ホワイトデー？
誰だいそんなこと言っているやつは？
そんなリア充イベントにうつつ抜かしてるんじゃねー！
俺にはホワイトデーなんてかんけーねーんだよ！

3月14日は
3.14ということで円周率の日なんですよ！

ということで、円周率を求めるプログラムを書いてみます。

ライプニッツの公式

円周率をもとめるのには、ライプニッツの公式というのがあります。
詳しい説明はwkipedia見てね!
http://ja.wikipedia.org/wiki/ライプニッツの公式

ということで、以下に書いてみます。

pi.c

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

double leibniz(int n){
  int i=0;
  double ret = 0.0;
  for (i = 0; i < n;i++) {
    ret += pow(-1, i) / (2 * i + 1);
  }
  return 4 * ret;
}

int main(void){
  printf("%1.9f\n", leibniz(10));
  printf("%1.9f\n", leibniz(100));
  printf("%1.9f\n", leibniz(1000));
  printf("%1.9f\n", leibniz(10000));
  printf("%1.9f\n", leibniz(100000));
  printf("%1.9f\n", leibniz(1000000));
  printf("%1.9f\n", leibniz(10000000));
  printf("%1.9f\n", leibniz(100000000));
  printf("%1.9f\n", leibniz(1000000000));
  return 0;
}

gccでコンパイルするときは、-lmオプション忘れないでね!
結果は以下のようになります。

$ ./a.out
3.041839619
3.131592904
3.140592654
3.141492654
3.141582654
3.141591654
3.141592554
3.141592644
3.141592653

モンテカルロ法

モンテカルロ法というのはシミュレーションを乱数を用いて行う手法で、これを使って円周率を求められるようです。

長さ1の正方形の中をランダムな値でn個プロットしていき、それが正方形の一点を中心として半径1の扇形に入る点をr個としたときに、

n : r = 1 : π/4

というカタチが成り立つのを利用します。

pi2.c

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

/* n : r = 1 : pi/4 */
double montecarlo(int n) {
  int i, r;
  double x, y;
  for (i = 0, r = 0; i < n; i++) {
    x = (double)rand()/(RAND_MAX+1.0);
    y = (double)rand()/(RAND_MAX+1.0);
    if ((x*x + y*y) < 1.0)
      r++;
  }
  return 4.0* r / n;
}

int main(void){
  printf("%1.9f\n", montecarlo(10));
  printf("%1.9f\n", montecarlo(100));
  printf("%1.9f\n", montecarlo(1000));
  printf("%1.9f\n", montecarlo(10000));
  printf("%1.9f\n", montecarlo(100000));
  printf("%1.9f\n", montecarlo(1000000));
  printf("%1.9f\n", montecarlo(10000000));
  printf("%1.9f\n", montecarlo(100000000));
  printf("%1.9f\n", montecarlo(1000000000));
  return 0;
}

とすると、

$ ./a.out
3.600000000
2.720000000
3.216000000
3.145600000
3.140280000
3.141728000
3.141824800
3.141367920
3.141588128

になりました。

ガウス＝ルジャンドルのアルゴリズム

上記の２つ書いたら、
「そんな収束の遅いアルゴリズムなんて使ってんじゃねー！」的に指摘されてしまったので、
ガウス＝ルジャンドルのアルゴリズムを使って計算してみます。
私、このアルゴリズム詳しくないので、Wikipediaのをそのままコード化してしまいました><
http://ja.wikipedia.org/wiki/ガウス＝ルジャンドルのアルゴリズム

pi3.c

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

double GaussLegendre(int n) {
  int i;
  double a = 1.0, b = 1.0 / sqrt(2.0), t = 1.0 / 4, p = 1.0, tmp = 0;
  double ret;
  for (i = 0; i < n; i++) {
    tmp = a;
    a = (tmp + b) / 2;
    b = sqrt(tmp * b);
    t = t - (p * (a - tmp) * (a - tmp));
    p = 2 * p;
  }
  return (a + b) * (a + b) / (4 * t);
}

int main(void){
  printf("%1.20f\n", GaussLegendre(1));
  printf("%1.20f\n", GaussLegendre(2));
  printf("%1.20f\n", GaussLegendre(3));
  printf("%1.20f\n", GaussLegendre(4));
  printf("%1.20f\n", GaussLegendre(5));
  return 0;
}

結果は以下です。
確かに収束速い…。

$ ./a.out
3.14057925052216857509
3.14159264621354283875
3.14159265358979400418
3.14159265358979400418
3.14159265358979400418

その他

たぶんその他にも円周率を求める方法あると思うので、教えてくれたらうれしいです。
またなにか不備がありましたらコメントください。