Pi Chudnovsky

[srvaldez]

this is a translation of the Python program pi_chudnovsky_bs_gmpy.py found at http://www.craig-wood.com/nick/articles/pi-chudnovsky/
it's about 3+ times faster than the Brent-Salamin algorithm.

Code: Select all

/'
python3 program to calculate pi using python long integers, binary
splitting and the chudnovsky algorithm

see: https://www.craig-wood.com/nick/articles/pi-chudnovsky/ for more info

nick craig-wood <nick@craig-wood.com>
'/

#Include Once "gmp.bi"

Dim Shared As __mpz_struct temp, c3_over_24

mpz_init(@temp)
mpz_init(@c3_over_24)

mpz_set_ui(@c3_over_24, 640320)
mpz_mul(@temp, @c3_over_24, @c3_over_24)
mpz_mul(@c3_over_24, @temp, @c3_over_24)
mpz_tdiv_q_ui(@c3_over_24,@c3_over_24,24)

Sub bs(Byval a As Uinteger, Byval b As Uinteger, Byref pab As __mpz_struct, Byref qab As __mpz_struct, Byref tab_ As __mpz_struct)
    /'
    computes the terms for binary splitting the chudnovsky infinite series

    a(a) = +/- (13591409 + 545140134*a)
    p(a) = (6*a-5)*(2*a-1)*(6*a-1)
    b(a) = 1
    q(a) = a*a*a*c3_over_24

    returns p(a,b), q(a,b) and t(a,b) '/
    
    Dim As Uinteger m, t
    
    If b - a = 1 Then
        ' directly compute p(a,a+1), q(a,a+1) and t(a,a+1)
        If a = 0 Then
            mpz_set_ui(@pab, 1)
            mpz_set_ui(@qab, 1)
        Else
            t=6*a
            mpz_set_ui(@pab, t-1)
            mpz_mul_ui(@pab,@pab,a+a-1)
            mpz_mul_ui(@pab,@pab,t-5)
            mpz_set_ui(@temp,a)
            mpz_mul(@qab,@temp,@temp)
            mpz_mul(@qab,@qab,@temp)
            mpz_mul(@qab,@qab,@c3_over_24)
        End If
        mpz_set_ui(@tab_,545140134)
        mpz_mul_ui(@tab_,@tab_,a)
        mpz_add_ui(@tab_,@tab_,13591409)
        mpz_mul(@tab_,@tab_,@pab) ' a(a) * p(a)
        If a And 1 Then
            mpz_neg(@tab_,@tab_)' = -tab
        End If
    Else
        Dim As __mpz_struct pam, qam, tam, pmb, qmb, tmb
        mpz_init(@pam)
        mpz_init(@qam)
        mpz_init(@tam)
        mpz_init(@pmb)
        mpz_init(@qmb)
        mpz_init(@tmb)
        ' recursively compute p(a,b), q(a,b) and t(a,b)
        ' m is the midpoint of a and b
        m = (a + b) / 2
        ' recursively calculate p(a,m), q(a,m) and t(a,m)
        bs(a, m, pam, qam, tam)
        ' recursively calculate p(m,b), q(m,b) and t(m,b)
        bs(m, b, pmb, qmb, tmb)
        ' now combine
        mpz_mul(@pab, @pam, @pmb)
        mpz_mul(@qab, @qam, @qmb)
        mpz_mul(@temp, @qmb, @tam)
        mpz_mul(@tab_, @pam, @tmb)
        mpz_add(@tab_, @tab_, @temp)
        mpz_clear(@tmb)
        mpz_clear(@qmb)
        mpz_clear(@pmb)
        mpz_clear(@tam)
        mpz_clear(@qam)
        mpz_clear(@pam)
    End If

End Sub
        
Sub pi_chudnovsky_bs(Byref pi As __mpz_struct, Byval digits As Uinteger)
    /'
    compute int(pi * 10**digits)

    this is done using chudnovsky's series with binary splitting
    '/
    Dim As __mpz_struct one_squared, p, q, sqrtc, t
    Dim As Uinteger n
    Dim As Double digits_per_term

    mpz_init(@one_squared)
    mpz_init(@p)
    mpz_init(@q)
    mpz_init(@sqrtc)
    mpz_init(@t)
    
