Original Scheme implementation for Continued fraction/Arithmetic bivariate

(cond-expand
  (r7rs)
  (chicken (import r7rs)))

(define-library (continued-fraction)

  (export define-continued-fraction
          make-continued-fraction
          continued-fraction?
          continued-fraction-ref
          continued-fraction->stream
          stream->continued-fraction
          convergent-stream
          convergent
          number->continued-fraction

          *continued-fraction-max-terms*
          continued-fraction->string

          *ng8-overflow-threshold*
          apply-ng8
          continued-fraction+
          continued-fraction-
          continued-fraction*
          continued-fraction/

          golden-ratio
          silver-ratio
          sqrt2)

  (import (scheme base)
          (scheme case-lambda)
          (scheme inexact))
  (import (srfi 41))

  (begin

    (define-record-type <continued-fraction>
      (%%continued-fraction terminated? ;; No more terms?
                            m           ;; No. of terms memoized.
                            memo        ;; Memoized terms.
                            gen)        ;; Term-generating thunk.
      continued-fraction?
      (terminated? cf-terminated? set-cf-terminated?!)
      (m cf-m set-cf-m!)
      (memo cf-memo set-cf-memo!)
      (gen cf-gen set-cf-gen!))

    (define-syntax define-continued-fraction
      (syntax-rules ()
        ((_ cf gen)
         ;; Make a continued fraction from a generator thunk.
         (define cf (make-continued-fraction gen)))
        ((_ (make-cf args ...) body body* ...)
         ;; Make a procedure that makes a continued fraction from a
         ;; procedure body that returns a generator thunk.
         (define (make-cf args ...)
           (make-continued-fraction
            (lambda (args ...)
              body body* ...))))))

    (define (make-continued-fraction gen)
      (%%continued-fraction #f 0 (make-vector 32) gen))

    (define (%%cf-get-more-terms! cf needed)
      (do () ((or (cf-terminated? cf) (= (cf-m cf) needed)))
        (let ((term ((cf-gen cf))))
          (if term
              (begin
                (vector-set! (cf-memo cf) (cf-m cf) term)
                (set-cf-m! cf (+ (cf-m cf) 1)))
              (set-cf-terminated?! cf #t)))))

    (define (%%cf-update! cf needed)
      (cond ((cf-terminated? cf) (begin))
            ((<= needed (cf-m cf)) (begin))
            ((<= needed (vector-length (cf-memo cf)))
             (%%cf-get-more-terms! cf needed))
            (else ;; Increase the storage space for memoization.
             (let* ((n1 (+ needed needed))
                    (memo1 (make-vector n1)))
               (vector-copy! memo1 0 (cf-memo cf) 0 (cf-m cf))
               (set-cf-memo! cf memo1))
             (%%cf-get-more-terms! cf needed))))

    (define (continued-fraction-ref cf i)
      (%%cf-update! cf (+ i 1))
      (and (< i (cf-m cf))
           (vector-ref (cf-memo cf) i)))

    (define *continued-fraction-max-terms* (make-parameter 20))

    (define continued-fraction->string
      (case-lambda
        ((cf) (continued-fraction->string
               cf (*continued-fraction-max-terms*)))
        ((cf max-terms)
         (let loop ((i 0)
                    (s "["))
           (let ((term (continued-fraction-ref cf i)))
             (cond ((not term) (string-append s "]"))
                   ((= i max-terms) (string-append s ",...]"))
                   (else
                    (let ((separator (case i
                                       ((0) "")
                                       ((1) ";")
                                       (else ",")))
                          (term-str (number->string term)))
                      (loop (+ i 1) (string-append s separator
                                                   term-str))))))))))

    ;;
    ;; I do not demonstrate the "continued-fraction->stream" procedure
    ;; in this program, but I wrote it and here it is.
    ;;
    (define (continued-fraction->stream cf)
      (stream-take-while
       (lambda (term) term)
       (stream-of (continued-fraction-ref cf i)
         (i in (stream-from 0)))))

    ;;
    ;; I do not demonstrate the "stream->continued-fraction" procedure
    ;; in this program, but I wrote it and here it is.
    ;;
    (define-continued-fraction (stream->continued-fraction lazylst)
      (lambda ()
        (stream-match lazylst
          (() #f)
          ((term . rest)
           (set! lazylst rest)
           term))))

    (define (convergent-stream cf)
      ;; Return the convergents as a SRFI-41 stream.
      (let ((cf^ (continued-fraction->stream cf))
            (p0 0)
            (p1 1)
            (q0 1)
            (q1 0))
        (stream-of (let ((numer (+ p0 (* p1 term)))
                         (denom (+ q0 (* q1 term))))
                     (set! p0 p1)
                     (set! p1 numer)
                     (set! q0 q1)
                     (set! q1 denom)
                     (/ numer denom))
          (term in cf^))))

    ;;
    ;; I do not demonstrate the "convergent" procedure in this
    ;; program, but I wrote it and here it is.
    ;;
    (define (convergent cf i)
      ;; Return the ith convergent, where i starts at zero, and the
      ;; continued fraction is (if necessary) extended infinitely with
      ;; implicit infinities.
      (when (negative? i)
        (error "convergent: index negative" i))
      (let loop ((j 0)
                 (p0 0)
                 (p1 1)
                 (q0 1)
                 (q1 0))
        (if (= j (+ i 1)) (/ p1 q1)
            (let ((term (continued-fraction-ref cf j)))
              (if (not term)
                  (if (zero? q1) +inf.0 (/ p1 q1))
                  (let ((numer (+ p0 (* p1 term)))
                        (denom (+ q0 (* q1 term))))
                    (loop (+ j 1) p1 numer q1 denom)))))))

    (define-continued-fraction (number->continued-fraction num)
      (lambda ()
        (and (not (infinite? num))
             (let ((num^ (exact num)))
               (let ((n (numerator num^))
                     (d (denominator num^)))
                 (let-values (((q r) (floor/ n d)))
                   (set! num (if (zero? r) +inf.0 (/ d r)))
                   q))))))

    (define *ng8-overflow-threshold*
      ;; Some very large positive number. Probably it should not be as
      ;; large as this.
      (make-parameter (expt 2 128)))

    (define (%%too-big? arg)
      (let ((threshold (*ng8-overflow-threshold*)))
        (let too? ((arg arg))
          (cond ((null? arg) #f)
                ((pair? arg) (or (too? (car arg))
                                 (too? (cdr arg))))
                (else (>= (abs arg) threshold))))))

    (define-continued-fraction (apply-ng8 ng8-as-list x y)
      (let ((x (if (number? x) (number->continued-fraction x) x))
            (y (if (number? y) (number->continued-fraction y) y))
            (state-ref `(,@ng8-as-list 0 0 #f)))
        (lambda ()
          (define (loop state)
            (let-values (((a12 a1 a2 a b12 b1 b2 b i j overflow?)
                          (apply values state)))
              (cond ((and (zero? b) (zero? b2))
                     (if (and (zero? b1) (zero? b12))
                         #f ;; b12 = b1 = b2 = b = 0
                         (absorb-x-term state)))
                    ((or (zero? b) (zero? b2))
                     (absorb-y-term state))
                    ((zero? b1)
                     (absorb-x-term state))
                    (else
                     (let-values
                         (((q r) (floor/ a b))
                          ((q1 r1) (floor/ a1 b1))
                          ((q2 r2) (floor/ a2 b2))
                          ((q12 r12) (if (zero? b12)
                                         (values +inf.0 0)
                                         (floor/ a12 b12))))
                       (let ((q1-diff (abs (- q1 q)))
                             (q2-diff (abs (- q2 q))))
                         (cond ((> q1-diff q2-diff)
                                (absorb-x-term state))
                               ((< q1-diff q2-diff)
                                (absorb-y-term state))
                               ((not (= q q12))
                                ;;
                                ;; Because q = q1 = q2, the following
                                ;; also are true:
                                ;;
                                ;;   a1/b1 - a/b = r1/b1 - r/b
                                ;;   a2/b2 - a/b = r2/b2 - r/b
                                ;;
                                ;; Rather than compare fractions, we
                                ;; will put the numerators over a
                                ;; common denominator of b*b1*b2, and
                                ;; then compare the new
                                ;; numerators. Doing this might keep
                                ;; the calculations entirely in
                                ;; fixnums.
                                ;;
                                (let ((n (* r b1 b2))
                                      (n1 (* r1 b b2))
                                      (n2 (* r2 b b1)))
                                  (if (> (abs (- n1 n))
                                         (abs (- n2 n)))
                                      (absorb-x-term state)
                                      (absorb-y-term state))))
                               (else ;; q = q1 = q2 = q12
                                (set! state-ref
                                  (list b12 b1 b2 b r12 r1 r2 r
                                        i j overflow?))
                                q))))))))

          (define (absorb-x-term state)
            (let-values (((a12 a1 a2 a b12 b1 b2 b i j overflow?)
                          (apply values state)))
              (let ((term (and (not overflow?)
                               (continued-fraction-ref x i))))
                (if term
                    (let ((ng (list (+ a2 (* a12 term))
                                    (+ a (* a1 term)) a12 a1
                                    (+ b2 (* b12 term))
                                    (+ b (* b1 term)) b12 b1)))
                      (if (not (%%too-big? ng))
                          (loop `(,@ng ,(+ i 1) ,j ,overflow?))
                          (loop (list a12 a1 a12 a1 b12 b1 b12 b1
                                      i j #t))))
                    (loop (list a12 a1 a12 a1 b12 b1 b12 b1
                                (+ i 1) j overflow?))))))

          (define (absorb-y-term state)
            (let-values (((a12 a1 a2 a b12 b1 b2 b i j overflow?)
                          (apply values state)))
              (let ((term (and (not overflow?)
                               (continued-fraction-ref y j))))
                (if term
                    (let ((ng (list (+ a1 (* a12 term)) a12
                                    (+ a (* a2 term)) a2
                                    (+ b1 (* b12 term)) b12
                                    (+ b (* b2 term)) b2)))
                      (if (not (%%too-big? ng))
                          (loop `(,@ng ,i ,(+ j 1) ,overflow?))
                          (loop (list a12 a12 a2 a2 b12 b12 b2 b2
                                      i j #t))))
                    (loop (list a12 a12 a2 a2 b12 b12 b2 b2
                                i (+ j 1) overflow?))))))

          (loop state-ref))))

    (define continued-fraction+
      ;; Not an efficient implementation!
      (case-lambda
        (() (number->continued-fraction 0))
        ((x . rest)
         (let loop ((z (if (number? x)
                           (number->continued-fraction x)
                           x))
                    (lst rest))
           (if (null? lst)
               z
               (loop (apply-ng8 '(0 1 1 0 0 0 0 1) z (car lst))
                     (cdr lst)))))))

    (define continued-fraction-
      ;; Not an efficient implementation!
      (case-lambda
        ((x) (apply-ng8 '(0 -1 0 0 0 0 0 1) x 0))
        ((x . rest)
         (let loop ((z x)
                    (lst rest))
           (if (null? lst)
               z
               (loop (apply-ng8 '(0 1 -1 0 0 0 0 1) z (car lst))
                     (cdr lst)))))))

    (define continued-fraction*
      ;; Not an efficient implementation!
      (case-lambda
        (() (number->continued-fraction 1))
        ((x . rest)
         (let loop ((z (if (number? x)
                           (number->continued-fraction x)
                           x))
                    (lst rest))
           (if (null? lst)
               z
               (loop (apply-ng8 '(1 0 0 0 0 0 0 1) z (car lst))
                     (cdr lst)))))))

    (define continued-fraction/
      ;; Not an efficient implementation!
      (case-lambda
        ((x) (apply-ng8 '(0 0 0 1 0 0 1 0) x x))
        ((x . rest)
         (let loop ((z (if (number? x)
                           (number->continued-fraction x)
                           x))
                    (lst rest))
           (if (null? lst)
               z
               (loop (apply-ng8 '(0 1 0 0 0 0 1 0) z (car lst))
                     (cdr lst)))))))

    (define-continued-fraction golden-ratio
      ;; The golden ratio, (1+sqrt(5))/2.
      (lambda () 1))

    (define-continued-fraction silver-ratio
      ;; The silver ratio, 1+sqrt(2).
      (lambda () 2))

    (define-continued-fraction sqrt2
      ;; The square root of two.
      (let ((term 0))
        (lambda ()
          (set! term (min 2 (+ term 1)))
          term)))

    )) ;; end library (continued-fraction)

(import (scheme base))
(import (scheme write))
(import (continued-fraction))

(import (scheme inexact))
(import (srfi 41))

(newline)

(let* ((x (continued-fraction+ 13/11 1/2))
       (convlst (stream->list (convergent-stream x))))
  (display "13/11 + 1/2 => ")
  (display (continued-fraction->string x))
  (display "\nIts convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(let* ((x (continued-fraction+ 22/7 1/2))
       (convlst (stream->list (convergent-stream x))))
  (display "22/7 + 1/2 => ")
  (display (continued-fraction->string x))
  (display "\nIts convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(let* ((x (continued-fraction/ 22/7 4))
       (convlst (stream->list (convergent-stream x))))
  (display "(22/7)/4 => ")
  (display (continued-fraction->string x))
  (display "\nIts convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(let* ((x (continued-fraction/ sqrt2))
       (convlst (stream->list (stream-take 7 (convergent-stream x)))))
  (display "1/sqrt(2) => ")
  (display (continued-fraction->string x))
  (display "\nIts first 7 convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(let* ((x (continued-fraction+ 1/2 (continued-fraction/ sqrt2 4)))
       (convlst (stream->list (stream-take 7 (convergent-stream x)))))
  (display "(1 + 1/sqrt(2))/2 => ")
  (display (continued-fraction->string x))
  (display "\nIts first 7 convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(let* ((x (continued-fraction* sqrt2 sqrt2))
       (convlst (stream->list (convergent-stream x))))
  (display "sqrt(2)*sqrt(2) => ")
  (display (continued-fraction->string x))
  (display "\nIts convergents:\n  ")
  (display convlst)
  (display "\n  ")
  (display (map inexact convlst))
  (newline))

(newline)

(display "Obviously the last term of the previous continued fraction SHOULD be\n")
(display "'infinity'. We could invent some heuristic to elide the term. For\n")
(display "instance, stopping before any term that exceeds a threshold in size.\n")
(display "In any case, the current implementation DOES terminate, with a very\n")
(display "close answer.\n")

(newline)

(display "Indeed, the answer can be improved by changing a parameter:\n")

(newline)

(parameterize ((*ng8-overflow-threshold* (expt 2 1024)))
  (let* ((x (continued-fraction* sqrt2 sqrt2))
         (convlst (stream->list (convergent-stream x))))
    (display "sqrt(2)*sqrt(2) => ")
    (display (continued-fraction->string x))
    (display "\nIts convergents:\n  ")
    (display convlst)
    (display "\n  ")
    (display (map inexact convlst))
    (newline)))

(newline)