多倍長演算で・・・

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら　　※コードを貼り付ける場合は [code] と [/code] で囲って下さい。詳しくはこちら
閉鎖
アバター
lbfuvab
記事: 72
登録日時: 13年前

多倍長演算で・・・

#1

投稿記事 by lbfuvab » 13年前

ネイピア数を100万桁求めてみようと思い立って、プログラムを組んでみたのですが、うまく動きません。
具体的には
LongInt_t DivLongInt(LongInt_t x,uint32 y)がどこかおかしいようで、小数点以下240ケタくらいからゴミが入ります。
変な場所を参照しているのかと思い、色々修正してみたのですが、手詰まりになりました・・・

環境はVC++ 2010 Express Editon を Win7 Proで動かしています。

コードは以下です(main関数のコメントアウト部分が本来のプログラムであり、下はデバッグ用です)

コード: 

#include<stdio.h>
#include<string.h>
#include<stdlib.h>

#define N_BUF	(100)				//10の倍数を推奨
#define S_PLUS	(1)					//0の符号は未定(どちらでもおｋだが、+を推奨)
#define S_MINUS	(-1)

#define ABS(x)	((x)<0?(-(x)):(x))

typedef unsigned short uint16;
typedef unsigned long uint32;
typedef unsigned __int64 uint64;

typedef struct{
	int sign;
	uint16 buf[N_BUF];
}LongInt_t;

//常にxは左辺
LongInt_t	ConvLongInt(int src);
int			CompLongInt(LongInt_t x,LongInt_t y);		//x>yで＋、x==yで0、x<yで-を返す
void		PrintLongInt(LongInt_t x,FILE *fp);

LongInt_t	AddLongInt(LongInt_t x,LongInt_t y);
LongInt_t	SubLongInt(LongInt_t x,LongInt_t y);
LongInt_t	MulLongInt(LongInt_t x,uint16 y);
LongInt_t	DivLongInt(LongInt_t x,uint16 y);
LongInt_t	DivLongInt(LongInt_t x,uint32 y);


LongInt_t ConvLongInt(int src){
	LongInt_t Ret={0,{0}};

	Ret.sign = (src<0)?S_MINUS:S_PLUS;
	Ret.buf[0]=(uint16)ABS(src);

	return Ret;
}

int	CompLongInt(LongInt_t x,LongInt_t y){
	int i;

	if(x.sign != y.sign)
		return x.sign;

	if(x.sign == S_MINUS){
		x.sign = y.sign = S_PLUS;
		return CompLongInt(y,x);
	}

	for(i=0;i<N_BUF;i++)
		if(x.buf[i] ^ y.buf[i])		//xとyが違っていたら・・・
			return (int)x.buf[i] - y.buf[i];

	return 0;
}

void PrintLongInt(LongInt_t x,FILE *fp){
	int i,nDigit=N_BUF*16*3/10;						//全体の桁数(0.3はlog2の事)
	char StrBuf[(N_BUF+1)*16*3/10];						//全体の桁数+

	sprintf(StrBuf,"%s%d.",(x.sign==S_MINUS)?"-":"",x.buf[0]);
	
	for(i=strlen(StrBuf);i<nDigit;i++){
		x.buf[0]=0;
		x=MulLongInt(x,10);
		StrBuf[i]=x.buf[0]+'0';
	}
	StrBuf[i]='\0';
	while(i--){
		if(StrBuf[i]=='0' && StrBuf[i-1]!='.')
			StrBuf[i]='\0';
		else
			break;
	}
	fputs(StrBuf,fp);
}

LongInt_t AddLongInt(LongInt_t x,LongInt_t y){		//ret = x+y
	LongInt_t ret;
	int i;
	uint32 tmp=0;

	//xとyの符号が違うならyの符号を反転して引き算へ
	if(x.sign != y.sign){
		y.sign *= -1;
		return SubLongInt(x,y);
	}
	//加算処理ここから
	ret.sign = x.sign;

	i=N_BUF;
	while(i--){
		tmp += x.buf[i] + y.buf[i];
		ret.buf[i] = tmp & 0xffff;
		tmp >>= 16;
	}
	if(tmp)					//tmpが残っているという事は桁あふれなので
		puts("加算オーバーフロー");

	return ret;
}

LongInt_t SubLongInt(LongInt_t x,LongInt_t y){
	LongInt_t ret={S_PLUS,{0}};
	int i,carry=0,tmp=0;

	//xとyの符号が違うならyの符号を反転して足し算へ
	if(x.sign != y.sign){
		y.sign *= -1;
		return AddLongInt(x,y);
	}
	// x=-X,y=-Yならx-y=Y-X
	if(x.sign == S_MINUS){
		x.sign = y.sign = S_PLUS;
		return SubLongInt(y,x);
	}

	if(CompLongInt(x,y)<0){
		ret = SubLongInt(y,x);
		ret.sign = S_MINUS;
		return ret;
	}

	i=N_BUF;
	while(i--){
		tmp = (int)x.buf[i] - y.buf[i] - carry;
		if(tmp<0){
			carry=1;
			tmp += 0x10000;
		}
		else{
			carry=0;
		}
		ret.buf[i] = tmp;
	}
	return ret;
}

