円周率

マチンの公式を用いて円周率を求めるプログラムを作ったのですが、微妙に（３．１４２３６６６...）間違っています。
たぶん、配列repynの初期設定が違っているのだと思いますが、どうすればいいのでしょう？
教えてください。

関数：
add　配列型のデータ同士を加算する。
subtract　配列型のデータ同士を減算する。
divide　配列型のデータとint型の変数を除算する。

確認してあるので、関数に間違いはないと思います。

int main(void)
{
int *repyn,*tn,*pi;
int i;

repyn=(int *)calloc(DIGIT,sizeof(int));
tn=(int *)calloc(DIGIT,sizeof(int));
pi=(int *)calloc(DIGIT,sizeof(int));

for(i=0;i<DIGIT;i++) *(pi+i)=0;
for(i=0;i<DIGIT;i++) *(repyn+i)=0;
*(repyn+1)=8;
*(repyn+2)=0;

i=1;
do{
divide(repyn,5*5,repyn);
divide(repyn,i,tn);
add(pi,tn,pi);
i+=2;
divide(repyn,5*5,repyn);
divide(repyn,i,tn);
subtract(pi,tn,pi);
i+=2;
}while(i<=CN1);

for(i=0;i<DIGIT;i++) *(repyn+i)=0;
*(repyn+0)=9;
*(repyn+1)=5;
*(repyn+2)=6;

i=1;
do{
divide(repyn,239,repyn);
divide(repyn,239,repyn);
divide(repyn,i,tn);
subtract(pi,tn,pi);
i+=2;
divide(repyn,239,repyn);
divide(repyn,239,repyn);
divide(repyn,i,tn);
add(pi,tn,pi);
i+=2;
}while(i<=CN2);

for(i=0;i<DIGIT;i++){
printf("%d",*(pi+i));
}

return 0;

}

No.2 ベストアンサー 回答者： kacchann 回答日時： 2003/10/26 11:14

add, subtract, divideが間違ってるかも。

-------------------------------
const int kK = 10;

void add(int a[] , int b[] , int c[])
{
int i,cy=0;
for(i=DIGIT-1 ; i>=0 ; i--){
c[i]=(int)(a[i]+b[i]+cy);
if(c[i]<kK)
cy=0;
else{
c[i]=(int)(c[i]-kK);
cy=1;
}
}
}

void subtract(int a[] , int b[] , int c[])
{
int i,cy=0;
for(i=DIGIT-1 ; i>=0 ; i--){
c[i]=(int)(a[i]-b[i]+cy);
if(c[i]>0)
cy=0;
else{
c[i]=(int)(kK+c[i]);
cy=-1;
}
}

}
void divide(int a[] , int b, int c[])
{
int i;
long d,rem=0;
for(i=0;i<DIGIT;i++){
d=a[i];
c[i]=(int)((d+rem)/b);
rem=((d+rem)%b)*kK;
}
}

参考URL：http://www.ccad.sccs.chukyo-u.ac.jp/~mito/ss/pro …

この回答へのお礼

回答ありがとうございます。
ご指摘のとおり関数が間違っていました。
一応載せておきます。

/*多数桁加算*/
void add(int *a,int *b,int *result)
{
/*carryは繰り上げ、iは制御変数*/
int carry,i;

/*carryの初期化*/
carry=0;

/*result=a+bの演算*/
for(i=DIGIT1-1;i>=0;i--){
*(result+i)=*(a+i)+*(b+i)+carry;
if(*(result+i)>=10){
*(result+i)-=10;
carry=1;
}
else carry=0;
}
}

/*多数桁減算*/
void subtract(int *a,int *b,int *result)
{
/*borrowは繰り下げ、iは制御変数*/
int borrow,i;

/*borrowの初期化*/
borrow=0;

/*result=a-bの演算*/
for(i=DIGIT1-1;i>=0;i--){
*(result+i)=*(a+i)-*(b+i)-borrow;
if(*(result+i)<0){
*(result+i)+=10;
borrow=1;
}
else borrow=0;
}
}

/*多数桁除算*/
void divide(int *a,int b,int *result)
{
/*remは余り、iは制御変数*/
int rem1,rem2,i;

/*result=a/bの演算*/
rem2=0;
for(i=0;i<DIGIT1;i++){
rem1=(*(a+i)+rem2*10)%b;
*(result+i)=(*(a+i)+rem2*10-rem1)/b;
rem2=rem1;
}
}

お礼日時：2003/10/27 23:08

No.1 回答者： paspas 回答日時： 2003/10/25 21:38

どなたもアドバイスがないようなので少々コメントさせていただきます。

マチンの公式
（π／４）＝４×arctan(1/5)-arctan(1/239)
ですね。
ということは、
π= 16 * arctan(1/5) - 4 * arctan(1/239)

arctan x = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 ....

とりあえず項ごとの収束値を算出し、間違いを探ってみてはいかがでしょう。

arctan(1/5) ≒ 0.1973955...
arctan(1/239) ≒ 0.0041840...

16*arctan(1/5) ≒ 3.15832895....
4*arctan(1/239) ≒ 0.01673630...

16*arctan(1/5) - 4*arctan(1/239) ≒ 3.14159265...