2015年9月30日水曜日

テーラー展開で円周率を計算する[C言語]

思い立ったが吉日。任意の桁数の円周率の計算をしようと思う。ただ計算速度とかはとりあえず気にしない方向でいきます。

1.計算方法

計算方法はarcsinのテーラー展開したものを用いる。しかしちょっと工夫は加える。とりあえずarcsinのよくある公式としては,
である。これにx = 0.5を代入すれば π/6 となります。ここでちょっと工夫を加える。ホーナー法というやつを使う。
 コンピューターで数式を計算するときによくやる手法ですが,x^3+x^2+x+1のようなx^nを多く含んだ式を計算したい時,x*(1+x*(1+x)))+1という書き方をするほうが掛け算の回数が少なくなり,計算速度が上がるのだ。
 これを上の式に適用すると次のようになる。

2.double型で実装

 任意の桁数で計算を行うプログラムの前に,arcsinのプログラムをdouble型の桁数の範囲で実装してみる。式があってるかの確認のため。
/**
 * main.c
 * 2015/09/30
 * テーラー展開とホーナー法で円周率の計算を行うプログラム1
 */
#include <stdio.h>

/**
 * arcsin(x)を計算する関数 
 * n テイラー展開をn次の項まで計算する
 */
double arcsin(int n,double x){
   double pi = 0;
   double x2 = x * x;
   
   do{
      pi = pi + 1.0 / (2.0*n-1.0);
      pi = pi * (x2);
      pi = pi * (2.0*n-3.0) / (2.0*n-2.0);
      
      n--;
   }while(n >= 2);
   
   pi = pi + 1;
   pi = pi * x;
   
   return pi;
}

int main(void){
   printf("pi = %lf\n",6.0*arcsin(10,0.5));
   return 0;
}

 計算結果はpi = 3.141593でまぁ大丈夫そうだ。

3.任意の桁数で計算する

 上のプログラムを任意の桁数で計算できるように実装し直す。具体的にはN個の配列を作り,3.14ならpi[0]=3,pi[1]=1,pi[2]=4…というように,配列の要素1つに1桁の数字を割り当てる。多倍長整数でググると似たような情報が得られるはず。興味のある方はぜひググってください。
 とりあえずソースコードを貼る。

/**
 * main.c
 * 2015/09/30
 * テーラー展開とホーナー法で円周率の計算を行うプログラム1
 */
#include <stdio.h>
#include <math.h>

/* 計算する桁数 */
#define DIGIT 1000
/* 余分に計算する桁数 */
#define EXTRA_DIGIT 10
/* 配列の大きさ */
#define N (DIGIT+EXTRA_DIGIT)

/* i番目の要素に円周率の小数第i位の桁の数字が代入されている */
int pi[N] = { 0 };
/* 計算用の一時変数 */
int temp[N];

void assign(int*,int,int);
void divide(int*,int,int);
void multi(int*,int,int);
void add(int*,int*,int);
void arrange(int*,int);
void calculate();
void show();

int main(void){
   calculate();
   printf("pi = 3.\n");
   show();
   return 0;
}

/* 小数点以下を表示する関数 */
void show(){
   int i;
   for( i = 1 ; i < DIGIT + 1 ; i++ ){
      printf("%d",pi[i]);
      if( i % 10 == 0 ) printf(" ");
      if( i % 50 == 0 ) printf("\n");
   }
}

/* 配列arrに数字numを代入する */
void assign(int *arr,int length,int num){
   int i;
   
   arr[0] = num;
   for(i = 1 ; i < length ; i++ )
      arr[i] = 0;
}

/* 配列arrに数字num掛ける */
void multi(int *arr,int length,int num){
   int i;
   
   for( i=0 ; i < length ; i++ ){
      arr[i] *= num;
   }
   
   arrange(arr,length);
}

/* 配列leftに配列rightを加える */
void add(int *left,int *right,int length){
   int i;
   
   for( i = 0 ; i < length ; i++ ){
      left[i] += right[i];
   }
   
   arrange(left,length);
}

/* 配列arrをnumで割る */
void divide(int *arr,int length,int num){
   int i;
   
   for( i = 0 ; i < N - 1 ; i++ ){
      arr[i+1] += 10 * ( arr[i] % num );
      arr[i] /= num;
   }
}

/* 正しい10進数の小数にする */
void arrange(int *arr,int length){
   int i;
   
   for( i = N - 1 ; i >= 1 ; i-- ){
      arr[i-1] += arr[i] / 10;
      arr[i] %= 10;
   }
}

