一般法によるPIの求め,および超長数演算による長PIの求め(まとめ)

一般法PIを求める
簡単ですが、ちょっと例を挙げます.Leibnizの定理によると:
1/1-1/3+1/5-1/7+1/9...=PI/4
コードは次のように実装されます.

//   
int i,j=0;
double pa=0;
for(i=1;i<40000;i=i+2)//    i            
{
    if(j==0)
    {
        pa=pa+1.0/i;
        j=1;
    }
    else
    {
        pa=pa-1.0/i;
        j=0;
    }
}
pa=pa*4;
printf("%f
",pa);

注意:同じ理屈でwallis式でPIを求めることができます.
wallis式:
(2/1)*(2/3)*(4/3)*(4/5)*(6/5)*(6/7)*(8/7)*(8/9)...=PI/2
(コード略)
超長数演算を用いて長PI原理を求める
円周率の後の小数位数は無限で、どのようにコンピュータを使ってこの無限の小数を計算するかはいくつかの数学者とプログラムデザイナーが興味を持っているので、ここで1つの公式を紹介して大数の演算を合わせて、指定の桁数の円周率を計算することができます.
まずJ.Marchinの円周率の公式を紹介します.
PI = [16/5 - 16/(3*53) + 16/(5*55) - 16/(7*57) +......] -
      [4/239 - 4/(3*2393) +4/(5*2395)- 4/(7*2397)+ ......]
この式は次のように整理できます.
PI = [16/5 - 4/239] - [16/(53) - 4/(2393)]/3+[16/(55)- 4/(2395)]/5+ ......
すなわち、n番目の項は、奇数であれば正数、偶数であれば負数であり、項数の表現は、
[16/52*n-1 - 4/2392*n-1]/(2*n-1)
円周率が10までの負のL次数を計算する場合、[16/52*n-1-4/2392*n-1]の16/52*n-1は4/2392*n-1より大きく、決定的であるため、少なくともn番目の項に計算しなければならないことを示す.
[16/52*n-1 ]/(2*n-1) = 10-L
上記の式をlogとして簡略化し、求めることができます.
n = L/(2log5) = L/1.39794
したがって、小数以下L桁までの精度が要求される場合、式のn番目の項にのみ要求され、nは次のようになります.
n = [L/1.39794] + 1
上式では[]はガウス記号、すなわち整数(L/1.39794以下の整数)とし、簡略化の便宜上、以下の式を用いてn項目を簡略化することができる.
[Wn-1/52- Vn-1/(2392)]/(2*n-1)
この式の演算アルゴリズムと大数演算関数式を組み合わせた演算アルゴリズムはdiv(w,25,w);
div(v, 239, v);
div(v, 239, v);
sub(w,v, q);
div(q,2*k-1, q)
大数演算の演算アルゴリズムについては、前の文章を参考にして、出力時に注意しなければならないのは、出力アレイの整数値であるため、アレイの整数ビット数が4ビット未満であれば、0を補う必要があり、C言語ではフォーマットで字%04 dを指定すれば、不足ビット数部分が自動的に0を補って出力され、Javaの部分については、フォーマットにはNumberFormatを使用します.
超長整数演算による長PI実現

#define L 2000  // L   
#define N L/4+1 // N array  （               4 ，   1  ）
void add(int*, int*, int*);
void sub(int*, int*, int*);
void divv(int*, int, int*);

//   
int s[N+3] = {0};
int w[N+3] = {0};
int v[N+3] = {0};
int q[N+3] = {0};
int n = (int)(L/1.39793 + 1);//           
int k;

w[0] = 16*5;
v[0] = 4*239;
for(k = 1; k <= n; k++) {
    //     
    divv(w, 25, w);//   ，  ， 
    divv(v, 239, v);
    divv(v, 239, v);
    sub(w, v, q);
    divv(q, 2*k-1, q);
    if(k%2) //    
        add(s, q, s);
    else    //    
        sub(s, q, s); 
} 
printf("%d.", s[0]); 
for(k = 1; k < N; k++) 
    printf("%04d", s[k]); //     s    
printf("
");

//        
void add(int *a, int *b, int *c) {
    int i, carry = 0;
    for(i = N+1; i >= 0; i--) {
        c[i] = a[i] + b[i] + carry;
        if(c[i] < 10000)
            carry = 0;
        else { //   
            c[i] = c[i] - 10000;  
            carry = 1;
        }
    }
}

//        
void sub(int *a, int *b, int *c) {
    int i, borrow = 0;
    for(i = N+1; i >= 0; i--) {
        c[i] = a[i] - b[i] - borrow;
        if(c[i] >= 0)
            borrow = 0;
        else { //   
            c[i] = c[i] + 10000;
            borrow = 1;
        }
    }
}

//        
void divv(int *a, int b, int *c) {//        
    int i, tmp, remain = 0; 
    for(i = 0; i <= N+1; i++) { 
        tmp = a[i] + remain; 
        c[i] = tmp / b; 
        remain = (tmp % b) * 10000; 
    } 
    
}