Pi to 1 million

[srvaldez]

here in America they have what is called Pi day which is March 14, march being the third month
it's a bit early but this program will compute Pi to a million digits in mere seconds, or less if you have a fast PC
formula from https://en.wikipedia.org/wiki/Borwein%2 ... nce_(1985)
in case you want to check your result https://stuff.mit.edu/afs/sipb/contrib/pi/
my time about 1.33 seconds

Code: Select all

#include once "gmp.bi"

dim as __mpf_struct a, y, y2, y4, pi
dim as long prec, p2, n, n2=1000000

dim as double lt=3.33333, t=timer

prec=n2*lt
mpf_init2(@a, prec)
mpf_init2(@y, prec)
mpf_init2(@y2, prec)
mpf_init2(@y4, prec)
mpf_init2(@pi, prec)

mpf_set_ui(@a, 2)		'a=2
mpf_sqrt(@a, @a)		'a=sqr(2)
mpf_set(@y, @a)			'y=sqr(2)
mpf_sub_ui(@y, @y, 1)	'y=sqr(2)-1
mpf_mul_ui(@a, @a, 4)	'a=sqr(2)*4
mpf_ui_sub(@a, 6, @a)	'a=6-sqr(2)*4

p2=2

n=log(n2)*log(2)

for i as long=1 to n
	p2*=4
	mpf_mul(@y2, @y, @y)
	mpf_mul(@y4, @y2, @y2)
	mpf_ui_sub(@y, 1, @y4)
	mpf_sqrt(@y, @y)
	mpf_sqrt(@y, @y)
	
	mpf_ui_sub(@y2, 1, @y)	'y2=1-y
	mpf_add_ui(@y4, @y, 1)	'y4=1+y
	mpf_div(@y, @y2, @y4)	'y=(1-y)/(1+y)
	
	mpf_add_ui(@y2, @y, 1)	'y2=1+y
	mpf_mul(@y2, @y2, @y2)	'y2=y2*y2
	mpf_mul(@y4, @y2, @y2)	'y4=y2*y2
	
	mpf_mul(@y2, @y, @y)	'y2=y*y
	mpf_add(@y2, @y2, @y)	'y2=y+y*y
	mpf_add_ui(@y2, @y2, 1)	'y2=1+y+y*y
	mpf_mul(@y2, @y2, @y)
	mpf_mul_ui(@y2, @y2, p2)
	mpf_mul(@a, @a, @y4)
	mpf_sub(@a, @a, @y2)
	'a=a*y4-p2*y*(1+y+y*y)'/
next
mpf_ui_div(@pi, 1, @a) 'pi=1/a
t=timer-t

print "it took ";n;" loops in ";t;" seconds"
print "press return to continue"
sleep

dim frmt as string
dim as long length = mpf_get_prec(@pi)*0.3010299956639812
if length>n2 then length=n2
frmt="%"+str(length+4)+"."+str(length-1)+"f"
frmt = left(frmt,len(frmt)-1)+"F"+right(frmt,1)
length = gmp_snprintf(0, 0, frmt, @pi)
dim as zstring ptr rstring=allocate(length+20)
gmp_sprintf(rstring, frmt, @pi)

print *rstring
print length
print "press return to end"
sleep
deallocate(rstring)

mpf_clear(@pi)
mpf_clear(@y4)
mpf_clear(@y2)
mpf_clear(@y)
mpf_clear(@a)

[badidea]

Here 32 vs 64-bit fbc makes big difference: 8.5 vs 2.7 second.
I check the first 3 digits, those looked familiar.

[srvaldez]

I check the first 3 digits, those looked familiar.

LOL
you could use mpfr, a bit faster about 0.76 seconds

Code: Select all

#include once "gmp.bi"
#include once "mpfr.bi"

dim as double t=timer
dim as __mpfr_struct pi
mpfr_init2(@pi, 3321928)
mpfr_const_pi(@pi, MPFR_RNDN)
dim as zstring ptr s=allocate(1000020)
mpfr_sprintf(s, "%1000005.1000000R*f", MPFR_RNDN, @pi)
t=timer-t
? t
? "press return to continue"
sleep
? *s
sleep
deallocate(s)
mpfr_clear(@pi)

