精确算出800位pi的程序!

Wuxuehui 2001-01-20 10:59:00
#include <stdio.h>

long a=10000,b,c=2800,d,e,f[2801],g;

main(){

for(;b-c;)f[b++]=a/5;

for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)

for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);

}

在gcc,C++ Builder 5编译通过.no error and no warnning.

中华技术网 http://www.asfocus.com
Wuxuehui oicq:178927,icq 2759606
...全文
357 点赞 收藏 22
写回复
22 条回复
yug 2001年06月05日
Wuxuehui给出的程序为什么这么恶心???
我觉得,恐怕其唯一的目的就是不想让人看懂,其心可诛!
我把程序整理了一下,供大家参考:

#define MAX 10000
#define NUM 2800
#define STEP 14
long b,c,d,e,g
long f[NUM+1];
main(){
int i;
for(i=0;i<NUM;i++)f[i]=MAX/5;
f[MAX]=0;
for(c=NUM;c>0;c-=STEP);
g=c*2;
d=0;
for(b=c;b>0;b--){
d=d*b+f[b]*MAX;
g=g-1;
f[b]=d%g;
d=d/g;
g=g-1;
}
printf("%04d"),e+d/a);
e=d%MAX;
}
}

我用MathMatic算过,结果是正确的,速度也不错,800位不到十秒(MathMatic算50000位用了一分钟),很有意思,但这个算法的本质我还没看明白,有明白的不妨说一声。

回复 点赞
Go_Rush 2001年06月02日
同意楼上!
程序没有一点注释,
变量命名不规范,
可读性奇差, 就算能解决问题,也没有什么用!

TO: Wuxuehui
你贴出来的程序是用来帮助别人解决问题的,而不是用来炫耀你的能力的


回复 点赞
starfish 2001年06月02日
这个公式
PI=16arctg(1/5)-4arctg(1/239)
的展开形式是我所知道的最好的计算方法了,收敛速度很快,而且可以算到任意的精度。

TO: Wuxuehui
这种程序就不要贴出来了,这么恶心的程序怎么好意思贴出来的呀,完全是卖弄技巧,可读性非常差!C语言的名声就是给这种程序弄坏了。
回复 点赞
yug 2001年06月02日
这个算法很好
回复 点赞
kiss2 2001年05月09日


我要DELPHI代码,谁能翻译下???谢谢!


回复 点赞
duz 2001年05月09日
Pi c!
--- = Sum{ --------- , 0<=c<=+oo}
2 (2c+1)!!
回复 点赞
duz 2001年05月09日
哈哈,很早就看到这个东东了。
看懂这个计算Pi的c程序我足足花了一个小时,主要使用在数学推导上,用的公式挺巧的,我记得用这个程序计算到2000为左右还是没有问题的,不过要继续提高精度就要求能够支持更高位数的整数了,计算时间复杂度位O(n^2),n为位数,计算时间复杂度并不低,不过前面的常系数非常小。在对Pi的位数要求不是太高时,这应该是可以找到的最好的方法了。
回复 点赞
HuaRoaD 2001年05月09日
我不太懂C ,有哪位朋友能够帮个忙把上面这段C程序翻译成Pascal。
(我改了半天都没改好)
回复 点赞
HuaRoaD 2001年05月03日
简直不明白,这样的程序居然能够算出Pi,基于什么数学原理,根本不可能!
回复 点赞
iamfool 2001年03月04日
为什么要把三重的FOR这样写,好看么,?????
回复 点赞
iamfool 2001年03月04日
为什么要把三重的FOR这样写,好看么,?????
回复 点赞
vonvon 2001年02月27日
Nowcan你是个酷哥,我喜欢
这年头执著于汇编的人可不多了,我也算是一个
回复 点赞
BCB 2001年02月27日
include <stdio.h>

long a=10000,b,c=2800,d,e,f[2801],g;

main(){

for(;b-c;)f[b++]=a/5;

for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)

for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);

}

这个程序好象是对的,算法水平很高,可编出的程序
却晦湿却懂,居然是对的,我就想不通了,是不是水
平高得太离奇了!



