多倍長で円周率を計算する(Machinの公式を使って)

ひろゆき
記事: 0
登録日時: 15年前
住所: 大阪市

多倍長で円周率を計算する(Machinの公式を使って)

投稿記事 by ひろゆき » 14年前

1. 自分は今何がしたくて

円周率の小数第600桁までの計算が出来るようにしないといけないんですけど
出来ないんです。
どこが間違っているか教えてくれませんか?

2. どう取り組んで(作ったプログラムはどれで
ここで使ってるdoはdo~while文は条件式を後判定して反復制御を行うために使っています。
まずは浮動小数点型を使って円周率を求めるプログラムを書いてみました。

#include

double arctan(double x){
double a=x;//乗数
double b=1.0;//分母
int flag=0;//引き算/足し算を交互に行うフラグ
double pos;
double ans=x;

do{
a=a*x*x;
b+=2;
pos=a/b;

if(flag=(!flag)){
ans-=pos;
}else{
ans+=pos;
}
}while(pos>0.00001);//小数点以下5桁の精度であきらめる
return ans;
}


int main(){
printf("%.13lf\n", 4.0*(4.0 * arctan(1.0/5.0) - arctan( 1.0 / 239.0)));
return 0;
}

このやり方では小数点型の計算ができなかったので、
浮動小数点型のプログラムを整数型に修正してみした。
#include

int test(int x){
int a=x;//乗数
int b=1;//分母
int flag=0;//引き算/足し算を交互に行うフラグ
int pos;
int ans;

ans=10000000/x;
do{
a=a*x*x;
b+=2;
pos=10000000/(b*a);

if(flag=(!flag)){
ans-=pos;
}else{
ans+=pos;
}
}while(pos>1);

return ans;
}


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

次に、多倍長の計算方法のメイン関数部分を上のプログラムと同じように作りました。

#define SIZE 304 //600桁求める

short a[SIZE],x[SIZE]//乗数
,b[SIZE]//分母
,c[SIZE]//分子
,ans[SIZE]
,ans2[SIZE]
,pos[SIZE]
,pos2[SIZE]
,val_2[SIZE]={2}// 2
,val_4[SIZE]={4};// 4

int main(){
int i;
int flag=0;//引き算/足し算を交互に行うフラグ

b[0]=1;//分母
c[SIZE-1]=1;//分子
x[0]=5;//乗数
for(i=0;i


void add(const short*ar,const short*ar2,short*ar3,const int size);

void sub(const short*ar,short*ar2,short*ar3,const int size);

void mul(const short*ar,const short*ar2,short*ar3,int size);

void div(short*ar,short*ar2,short*ar3,const int size);

void ar_print(short*ar,const int size);

void dp_print(short*ar,const int size,const int low);

void turn(short*ar,const int size);
int cmp(const short*ar,const short*ar1,const int size);
int cmp0(const short*ar,const int size);
void r_shift(short*ar,const int size);
void l_shift(short*ar,const int size);




#define SIZE 304

short a[SIZE],x[SIZE]
,b[SIZE]
,c[SIZE]
,ans[SIZE]
,ans2[SIZE]
,pos[SIZE]
,pos2[SIZE]
,val_2[SIZE]={2}
,val_4[SIZE]={4};

int main(){
int i;
int flag=0;

b[0]=1;
c[SIZE-1]=1;
x[0]=5;
for(i=0;i99){
ar3-=100;
ar3[i+1]+=1;
}else if(ar30){
if(0>ar3){
ar3[i+1]-=1;
ar3+=100;
}
}else{
if(099 || ar3[j+i]=0){
if(cmp(ar,ar2,size)>=0){
sub(ar,ar2,ar,size);
if(pos%2){
ar3[pos/2]+=10;
}else{
ar3[pos/2]++;
}
}else{
if(pos) r_shift(ar2,size);
pos--;
}
}

if(-1==div_flag) turn(ar,size);
if(-1==met_flag) turn(ar2,size);
if(div_flag!=met_flag) turn(ar3,size);
}

void ar_print(short*ar,const int size){
int i;
int start_flag=0;
int p_flag=cmp0(ar,size);
if(p_flag==0){
printf("0");
return;
}
if(p_flag==-1){
printf("-");
turn(ar,size);
}
for(i=size-1;i>=0;i--){
if(ar || start_flag){
if(!start_flag) printf("%d",ar);
else printf("%02d",ar);
start_flag=1;
}
}
if(p_flag==-1) turn(ar,size);
}

void dp_print(short*ar,const int size,const int low){
int i,j=0,k=0;
int p_flag=cmp0(ar,size);
if(p_flag==0){
printf("0");
return;
}
if(p_flag==-1){
printf("-");
turn(ar,size);
}

printf("%d.",ar[size-1]);

for(i=size-2;i>=low;i--){
if(j++%5==0) printf(" ");
if(k++%25==0) printf("\n");
printf("%02d",ar);

}
if(p_flag==-1) turn(ar,size);
}

void turn(short*ar,const int size){
int i;
for(i=0;i=0;i--){
if(arar1) return 1;
}
return 0;
}


int cmp0(const short*ar,const int size){
int i;
for(i=size-1;i>=0;i--){
if(ar[i]0) return 1;
}
return 0;
}

3. どのようなエラーやトラブルで困っていて
このような実行エラーが出まして
workf1.c: In function ‘main’:
workf1.c:57: 警告: 真偽値として使われる代入のまわりでは、丸括弧の使用をお勧めします
workf1.c:84: 警告: 真偽値として使われる代入のまわりでは、丸括弧の使用をお勧めします
/tmp/ccAcaEZ5.o: In function `div':
workf1.c:(.text+0x8cb): undefined reference to `l_shift'
workf1.c:(.text+0x95f): undefined reference to `r_shift'
collect2: ld はステータス 1 で終了しました

このように考え直したら
#include

double arctan(double x){
double a=x;
double b=1.0;
int flag=0;
double pos;
double ans=x;

do{
a=a*x*x;
b+=2;
pos=a/b;

if(flag!=flag){
ans-=pos;
}else{
ans+=pos;
}
}while(pos>0.00001);
return ans;
}

int main(void){
double c;
c=4.0*(4.0 * arctan(1.0/5.0) - arctan( 1.0 / 239.0));
printf("%.13lf\n", c);
return 0;
}

#include

// ar+ar2=ar3
void add(const short*ar,const short*ar2,short*ar3,const int size);
// ar-ar2=ar3
void sub(const short*ar,short*ar2,short*ar3,const int size);
// ar*ar2=ar3
void mul(const short*ar,const short*ar2,short*ar3,int size);
// ar/ar2=ar3...ar
void div(short*ar,short*ar2,short*ar3,const int size);
// printf(ar)
void ar_print(short*ar,const int size);
//配列を固定小数点型として画面表示する ar[size-1].ar[size-2]~ar[low]
void dp_print(short*ar,const int size,const int low);

void turn(short*ar,const int size);
int cmp(const short*ar,const short*ar1,const int size);
int cmp0(const short*ar,const int size);
void r_shift(short*ar,const int size);
void l_shift(short*ar,const int size);




#define SIZE 304 //600桁求める

short a[SIZE],x[SIZE]//乗数
,b[SIZE]//分母
,c[SIZE]//分子
,ans[SIZE]
,ans2[SIZE]
,pos[SIZE]
,pos2[SIZE]
,val_2[SIZE]={2}// 2
,val_4[SIZE]={4};// 4

int main(){
int i;
int flag=0;//引き算/足し算を交互に行うフラグ

b[0]=1;//分母
c[SIZE-1]=1;//分子
x[0]=5;//乗数
for(i=0;i99){
ar3[i]-=100;
ar3[i+1]+=1;
}else if(ar3[i]0){
if(0>ar3[i]){
ar3[i+1]-=1;
ar3[i]+=100;
}
}else{
if(099 || ar3[j+i]=0){
if(cmp(ar,ar2,size)>=0){
sub(ar,ar2,ar,size);
if(pos%2){
ar3[pos/2]+=10;
}else{
ar3[pos/2]++;
}
}else{
if(pos) r_shift(ar2,size);
pos--;
}
}

if(-1==div_flag) turn(ar,size);/*被除数の符号を戻す*/
if(-1==met_flag) turn(ar2,size);/*法の符号を戻す*/
/*被除数と法の符号が違うため、答えをマイナスにする*/
if(div_flag!=met_flag) turn(ar3,size);
}

/*画面表示*/
void ar_print(short*ar,const int size){
int i;
int start_flag=0;
int p_flag=cmp0(ar,size);
if(p_flag==0){
printf("0");
return;
}
if(p_flag==-1){
printf("-");
turn(ar,size);
}
for(i=size-1;i>=0;i--){
if(ar[i] || start_flag){
if(!start_flag) printf("%d",ar[i]);
else printf("%02d",ar[i]);
start_flag=1;
}
}
if(p_flag==-1) turn(ar,size);
}

/*配列を固定小数点型として画面表示する ar[size-1].ar[size-2]~ar[low]*/
void dp_print(short*ar,const int size,const int low){
int i,j=0,k=0;
int p_flag=cmp0(ar,size);
if(p_flag==0){
printf("0");
return;
}
if(p_flag==-1){
printf("-");
turn(ar,size);
}

printf("%d.",ar[size-1]);

for(i=size-2;i>=low;i--){
if(j++%5==0) printf(" ");
if(k++%25==0) printf("\n");
printf("%02d",ar[i]);

}
if(p_flag==-1) turn(ar,size);
}

/*符号の反転*/
void turn(short*ar,const int size){
int i;
for(i=0;iar1 の場合は+1を返す*/
/*ar==ar1 の場合は0を返す*/
int cmp(const short*ar,const short*ar1,const int size){
int i;
for(i=size-1;i>=0;i--){
if(ar[i]ar1[i]) return 1;
}
return 0;
}

/*arがプラスの値をとる場合は+1を返す*/
/*arがマイナスの値をとる場合は-1を返す*/
/*arがゼロの値をとる場合は0を返す*/
int cmp0(const short*ar,const int size){
int i;
for(i=size-1;i>=0;i--){
if(ar[i]0) return 1;
}
return 0;
}


/*右1シフト*/
void r_shift(short*ar,const int size){
int i,pos;
ar[0]=ar[0]/10;
for(i=1;i=0;i--){
ar[i]=ar[i]*10;
pos=ar[i]/100;
ar[i]-=(pos*100);
if(pos) ar[i+1]+=pos;
}

}

このような実行エラーがいくつもでます。
workf.c:196: error: expected identifier or ‘(’ before ‘/’ token
workf.c:197: error: stray ‘\213’ in program
workf.c:334:1: error: 終端されていないコメント

コメントはまだありません。