[dodicat]

Here was my version way back.

Code: Select all

#Include once "gmp.bi"
type mpf_t as  __mpf_struct 
Dim Shared As Zstring * 100000000 outtext

Function Pi_ui Overload(places As UInteger) As mpf_t
    Dim As __mpf_struct a, b, t, aa, bb, tt, pi
    mpf_init2(@a, 4*places)
    mpf_init2(@b, 4*places)
    mpf_init2(@t, 4*places)
    Dim As UInteger p
    mpf_init2(@aa,4*places)
    mpf_init2(@bb,4*places)
    mpf_init2(@tt,4*places)
    mpf_init2(@pi,4*places)
    mpf_set_ui(@a, 1)
    mpf_set_ui(@b, 2) : mpf_sqrt(@b, @b)
    mpf_set_str(@t,".25",10)
    mpf_ui_div(@b,1,@b)
    Do
        mpf_add(@aa, @a, @b)
        mpf_div_2exp(@aa, @aa, 1)
        mpf_mul(@bb, @a, @b)
        mpf_sqrt(@bb, @bb)
        mpf_sub(@tt, @a, @aa)
        mpf_mul(@tt,@tt,@tt)
        mpf_mul_2exp(@tt, @tt, p)
        p += 1
        mpf_sub(@tt, @t, @tt)
        mpf_swap(@a, @aa)
        mpf_swap(@b, @bb)
        mpf_swap(@t, @tt)
    Loop Until  Mpf_cmp(@a, @aa) = 0
    mpf_add(@pi, @a, @b)
    mpf_mul(@pi, @pi, @pi)
    mpf_div_2exp(@pi, @pi, 2)
    mpf_div(@pi, @pi, @t)
    mpf_clear(@a) : mpf_clear(@aa)
    mpf_clear(@b) : mpf_clear(@bb)
    mpf_clear(@t) : mpf_clear(@tt)
    Return pi
End Function

Function Pi_ui Overload(places As Integer) As String
    Dim As Mpf_t ans
    Var pl=CUInt(places)
    ans=Pi_ui(pl)
    gmp_sprintf(@outtext, "%." & pl & "Ff", @ans )
    mpf_clear(@ans)
    Var outtxt=Trim(outtext)
    If InStr(outtxt,".") Then outtxt= RTrim(outtxt,"0"):outtxt=RTrim(outtxt,".")
    Return Trim(outtxt)
End Function

dim as double t=timer
var pi= pi_ui(1000000)
var t2= timer-t
print pi,"length  ";len(pi)
print t2; " seconds"
sleep
 

My 64 bit result:

Code: Select all

............606242994120327362558244174983450947309453436615907284163193683075719798068231535737155571816122156787936425013887117023275555779302266785803199930810830576307652332050740013939095807901637717629259283764874790177274125678190555562180504876746991140839977919376542320623374717324703369763357925891515260315614033321272849194418437150696552087542450598956787961303311646283996346460422090106105779458151              length   1000002
 2.428122299999814 seconds 

But there is a faster algorithm, but I cannot find it at the moment.

[srvaldez]

hi dodicat
your version is a tad faster than the one I posted
perhaps you were thinking of Pi Chudnovsky?
I found it and it's well over 3 times faster except for the printing to screen viewtopic.php?f=7&t=20398
the print to screen is extremely slow it's a 1000 times faster to print to a string and then print that, replace the statement

Code: Select all

gmp_printf(!"%Zd\n",@pi)

with

Code: Select all

dim as long length = gmp_snprintf(0, 0, "%Zd", @pi)
dim as zstring ptr rstring=allocate(length+20)
gmp_sprintf(rstring, "%Zd", @pi)
print *rstring
deallocate(rstring)

[dafhi]

new Asus G15 .. 0.637 (mpfr)