用Machin公式计算圆周率的源程序

类别:编程语言 点击:0 评论:0 推荐:
用Machin公式计算圆周率的源程序

 

/* Program to compute PI, by Jason Chen, May 1999 ** ** Open VC++ IDE, new a win32 console application project, ** add this file to the project and build it. ** ** Homepage : http://www.jason314.com ** Email : [email protected] */

#include <stdlib.h> #include <stdio.h> #include <string.h> #include <malloc.h> #include <math.h> #include <time.h> #include <conio.h> #include <io.h>

int *arctg5, *arctg239, *tmp; int DecLen, BinLen; int x, n, sign, NonZeroPtr; int Step; char fn_status[] = "~pi_sts.___"; char fn_arctg5[] = "~atg5.___"; char fn_arctg239[] = "~atg239.___"; char fn_tmp[] = "~tmp.___"; time_t TotalTime, time1, time2;

void __cdecl FirstDiv(int *arctg_array) { __asm { mov esi, arctg_array mov dword ptr [esi], 1 xor edx, edx mov ebx, x mov ecx, BinLen fd1: mov eax, [esi] div ebx mov [esi], eax add esi, 4 loop fd1 mov esi, arctg_array mov edi, tmp mov ecx, BinLen pushf cld rep movsd popf } }

void __cdecl arctgx(int *arctg_array) { int NonZeroBytePtr; int key; FILE *fp;

for (;NonZeroPtr<BinLen;) { NonZeroBytePtr = NonZeroPtr * 4;

if (_kbhit()) { key=_getch(); if (key==0 || key==0xe0) _getch(); if (key=='p') { printf("Swap data to disk ...\n"); if (x==25) Step = 1; else Step = 2; time(&time2); TotalTime += time2 - time1; if ((fp=fopen(fn_status,"wt"))==NULL) { printf("Create file %s error!!\n", fn_status); exit(0); } fprintf(fp, "%d %d %d %d %d %d %d\n", Step, DecLen, BinLen, n, sign, NonZeroPtr, TotalTime); fclose(fp); if ((fp=fopen(fn_arctg5,"wb"))==NULL) { printf("Create file %s error!!\n", fn_arctg5); exit(0); } if (fwrite(arctg5, 4, BinLen, fp) != (unsigned)BinLen) { printf("Write file %s error!!\n", fn_arctg5); exit(0); } fclose(fp); if ((fp=fopen(fn_arctg239,"wb"))==NULL) { printf("Create file %s error!!\n", fn_arctg239); exit(0); } if (fwrite(arctg239, 4, BinLen, fp) != (unsigned)BinLen) { printf("Write file %s error!!\n", fn_arctg239); exit(0); } fclose(fp); if ((fp=fopen(fn_tmp,"wb"))==NULL) { printf("Create file %s error!!\n", fn_tmp); exit(0); } if (fwrite(tmp, 4, BinLen, fp) != (unsigned)BinLen) { printf("Write file %s error!!\n", fn_tmp); exit(0); } fclose(fp); printf("Exit.\n"); exit(0); } }

__asm { mov esi, tmp add esi, NonZeroBytePtr xor edx, edx mov ebx, x mov ecx, BinLen sub ecx, NonZeroPtr arctg1: mov eax, [esi] div ebx mov [esi], eax add esi, 4 loop arctg1

cmp sign, 1 jne sub_

mov esi, tmp add esi, NonZeroBytePtr mov edi, arctg_array add edi, NonZeroBytePtr xor edx, edx mov ebx, n mov ecx, BinLen sub ecx, NonZeroPtr add_1: mov eax, [esi] div ebx add [edi], eax adc dword ptr [edi-4], 0 jnc add_3 push edi sub edi, 4 add_2: sub edi, 4 add dword ptr [edi], 1 jc add_2 pop edi add_3: add esi, 4 add edi, 4 loop add_1 jmp adj_var

sub_: mov esi, tmp add esi, NonZeroBytePtr mov edi, arctg_array add edi, NonZeroBytePtr xor edx, edx mov ebx, n mov ecx, BinLen sub ecx, NonZeroPtr sub_1: mov eax, [esi] div ebx sub [edi], eax sbb dword ptr [edi-4], 0 jnc sub_3 push edi sub edi, 4 sub_2: sub edi, 4 sub dword ptr [edi], 1 jc sub_2 pop edi sub_3: add esi, 4 add edi, 4 loop sub_1

adj_var: add n, 2 neg sign mov esi, tmp add esi, NonZeroBytePtr cmp dword ptr [esi], 0 jne adj_var_ok inc NonZeroPtr adj_var_ok: } } }

void __cdecl mul_array(int *array, int multiplicator) { __asm { mov esi, BinLen dec esi shl esi, 2 add esi, array mov ecx, BinLen mov ebx, multiplicator xor edi, edi mul1: mov eax, [esi] mul ebx add eax, edi adc edx, 0 mov [esi], eax mov edi, edx sub esi, 4 loop mul1 mov [esi], edx } }