LongInt_t MulLongInt(LongInt_t x,uint16 y){
	LongInt_t ret;
	uint32 tmp=0;
	int i=N_BUF;

	ret.sign=x.sign;
	while(i--){
		tmp += (uint32)x.buf[i]*y;
		ret.buf[i] = tmp & 0xffff;
		tmp >>= 16;
	}
	if(tmp)
		puts("乗算オーバーフロー");
	return ret;
}

LongInt_t DivLongInt(LongInt_t x,uint16 y){
	LongInt_t ret;
	uint32 tmp=0;

	if(!y)
		return x;

	ret.sign = x.sign;

	for(int i=0;i<N_BUF;i++){
		tmp += x.buf[i];
		ret.buf[i]=tmp / y;
		tmp %= y;
		tmp <<= 16;
	}
	return ret;
}

LongInt_t DivLongInt(LongInt_t x,uint32 y){
	LongInt_t ret;
	uint64 tmp=0;

	if(!y)
		return x;

	ret.sign = x.sign;

	for(uint32 i=0;i<N_BUF/2;i += 2){
		tmp += x.buf[i]<<16 | x.buf[i+1];
		ret.buf[i]=((tmp / y)>>16) & 0xffff;
		ret.buf[i+1]=(tmp /y) & 0xffff;
		tmp %= y;
		tmp <<= 32;
	}
	return ret;
}

int main(){
	LongInt_t e,tmp;
	uint32 i;
	FILE *fp;

	//system("cd");
	//if((fp=fopen("Napier.txt","w"))==NULL)
	//	return -1;

	e=tmp=ConvLongInt(1);

	for(i=1;i<20;i++){
		tmp=DivLongInt(e,i);
		//e=AddLongInt(e,tmp);
		PrintLongInt(tmp,stdout);
		putchar('\n');
	}
	//PrintLongInt(e,fp);
	//fclose(fp);
	return 0;
}
実行例は以下の通りです。

1.000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000002228018326060731476931755622220920839659146703827144034231019965588568740
10212907261538280864933615857422729622297613931439156398521574154113396854520809
58952311155483226586851536505191374526896887441135276579672413529718296554336239

0.500000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000018
50054219208735440092783922153277948415904559524257908962280477347280915748313292
84514197200287959920918214676609867494011223208501564297295661034882014630472924
21695745651570773474214944797136836187481892642367499595838469382426494492422711

0.333333333333333333333333333333333333333333333333333333333333333333333333333333
33333333333333333333333333333333333333333333333333333333333333333333333333333333
33333333333333333333333333333333333333333333333333333333333333333333333333333319
07249872693266431595145726673514914762740235366717861841575465378137627444008503
43186972991444697560958876186779893806699682110113377520834594618945113722343787
58276196060247528780261965980912235236336985481284329058042844767780922760075913

0.250000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000019
65682607909281405098582917287857820191898594494524028272423007181485972982582873
64796334525305957415975603093897984212386924659032912065876639849562140544877481
98051729754793946816353373199849798247500195413585023420510951837530253872009942

0.199999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999998
84371611299454034994201004865420128224005965029733880689857470165794942765730419
19717862674982002504942611582711883281624298549468652231419021185319874085595442
23644015896776826657843157507402075616795042358583495457701344053318539163110121

0.166666666666666666666666666666666666666666666666666666666666666666666666666666
66666666666666666666666666666666666666666666666666666666666666666666666666666666
66666666666666666666666666666666666666666666666666666666666666666666666666666623
11330692279435318114904516597491496437558016119976172651298042911609510842512456
42706160757655427686171702948814270274515245363319234050116464647048590557428324
23924598778593804112724272372826520465246552995716940561251486579034148912078257

0.142857142857142857142857142857142857142857142857142857142857142857142857142857
14285714285714285714285714285714285714285714285714285714285714285714285714285714
28571428571428571428571428571428571428571428571428571428571428571428571428571421
30335842453711077106406316296926520265180351615470107193389812470711068813162634
95369422528458301459639272805617552055924162310945814026062418879153494252314208
3433381420831148184938924833151310234939406089465744042496454463163775930379608
0.125000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000021
57496178268385520300582486910583169542981792191903875596605594834599199171986315
76248627296291378941374500755110856622656412476035269825542961737088452914773383
05288724484775803213211736670233186119681263034645807487992540179908066892721323

0.111111111111111111111111111111111111111111111111111111111111111111111111111111
11111111111111111111111111111111111111111111111111111111111111111111111111111111
11111111111111111111111111111111111111111111111111111111111111111111111111111053
16844077339307753598294799367164203225187803156664465680635449419279909704935450
88084007379653681081013091533673262223617627312262461818886506064362579178171605
40383463271816535854991327559279433968143703836394629894884148167659972272192399

0.099999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999959
91549191714406546465634835334564445098873454364107863915058965747558015878654532
16885906066042753504677201534011953762975683048246610689192734424422301633975330
86325884421596657472479059062170208495427037420711176255349642471987248688371208

