精确算出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
...全文
407 22 打赏 收藏 举报
写回复
22 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
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;

}
  • 打赏
  • 举报
回复
加载更多回复(2)
相关推荐
发帖
数据结构与算法

3.2w+

社区成员

数据结构与算法相关内容讨论专区
社区管理员
  • 数据结构与算法社区
加入社区
帖子事件
创建了帖子
2001-01-20 10:59
社区公告
暂无公告