テーラー展開で円周率を計算する[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がそこそこ大きければこの式で十分求められると思う。ただ厳密に計算したほうが計算回数は減って早いプログラムになる。その気があったらまたやります。その時はラマヌジャンの公式とかいい感じのやつ使います。