用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);
}
}
[url=http://www.jason314.com/palgorithm.htm][/url]
