[gambit-list] Computing a billion decimal digits of pi

Bakul Shah bakul at bitblocks.com
Wed Jul 10 15:33:04 EDT 2013


> On Wed, Jul 10, 2013 at 11:47:48AM -0400, Bradley Lucier wrote:
> > 
> > This is the same program, run on my new desktop:
> > 
> > model name    : Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz
> > 
> > with 16GB of RAM, which seems to run 2-3 times as fast as the server in 
> > the math department.  (Luckily, a new server will be up and running soon.)
	...
> > Chudnovsky's algorithm using binary splitting in Gambit Scheme: digits 
> > 1000000000, CPU time: 14459.672.
> > Last 5 digits 45519.

Neat!

On Wed, 10 Jul 2013 14:28:48 EDT Hendrik Boom <hendrik at topoi.pooq.com> wrote:
> Is it now faster to calculate them than to download them?

A C version of the same algorithm takes about 40 minutes (on
an 3.6Ghz AMD FX with 16GB).  If you have an internet link
that can sustain 417KB/s you can download faster :-)

>From earlier in this thread:
On Wed, 12 Jun 2013 22:25:15 PDT Bakul Shah <bakul at bitblocks.com> wrote:
> On Wed, 12 Jun 2013 22:05:20 CDT Phil Bewig <pbewig at gmail.com> wrote:
> >
> > What is the time complexity of the Chudnovsky algorithm?
>       ...
> It probably wouldn't be hard to figure out the number of
> additions, multiplications etc. as a function of how many
> terms are being evaluated.

I had measured these a while ago:

> (test 1000000)
1000000: 211545 +, 356108 -, 775665 *, 141030 quotients

Roughly:

21% adds
36% subtracts
77% multiplies
14% quotients

These counts scale linearly.

-- Bakul

chud-test.scm:

(define _* *)
(define _+ +)
(define _- -)
(define _quotient quotient)
(define +count 0)
(define -count 0)
(define *count 0)
(define qcount 0)
(set! + (lambda x (set! +count (_+ 1 +count)) (apply _+ x)))
(set! - (lambda x (set! -count (_+ 1 -count)) (apply _- x)))
(set! * (lambda x (set! *count (_+ 1 *count)) (apply _* x)))
(set! quotient (lambda x (set! qcount (_+ 1 qcount)) (apply _quotient x)))

(define (pi digits)
  (let* ((A 13591409)
	 (B 545140134)
	 (C 640320)
	 (C^3 (expt 640320 3))
	 (D 12))
    (define (split m n)
      (if (= 1 (- n m))
	(let ((g (* (- (* 6 n) 5) (- (* 2 n) 1) (- (* 6 n) 1))))
	  (list g (quotient (* C^3 (expt n 3)) 24) (* (expt -1 n) g (+ (* n B) A))))
	(let* ((mid (quotient (+ m n) 2))
	       (gpq1 (split m mid))
	       (gpq2 (split mid n))
	       (g1 (car gpq1)) (p1 (cadr gpq1)) (q1 (caddr gpq1))
	       (g2 (car gpq2)) (p2 (cadr gpq2)) (q2 (caddr gpq2)))
	  (list (* g1 g2) (* p1 p2) (+ (* q1 p2) (* q2 g1))))))
    (let* ((num-terms (inexact->exact (floor (+ 2 (/ digits 14.181647462)))))
	   (sqrt-C (integer-sqrt (* C (expt 100 digits))))
	   (gpq (split 0 num-terms))
	   (g (car gpq)) (p (cadr gpq)) (q (caddr gpq)))
      (quotient (* p C sqrt-C) (* D (+ q (* p A)))))))

(define (print . x) (for-each display x))

(define (test count)
  (set! +count 0)
  (set! -count 0)
  (set! *count 0)
  (set! qcount 0)
  (pi count)
  (print count ": " +count " +, " -count " -, " *count " *, " qcount " quotients\n"))


More information about the Gambit-list mailing list