The MASM Forum
General => The Workshop => Topic started by: MichaelW on July 20, 2012, 01:01:47 AM
-
Pi to 1000000 digits
This combines the Brent-Salamin algorithm (http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm) with GMP (http://gmplib.org/index.html#dir) (using the gmp-dynamic-vc-4.1.2 library from here (http://cs.nyu.edu/exact/core/gmp/)) to calculate the first 1000, 100000, or 1000000 digits of pi.
;==============================================================================
include \masm32\include\masm32rt.inc
includelib gmp.lib
;==============================================================================
PIDIGITS = 1000
;==============================================================================
;-----------------------------------------------
; Sorry, no time to translate the full headers.
;-----------------------------------------------
mp_exp_t TYPEDEF DWORD
mp_limb_t TYPEDEF DWORD
__mpf_struct STRUCT
_mp_prec DWORD ?
_mp_size DWORD ?
_mp_exp mp_exp_t ?
_mp_d mp_limb_t ?
__mpf_struct ENDS
__gmpf_init PROTO C :DWORD
__gmpf_init2 PROTO C :DWORD,:DWORD
__gmpf_set_str PROTO C :DWORD,:DWORD,:DWORD
__gmpf_init_set_str PROTO C :DWORD,:DWORD,:DWORD
__gmpf_set_prec PROTO C :DWORD,:DWORD
__gmpf_clear PROTO C :DWORD
__gmpf_add PROTO C :DWORD,:DWORD,:DWORD
__gmpf_sub PROTO C :DWORD,:DWORD,:DWORD
__gmpf_mul PROTO C :DWORD,:DWORD,:DWORD
__gmpf_mul_ui PROTO C :DWORD,:DWORD,:DWORD
__gmpf_div PROTO C :DWORD,:DWORD,:DWORD
__gmpf_div_ui PROTO C :DWORD,:DWORD,:DWORD
__gmpf_ui_div PROTO C :DWORD,:DWORD,:DWORD
__gmpf_sqrt PROTO C :DWORD,:DWORD
__gmpf_pow_ui PROTO C :DWORD,:DWORD,:DWORD
__gmpf_swap PROTO C :DWORD,:DWORD
__gmp_printf PROTO C :DWORD,:VARARG
;==============================================================================
.data
a __mpf_struct <>
b __mpf_struct <>
t __mpf_struct <>
p __mpf_struct <>
A __mpf_struct <>
B __mpf_struct <>
T __mpf_struct <>
P __mpf_struct <>
Pi __mpf_struct <>
.code
;==============================================================================
start:
;==============================================================================
invoke GetTickCount
push eax
IF PIDIGITS EQ 1000
;PRECISION = 3297 ; minimum for correct value of Pi
PRECISION = 3297+10
ELSEIF PIDIGITS EQ 100000
;PRECISION = 332193 ; minimum for correct value of Pi
PRECISION = 332193+10
ELSEIF PIDIGITS EQ 1000000
;PRECISION = 3321921 ; minimum for correct value of Pi
PRECISION = 3321921+10
ELSE
echo ------------------------------------------------
echo ERROR PIDIGITS MUST BE 1000, 100000, OR 1000000
echo ------------------------------------------------
.err
ENDIF
FOR var,<a,b,t,p,A,B,T,P,Pi>
invoke __gmpf_init2, ADDR var, PRECISION
ENDM
invoke __gmpf_set_str, ADDR a, chr$("1.0"),10
invoke __gmpf_set_str, ADDR b, chr$("2.0"),10
invoke __gmpf_set_str, ADDR t, chr$("0.25"),10
invoke __gmpf_set_str, ADDR p, chr$("1.0"),10
invoke __gmpf_sqrt, ADDR b, ADDR b
invoke __gmpf_ui_div, ADDR b, 1, ADDR b
;invoke __gmp_printf, cfm$("%.20Ff\n"), ADDR b ; should be 0.707106781...
IF PIDIGITS EQ 1000
mov ebx, 10
ELSEIF PIDIGITS EQ 100000
mov ebx, 17
ELSEIF PIDIGITS EQ 1000000
mov ebx, 20
ENDIF
@@:
;-------------
; A = (a+b)/2
;-------------
invoke __gmpf_add, ADDR A, ADDR a, ADDR b
invoke __gmpf_div_ui, ADDR A, ADDR A, 2
;---------------
; B = sqrt(a*b)
;---------------
invoke __gmpf_mul, ADDR B, ADDR a, ADDR b
invoke __gmpf_sqrt, ADDR B, ADDR B
;-----------------
; T = t-p*(a-A)^2
;-----------------
invoke __gmpf_sub, ADDR T, ADDR a, ADDR A
invoke __gmpf_pow_ui, ADDR T, ADDR T, 2
invoke __gmpf_mul, ADDR T, ADDR T, ADDR p
invoke __gmpf_sub, ADDR T, ADDR t, ADDR T
;---------
; P = 2*p
;---------
invoke __gmpf_mul_ui, ADDR P, ADDR p, 2
;----------------------------
; a = A, b = B, t = T, p = P
;----------------------------
invoke __gmpf_swap, ADDR a, ADDR A
invoke __gmpf_swap, ADDR b, ADDR B
invoke __gmpf_swap, ADDR t, ADDR T
invoke __gmpf_swap, ADDR p, ADDR P
sub ebx, 1
jnz @B
;-------------------
; Pi ~= (a+b)^2/4*t
;-------------------
invoke __gmpf_add, ADDR Pi, ADDR a, ADDR b
invoke __gmpf_pow_ui, ADDR Pi, ADDR Pi, 2
invoke __gmpf_div_ui, ADDR Pi, ADDR Pi, 4
invoke __gmpf_div, ADDR Pi, ADDR Pi, ADDR t
IF PIDIGITS EQ 1000
invoke __gmp_printf, cfm$("%.1000Ff\n\n"), ADDR Pi
ELSEIF PIDIGITS EQ 100000
invoke __gmp_printf, cfm$("%.100000Ff\n\n"), ADDR Pi
ELSEIF PIDIGITS EQ 1000000
invoke __gmp_printf, cfm$("%.1000000Ff\n\n"), ADDR Pi
ENDIF
invoke GetTickCount
pop edx
sub eax, edx
printf("%dms\n\n", eax)
inkey
exit
;==============================================================================
end start
These from http://gc3.net84.net/pi.htm
Last 40 digits of 1000:
1712268066 1300192787 6611195909 2164201989
Last 40 digits of 100000:
7728658035 7127909137 6742080565 5493624646
Last 40 digits of 1000000:
3311646283 9963464604 2209010610 5779458151
I determined the minimum bits of precision that will produce the correct values of Pi experimentally.
I derived the number of iterations required from the quadratic convergence of the algorithm.
On my 3.0GHz P4 the 1000000-digit calculation runs in about 277 seconds.
For some reason I currently cannot add an attachment, so here is the batch file that I used to build the app with a larger stack:
set file="test"
set subsys=CONSOLE
set path=masm32\bin;%PATH%
set include=masm32\include
set lib=masm32\lib
: rc /v rsrc.rc
: cvtres /machine:ix86 rsrc.res
: pause
if exist %file%.obj del %file%.obj
if exist %file%.exe del %file%.exe
ml /c /coff %file%.asm
pause
: Link /SUBSYSTEM:%subsys% %file%.obj rsrc.res
:
: A larger stack is necessary to avoid a
: stack overflow with 1000000 digits
: /STACK:reserve[,commit]
Link /STACK:0x400000 /SUBSYSTEM:%subsys% %file%.obj
pause
-
> For some reason I currently cannot add an attachment
Strange, its working perfectly from here.
-
It works for me now too, but when I posted, on the same system that I used just now to add the attachment, I tried multiple times with no success. Concurrent with this, every time I clicked the Preview button there would be a 10-20 second delay before “Fetching preview” appeared, and I’ve seen that problem multiple times now.
-
Interesting stuff :t, but apparently it needs a full C installation which I don't have. After including gmp.lib, the linker complains:
POLINK: fatal error: File not found: 'LIBC.lib'
-
And past LIBC.lib you'll find other things missing. That is why I specified the dynamic library instead of the static library that I would have preferred.
-
I was able to do 10 million digits OK, in 252 minutes. But running on a system with 512MB of memory, in the time I had to spend on it, for 100 million digits I could not find any way to avoid an eventual stack overflow.
-
So happend, that I can't prove this result :dazzled:
I'am using own version of GMPXX 5 - and testing you algo linked to it give 3.14159265 4258268 - 8 numbers proved. And to me really not undastandable how on arbitrary precision algos of GMP (that out of 2 limbs begin approximating weight of number bits) possible count with such exactness.
I attached example with exe and source, that counting Pi with 40 numbers proved (more not giving my GMPXX) - maybe you could complie it on you version and post for checking perposes :eusa_clap:
-
not working due to missing GMP.dll. But here it is, i attached it:
-
cannot compile due to lack of libc.lib
how can we install your gmp package MichaelW?
-
how can we install your gmp package MichaelW?
Sorry about being so slow to respond. Copy test.asm and makeit.bat from my attachment to your working directory, get the gmp-dynamic-vc-4.1.2.zip that I provided a link to and copy gmp.dll and the matching import library gmp.lib from it to your working directory, then run makeit.bat. As I stated, I would have preferred a static library, but in the time I had available I could not make it work. Also, I learned on the FreeBASIC forum that using a more recent version of GMP can speed the calculations up substantially (~5x IIRC), but I could not build the required binary in the time I had available.
-
this line has a problem about syntax, says compiler?
printf("%dms\n\n", eax)
-
I’m guessing here that the problem is your version of the MASM32 package. The printf macro was first included in the version 11 macros.asm.
-
you are right i have then 10th version now i am downloading the 11th one :icon_redface:
-
277 seconds? :dazzled: that is super fast, did not some folks tested it for a few hours or days?