;; Program to calulate the decimal expansion of pi

;; This function calculates the arc cotan of a value scaled by the scale 
;; factor. It uses scaled integer arithmetic to get more precision.
(define (scaled-acot x scale n value)
    (define (calc-term x n)
        (define two-n-minus-1 (- (* n 2) 1))
        (quotient
            (* (if (odd? n) 1 -1) scale)
            (* (expt x two-n-minus-1) two-n-minus-1)))
    (let ((term (calc-term x n)))
        (if (< (abs term) 1)
            value
            (scaled-acot x scale (+ n 1) (+ value term)))))

;; This function sums up the arc cotan function values of the odd terms 
;; in the Fibonacci sequence to produce an approximation of pi/4
(define (sum-acot-even-fibonacci F_even F_odd scale sum)
    (let ((term (scaled-acot F_odd scale 1 0)))
        (if (< (abs term) 1)
            sum
            (sum-acot-even-fibonacci 
                (+ F_even F_odd)        ;; i.e. Fn+1
                (+ F_odd F_even F_odd)  ;; equivalent to (+ Fn+1 Fn), i.e., Fn+2
                scale 
                (+ sum term)))))


(define (calc-extra-digits precision)
    (define (log10 x)
        (/ (log x) (log 10)))
    (inexact->exact 
        (+  2
            (ceiling
                (log10 (* 3 precision))))))


(define (calc-pi precision)
    (let* 
        ((extra-digits (calc-extra-digits precision))
         (scale (expt 10 (+ precision extra-digits))))
        ;; Add a couple of extra decimal places to the precision to take 
        ;; accumulation of error in lower order terms into account. Scale 
        ;; back down for the output
        (quotient
            (* 4 (sum-acot-even-fibonacci 1 2 scale 0))
            (expt 10 extra-digits))))


(begin
    (display (calc-pi 1000))
    (display "\n"))