Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

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

になりました。

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

上記の2つ書いたら、
「そんな収束の遅いアルゴリズムなんて使ってんじゃねー!」的に指摘されてしまったので、
ガウス=ルジャンドルのアルゴリズムを使って計算してみます。
私、このアルゴリズム詳しくないので、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

その他

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

sassy_watson
Androidアプリやフロントエンドのエンジニア。
http://sassy.github.io/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした