vier Gourdon坛子上的那个算是好的开始的计算pi程序，piclassic.c的程序，我下载了，/o2编译，发现与我在没有改动思路的情形下，仍然比其快一点，当然不能和上面的程序相比。它用时770ms ，我重写了一下计算10000位用时465ms。是重写而不是改写，可见/o2有时也不如人工调优，呵呵。
我写这个程序，其效率不是主要学习的东东，主要是体会一下多精度计算的一种思想，用多字节模拟大数据，当然有些3两行的代码或许也可计算pi
但很难体会不到这一点。
像3.14159265......可以定义成以下方式的数据
X = x(0) + x(1)/B^1 + ... + x(n-1)/B^(n-1) ，基为B (可见，x(0)是唯一的正整数部分)
比如当B=10的时候，有
3 1 4 1 5 9 2 6 5 ...分别对应
x(0) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8)
当B=100的时候，有
3 14 15 92 65...分别对应
x(0) x(1) x(2) x(3) x(4)...
...
当B=10000的时候，有
3 1415 9265...分别对应
x(0) x(1) x(2).....
可以看出当x(i)同时为dword型的时候，对于pi的数据存储，基10000比基10需要更少的存储空间。
对于pi的计算，有很多公式可以实现，我们选一个简单而效率相对较高的Guass公式：
Pi/4=12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
更多的公式可以参考：http://en.wikipedia.org/wiki/Pi
对于效率而言，当前Chudnovsky+FFT用来计算pi是无敌的，但若没有艰深的数学功底，FFT是无法透测理解的，而高效的大数库通常有快速傅里叶变换(FFT)的大数乘。所以不谈效率，若再经过简单的优化仍然很快，比如有些除法是不需要的，但作为除法的学习也就用上了，或用mpich2实现多核计算。
哦，不说多了.......不能浮躁.......
效果图：
以下程序将pi计算到小数点后10000位（当然，稍加修改可以得到更多位）。
.....................................................................
;************************************************************
;*--==--* Simple asm program to compute Pi with many digits.
;*--==--* By G-Spider 2010.11.27
;*--==--* Web:http://blog.youkuaiyun.com/G_Spider
;*--==--* ---------------------------------------------------
;*--==--* ml /c /coff pi.asm
;*--==--* link /subsystem:console pi.obj
;*--==--* 请保持以上完整性!!
;************************************************************
.386
.model flat,stdcall
option casemap:none
include windows.inc
include user32.inc
include kernel32.inc
include msvcrt.inc
includelib user32.lib
includelib kernel32.lib
includelib msvcrt.lib
;***********************************************************
.data
num dd 0
flag dd 0
fmt db '%.4d ',0
fmt1 db '%d.',0dh,0ah,0
fmt2 db ' :%d',0dh,0ah,0
fmt3 db 0dh,0ah,0
szPause db 'Pause',0
;NbDigits dd 10000
B dd 10000
LB dd 4
MaxDiv dd 450
_size dd 2502 ;2+NbDigits/LB(决定精度)
p dd 18,57,239
m dd 12,8,-5
dwlm dd 1000 ;1000是毫秒为单位,1000000则是微秒为单位
fmt6 db /
'=================================================',0dh,0ah,/
'Simple asm program to compute Pi with many digits.',0dh,0ah,/
'Pi/4=12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)',/
0dh,0ah,'By G-Spider @2010',0dh,0ah,/
'Computation of 10000 digits of Pi',0dh,0ah,/
'Total computation time: ',/
'%6lld ms',0dh,0ah,/
'==================================================',0dh,0ah,0
.data?
x dd ?
y dd ?
_sizeChar dd ?
lpPi dd ?
lpArctan dd ?
lpBuffer1 dd ?
lpBuffer2 dd ?
;-------------------------------------
dqTickCounter1 dq ?
dqTickCounter2 dq ?
dqFreq dq ?
dqTime dq ?
;***********************************************************
.code
;-----------------------------------------------------------
_SetToInteger proc n:dword,lpx:dword,integer:dword
mov eax,1
mov edx,n
mov esi,lpx
@@:
mov dword ptr [esi+eax*4],0
inc eax
cmp eax,edx
jl @B
mov eax,integer
mov [esi],eax
xor eax,eax
ret
_SetToInteger endp
;-----------------------------------------------------------
_IsZero proc n:dword,lpx:dword
;全0则返回1
xor eax,eax
mov esi,lpx
mov ecx,n
@@:
mov edx,[esi+eax*4]
test edx,edx
jnz exit
inc eax
cmp eax,ecx
jl @B
mov eax,1
ret
exit:
xor eax,eax
ret
_IsZero endp
;-----------------------------------------------------------
_Add proc n:dword,lpx:dword,lpy:dword
;x+=y
mov edi,lpx
mov esi,lpy
xor ecx,ecx ;carry
mov eax,n
dec eax ;n-1 下标
@@:
mov edx,[esi+eax*4]
add edx,ecx
add edx,[edi+eax*4]
cmp edx,B
jl A111
mov ecx,1
sub edx,B
mov [edi+eax*4],edx
dec eax
jns @B ;不为负则跳转
A111:
xor ecx,ecx
mov [edi+eax*4],edx
dec eax
jns @B ;不为负则跳转
xor eax,eax
ret
_Add endp
;-----------------------------------------------------------
_Sub proc n:dword,lpx:dword,lpy:dword
;x-=y
mov edi,lpx
mov esi,lpy
mov eax,n
@@:
dec eax ;n-1
mov ebx,[edi+eax*4]
sub ebx,[esi+eax*4]
mov [edi+eax*4],ebx ;x[i]-=y[i]
cmp ebx,0 ;x[i]<0
jge A000
test eax,eax
jz A001
add ebx,B
mov [edi+eax*4],ebx ;x[i]+=B
dec eax
sub dword ptr [edi+eax*4],1 ;x[i-1]--
inc eax
A000:
cmp eax,0
jg @B
A001:
xor eax,eax
ret
_Sub endp
;-----------------------------------------------------------
_Mul proc n:dword,lpx:dword,q:dword
;x*=q
xor ecx,ecx ;carry=0
mov edi,lpx
mov eax,n
@@:
dec eax ;n-1
mov ebx,[edi+eax*4]
imul ebx,q ;xi=xi*q
push eax
pop eax
add ebx,ecx
cmp ebx,B
jl A000
;******************************
push eax
mov eax,ebx
cdq
idiv B ;eax=eax/B=xi/B
mov ecx,eax ;ecx=carry
imul eax,B ;carry*B
sub ebx,eax ;xi-=(carry*B)
pop eax
mov [edi+eax*4],ebx
test eax,eax
jnz @B
;******************************
jmp exit
A000:
mov ecx,0
mov [edi+eax*4],ebx
test eax,eax
jnz @B
exit:
xor eax,eax
ret
_Mul endp
;-----------------------------------------------------------
_Div proc n:dword,lpx:dword,d:dword,lpy:dword
;y=x/d
xor ecx,ecx ;carry=0
xor eax,eax
mov esi,lpx
mov edi,lpy
@@:
mov ebx,ecx
imul ebx,B
add ebx,[esi+eax*4] ;ebx=x[i]+carry*B
push eax ;*************
mov eax,ebx
cdq
idiv d ;eax=ebx/d
mov edx,eax ;save eax to edx 商
imul eax,d ;eax=eax*d
mov ecx,ebx
sub ecx,eax ;carry=xi-eax*d
pop eax ;*************
mov [edi+eax*4],edx
inc eax
cmp eax,n
jl @B
xor eax,eax
ret
_Div endp
_ArcCot proc _p:dword,n:dword,lpx:dword,buf1:dword,buf2:dword
local @sign,p2,k
mov @sign,0
mov eax,_p
imul eax,_p
mov p2,eax
mov k,3
invoke _SetToInteger,n,lpx,0
invoke _SetToInteger,n,buf1,1
invoke _Div,n,buf1,_p,buf1
invoke _Add,n,lpx,buf1
@@:
invoke _IsZero,n,buf1
test eax,eax
jnz exit
mov eax,MaxDiv
.if p<eax
invoke _Div,n,buf1,p2,buf1
.else
invoke _Div,n,buf1,_p,buf1
invoke _Div,n,buf1,_p,buf1
.endif
invoke _Div,n,buf1,k,buf2
.if @sign
invoke _Add,n,lpx,buf2
.else
invoke _Sub,n,lpx,buf2
.endif
add k,2
neg @sign
add @sign,1
jmp @B
exit:
xor eax,eax
ret
_ArcCot endp
_Print proc n:dword,lpx:dword
local @i
mov esi,lpx
dec n
mov @i,1
invoke crt_printf,offset fmt1,dword ptr[esi]
@@:
mov eax,@i
invoke crt_printf,offset fmt,dword ptr[esi+eax*4]
mov eax,@i
mov ecx,10
cdq
idiv ecx
test edx,edx
jnz A000
mov edx,@i
shl edx,2
invoke crt_printf,offset fmt2 ,edx
A000:
inc dword ptr @i
mov eax,@i
cmp eax,n
jl @B
invoke crt_printf,offset fmt3
xor eax,eax
ret
_Print endp
_PiCalc proc
local @NbArctan,@flag
mov @NbArctan,3
mov eax,_size
shl eax,2
mov _sizeChar,eax
;-----------------------------------
invoke crt_malloc,_sizeChar
.if eax
mov lpPi,eax
.else
invoke MessageBox,NULL,0,0,0
.endif
invoke crt_malloc,_sizeChar
.if eax
mov lpArctan,eax
.else
invoke MessageBox,NULL,0,0,0
.endif
invoke crt_malloc,_sizeChar
.if eax
mov lpBuffer1,eax
.else
invoke MessageBox,NULL,0,0,0
.endif
invoke crt_malloc,_sizeChar
.if eax
mov lpBuffer2,eax
.else
invoke MessageBox,NULL,0,0,0
.endif
;-----------------------------------
;时间测试开始
invoke QueryPerformanceCounter,addr dqTickCounter1
;//////////////////////////////////////////////////
invoke _SetToInteger,_size,lpPi,0
xor ecx,ecx
@@:
push ecx
invoke _ArcCot,p[ecx*4],_size,lpArctan,lpBuffer1,lpBuffer2
pop ecx
mov eax, m[ecx*4]
and eax,80000000
.if eax!=0
neg m[ecx*4]
mov @flag,1
.else
mov @flag,0
.endif
push ecx
invoke _Mul,_size,lpArctan,m[ecx*4]
.if @flag
invoke _Sub,_size,lpPi,lpArctan
.else
invoke _Add,_size,lpPi,lpArctan
.endif
pop ecx
inc ecx
cmp ecx,@NbArctan
jl @B
invoke _Mul,_size,lpPi,4
;//////////////////////////////////////////////////
;时间测试结束
invoke QueryPerformanceCounter,addr dqTickCounter2
invoke QueryPerformanceFrequency,addr dqFreq
mov eax,dword ptr dqTickCounter1
mov edx,dword ptr dqTickCounter1[4]
sub dword ptr dqTickCounter2,eax
sub dword ptr dqTickCounter2[4],edx
finit
fild dqFreq
fild dqTickCounter2
fimul dwlm
fdivr
fistp dqTime ;dqTime中的64位值就是时间间隔(以毫秒为单位)
invoke crt_printf,offset fmt6,dqTime
;---------------------------------------------------
invoke _Print,_size,lpPi ;输出结果到控制台
;----------------------------------------------------
invoke crt_free,lpPi
invoke crt_free,lpArctan
invoke crt_free,lpBuffer1
invoke crt_free,lpBuffer2
xor eax,eax
ret
_PiCalc endp
;***********************************************************
start:
invoke _PiCalc
invoke crt_system,offset szPause
invoke ExitProcess,0
end start
扫一扫
7
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
