Aaron S. Hawley (aaronhawley) wrote,
Aaron S. Hawley
aaronhawley

Frobbing pi in Emacs Lisp

What follows is 500 years of technique used to calculate pi -- in Emacs -- up to 1997. No, Leibniz didn't write Emacs. However, the Lisp available in Emacs is pretty expressive although it has limited floating point support (it lacks bignum). For the traditional mathematical notation, see the article on pi at Wikipedia which this article is based on.

This is 17th century technology for calculating pi in 3 Lisp functions:

(defun pi-Gregory-Leibniz-term (i)
  (/ (* (expt -1 i) 4)
     (1+ (* 2.0 i))))

(defun pi-Gregory-Leibniz-series (n)
  (mapcar (lambda (i) (pi-Gregory-Leibniz-term i))
          (number-sequence 0 n)))

(defun pi-Gregory-Leibniz-sum (n)
  (apply '+ (pi-Gregory-Leibniz-series n)))

After summing 3 million terms (wait for it!) you get the 6th digit of pi.

(pi-Gregory-Leibniz-sum 2900000)
;=> 3.1415929984172606

The numbers after the 6th decimal digit are incorrect.

Here's an iterative version of the above defined in a single function:

(defun pi-Gregory-Leibniz-sum (n) ;; Now 20% faster!
  (let ((sum 0.0))
    (dotimes (i n)
      (setq sum (+ sum
                   (/ (* (expt -1 i) 4)
                      (1+ (* 2.0 i))))))
    sum))

Using the summed series to calculate pi, discovered at about the time calculus was invented, was a simplification of previous formula using infinite products discovered in the 16th century:

(defun pi-Viete (n)
  (let ((prod 1.0)
        (sqrt2 (sqrt 2))
        (sqrt-sum 0))
    (dotimes (i n)
      (setq prod (* prod
                    (sqrt (+ 2 sqrt-sum)))
            sqrt-sum (sqrt (+ 2 sqrt-sum))))
    (/ (expt 2 (1+ n)) prod)))

After 25 iterations, the product series finds the 15th digit.

(pi-Viete 25)
;=> 3.141592653589792

This is early 19th century technology:

(defun pi-Gauss-Legendre (n)
  (let ((a 1.0)
        (b (/ 1 (sqrt 2)))
        (z (/ 1.0 4))
        (p 1.0))
    (dotimes (i n)
      (let* ((a_n+1 (/ (+ a b) 2.0))
             (b_n+1 (sqrt (* a b)))
             (z_n+1 (- z (* p (expt (- a a_n+1) 2))))
             (p_n+1 (* 2 p)))
        (setq a a_n+1
              b b_n+1
              z z_n+1
              p p_n+1)))
    (/ (expt (+ a b) 2)
       (* 4 z))))

In 3 iterations it gets 14 digits of pi.

(pi-Gauss-Legendre 3)
;=> 3.141592653589794

This is 1914 technology (fac is the factorial function):

(defun pi-Ramanujan (n)
  (let ((sum 0.0))
    (dotimes (k n)
      (let ((4k (* 4 k)))
        (setq sum (+ sum
                     (/ (* (fac 4k) (+ 1103 (* 26390.0 k)))
                        (* (expt (fac k) 4) (expt 396.0 4k))))))
    (/ 9801 (sqrt 2) 2 sum))))

It gets 15 digits in 2 iterations.

(pi-Ramanujan 2)
;=> 3.1415926535897936

This route was found in 1987:

(defun pi-Chudnovsky (n)
  (let ((sum 0.0))
    (dotimes (k n)
      (let ((3k (* 3 k))
            (6k (* 6 k))
            (k! (fac k)))
      (setq sum (+ sum
                   (/ (* (fac 6k) (+ 13591409 (* 545140134.0 k)))
                      (* (fac 3k) (expt k! 3) (expt -640320.0 3k))))))
    (/ (* 426880 (sqrt 10005)) sum)))

It gets 14 digits in 1 iteration.

(pi-Chudnovsky 1)
;=> 3.141592653589734

This is 1995 technology:

(defun pi-Ploufe (n)
  (let ((sum 0))
    (dotimes (i n)
      (let ((8i (* 8 i)))
        (setq sum (+ sum
                     (/ (- (/ 4.0 (+ 8i 1))
                           (/ 2.0 (+ 8i 4))
                           (/ 1.0 (+ 8i 5))
                           (/ 1.0 (+ 8i 6)))
                        (expt 16.0 i))))))
    sum))

It can get 14 digits of pi in 11 iterations.

(pi-Ploufe 11)
;=> 3.141592653589793

This is 1997 technology.

(defun pi-Bellard (n)
  (let ((2^2  4.0)
        (2^5  32.0)
        (2^6  64.0)
        (2^8  256.0)
        (2^10 1024.0)
        (sum 0))
    (dotimes (i n)
      (let ((4i (* 4 i))
            (10i (* 10 i)))
        (setq sum (+ sum
                     (*
                      (/ (expt -1 i)
                         (expt 2^10 i))
                      (+ (- (/ 2^5 (1+ 4i)))
                         (- (/ 1.0 (+  4i 3)))
                            (/ 2^8 (+ 10i 1))
                         (- (/ 2^6 (+ 10i 3)))
                         (- (/ 2^2 (+ 10i 5)))
                         (- (/ 2^2 (+ 10i 7)))
                            (/ 1.0 (+ 10i 9))))))))
    (* (/ 1 2^6) sum)))

In 5 iterations it gets 14 digits of precision of pi.

(pi-Bellard 5)
;=> 3.1415926535897927

This article is a follow-up to another Lisp math article, Frobbing primes with Emacs Lisp, which -- if you got this far -- you may also enjoy.

Tags: emacs, free software, programming
Subscribe
  • Post a new comment

    Error

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

    When you submit the form an invisible reCAPTCHA check will be performed.
    You must follow the Privacy Policy and Google Terms of use.
  • 2 comments
Thank you, very interesting.

There seems to be a paren to much in the first defun (pi-Gregory-Leibniz-term) and a missing paren in the pi-Ramanujan defun. Defining fac is left as an exercise to the reader?
Yes, fac is just the factorial function and needs to be user-defined.

Thanks for reporting the typos.