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

[ＨＩＲＯＹＵＫＩ]

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

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

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

#include <stdio.h>

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 <stdio.h>

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<SIZE;i++) a=x;
div(c,x,ans,SIZE);
for(i=0;i<SIZE;i++) c=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

if(flag=(!flag)){
sub(ans,pos,ans,SIZE);
}else{
add(ans,pos,ans,SIZE);
}
}while(cmp0(pos,SIZE)!=0);


flag=0;//引き算/足し算を交互に行うフラグ

for(i=0;i<SIZE;i++){ c=0;x=0;b=0;}
b[0]=1;//分母
c[SIZE-1]=1;//分子
x[0]=39;//乗数
x[1]=2;
for(i=0;i<SIZE;i++) a=x;
div(c,x,ans2,SIZE);
for(i=0;i<SIZE;i++) c=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

if(flag=(!flag)){
sub(ans2,pos,ans2,SIZE);
}else{
add(ans2,pos,ans2,SIZE);
}
}while(cmp0(pos,SIZE)!=0);

mul(ans,val_4,pos,SIZE);
sub(pos,ans2,pos,SIZE);
mul(pos,val_4,pos2,SIZE);


dp_print(pos2,SIZE,3);//誤差が生じるため下6桁を切り捨てて固定小数点形式で表示
printf("\n");

return 0;
}

これらをまとめた円周率600桁を計算するソースコードはこうなったんですが

#include <stdio.h>


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;i<SIZE;i++) a[i]=x[i];
div(c,x,ans,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;

do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;

if(flag=(!flag)){
sub(ans,pos,ans,SIZE);
}else{
add(ans,pos,ans,SIZE);
}
}while(cmp0(pos,SIZE)!=0);


flag=0;

for(i=0;i<SIZE;i++){ c[i]=0;x[i]=0;b[i]=0;}
b[0]=1;
c[SIZE-1]=1;
x[0]=39;
x[1]=2;
for(i=0;i<SIZE;i++) a[i]=x[i];
div(c,x,ans2,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;
do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;

if(flag=(!flag)){
sub(ans2,pos,ans2,SIZE);
}else{
add(ans2,pos,ans2,SIZE);
}
}while(cmp0(pos,SIZE)!=0);

mul(ans,val_4,pos,SIZE);
sub(pos,ans2,pos,SIZE);
mul(pos,val_4,pos2,SIZE);


dp_print(pos2,SIZE,3);
printf("\n");

return 0;
}


void add(const short*ar,const short*ar2,short*ar3,const int size){
int i,p_flag;
for(i=0;i<size;i++){
ar3[i]=ar[i]+ar2[i];
}
for(i=0;i<size;i++){
if(ar3[i]>99){
ar3[i]-=100;
ar3[i+1]+=1;
}else if(ar3[i]<-99){
ar3[i]-=-100;
ar3[i+1]+=-1;
}
}

p_flag=cmp0(ar3,size);
if(p_flag==0) return;

for(i=0;i<size;i++){
if(p_flag>0){
if(0>ar3[i]){
ar3[i+1]-=1;
ar3[i]+=100;
}
}else{
if(0<ar3[i]){
ar3[i+1]-=-1;
ar3[i]+=-100;
}
}
}
}

void sub(const short*ar,short*ar2,short*ar3,const int size){
turn(ar2,size);
add(ar,ar2,ar3,size);
turn(ar2,size);
}

void mul(const short*ar,const short*ar2,short*ar3,int size){
int i,j,pos;
for(i=0;i<size;i++){
ar3[i]=0;
}
for(i=0;i<size;i++){
if(ar2[i]){
for(j=0;j<size;j++){
ar3[j+i]+=ar[j]*ar2[i];

if(ar3[j+i]>99 || ar3[j+i]<-99){
pos=ar3[j+i]/100;
ar3[j+i]-=(pos*100);
ar3[j+1+i]+=pos;
}
}
}
}
}


void div(short*ar,short*ar2,short*ar3,const int size){
int i,pos=0;
int div_flag;
int met_flag;

for(i=0;i<size;i++){
ar3[i]=0;
}

div_flag=cmp0(ar,size);
met_flag=cmp0(ar2,size);

if(!div_flag || !met_flag) return;
if(-1==div_flag) turn(ar,size);
if(-1==met_flag) turn(ar2,size);

while(1){
if(cmp(ar,ar2,size)==1 && (ar2[size-1]/10)==0){
pos++;
l_shift(ar2,size);
}else break;
}

while(pos>=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);
}

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;i<size;i++){
ar[i]=-ar[i];
}
}


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;
if(ar[i]>ar1[i]) 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;
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 <stdio.h>

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 <stdio.h>

// 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;i<SIZE;i++) a[i]=x[i];
div(c,x,ans,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

if(flag=(!flag)){
sub(ans,pos,ans,SIZE);
}else{
add(ans,pos,ans,SIZE);
}
}while(cmp0(pos,SIZE)!=0);


