Pi to 1000000 digits with Brent-Salamin and GMP

Post your FreeBASIC tips and tricks here. Please don’t post your code without including an explanation.
MichaelW
Posts: 3500
Joined: May 16, 2006 22:34
Location: USA

Pi to 1000000 digits with Brent-Salamin and GMP

Postby MichaelW » Jul 22, 2012 6:03

This combines the Brent-Salamin algorithm with GMP (using the gmp-static-mingw-4.1 library from here) to calculate the first 1000000 digits of pi.

Code: Select all

''=============================================================================
#include "gmp.bi"
''=============================================================================

dim as double t1,t2
t1 = timer

#define PRECISION 3321921

dim as __mpf_struct a,b,t,p,aa,bb,tt,pp,pi

mpf_init2( @a, PRECISION )
mpf_init2( @b, PRECISION )
mpf_init2( @t, PRECISION )
mpf_init2( @p, PRECISION )
mpf_init2( @aa, PRECISION )
mpf_init2( @bb, PRECISION )
mpf_init2( @tt, PRECISION )
mpf_init2( @pp, PRECISION )
mpf_init2( @pi, PRECISION )

mpf_set_str( @a, "1.0", 10 )
mpf_set_str( @b, "2.0", 10 )
mpf_set_str( @t, "0.25", 10 )
mpf_set_str( @p, "1.0", 10 )
mpf_sqrt( @b, @b )
mpf_ui_div( @b, 1, @b )

for i as integer = 1 to 20

    ''--------------
    '' aa = (a+b)/2
    ''--------------

    mpf_add( @aa,  @a,  @b )
    mpf_div_ui( @aa, @aa , 2 )

    ''----------------
    '' bb = sqrt(a*b)
    ''----------------

    mpf_mul( @bb, @a, @b )
    mpf_sqrt( @bb, @bb )

    ''-------------------
    '' tt = t-p*(a-aa)^2
    ''-------------------

    mpf_sub( @tt, @a, @aa )
    mpf_pow_ui( @tt, @tt, 2 )
    mpf_mul( @tt, @tt, @p )
    mpf_sub( @tt, @t, @tt )

    ''----------
    '' pp = 2*p
    ''----------

    mpf_mul_ui( @pp, @p, 2 )

    ''--------------------------------
    '' a = aa, b = bb, t = tt, p = pp
    ''--------------------------------

    mpf_swap( @a, @aa )
    mpf_swap( @b, @bb )
    mpf_swap( @t, @tt )
    mpf_swap( @p, @pp )

next

''-------------------
'' pi ~= (a+b)^2/4*t
''-------------------

mpf_add( @pi, @a, @b )
mpf_pow_ui( @pi, @pi, 2 )
mpf_div_ui( @pi, @pi, 4 )
mpf_div( @pi, @pi, @t )

t2 = timer

gmp_printf( !"%.1000000Ff\n\n", @pi )

print t2-t1; " seconds"

sleep


To use the static library I copied it to my working directory and changed the name to differentiate it from the libgmp.dll.a import library in the FreeBASIC lib directory. I then copied gmp.bi to my working directory and modified the #inclib statement to match the name I selected for the static library, which was libgmp_static.a, so the statement was then:

#inclib "gmp_static"

I determined the minimum bits of precision that will produce the correct value of Pi experimentally. If the precision is too low then the digits at the end will be incorrect (usually but not necessarily zeros).

I derived the number of iterations required from the quadratic convergence of the algorithm.

Although the default stack size worked OK for much smaller numbers of digits, to calculate 1000000 digits I had to increase it by adding -t 3075 to the FBC command line. If the stack is too small, then depending on how much too small it is and/or the version of the library, the app may be terminated silently, or it may trigger an access-violation or stack-overflow exception.

Last 40 digits of 1000000:

3311646283 9963464604 2209010610 5779458151

Source: http://gc3.net84.net/pi.htm

On my 3.0GHz P4 the 1000000-digit calculation runs in ~59 seconds (where in a previous test with the VC dynamic library from the same page it ran in ~277 seconds).
dodicat
Posts: 6755
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby dodicat » Jul 22, 2012 8:44

Thanks MW.
I have the libgmp-10.dll from a Haskell distribution on this box.
(~3000 Mhz Pentium Xeon X 2 -- XP)
I simply copied the gmp.bi from the free basic files and went -- #inclib "gmp-10", and put it alongside the dll in a temp folder.
Here is the tail end of my result:
---

