## Pi Chudnovsky

srvaldez
Posts: 2152
Joined: Sep 25, 2005 21:54

### Pi Chudnovsky

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, binarysplitting and the chudnovsky algorithmsee: http://www.craig-wood.com/nick/articles/pi-chudnovsky/ for more infonick craig-wood <nick@craig-wood.com>'/#include once "gmp.bi"dim shared as mpz_t temp, c3_over_24mpz_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 ifend 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 subdim as double t1, t2dim as mpz_t pimpz_init(@pi)t1=timerpi_chudnovsky_bs(pi,1000000)t2=timergmp_printf(!"%Zd\n",@pi)print "elapsed time: ";t2-t1 print "press RETURN to end ";sleepmpz_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

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.