void __cdecl sub2array(int *array1, int *array2) { __asm { mov esi, array1 mov edi, array2 mov ecx, BinLen dec ecx sub1: mov eax, [edi+ecx*4] sbb [esi+ecx*4], eax loop sub1 } }

void main(void) { struct tm *ts; FILE *pi, *fp; int i, tail, p10tail;

printf("\nProgram to compute PI, by Jason Chen, May 1999.\n"); printf(" Dec Length Time(h:m:s)\n"); printf(" 20000 00:00:07\n"); printf(" 100000 00:02:54\n"); printf(" (Running on PII-233, 128MB, Win98 Dos mode)\n"); printf(" Homepage: http://www.jason314.com\n"); printf(" Email: [email protected]\n\n");

if ((fp=fopen(fn_status,"rt"))==NULL) { printf("Decimal length = "); scanf("%d", &DecLen); if (DecLen < 100) DecLen = 100; BinLen = (int)(DecLen / log10(2) / 32) + 2; Step = 0; TotalTime = 0; } else { printf("Reading previous data ...\n"); fscanf(fp, "%d %d %d %d %d %d %d", &Step, &DecLen, &BinLen, &n, &sign, &NonZeroPtr, &TotalTime); fclose(fp); if (Step*DecLen*BinLen*n*sign*NonZeroPtr*TotalTime == 0) { printf("File %s error !!\nExit!\n", fn_status); exit(0); } }

arctg5 = calloc(BinLen, 4); arctg239 = calloc(BinLen, 4); tmp = calloc(BinLen, 4); if (arctg5==NULL || arctg239==NULL || tmp==NULL) { printf("Not enough memory !!\n"); exit(0); }

if (Step == 0) { memset(arctg5, 0, BinLen*4); memset(arctg239, 0, BinLen*4); memset(tmp, 0, BinLen*4); } else { if ((fp=fopen(fn_arctg5,"rb"))==NULL) { printf("Open file %s error!!\n", fn_arctg5); exit(0); } if (fread(arctg5, 4, BinLen, fp) != (unsigned)BinLen) { printf("File %s error!!\n", fn_arctg5); exit(0); } fclose(fp); if ((fp=fopen(fn_arctg239,"rb"))==NULL) { printf("Open file %s error!!\n", fn_arctg239); exit(0); } if (fread(arctg239, 4, BinLen, fp) != (unsigned)BinLen) { printf("File %s error!!\n", fn_arctg239); exit(0); } fclose(fp); if ((fp=fopen(fn_tmp,"rb"))==NULL) { printf("Open file %s error!!\n", fn_tmp); exit(0); } if (fread(tmp, 4, BinLen, fp) != (unsigned)BinLen) { printf("File %s error!!\n", fn_tmp); exit(0); } fclose(fp); }

printf("Working ......\n"); printf("Press 'p' to pause & exit\n"); time(&time1);

if (Step == 0) { x = 5; FirstDiv(arctg5); x = x * x; n = 3; sign = -1; NonZeroPtr = 1; arctgx(arctg5); } else if (Step == 1) { x = 5 * 5; arctgx(arctg5); }

if (Step == 0 || Step == 1) { x = 239; FirstDiv(arctg239); x = x * x; n = 3; sign = -1; NonZeroPtr = 1; arctgx(arctg239); } else { x = 239 * 239; arctgx(arctg239); }

mul_array(arctg5, 16); mul_array(arctg239, 4); sub2array(arctg5, arctg239);

if ((pi=fopen("pi.txt","wt"))==NULL) { printf("Create file pi.txt error!!\n"); exit(0); }

printf("Writing result to file: pi.txt ...\n"); fprintf(pi, "%d", arctg5[0]); for (i=1; i<=DecLen/9; i++) { arctg5[0] = 0; mul_array(arctg5, 1000000000); fprintf(pi, "%09d", arctg5[0]); } tail = DecLen % 9; p10tail = 1; for (i=1;i<=tail;i++) p10tail *= 10; arctg5[0] = 0; mul_array(arctg5, p10tail); fprintf(pi, "%0*d", tail, arctg5[0]); fclose(pi);

time(&time2); TotalTime += time2 - time1; printf("Done !!\n"); ts = gmtime(&TotalTime); ts->tm_mon --; ts->tm_mday = ts->tm_mday - 1 + ts->tm_mon * 31; printf("Time : "); if (ts->tm_mday > 0) printf("%d Day(s) ", ts->tm_mday); printf("%02d:%02d:%02d\n", ts->tm_hour, ts->tm_min, ts->tm_sec);

if (_unlink(fn_status) == 0) { _unlink(fn_arctg5); _unlink(fn_arctg239); _unlink(fn_tmp); } }

 

本文地址:http://com.8s8s.com/it/it27739.htm