Pi with Machin's formula (Lisp)
In order to obtain n digits, we will use fixed-point arithmetic to compute π × 10n as a regular integer.
To calculate the arccot of an argument x, we start by dividing the number 1 (represented by 10n, which we provide as the argument unity) by x to obtain the first term. We then repeatedly divide by x2 and a counter value that runs over 3, 5, 7, ..., to obtain each successive term (to avoid recomputation we compute x2 only once). We use two mutually recursive functions, one to handle positive terms and one negative.
When our term value reaches zero, which in this fixed-point representation corresponds to a real value less than 10-n, we stop and return back up the call stack, adding and subtracting terms alternately to get the final sum.
<<arccot>>= (defun arccot (x unity) (let ((xpower (floor (/ unity x)))) (arccot-plus-helper (* x x) 1 xpower))) (defun arccot-plus-helper (xsq n xpower) (let ((term (floor (/ xpower n)))) (if (= term 0) 0 (+ (arccot-minus-helper xsq (+ n 2) (floor (/ xpower xsq))) term)))) (defun arccot-minus-helper (xsq n xpower) (let ((term (floor (/ xpower n)))) (if (= term 0) 0 (- (arccot-plus-helper xsq (+ n 2) (floor (/ xpower xsq))) term))))
Finally, the main function, which uses Machin's formula to compute π using the necessary level of precision:
<<pi.lisp>>= (defun pidigits (digits) (let* ((unity (expt 10 (+ digits 10))) (pi (* 4 (- (* 4 (arccot 5 unity)) (arccot 239 unity))))) (floor (/ pi (expt 10 10)))))
To avoid rounding errors in the result, we use 10 guard digits internally during the calculation. We may now reproduce Machin's result:
>>> pi(100) 31415926535897932384626433832795028841971693993751058209749445923078164062862089 986280348253421170679
The program can be used to compute tens of thousands of digits in just a few seconds on a modern computer. (More sophisticated techniques are necessary to calculate millions or more digits in reasonable time, although in principle this program will also work.)