リレー式計算機で円周率を計算する その1

サークルでの作業やら何やらで更新が遅れてしまいました。今回は使用する円周率計算公式の選定です。

円周率計算公式の候補

ライプニッツの公式

ライプニッツの公式(Wikipedia)

別名マータヴァ・ライプニッツ級数。円周率を求めるための公式の一つであるが、収束が非常に遅い(Wikipediaより)

マチンの公式

マチンの公式(Wikipedia)

グレゴリー級数を利用した、非常に収束の速い級数。1兆2411億桁まで計算したという記録もある。(Wikipediaより)

をグレゴリー級数

を利用することによって計算する。

実際に計算してみた

実際にテストプログラムを書き、計算してみました。固定小数点演算を前提とした計算のため、予め分子を10000倍し、項が1未満になった時点で計算を終了します。

ライプニッツの公式

プログラム

#include <stdio.h>
#define N 10000
#define MAX_STEP 200
int main(){
  int result=0;
  int step = 1;
  int u=1;
  int l=1;
  while(u >= 1 && l >= 1){
    u = (int)(N/(step*2-1));
    l = (int)(N/(step*2+1));
    result = u - l + result;
    step+=2;
    if(result*4 >= 31400 && result*4 < 31500){
      break;
    }
  }
  printf("Step:%d\n",(step+1)/2);
  printf("%d\n",result*4);
}

結果

Step:210
31400

目標値3.14を計算するまでに210ステップでした。

マチンの公式

プログラム

#include <stdio.h>
#define N 10000
#define MAX_STEP 200
int step=0;   
int test(int x){
    int a=x;
    int b=1;
    int pos=1;
    int ans;
    step++;
    ans=N/x;
    printf("Step:%d\n",step);
    int flag=1;
    while(pos>=1){
          a=a*x*x;
      b+=2;
      pos=N/(b*a);
      if(flag==0){
        ans+=pos;
            flag = 1;
      }else{
        ans-=pos;
            flag = 0;
      }
      step++;
      printf("Step:%d\n",step);
    };
    return ans;
}


int main(){
    printf("%d", 4*(4 * test(5) - test(239)));
    return 0;
}

結果

Step:5
31420

1つの項あたりの計算量がかなり異なるため一概には言えませんが、収束はかなり速いことが分かります。 計算誤差も桁数が増加した場合はマチンの公式のほうが有利であると考えられます。

まとめ

多倍長計算や乗算・除算の計算コストなどを全く考慮していないため、どちらの公式がいいとは言えませんが、桁数が増加した場合の誤差を考慮し、CPUエミュレータ作成まではマチンの公式を利用する方針とします。

（と言っても、利用する四則演算は変わらないんですけどね...)