/* 円周率を計算する関数 */
void calculate(){
   /* 誤差が10^(-N)未満になるようなnの値を計算し代入 */
   int n = (int)( N * log(10) / log(2) - 1 )/2;
   
   do{
      /* arcsin(0.5)を計算 */
      assign(temp,N,1); /* tempに1を代入 */
      divide(temp,N,2*n-1); /* tempを2n-1で割る */
      add(pi,temp,N); /* piにtempを加える */
      divide(pi,N,4); /* piを4で割る(0.5の二乗を掛けている) */
      multi(pi,N,2*n-3); /* piに2n-3を掛ける */
      divide(pi,N,2*n-2); /* piを2n-2で割る */
   
      n--;
   }while( n>= 2 );
   
   /* piに1を加え2で割る(0.5を掛ける) */
   assign(temp,N,1);
   add(pi,temp,N);
   divide(pi,N,2);
   
   /* piには今π/6が代入されているので6をかける */
   multi(pi,N,6);
}

 ひとつ目のdouble型でやった時のプログラムではテイラー展開を何次の項で打ち切るかをこちらで決めていたが,任意の桁数を計算するためにはそうはいかない。一応誤差を評価して適切な次数を求めにゃならん。
 本当はベルヌーイだかコーシーの剰余項を見てちゃんと計算しなきゃならんが面倒なので最終項の大きさで適当にやった。正しくないけど計算があってるから大丈夫だろう。まぁ大筋は一応書いとく。
 まぁ割りと誤差は大きめにとったから
 で計算してみた。Nがそこそこ大きければこの式で十分求められると思う。ただ厳密に計算したほうが計算回数は減って早いプログラムになる。その気があったらまたやります。その時はラマヌジャンの公式とかいい感じのやつ使います。
 

2015年9月29日火曜日

モンテカルロ法で円周率を計算[C言語]

まぁやりつくされたネタではあるけれども,暇つぶしには丁度いいんじゃないかなと思って書くことにします。

1.モンテカルロ法とは

 簡単に言えば何かをシミュレーションをして,結果から求めたい値を探すというもの。1番わかりやすい例としてはサイコロ。サイコロが2個あります。2つとも同じ目が出る確率はいくつですか,という問題に対して実際に100回ぐらいサイコロふってみればわかるじゃないか,という話。そこから発展して,確率を求めることで円周率を計算したり数値的に積分したりできる。

2.具体的な計算法

 モンテカルロ法のWikipediaに良い図があるけど,一応一通り説明します。
 一辺が1の正方形の中に半径1の円が描かれています。ここで正方形の中に適当に砂を1粒落とすと,砂が円の内部に入る確率はいくつでしょう。
 正方形全体の面積は1,円の面積は半径1の円の4分の1の面積が描かれているからπ/4。したがって(確率) = (円の面積) ÷ (正方形全体の面積) だから答えはπ/4になる。

 じゃあ実際に適当に砂を,n粒落としたらm粒円の中に入ったとすると,
 
 と,こういう方法で近似的に円周率を求めることができる。

3.プログラム

 まぁ説明はあまり得意じゃないからね。プログラムで(お察しください)ということです。アルゴリズムは比較的単純。砂粒を落とすという表現を,一辺1の正方形の中にランダムに点をおくということに言い換える。

 m=0から始めて,
 「x,yに0から1の乱数を代入→円の中ならmを1増やす」という作業をn回繰り返す。
 円の中なら原点からの距離が1以下(未満)であれば良い。

/**
 * main.c
 * 2015/09/29
 * モンテカルロ法で円周率を求めるプログラム
 */

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

/* 試行回数 */
#define N 100000

/* ランダムに点を決め,円の中なら1,そうでなければ0を返す */
int step(){
   double x,y;
   
   /* x,yに0から1の乱数を代入する */
   x = (double)rand() / ( RAND_MAX + 1);
   y = (double)rand() / ( RAND_MAX + 1);
   
   /* 原点からの距離の二乗が1未満なら円の中 */
   if( x * x + y * y < 1 )
      return 1;
   
   return 0;
}

int main(void){
   int i;
   int m = 0;
   
   /* 乱数の種をまく */
   srand(time(NULL));
   
   for( i = 0 ; i < N ; i++ ){
      /* 点が円の中なら */
      if( step() == 1 ){
         /* mを1増やす */
         m++;
      }
   }
   
   /* 結果を出力する */
   printf("pi = %lf\n", 4.0 * m / N );
      
   return 0;
}