60624299412032736255824417498345094730945343661590728416319368307571979806823153
57371555718161221567879364250138871170232755557793022667858031999308108305763076
52332050740013939095807901637717629259283764874790177274125678190555562180504876
74699114083997791937654232062337471732470336976335792589151526031561403332127284
91944184371506965520875424505989567879613033116462839963464604220901061057794581
51

31.18116026572386 seconds
srvaldez
Posts: 2574
Joined: Sep 25, 2005 21:54

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby srvaldez » Jul 22, 2012 14:30

my setup
VMware Fusion on Mac Pro with 2 6-Core Intel Xeon CPU's and 48 GB RAM
virtual machine has 4 cores and 3 GB RAM running windows 7 32-bit
FBIde compiler command: "<$fbc>" "<$file>" -t 30000 -exx

with the static lib that you pointed to it takes 25 seconds on my virtual machine
with gmp 5.05 dll it takes about 8 seconds.
I added the following line after #include "gmp.bi"

Code: Select all

extern gmp_version alias "__gmp_version" as zstring ptr

and

Code: Select all

print "gmp version ";*gmp_version

after the elapsed time print statement.

3311646283996346460422090106105779458151

8.340624646429205 seconds
gmp version 5.0.5


@dodicat
you need a matching import lib to your dll, else your result or time may not be what you would expect,
if you don't have the matching import lib you can build it by something like this at command prompt
pexports -o libgmp-10.dll >libgmp-10.def
dlltool -d libgmp-10.def -l libgmp-10.dll.a

and of course make sure the right import lib and dll are used. :)
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby frisian » Jul 23, 2012 22:16

MichaelW

A while back I played with Brent-Salamin algorithm and I found out that it takes 19 loops to calculate Pi with 1e6 digits.
Your program uses 20 loops (1 to 20).

So
for i as integer = 1 to 20
can become
for i as integer = 1 to 19

The next suggestions give a small speedup when calculating Pi with 1e7 digits or higher.
On the early computer's A=B+B was faster then A=B*2 and A=B*B was faster then A=B^2.

Replace
mpf_pow_ui( @tt, @tt, 2 ) with mpf_mul(@tt, @tt, @tt)
and
mpf_mul_ui( @pp, @p, 2 ) with mpf_add(@pp, @p, @p)

For 1e6 digits the difference is to small to see, but for bigger numbers the difference becomes clearer.
Perhaps it will work for you.

Timings:
GMP 5.0.4 MinGW ??? build dll (not optimized)
AMD Athlon II x4 645 (4 core) 3255Mhz
Win7-64bit
20 loops : 6.88 sec.
19 loops : 6.58 sec.

OOPS
After building GMP 5.0.5 found out that with --build=none-pc-mingw32 the timing is 25 sec.
After a new build for a i586 the timing was 7.5 sec. Not a good idea to have different builds on the same disk, should stick to one. Since the only other build of GMP 5.0.4 I made is for the k10 the timing is for a k10 build.
Last edited by frisian on Jul 24, 2012 8:58, edited 1 time in total.
srvaldez
Posts: 2574
Joined: Sep 25, 2005 21:54

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby srvaldez » Jul 23, 2012 22:22

frisian wrote:Timings:
GMP 5.0.4 MinGW generic build dll (not optimized)
AMD Athlon II x4 645 (4 core) 3255Mhz
Win7-64bit
20 loops : 6.88 sec.
19 loops : 6.58 sec.

that's why I made a comment to dodicat, his timing should be close to yours.
@MichaelW, in your code you forgot to clear the mpf variables at the end, also wonder how hard it would be to implement the pi_chudnovsky in FB
dodicat
Posts: 6755
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby dodicat » Jul 24, 2012 0:45

Hi srvaldez.
I got a libgmp-10.def from tcc compiler, and the functions are shown.
(pexports is unrecognozed by the command interperator here, and I get an empty .def file)
I ran the dlltool from command prompt and got libgmp-10.dll.a -- ok.

But I still get the same speed for Pi.
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby frisian » Jul 24, 2012 9:18

After my recent attempt to build GMP 5.0.5 (see my previous post) it clear to me what going on.

You can build a generic version but it's slow, something the GMP PDF confirms but it works on every computer.

To have a faster version it needs to be build for a I486/I586 type processor or better build your one version.
greenllll
Posts: 57
Joined: Jul 02, 2007 19:33

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby greenllll » Sep 10, 2012 2:08

I'd like to create a txt file with pi to a 100 million digits, accurately. Can someone please show me what variables to change in this program to generate this? thanks.

