Another slice of π ?

Mon, Aug 10, 2009

The previous Common Lisp and Haskell functions to generate the digits of PI where only accurate between 10000 and 20000 digits. The algorithm uses an approximation where we discard a certain number of ‘guard’ digits to get an accurate result. Some background regarding how the number of guard digits is calculated : There is ‘wobble’ in the number of contaminated digits depending on the number of input digits. When only the order of magnitude of rounding error is of interest, ulps(units in the last place) and ε (machine precision) may be used interchangeably, since they differ by at most a factor of β (the ‘wobble’ in the number of contaminated digits). For example, when a floating-point number is in error by n ulps, that means that the number of contaminated digits is about log_βn. _{[1] [2]}
Here are new versions of the functions using a guard digits calculation :
Common Lisp :

  "Accurately calculates PI digits using Machin's formula
   with fixed point arithmetic and guard digits."
  (labels
      ((arccot-recur (xsq n xpower op)
         (let ((term (floor xpower n))
               (opfun (if op #'+ #'-)))
           (if (= term 0)
               0
               (funcall opfun
                        (arccot-recur xsq (+ n 2) (floor xpower xsq) (not op))
                        term))))
       (arccot (x unity)
         (let ((xpower (floor (/ unity x))))
           (arccot-recur (* x x) 1 xpower t))))
    (if (> digits 0)
        (let* ((guard (+ 10 (ceiling (log digits 10))))
               (unity (expt 10 (+ digits guard)))
               (thispi (* 4 (- (* 4 (arccot 5 unity)) (arccot 239 unity)))))
          (floor thispi (expt 10 guard))))))

Haskell :

{-
Accurately calculates PI digits using Machin's formula
with fixed point arithmetic and variable guards digits.
-}
arccot :: Integer -> Integer -> Integer
arccot x unity = arccot' 0 (unity `div` x) 1 1
    where
      arccot' sum xpower n sign
          | xpower `div` n == 0 = sum
          | otherwise           =
              arccot' (sum + sign * (xpower `div` n))
                      (xpower `div` (x * x)) (n + 2) (- sign)
machin_pi :: Integer -> Integer
machin_pi digits =
    if digits < 1 then 0 else
        pi' `div` (10 ^ guard)
            where
              guard = 10 + (ceiling (logBase 10 (fromInteger digits)))
              unity = 10 ^ (digits + guard)
              pi' = 4 * (4 * arccot 5 unity - arccot 239 unity)

“What Every Computer Scientist Should Know About Floating-Point Arithmetic”, by David Goldberg, March 1991
1. Relative error and Ulps
2. Guard digits