# Another slice of π ?

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