回复 点赞
yudi1226 2001年02月23日
请问:它们是外星人吗?
回复 点赞
hyqryq 2001年02月02日
可以贴到文档中心!
回复 点赞
NowCan 2001年02月02日
概述:
用汇编语言编制计算程序并不是强项,特别是在涉及到浮点计算时,但汇编的一个好处就是速度快,所以在整数计算时可以试一下。本文的理论基础来自是电脑杂志1996年第10期,作者郭继展发表的一篇文章,作者提出一个公式:
PI=16arctg(1/5)-4arctg(1/239)
在展开成两个级数之和,然后整理得到:
PI=16x(1/5-1/(5^3/3)+1/(5^5/5)-1/(5^7/7)+...)-4x(1/239-1/(239^3/3)+1/(239^5/5)-1/(239^7/7)+...)
=4x(4x5/25-239/57121)/1-4x(4x5/25^2-239/57121^2)/3+4x(4x5/25^3-239/57121^3)/5-...
我对以上公式和推导一看就头疼,但根据它编出的程序却可以在4分钟内算出圆周率的小数点下8万位!(在P5/200上)想当年祖冲之算了一生才算到3.14159265,十九世纪英国人香克思用了一生才算到小数点下707位。
本程序的难点就是如何达到小数点下这么多位的精度,这个办法就是:在计算机中一个 WORD 可以表示0到65535,我们可以在内存定义一个字来表示五位数,如果要算到小数点下10000位,则定义2000个字来表示它,如计算239/57121时,可以用23900000/57121,得到小于五位的结果存到第一个字中,然后用余数乘以100000再除57121,得到小于五位的结果存到第二个字中,依此类推。为了计算时不至于溢出,本程序动用一个双字来表示五位数,再用一个段64K来表示一个高精度数,共可以表示(65536/4)*5 共有小数点下 81920 位。一共用到三个段,第一个段存储(4*4*5/25^n),第二个段存储(4*239/57121^n),第三个段存储最后的结果即 PI。
本程序的意义就在于提出了一个表示高精度数的办法,如果内寸足够的话,理论上可以进行任何精度的计算。这里是编译好的可执行文件:pi.com,计算结果8万位的PI值:pi.txt,还有本程序要用到的两个公用子程序 scanf.asm 和 printf.asm,用于读入键盘数字输入和屏幕数字输出。(引用本程序时请注明出处,并请勿改动版权信息)

源程序:

.386
CODE SEGMENT USE16
ASSUME CS:CODE,DS:CODE
ORG 100H
start:
jmp install

HANDLE DW ?
_MIN DW ?
_SEC DW ?
_SEC1 DW ?
_A DD ?
_B DD ?
_FS_POINT DW 0
_GS_POINT DW 0
_DIV DD 1
FLAG DB ? ;=1 +
DIGITAL DD 5000 ;how many points want to calculate
POINT DW ? ;total point /5
NUM_POINT DW ? ;total point /5 * 4
_COUNT DD ?
TMP_NUM0 DD ?
TMP_NUM DD 10 dup (?)
KEY_BUFF DB 6,0,6 dup (0)

D_SAVE DB 'Saving result to file %c ...',0
DW 80h
D_SAVE_OK DB 8,8,8,8,', OK !',0dh,0ah,0
D_SAVE_ERR DB 8,8,8,8,', error !',0dh,0ah,0

D_TIME DB 0dh,0ah,'Total %ld points, time used %d:%02d.%03d, '
db 'calculate %ld times.',0dh,0ah,0
dw digital,_min,_sec,_sec1,_div
D_SCAN DB '<<PI calculater>> Dec 18, 1996',0dh,0ah
DB 'Copyright(C) by Luo Yun Bin, phone 0576-4114689',0dh,0ah,0ah
DB 'How many points (10-80000): ',0
D_ABORT DB 'User pressed Esc, calculate aborted !%20r ',0dh,0ah,0
D_CAL DB 'Calculating, please waiting ... (Esc to cancel)',0dh,0
D_CAL1 DB '%d %% calculated, please waiting ... (Esc to cancel)',0dh,0
DW PERCENT
PERCENT DW ?
D_STR1 DB ' PI = %1ld.%c',0dh,0ah,0
DW tmp_num0,d_sub_str
D_STR2 DB '%5ld : %c'
DB 0dh,0ah,0
DW _count,d_sub_str
D_SUB_STR DB '%05ld %05ld %05ld %05ld %05ld '
DB '%05ld %05ld %05ld %05ld %05ld',0
DW tmp_num,tmp_num+4,tmp_num+8,tmp_num+12,tmp_num+16
DW tmp_num+20,tmp_num+24,tmp_num+28,tmp_num+32,tmp_num+36

