求教PI小数点后8000位的程序！
PI=3.1415926.........
哪位能给出源程序，最好是c写的或者告诉我算法！
在此谢过。
问题点数：50、回复次数：18Top
1 楼PRCLionel（lionel）回复于 2002-04-05 21:33:40 得分 15
// This is the beginning of part 1
/*这里有个程序你可以试试，原来有点问题，我修改了一下，在Microsoft
**Virual C++ 6.0下编译通过，由于程序过长，我将它分成几个部分
** Program to compute PI, by Jason Chen, May 1999
** checked by PRCLionel,April 2002
** Compiled & Linked by Microsoft Virual C++ 6.0
**
** Open VC++ IDE, new a win32 console application project,
** add this file to the project and build it.
**
*/
#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:
}
}
}
/* this is the end of part 1 */Top
2 楼PRCLionel（lionel）回复于 2002-04-05 21:34:44 得分 0
// this is the beginning of part 2
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(" checked by PRCLionel,April 2002.\n");
printf("Compiled & Linked by Microsoft Virual C++ 6.0\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");
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);
}
}
/*
a bit changed by PRCLionel april 2002
-------------------------------------------------------
*/
arctg5 = (int*)calloc(BinLen, 4);
arctg239 = (int*)calloc(BinLen, 4);
tmp = (int*)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);
}
}
//this is the end of part 2
//this is the end of fileTop
3 楼NowCan（城市浪人）回复于 2002-04-07 13:48:04 得分 0
天，汇编的啊。Top
4 楼sxbobo2002（五月雪）回复于 2002-04-07 14:11:14 得分 0
呵呵~~汇编不会！
To:PRCLionel(Lionel)
谢谢你的程序！但你能否把程序的设计算法说一下，多谢！Top
5 楼sky_blue（蓝天2007）回复于 2002-04-07 14:15:17 得分 0
用级数可以实现的吧。Top
6 楼sxbobo2002（五月雪）回复于 2002-04-07 14:40:10 得分 0
级数?
pi/4=1-1/3+1/5-1/7....
是这个吗？
说详细一些，多谢！Top
7 楼mathe（）回复于 2002-04-07 14:40:48 得分 0
上面的代码是使用：
Pi=16*arctgn(1/5)-4*arctgn(1/239)
然后对arctgn作泰勒展开。Top
8 楼PRCLionel（lionel）回复于 2002-04-08 15:52:27 得分 0
嘿嘿，已经有人回答了。Top
9 楼qsyang（二拳映月）回复于 2002-04-08 16:09:12 得分 10
同意mathe，那是著名的Machin公式。
你只要写个反正切的函数就行了。函数如下：
double arctan(double x)
{
int I;
double r, e, f, sqr;
sqr=x*x;
r =0;
e =x;
I =1;
while ((double)(e/I) >exp(-？？))
// ??：位数，不过你要8000位，还要加个移位器
{
f=e/I;
r=(I%4==1)?(r+f):(r-f);
e=e*sqr;
I+=2;
}
return r;
}Top
10 楼Viper（）回复于 2002-04-08 16:49:54 得分 0
kao，关注一下。Top
11 楼mengxihe（濛溪河）回复于 2002-04-10 12:45:35 得分 0
好像有个经典的算法就几行呀，怎么你们的程序那么长？
Top
12 楼microblue（）回复于 2002-04-10 13:26:09 得分 0
目前 pi 值已计算到小数点后 206,158,430,000 位。
即便是是在 PC 机上也计算到 12,884,901,372 位。
具体情况请访问以下网站：
http://www.jason314.com/
Top
13 楼steedhorse（晨星）回复于 2002-04-10 14:39:32 得分 0
收藏。Top
14 楼Nizvoo（）回复于 2002-04-10 15:06:08 得分 0
sfi yn .Top
15 楼microblue（）回复于 2002-04-10 18:45:26 得分 15
以下是我写的一个计算 pi 值的程序，可以在 Unix 和 Windows 用编译和运行。
#include <stdio.h>
#include <stdlib.h>
main(int argc, char * argv[])
{
long * pi, * t, m, n, r, s;
int t0[][3] = {48, 32, 20, 24, 8, 4}, k0[][3] = {1, 1, 0, 1, 1, 1};
int n0[][3] = {18, 57, 239, 8, 57, 239}, d, i, j, k, p, q;
d = (argc > 1) ? (((i = atoi(argv[1])) < 0) ? 0 : i) : 1000;
q = (argc > 2) ? 1 : 0;
printf("%s\n\n", "microblue (R) Pi value compute Program (C) Tue 1999.11.30");
printf("pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) [%s]\n",
k0[q][0] ? "" : "-", t0[q][0], n0[q][0], k0[q][1] ? "+" : "-", t0[q][1],
n0[q][1], k0[q][2] ? "+" : "-", t0[q][2], n0[q][2], q ? "Stomer" : "Gauss");
if ((t = (long *)calloc((d += 5) + 1, sizeof(long))) == NULL) return 1;
if ((pi = (long *)calloc(d + 1, sizeof(long))) == NULL) return 2;
for (i = d; i >= 0; i--) pi[i] = 0;
for (p = 0; p < 3; p++) {
for (k=k0[q][p], n=n0[q][p], t[i=j=d]=t0[q][p], i--; i >= 0; i--) t[i] = 0;
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % n;
t[i] = m / n;
k ? (pi[i] += t[i]) : (pi[i] -= t[i]);
}
while (j > 0 && t[j] == 0) j--;
for (k = !k, s = 3, n *= n; j > 0; k = !k, s += 2) {
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % n;
t[i] = m / n;
}
while (j > 0 && t[j] == 0) j--;
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % s;
m /= s;
k ? (pi[i] += m) : (pi[i] -= m);
}
}
}
for (n = i = 0; i <= d; pi[i++] = r) {
n = (m = pi[i] + n) / 10;
if ((r = m % 10) < 0) r += 10, n--;
}
printf("pi= %ld.", pi[d]);
for (i = d - 1; i >= 5; i--)
printf("%ld%s", pi[i], ((m = d - i + 5) % 65) ? ((m % 5) ? "" : " ") : "\n");
printf("%sDIGITS: %d\n", (m % 65) ? "\n" : "", d - 5);
return 0;
}
Top
16 楼microblue（）回复于 2002-04-10 19:05:06 得分 5
据 http://www.jason314.com 网站，目前PC机上流行的最快的圆周率计算程序是PiFast(该程序可能从 http://numbers.computation.free.fr/Constants/PiProgram/download.html 下载)。它除了计算圆周率，还可以计算e和sqrt(2)。PiFast可以利用磁盘缓存，突破物理内存的限制进行超高精度的计算，最高计算位数可达240亿位，并提供基于Fabrice Bellard公式的验算功能。
PC机上的最高计算记录
最高记录 12,884,901,372位
时间 2000年10月10日
记录创造者 Shigeru Kondo
所用程序 PiFast ver3.3
机器配置 Pentium III 1G, 1792M RAM，WindowsNT4.0，40GBx2 (IDE,FastTrak66)
计算时间 1,884,375秒 (21.8天)
验算时间 29小时
而我使用 PiFast version 4.1, 计算 100,000,000 位 pi 值，所用时间为 3820.52 seconds (~ 1.06 hours)。机器配置: Pentium 4 1.6G, 256M RAM, WindowsXP, 40GB(IDE)。
不知那位高手有 PiFast 的源程序？
Top
17 楼PRCLionel（lionel）回复于 2002-04-12 18:56:42 得分 5
To: sxbobo2002 (五月雪)
明白了么？
去
http://www.jason314.com/
好好看看吧。Top
18 楼willa（chocolate boy）回复于 2002-04-12 20:05:43 得分 0
关注Top
