Pi Chudnovsky

Post your FreeBASIC tips and tricks here. Please don’t post your code without including an explanation.
srvaldez
Posts: 2004
Joined: Sep 25, 2005 21:54

Pi Chudnovsky

Postby srvaldez » Sep 11, 2012 19:55

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: http://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_t 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_t, byref qab as mpz_t, byref tab_ as mpz_t)
    /'
    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_t 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_t, byval digits as uinteger)
    /'
    compute int(pi * 10**digits)

    this is done using chudnovsky's series with binary splitting
    '/
    dim as mpz_t 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_t 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
Posts: 88
Joined: Apr 03, 2011 3:44
Location: Inside the bomb
Contact:

Re: Pi Chudnovsky

Postby Destructosoft » Sep 12, 2012 0:17

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.

Return to “Tips and Tricks”

Who is online

Users browsing this forum: No registered users and 2 guests