Or, I'd like to download a file that has this. I found 1 with 1 million digits, but I'd like more.
I found 1 w/ 100 million digits, but it's formatted weird. I just want the numbers, nothing else.
srvaldez
Posts: 2574
Joined: Sep 25, 2005 21:54

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby srvaldez » Sep 10, 2012 5:02

try this, it took about 30 minutes on my machine.

Code: Select all

''=============================================================================
#include "gmp.bi"
''=============================================================================

dim as double t1,t2
t1 = timer

extern gmp_version alias "__gmp_version" as zstring ptr

#define PI_DIGITS 100000000
#define LOG2_10 3.321928094887362
#define PRECISION PI_DIGITS*LOG2_10

dim as __mpf_struct a,b,t,p,aa,bb,tt,pp,pi

mpf_init2( @a, PRECISION )
mpf_init2( @b, PRECISION )
mpf_init2( @t, PRECISION )
mpf_init2( @p, PRECISION )
mpf_init2( @aa, PRECISION )
mpf_init2( @bb, PRECISION )
mpf_init2( @tt, PRECISION )
mpf_init2( @pp, PRECISION )
mpf_init2( @pi, PRECISION )

mpf_set_str( @a, "1.0", 10 )
mpf_set_str( @b, "0.5", 10 )
mpf_set_str( @t, "0.25", 10 )
mpf_set_str( @p, "1.0", 10 )
mpf_sqrt( @b, @b )

for i as integer = 1 to int(log(PI_DIGITS)*1.44269504088896)

    ''--------------
    '' aa = (a+b)/2
    ''--------------

    mpf_add( @aa,  @a,  @b )
    mpf_div_ui( @aa, @aa , 2 )

    ''----------------
    '' bb = sqrt(a*b)
    ''----------------

    mpf_mul( @bb, @a, @b )
    mpf_sqrt( @bb, @bb )

    ''-------------------
    '' tt = t-p*(a-aa)^2
    ''-------------------

    mpf_sub( @tt, @a, @aa )
    mpf_pow_ui( @tt, @tt, 2 )
    'mpf_mul( @tt, @tt, @tt )
    mpf_mul( @tt, @tt, @p )
    mpf_sub( @tt, @t, @tt )

    ''----------
    '' pp = 2*p
    ''----------

    mpf_mul_ui( @pp, @p, 2 )

    ''--------------------------------
    '' a = aa, b = bb, t = tt, p = pp
    ''--------------------------------

    mpf_swap( @a, @aa )
    mpf_swap( @b, @bb )
    mpf_swap( @t, @tt )
    mpf_swap( @p, @pp )

next

''-------------------
'' pi ~= (a+b)^2/4*t
''-------------------

mpf_add( @pi, @a, @b )
mpf_pow_ui( @pi, @pi, 2 )
mpf_div_ui( @pi, @pi, 4 )
mpf_div( @pi, @pi, @t )

t2 = timer
dim as zstring ptr pie=allocate(PI_DIGITS+10)

gmp_sprintf(pie, "%."+str(PI_DIGITS)+!"Ff\n\n", @pi )

open "pi_100_million.txt" for output as 1
print #1,*pie
close 1
print "pi written to the file pi_100_million.txt"
print t2-t1; " seconds"
print "gmp version ";*gmp_version

mpf_clear( @pi )
mpf_clear( @pp )
mpf_clear( @tt )
mpf_clear( @bb )
mpf_clear( @aa )
mpf_clear( @p )
mpf_clear( @t )
mpf_clear( @b )
mpf_clear( @a )

deallocate(pie)

sleep
integer
Posts: 392
Joined: Feb 01, 2007 16:54
Location: usa

Re: Pi to 1000000 digits with Brent-Salamin and GMP

Postby integer » Jan 03, 2014 17:52

srvaldez wrote:try this, it took about 30 minutes on my machine.

and that for 100 million digits.

Here is a program in c that calculates individual digits of pi
Fabrice Bellard wrote both of these programs in c.
http://bellard.org/pi/pi.c

And a second one.
http://bellard.org/pi/pi1.c
This one I ported to FreeBasic.
All in FB code using LONGINT variables.

The program required almost half an hour for the 100,000th digit of pi.
Whereas srvaldez computed 1000 times MORE in the same time.

I have not posted my FB implementation of the code.
It's unclear as to what restrictions/prohibitions/liabilities the DRM imposes for derived works.

The only advantage: the individual digit extraction requires a small amount of computer memory.

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 4 guests