Pi to 1000000 digits

[MichaelW]

Pi to 1000000 digits

This combines the Brent-Salamin algorithm with GMP (using the gmp-dynamic-vc-4.1.2 library from here) to calculate the first 1000, 100000, or 1000000 digits of pi.

Code Select Expand


;==============================================================================
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:

Code Select Expand


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


[hutch--]

> For some reason I currently cannot add an attachment

Strange, its working perfectly from here.

[MichaelW]

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.

[jj2007]

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'

[MichaelW]

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.

[MichaelW]

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.

[Adamanteus]

So happend, that I can't prove this result
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 

[bluedevil]

not working due to missing GMP.dll. But here it is, i attached it:

[bluedevil]

cannot compile due to lack of libc.lib
how can we install your gmp package MichaelW?

[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.

[bluedevil]

this line has a problem about syntax, says compiler?

Code Select Expand

printf("%dms\n\n", eax)

[MichaelW]

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.

[bluedevil]

you are right i have then 10th version now i am downloading the 11th one :icon_redface:

[Farabi]

277 seconds?  that is super fast, did not some folks tested it for a few hours or days?