install:
mov si,offset d_scan
call printf
mov ah,0Ah
mov dx,offset key_buff
int 21h
mov si,offset key_buff+2
call scanf
mov eax,dword ptr scan_num
mov digital,eax
mov si,offset d_cal
call printf

xor ax,ax
mov ds,ax
mov ax,ds:[046ch]
push cs
pop ds
mov _sec,ax

mov ax,cs
add ax,1000h ;result of 4*4*5/25^n
mov fs,ax
add ax,1000h ;result of 4*239/57121^n
mov gs,ax
add ax,1000h ;total result
mov bp,ax

mov ax,fs
call init_num
mov dword ptr fs:[4],4*239*100000

mov ax,gs
call init_num
mov dword ptr gs:[4],4*4*5*100000

mov ax,bp
call init_num

call pre
call calc

xor ax,ax
mov ds,ax
mov ax,ds:[046ch]
push cs
pop ds
mov _sec1,ax

push point
call num_out
pop point

mov ax,_sec1
sub ax,_sec
mov cx,55
mul cx
mov cx,1000
div cx
mov _sec1,dx
mov cx,60
xor dx,dx
div cx
mov _min,ax
mov _sec,dx
mov si,offset d_time
call printf

mov si,81h
mov di,80h
cmd_lop:
lodsb
cmp al,0dh
jz cmd_end
cmp al,20h
jbe cmd_lop
cmp al,'a'
jb cmd_store
cmp al,'z'
ja cmd_store
sub al,20h
cmd_store:
stosb
jmp short cmd_lop
cmd_end:
xor al,al
stosb

mov si,80h
cmp byte ptr ds:[si],0
jz quit

mov si,offset d_save
call printf

mov ah,3ch
xor cx,cx
mov dx,80h
int 21h
jb file_err

mov handle,ax
mov std_out,offset write_file
call num_out

mov std_out,offset prt_to_scr
mov si,offset d_save_ok
call printf

mov ah,3eh
mov bx,handle
int 21h
int 20h
file_err:
mov si,offset d_save_err
call printf
int 20h
quit:
int 20h

WRITE_FILE PROC

pusha
mov ah,40h
mov flag,al
mov bx,handle
mov cx,1
mov dx,offset flag
int 21h
popa
ret

WRITE_FILE ENDP

PRE PROC

mov eax,digital ;total 65536*5/4 points
cmp eax,(65536/4-3)*5 ;comp max points
jbe pre_1
mov eax,(65536/4-3)*5
pre_1:
cmp eax,10 ;comp min points
jae pre_2 ;must > 10 and < 81915
mov eax,10
pre_2:
xor edx,edx
mov ecx,5
div ecx
mov ebx,eax
inc ebx
mov point,bx ;point for print control

mul ecx
mov digital,eax ;

mov eax,ebx
inc eax
mov ecx,4
mul ecx
mov num_point,ax ;max used memory

ret

PRE ENDP


CALC PROC

mov es,bp
c_lop0:
mov ah,1
int 16h
jz calc_0
xor ah,ah
int 16h
cmp al,1bh
jnz calc_00
push cs
pop es
mov si,offset d_abort
call printf
int 20h
calc_00:
xor eax,eax
mov ax,_gs_point
mov ecx,500
mul ecx
mov ecx,4
div ecx
xor edx,edx
div digital
mov percent,ax
mov si,offset d_cal1
push cs
pop es
call printf
mov es,bp
calc_0:
xor eax,eax
mov ecx,100000
mov _a,eax
mov _b,eax
xor flag,00000001b ;init part
;============================================================
mov ebx,57121
c_lop1:
mov si,_fs_point ;if 4*5/25^n = 0 skip
cmp si,num_point
jae calc_1
cmp dword ptr fs:[si],0
jnz c_lop2
add _fs_point,4
jmp short c_lop1
c_lop2:
mul ecx
add eax,fs:[si]
adc edx,0
div ebx
mov fs:[si],eax
mov eax,edx

