Generating π in CL

July 22, 2009

Update 2009-07-23 : Faster version in CL and a Haskell version.
——————————————————————————–
A trivial approximation using the Leibniz formula.

(defun leibniz-pi()
  (labels ((local-pi(sum n)
             (if (< n (expt 10 6))
                 (local-pi
                  (+ sum (/ (expt -1 n) (+ (* 2 n) 1)))
                  (+ n 1))
                 sum)))
    (* 4 (local-pi 0f0 0f0))))

And here’s a longer version(but faster and more precise) using Machin’s formula with fixed point arithmetic to x digits.

(defun machin-pi (digits)
  "Calculates PI digits using fixed point arithmetic and Machin's formula with double recursion"
  (labels
      ((arccot-minus (xsq n xpower)
         (let ((term (floor (/ xpower n))))
           (if (= term 0)
               0
               (- (arccot-plus xsq (+ n 2) (floor (/ xpower xsq)))
                  term))))
       (arccot-plus (xsq n xpower)
         (let ((term (floor (/ xpower n))))
           (if (= term 0)
               0
               (+ (arccot-minus xsq (+ n 2) (floor (/ xpower xsq)))
                  term))))
       (arccot (x unity)
         (let ((xpower (floor (/ unity x))))
           (arccot-plus (* x x) 1 xpower))))
    (let* ((unity (expt 10 (+ digits 10)))
           (thispi (* 4 (- (* 4 (arccot 5 unity)) (arccot 239 unity)))))
      (floor (/ thispi (expt 10 10))))))

The first 10000 digits.

* (time (pidigits 10000))
Evaluation took:
  2.496 seconds of real time
  1.149998 seconds of total run time (0.974647 user, 0.175351 system)
  [ Run times consist of 0.325 seconds GC time, and 0.825 seconds non-GC time. ]
  46.07% CPU
  5,626,335,988 processor cycles
  217,893,792 bytes consed
31415926535897932384626433832795028841971693993751058209749445923078164062862089 ...