;; This LISP program computes the continued fraction for pi using
;; W. Gosper's technique.  Update formulas:

;;              r+2+p   p  q            m is leading term ratio as (quo rem)
;;              r+2     r  s            l is last    term ratio as (quo rem)

;;              ar+bq   a  b
;;              cr+dq   c  d                    -S.McD 5/4/83

(defun gcd(a b)                                         ;; gcd(a,b)
        (prog (rem)
        loop    (setq rem (remainder a b))
                (if (zerop rem) (return b))
                (setq a b)(setq b rem)(go loop)))
(defun input()
        (setq tmp a)(setq a (plus (times a r) (times b q)))(setq b tmp)
        (setq tmp c)(setq c (plus (times c r) (times d q)))(setq d tmp)
        (setq tmp r)(setq r (plus r 2))(setq s tmp)
        (setq tmp p)(setq p (plus p r))(setq o q)(setq q tmp)
        (if (greaterp (setq tmp (gcd d o)) 1)
                (if (greaterp (setq tmp (gcd b tmp)) 1)
                (if (greaterp (setq tmp (gcd c tmp)) 1) (setq tmp (gcd a tmp)))))
        (if (greaterp tmp 1)
                (prog ()
                        (setq a (*quo a tmp))
                        (setq b (*quo b tmp))
                        (setq c (*quo c tmp))
                        (setq d (*quo d tmp))))
        (if (plusp d) (setq l (Divide b d))))
(defun output()
        (setq m (car m))
        (setq tmp c)(setq c (difference a (times c m)))(setq a tmp)
        (setq tmp d)(setq d (difference b (times d m)))(setq b tmp)
        (setq l (Divide b d))
        (setq n (sub1 n))
        (print m)(terpr))

(defun decide()
        (prog ()
                (if (zerop c) (return 0))
                (if (zerop d) (return 0))
                (setq m (Divide a c))
                (if (equal (car m) (car l)) (return 1))
                (if (equal (car m) (sub1 (car l)))
                        (if (zerop (cadr l)) (return 1)))
                (setq l m)
                (return 0)))
(defun init()
        (setq p 1)(setq q 1)    (setq a 0)(setq b 4)
        (setq r 1)              (setq c 1)(setq d 0)
        (setq l (Divide a c)))

(defun cfpi(nn)         (setq n nn)             ;; compute first nn coeffs
        (do () ((zerop n))
                (if (plusp (decide)) (output) (input)))
        (patom -99999)(terpr))
(init)(cfpi 65536)(exit)