こんにちはゲストさん。会員登録(無料)して質問・回答してみよう!

解決済みの質問

円周率

マチンの公式を用いて円周率を求めるプログラムを作ったのですが、微妙に(3.1423666...)間違っています。
たぶん、配列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;

}

投稿日時 - 2003-10-24 23:17:21

QNo.688295

すぐに回答ほしいです

質問者が選んだベストアンサー

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/progfun/pai/

投稿日時 - 2003-10-26 11:14:37

お礼

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

/*多数桁加算*/
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:05

このQ&Aは役に立ちましたか?

2人が「このQ&Aが役に立った」と投票しています

回答(2)

ANo.1

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

マチンの公式
(π/4)=4×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...

投稿日時 - 2003-10-25 21:38:02

あなたにオススメの質問