flag=0;//引き算/足し算を交互に行うフラグ

for(i=0;i<SIZE;i++){ c[i]=0;x[i]=0;b[i]=0;}
b[0]=1;//分母
c[SIZE-1]=1;//分子
x[0]=39;//乗数
x[1]=2;
for(i=0;i<SIZE;i++) a[i]=x[i];
div(c,x,ans2,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

do{
mul(a,x,pos2,SIZE);
mul(x,pos2,a,SIZE);
add(b,val_2,b,SIZE);
mul(b,a,pos2,SIZE);

div(c,pos2,pos,SIZE);
for(i=0;i<SIZE;i++) c[i]=0;c[SIZE-1]=1;//割り算により分子の内容が破壊されるため再設定

if(flag=(!flag)){
sub(ans2,pos,ans2,SIZE);
}else{
add(ans2,pos,ans2,SIZE);
}
}while(cmp0(pos,SIZE)!=0);

mul(ans,val_4,pos,SIZE);
sub(pos,ans2,pos,SIZE);
mul(pos,val_4,pos2,SIZE);


dp_print(pos2,SIZE,3);//誤差が生じるため下6桁を切り捨てて固定小数点形式で表示
printf("\n");

return 0;
}

//ar+ar2=ar3
void add(const short*ar,const short*ar2,short*ar3,const int size){
int i,p_flag;
for(i=0;i<size;i++){
ar3[i]=ar[i]+ar2[i];
}
for(i=0;i<size;i++){
if(ar3[i]>99){
ar3[i]-=100;
ar3[i+1]+=1;
}else if(ar3[i]<-99){
ar3[i]-=-100;
ar3[i+1]+=-1;
}
}

p_flag=cmp0(ar3,size);
if(p_flag==0) return;

for(i=0;i<size;i++){
if(p_flag>0){
if(0>ar3[i]){
ar3[i+1]-=1;
ar3[i]+=100;
}
}else{
if(0<ar3[i]){
ar3[i+1]-=-1;
ar3[i]+=-100;
}
}
}
}

//ar-ar2=ar3
void sub(const short*ar,short*ar2,short*ar3,const int size){
turn(ar2,size);
add(ar,ar2,ar3,size);
turn(ar2,size);
}

//ar*ar2=ar3
void mul(const short*ar,const short*ar2,short*ar3,int size){
int i,j,pos;
for(i=0;i<size;i++){
ar3[i]=0;
}
for(i=0;i<size;i++){
if(ar2[i]){
for(j=0;j<size;j++){
ar3[j+i]+=ar[j]*ar2[i];
//3桁目を繰り上げる
if(ar3[j+i]>99 || ar3[j+i]<-99){
pos=ar3[j+i]/100;
ar3[j+i]-=(pos*100);
ar3[j+1+i]+=pos;
}
}
}
}
}

// ar/ar2=ar3...ar
// arが書き換えられて余りになる
void div(short*ar,short*ar2,short*ar3,const int size){
int i,pos=0;
int div_flag;/*被除数の符号フラグ*/
int met_flag;/*法の符号フラグ*/

for(i=0;i<size;i++){
ar3[i]=0;
}

div_flag=cmp0(ar,size);
met_flag=cmp0(ar2,size);
/*法、被除数のどちらかが0であるため計算中止*/
if(!div_flag || !met_flag)
return;
if(-1==div_flag) turn(ar,size);/*被除数が負数であるため正数にする*/
if(-1==met_flag) turn(ar2,size);/*法が負数であるため正数にする*/

while(1){
if(cmp(ar,ar2,size)==1 && (ar2[size-1]/10)==0){
pos++;
l_shift(ar2,size);
}else break;
}

while(pos>=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;i<size;i++){
ar[i]=-ar[i];
}
}

/*ar<ar1 の場合は-1を返す*/
/*ar>ar1 の場合は+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;
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;
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<size;i++){
pos=ar[i];
ar[i]=ar[i]/10;
ar[i-1]+=((pos-(ar[i]*10))*10);
}

}

/*左1シフト* \
void l_shift(short*ar,const int size){
int i,pos;
for(i=size-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: 終端されていないコメント

[box]

実行時エラー？　いいえ、コンパイルエラーです。

それから、適切に字下げしてないソースコードは、全くといっていいほど読む気が起きません。

[めるぽん]

>workf.c:196: error: expected identifier or ‘(’ before ‘/’ token
>workf.c:197: error: stray ‘\213’ in program
これは、// で始まるコメントがC言語に無いからじゃないでしょうか。
/* */ の形式に書き直してみて下さい。

>workf.c:334:1: error: 終端されていないコメント
これは
/*左1シフト* \
の後ろ側がバックスラッシュになっているからじゃないでしょうか。
ちゃんと /* */ にしてください。

[HIROYUKI]

Name: box [URL] ぴよぴよさん、めるぽんさんご指摘ありがとうございます
ご指摘されたところ直してみます

[HIROYUKI]

いろいろ触ったら出来ました。
お世話にりました。