出力結果はpi = 3.144840と,まぁそこそこあってる。でも物足りない。計算した気がしない。円周率をゴリゴリ計算するプログラムはこっち。まぁこれも遅いけど一応任意の桁数を計算できます。

追記:Processingで可視化した

逐次近似法で非線形方程式を解く[C言語]

先日ブックオフにて古めかしい『情報処理入門コース7 数値計算』という本が200円で売られていた。こういう専門書とかは古くて汚れてても内容がしっかりしてるから,なかなか安くなっとる奴がないんですわ。

こういった経緯があって,せっかく衝動買いしたのでFortranで書かれたコードを自分なりにC言語で書きなおして読んでいこうと思う。

前置きはここまでにして早速本題へ。
逐次近似法 ( successive approximation method ) で非線形方程式を解くって言ったけれども,ざっくり言うと漸化式を計算するだけ(厳密な説明は省略)。高校とかでやってた「この数列の極限は~」の応用みたいなもの。でもちゃんと解けるからすごい。

今回は次の方程式を解いてみようと思います。
本の方程式とは違うやつにしたかっただけで,特に意味は無い。この式が0となるxを求めていきます。

1.漸化式を作る

すべての方程式で使えるというわけではないが,次のプロセスでだいたい漸化式ができる。
まずx=の形をとりあえず作る
で,それっぽくすれば漸化式の完成。
これが正しいかどうかは解との差が収束するかどうかを証明しなければならないのだが,なれない数式入力にギブアップ。こいつの極限が解に近づくという前提ですすめます。

2.グラフを書いて解の目星をつける

逐次近似法の悪いところは最初にどの辺に解があるかの目星をつけなければならないところ。場合によっては解に収束しないということがあります。
拡大して見ていただくと分かる通り,大体-0.3ぐらいでx軸に交わっている事がわかります。なので解はおおよそ-0.3だろうという前提ですすめていきます。

3.漸化式を計算する

数列Xnの初期値X1=-0.3として計算を始める。プログラムは以下の通り。
/**
 * main.c
 * 2015/09/29
 * モンテカルロ法で円周率を求めるプログラム
 */

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

/* 試行回数 */
#define N 100000

/* ランダムに点を決め,円の中なら1,そうでなければ0を返す */
int step(){
 double x,y;
 
 /* x,yに0から1の乱数を代入する */
 x = (double)rand() / ( RAND_MAX + 1);
 y = (double)rand() / ( RAND_MAX + 1);
 
 /* 原点からの距離の二乗が1未満なら円の中 */
 if( x * x + y * y < 1 )
  return 1;
 
 return 0;
}

int main(void){
 int i;
 int m = 0;
 
 /* 乱数の種をまく */
 srand(time(NULL));
 
 for( i = 0 ; i < N ; i++ ){
  /* 点が円の中なら */
  if( step() == 1 ){
   /* mを1増やす */
   m++;
  }
 }
 
 /* 結果を出力する */
 printf("pi = %lf\n", 4.0 * m / N );
  
 return 0;
}
このプログラムを実行するとx = -0.335420となる。THRESHOLDをもう少し小さくすればもうすこし精度は上がるが,丸め誤差があるからあまり小さくすると収束しなくなるのでお気をつけて。

実際にf(-0.335420)を計算すると-5.793E-6となりほぼ0に等しいことがわかります。
まぁ以上が大方の流れです。この後漸化式の証明と,誤差の評価とかをしなくちゃならないが,そこら辺は難しそうな本とか他のサイトに任せます。

ブログ,始めました。

どうも。プログラミングを勉強したりしているものです。

 今までは一人でやって来ましたが,いいとも悪いとも誰も言ってくれないのでネットに公開してみようと思い至りました。プログラミングを中心に,自分が好きなことを適当に書いていくつもりです。
 まぁあとちまちましたプログラムが邪魔なのでネットに公開すれば保存しなくて済むし,簡単に検索できそうだったからこうします。

ただ必ずしも正しいことを書くとは限らないのでご了承ください。とりあえず定型文として,

このサイトに書かれている内容によって生じた一切の損害・不利益等の責任を負いかねます。
また当サイトはリンクフリーです。

とだけ書いておきます。(かきゃあええっちゅうもんじゃにゃーがの)