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

## 用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    : jason@szonline.net
*/
#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: jason@szonline.net/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);
}
}

• 本文已收录于以下专栏：

举报原因： 您举报文章：用Machin公式计算圆周率的源程序 色情 政治 抄袭 广告 招聘 骂人 其他 (最多只允许输入30个字)