'    mpz_set_ui (@c, 640320)
'    mpz_mul (@c3_over_24, @c, @c)
'    mpz_mul (@c3_over_24, @c3_over_24, @c)
'    mpz_cdiv_q_ui (@c3_over_24, @c3_over_24, 24)

    ' how many terms to compute
    digits_per_term = Log(mpz_get_d (@c3_over_24)/6/2/6)*0.4342944819032518
    n = Int(digits/digits_per_term + 1)
    ' calclate p(0,n) and q(0,n)
    bs(0, n, p, q, t)
    mpz_set_ui(@one_squared,10)
    mpz_pow_ui(@one_squared,@one_squared,2*digits)
    mpz_mul_ui(@sqrtc,@one_squared,10005)
    mpz_sqrt(@sqrtc,@sqrtc)
    mpz_mul_ui(@pi,@sqrtc,426880)
    mpz_mul(@pi,@pi,@q)
    mpz_tdiv_q(@pi,@pi,@t)

    mpz_clear(@t)
    mpz_clear(@sqrtc)
    mpz_clear(@q)
    mpz_clear(@p)
    mpz_clear(@one_squared)
    
End Sub

Dim As Double t1, t2
Dim As __mpz_struct pi

mpz_init(@pi)

t1=Timer
pi_chudnovsky_bs(pi,1000000)
t2=Timer

gmp_printf(!"%Zd\n",@pi)
Print "elapsed time: ";t2-t1 
Print "press RETURN to end ";
Sleep

mpz_clear(@pi)
mpz_clear(@c3_over_24)
mpz_clear(@temp)

[Destructosoft]

If we're gonna play "dueling πs", why not investigate the algorithms of Borwein and Borwein? In a Scientific American article on the Chudnovskys (or maybe it was just π), I seem to remember they printed something by B and B that supposedly outperformed C and C.

EDIT: I got it wrong, I now remember it was an article on Ramanujan, sorry about that.

[Coolman]

this code generates many errors when compiling with freebasic 1.09 under linux !

[SARG]

Replacing 'as mpz_t' by 'as __mpz_struct' compillation is ok but I can't execute because lgmp lib is missing.

[srvaldez]

hello Coolman
I corrected the code as per SARG's suggestion and tested it on Fedora without problems
@SARG
if you are on Windows and using FB with winlibs gcc then you can find libgmp-10.dll in the topmost folder of bin/libexec
you could just copy and rename the dll in the place where your program resides

[Coolman]

hello Coolman
I corrected the code as per SARG's suggestion and tested it on Fedora without problems
@SARG
if you are on Windows and using FB with winlibs gcc then you can find libgmp-10.dll in the topmost folder of bin/libexec
you could just copy and rename the dll in the place where your program resides

thank you. it works now and it's very fast : 0.4x

[jdebord]

My contributions (using GMP + MPFR) :

Code: Select all

' -------------------------------------------------------------------
' Compute Pi by Chudnovsky's formula 
' https://en.wikipedia.org/wiki/Chudnovsky_algorithm
' -------------------------------------------------------------------

' Uses GMP and MPFR 
' https://www.freebasic.net/forum/viewtopic.php?t=24110&start=45

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

const Digits = 200

'we need to set mpfr_digits before including the overloded operators
dim shared as long mpfr_digits_precision = Digits '<- set mpfr digits-precision

dim shared as mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
' rounding options are
'  MPFR_RNDN    round to nearest, with ties to even
'  MPFR_RNDZ    round toward zero
'  MPFR_RNDU    round toward +Inf
'  MPFR_RNDD    round toward -Inf
'  MPFR_RNDA    round away from zero
'  MPFR_RNDF    faithful rounding (not implemented yet)
'  MPFR_RNDNA   round to nearest, with ties away from zero (mpfr_round)

#include once "mpfr-fb.bi" ' mpfr overloaded operators
''=============================================================================

' The formula adds 14 digits at each iteration
 
const MaxIter = Digits \ 14 + 1

dim as mpfr C, L, DL, X, FX, M, S
dim as long iter
    
C   = mpfr(426880) * sqr(mpfr(10005))
L   = mpfr(13591409)
DL  = mpfr(545140134)
X   = 1
FX  = mpfr(-262537412640768000)
S   = L
    
for iter = 1 to MaxIter
  M = facui(6 * iter) / (facui(3 * iter) * facui(iter) ^ 3)
  L = L + DL
  X = X * FX
  S = S + M * L / X
next iter
    
print "Pi Chudnovsky : " : print C / S
print
print "Reference : " : print pi_const()

sleep

Code: Select all

' -------------------------------------------------------------------
' Compute Pi by Borwein's method 
' https://carma.edu.au/resources/jon/RAMA125f.pdf
' Each iteration multiplies the number of digits  by 4
' -------------------------------------------------------------------