0.090909090909090909090909090909090909090909090909090909090909090909090909090909
09090909090909090909090909090909090909090909090909090909090909090909090909090909
09090909090909090909090909090909090909090909090909090909090909090909090909090886
38569821880188323522492459175522517853208040583865293546292141437427967036160958
78096212526919321915236736533251526257349862425929898362409870548099345680783229
37009766700344959827057986164969448343247785335181393447716376516578438718269081

0.083333333333333333333333333333333333333333333333333333333333333333333333333333
33333333333333333333333333333333333333333333333333333333333333333333333333333333
33333333333333333333333333333333333333333333333333333333333333333333333333333322
52785291066882389277844669332492897346255483838174591920133686828213440502059235
16948379472432332098717920173072787346437186833880206650290593567094283770058574
98515090189723782191788475760225678279504168100228898563826959242147866661035089

0.076923076923076923076923076923076923076923076923076923076923076923076923076923
07692307692307692307692307692307692307692307692307692307692307692307692307692307
69230769230769230769230769230769230769230769230769230769230769230769230769230766
20664214167802860754064166571098796893554062385457841804242614279771388771910327
13108256226876006551388368754784925505786626975533398143711286176990439916172695
07992041576185546643616695557848151435168960638806076210339125089890213434556712

0.071428571428571428571428571428571428571428571428571428571428571428571428571428
57142857142857142857142857142857142857142857142857142857142857142857142857142857
14285714285714285714285714285714285714285714285714285714285714285714285714285722
20998368237056754163403530932202487048874263674856350184609198659990636724784147
26636636342281365404898888967870475824333008144974726665832964990422358797972594
12424138916256858583277946765423973857667222605077357820527413737813611784699077

0.066666666666666666666666666666666666666666666666666666666666666666666666666666
66666666666666666666666666666666666666666666666666666666666666666666666666666666
66666666666666666666666666666666666666666666666666666666666666666666666666666680
15216964521946064223493789918669162389992876547258921920689557914307405332007402
41348742423098913088206977809672706833660908916845338197544616369396480407263561
05651525491181981811783573646713698077400547439034896755277869734968989407819252

0.062500000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000020
81310996609827370104381912422437691967892629464790147582565537015691030216852454
45078471850323954911032991511186100930762626109564259834457618664242266459282039
74407713858017120158491801602562760307518360120607676425726406853561138177393091

0.058823529411764705882352941176470588235294117647058823529411764705882352941176
47058823529411764705882352941176470588235294117647058823529411764705882352941176
47058823529411764705882352941176470588235294117647058823529411764705882352941189
18971099235417379769671299421555060124169678790574371235097239940961511929906565
30162334104609737151513625531345754490368010073491884277920178726187267411391311
86974648664866671469399713241332494843162372506298374426585232602886576594495924

0.055555555555555555555555555555555555555555555555555555555555555555555555555555
55555555555555555555555555555555555555555555555555555555555555555555555555555555
55555555555555555555555555555555555555555555555555555555555555555555555555555489
38591105443603910738232867868287633923621624179710616219508863918022577339567248
13438786704423775518340172381871420242534352880665401521359545739615431668262125
61474238883986131166016941350402968618470002604434081825611909701947315598351557

0.052631578947368421052631578947368421052631578947368421052631578947368421052631
57894736842105263157894736842105263157894736842105263157894736842105263157894736
84210526315789473684210526315789473684210526315789473684210526315789473684210458
80503901834789972994753406879266433850734238929023770804309837048518763585102369
96860473374387795376298444688862237916590624081257449581798289556921771173011083
55917279025847002058922750751090651102859929322475177614188636610452645216220996

続行するには何かキーを押してください . . .
ページトップ
アバター
Justy
副管理人
記事: 122
登録日時: 13年前
住所: 神奈川県

Re: 多倍長演算で・・・

#2

投稿記事 by Justy » 13年前

for文の N_BUFを2で割っているのが原因では？
ページトップ
名前忘れた・・・・

Re: 多倍長演算で・・・

#3

投稿記事 by 名前忘れた・・・・ » 13年前

ゴミって絶対はいるのでは？関係ない話だったら申し訳ないが・・・・・

例えば10．0/3．0をプログラムで実行すると

3.33333325386047366330000000000000000

となるわけだし。丸め誤差だっけ？そんなものでないのか？
ページトップ
アバター
lbfuvab
記事: 72
登録日時: 13年前

Re: 多倍長演算で・・・

#4

投稿記事 by lbfuvab » 13年前

>>Justyさん
それが原因みたいです、お騒がせしました。
ありがとうございましたm(_ _)m

>>名前忘れた・・・さん
ケタ数を馬鹿みたいにデカくすると偶に下の方のケタが変なことになるのはそういう事ですかねえ。。。

後、
n!=10m
をmを定数として簡単に求める方法って、
スターリングの近似を使って
2 * √π * nn+1/2 / en = 10m
としてから対数を取るくらいしかないんですか？
ページトップ
閉鎖

“C言語何でも質問掲示板” へ戻る