add si,4
cmp si,num_point
jb c_lop2
;=======================================================
calc_1:
mov ebx,25
c_lop3:
mov si,_gs_point ;if 4*5/25^n = 0 skip
cmp si,num_point
jae calc_4
cmp dword ptr gs:[si],0
jnz c_lop4
add _gs_point,4
jmp short c_lop3
c_lop4:
mov eax,_a
mul ecx
add eax,gs:[si]
adc edx,0
div ebx
mov gs:[si],eax
mov _a,edx

mov eax,_b
mul ecx
add eax,gs:[si]
adc edx,0
sub eax,fs:[si]
sbb edx,0

cmp edx,0
jl _t1
div _div
mov _b,edx
jmp short _t2
_t1:
idiv _div
dec eax
add edx,_div
mov _b,edx
_t2:
test flag,00000001b
jnz calc_2
sub es:[si],eax
jmp short calc_3
calc_2:
add es:[si],eax
calc_3:
add si,4
cmp si,num_point
jb c_lop4
c_lop5:
add _div,2
jmp c_lop0
;====================================================
calc_4:
xor eax,eax
mov _a,eax
mov si,num_point
c_lop6:
sub si,4
mov eax,es:[si]
add eax,_a

cmp eax,ecx
jge calc_5
cmp eax,0
jle calc_6
mov _a,0
mov es:[si],eax
jmp short calc_8
calc_5:
xor edx,edx
div ecx
mov es:[si],edx
mov _a,eax
jmp short calc_8
calc_6:
mov edx,-1
idiv ecx
dec eax
add edx,ecx
mov es:[si],edx
mov _a,eax
calc_8:
cmp si,4
ja c_lop6

mov eax,_a
mov es:[0],eax

push cs
pop es
ret

CALC ENDP

INIT_NUM PROC

push es
mov es,ax
xor eax,eax
cld
xor di,di
mov cx,65536/4
cld
rep stosd
pop es
ret

INIT_NUM ENDP

NUM_OUT PROC

cld
xor esi,esi
mov _count,esi
mov di,offset tmp_num0
mov cx,11
cmp cx,point
jbe no_adjust
mov cx,point
no_adjust:
sub point,cx
mov ds,bp
rep movsd
push cs
pop ds

push si
mov si,offset d_str1
call printf
pop si
no_lop:
cmp point,0
jle no_quit
mov cx,10
mov di,offset tmp_num
xor eax,eax
rep stosd

mov cx,10
mov di,offset tmp_num
cmp point,10
jae no1
mov cx,point
no1:
sub point,cx
add _count,50
mov ds,bp ;bp = result segment
rep movsd
push cs
pop ds

push si
mov si,offset d_str2
call printf
pop si
jmp short no_lop
no_quit:
ret

NUM_OUT ENDP

INCLUDE SCANF.ASM
INCLUDE PRINTF.ASM

CODE ENDS
END START



--------------------------------------------------------------------------------




--------------------------------------------------------------------------------


(C) Copyright by LuoYunBin's Win32 ASM Page,http://asm.yeah.net/
回复 点赞
Wuxuehui 2001年02月02日
要啊,快贴出来看看 ,谢谢 :)
回复 点赞
leojay 2001年02月02日
我有一个4分钟以内算8万位的汇编程序你要不要?
回复 点赞
hyqryq 2001年01月21日
这个贴子有人贴过了!
回复 点赞
Wuxuehui 2001年01月20日
加上return,就不会有那个warnning!

#include <stdio.h>

long a=10000,b,c=2800,d,e,f[2801],g;

main(){

for(;b-c;)f[b++]=a/5;

for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)

for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);

return;

}
回复 点赞
发动态
发帖子
数据结构与算法
创建于2007-08-27

3.0w+

社区成员

3.4w+

社区内容

数据结构与算法相关内容讨论专区
社区公告
暂无公告