' Uses GMP and MPFR 
' https://www.freebasic.net/forum/viewtopic.php?t=24110&start=45

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

const Digits = 200

'we need to set mpfr_digits before including the overloded operators
dim shared as long mpfr_digits_precision = Digits '<- set mpfr digits-precision

dim shared as mpfr_rnd_t mpfr_rounding_mode = MPFR_RNDN '<- set the rounding mode
' rounding options are
'  MPFR_RNDN    round to nearest, with ties to even
'  MPFR_RNDZ    round toward zero
'  MPFR_RNDU    round toward +Inf
'  MPFR_RNDD    round toward -Inf
'  MPFR_RNDA    round away from zero
'  MPFR_RNDF    faithful rounding (not implemented yet)
'  MPFR_RNDNA   round to nearest, with ties away from zero (mpfr_round)

#include once "mpfr-fb.bi" ' mpfr overloaded operators
'=============================================================================

dim as mpfr p, r, d, a, y, y2, y4, z
dim as long iter
    
p = 0.25
r = sqr(mpfr(2))
d = 8
a = 6 - 4 * r
y = r - 1

print
print "Pi by Borwein's formula"
print

for iter = 1 to 4
  y2 = y * y
  y4 = y2 * y2
  z = (1 - y4)^p 
  y = (1 - z) / (1 + z)
  a = a * (1 + y)^4 - d * y * (1 + y + y * y)
  d = 4 * d
  print "Iteration"; iter: print 1 / a : print
next iter

print : print "Reference" : print pi_const

sleep

[dodicat]

My gmp method, alas, done so long ago I forget the method I used.
I used the 64 bit dll from libexec in the latest official build, and renamed it to libgmp.dll
(NOTE: I copied libgmp-10.dll to another folder and renamed it there )

Code: Select all

#cmdline "-gen gcc -Wc -O3"
#Include Once "gmp.bi"
Type mpf_t As __mpf_struct

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)
    ' remove big int's from memory
    mpf_clear(@a) : mpf_clear(@aa)
    mpf_clear(@b) : mpf_clear(@bb)
    mpf_clear(@t) : mpf_clear(@tt)
    Return pi
End Function

Dim Shared As ZString * 100000000 outtext

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 pistring= Pi_ui(1000000)
t=timer-t
print pistring
print "Time taken ";t
print len(pistring)
sleep 

Actually it makes no difference in speed by optimizing
#cmdline "-gen gcc -Wc -O3"
I am sure the dll is already optimised.
But my question
Why is the gmp dll in along with the fb package?

[srvaldez]

But my question
Why is the gmp dll in along with the fb package?

dodicat, I recommend that you install the NTCore Explorer Suite
launch it and drag&drop cc1.exe onto the explorer window, on the left pane click on imports, you will see a list of dependencies

[Coolman]

My gmp method, alas, done so long ago I forget the method I used.
I used the 64 bit dll from libexec in the latest official build, and renamed it to libgmp.dll

Code: Select all

#cmdline "-gen gcc -Wc -O3"
#Include Once "gmp.bi"
Type mpf_t As __mpf_struct

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)
    ' remove big int's from memory
    mpf_clear(@a) : mpf_clear(@aa)
    mpf_clear(@b) : mpf_clear(@bb)
    mpf_clear(@t) : mpf_clear(@tt)
    Return pi
End Function

Dim Shared As ZString * 100000000 outtext

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 pistring= Pi_ui(1000000)
t=timer-t
print pistring
print "Time taken ";t
print len(pistring)
sleep 

Actually it makes no difference in speed by optimizing
#cmdline "-gen gcc -Wc -O3"
I am sure the dll is already optimised.
But my question
Why is the gmp dll in along with the fb package?

1.5x

I confirm, the use of -O3 or -Ofast has no impact on the processing speed. it is likely that the library is already optimized...

[dodicat]

But my question
Why is the gmp dll in along with the fb package?

dodicat, I recommend that you install the NTCore Explorer Suite
launch it and drag&drop cc1.exe onto the explorer window, on the left pane click on imports, you will see a list of dependencies

Hi srvaldez.
This might do instead.

Code: Select all

#include "file.bi"
Function tally(somestring As String,partstring As String,arr() As Long) As Long
      Dim As Long i,j,ln,lnp,count
      If Instr(somestring,partstring)=0 Then Return 0
      ln=Len(somestring)
      lnp=Len(partstring)
      Redim arr(1 To ln)'start with a big enough array to handle all circumstances
      count=0
      i=-1
      Do
            i+=1
            If somestring[i] <> partstring[0] Then Continue Do
            If somestring[i] = partstring[0] Then 'might be!
                  For j=0 To lnp-1 
                        If somestring[j+i]<>partstring[j] Then Continue Do 'not this time!
                  Next j
            End If
            ' if got here then a partstring has been found
            count+=1
            arr(count)=i+1 
            i=i+lnp-1
      Loop Until i>=ln-1 
      Redim Preserve arr(1 To count)'resize the output array
      Return count
End Function

Sub showabit(s As String,position As Long,delim As String) ' only to view some characters.
      Dim As String  g
      Locate Csrlin-1,50
      For n As Long=position-20 To position+20 +Len(delim)
            If n=position-1 Then Color 5
            If n=position+Len(delim)-1 Then Color 15
           if n<=len(s)-1 then 
                 if s[n]=0 then s[n]=asc("-")
           Print Chr(s[n]);
           end if
      Next n
      Print
End Sub


Function loadfile(file as string) as String
	If FileExists(file)=0 Then Print file;" not found":Sleep:end
   dim as long  f=freefile
    Open file For Binary Access Read As #f
    Dim As String text
    If Lof(f) > 0 Then
      text = String(Lof(f), 0)
      Get #f, , text
    End If
    Close #f
    return text
end Function


dim as string file="C:\Users\Computer\Desktop\fb\FreeBASIC-1.09.0-winlibs-gcc-9.3.0\bin\libexec\gcc\x86_64-w64-mingw32\9.3.0\cc1.exe"

Dim As String s=loadfile(file)
dim as double t
Randomize

      t=timer
      Color 15
      
      
      Dim As String LookFor=".dll"
      Print "Look for ";LookFor
      
      Redim As Long i()
      Var num=tally(s,LookFor,i())
      
      If num Then
            Print "position ","    found";Tab(59);"string near sought string"
            For n As Long=Lbound(i) To Ubound(i)
                  Print i(n);Tab(20);Mid(s,i(n),Len(LookFor)):showabit(s,i(n),LookFor)
            Next
            Print "tally = ";num
      Else
            Print "not found"
      End If
      Print
      print "Total time taken ";timer-t;"  seconds"
      Print "Press any key to end"
      Sleep
     
     

My result

Code: Select all

Look for .dll
position          found                                   string near sought string
 15453172          .dll                          get_named_event_id-.dll-PLUGIN_START_PARSE_FU
 21249225          .dll                          d----------kernel32.dll-GlobalMemoryStatusEx-
 23286590          .dll                          x☺-ðx☺-ðx☺libgmp-10.dll---¶ðx☺¶ðx☺¶ðx☺libicon
 23286619          .dll                          ☺¶ðx☺¶ðx☺libiconv-2.dll--(ðx☺(ðx☺(ðx☺(ðx☺(ðx☺
 23287522          .dll                          x☺(ðx☺(ðx☺libisl-22.dll---<ðx☺<ðx☺<ðx☺<ðx☺<ðx
 23287621          .dll                          ðx☺<ðx☺<ðx☺libmpc-3.dll----Pðx☺Pðx☺Pðx☺Pðx☺Pð
 23287914          .dll                          x☺Pðx☺Pðx☺libmpfr-6.dll---dðx☺dðx☺dðx☺dðx☺dðx
 23287970          .dll                          x☺dðx☺dðx☺dðx☺zlib1.dll---xðx☺xðx☺xðx☺xðx☺xðx
 23288189          .dll                          ðx☺xðx☺xðx☺KERNEL32.dll----îðx☺îðx☺îðx☺îðx☺îð
 23288671          .dll                          ☺îðx☺îðx☺îðx☺msvcrt.dll--áðx☺áðx☺áðx☺áðx☺áðx☺
 23288732          .dll                          áðx☺libwinpthread-1.dll-┤ðx☺USER32.dll-------
 23288747          .dll                          ad-1.dll-┤ðx☺USER32.dll----------------------
tally =  12

Total time taken  0.1662011999625292  seconds
Press any key to end

 

Remember to put in your own path to cc1.exe
(a bit rough round the edges, but it shows the dll's used in cc1.x, or near enough)

[srvaldez]

ok dodicat, now show me the exports
strangely enough, cc1.exe exports a ton of functions just as if it were a dll

[dodicat]

Code: Select all

dim as string file="C:\Users\Computer\Desktop\fb\FreeBASIC-1.09.0-winlibs-gcc-9.3.0\bin\libexec\gcc\x86_64-w64-mingw32\9.3.0\cc1.exe"
shell "gendef "+file
shell "notepad cc1.def"