3
\$\begingroup\$

I'm writing a program that will compute a continued fraction expansion of a number. I've been learning Racket/Scheme for a few weeks and I wanted to see what I can do better.

Usage: This program takes computes an arbitrary number of terms of the continued fraction of an algebraic number (the root of a polynomial with rational coefficients), starting only with an initial decimal approximation.

Example: To compute 100 terms of the continued fraction of the cube root of two, you would enter

(cfrac 100 #e1.2599 (poly -2 0 0 1))

The last argument is read "the polynomial with list of coefficients (-2 0 0 1)" (the polynomial \$x^3 - 2\$, which is the minimal polynomial of cube root of two.)

Remarks: There are some things I haven't gotten around to yet, for instance there is a bug when you try to compute the golden ratio (which should be trivial to fix, the problem is that this is c.fraction with all ones). Also, I will eventually extend this to numbers which are zeros of differentiable functions, like pi. Finally, the algorithm itself is documented very well in the paper Shiu, "Computation of continued fractions without input values".

Should I include comments for each step in the algorithm anyways? The main procedure is pretty big and it seems like it would be unreadable if you weren't already familiar with it.

#lang racket

(define poly list)
(define (sq x) (* x x))

; Evaluate a polynomial (a degree-ascending list of coefficients) using
; Horner's rule
(define (horner-eval x coeffs)
  (foldr (lambda (coeff h.o.t.)
           (+ coeff (* x h.o.t.)))
         0
         coeffs))

(define (deriv poly)
  (define (iter coeff li)
    (if (null? li)
        null
        (cons (* coeff (car li)) (iter (+ coeff 1) (cdr li)))))
  (iter 1 (cdr poly)))

; Returns a function which corresponds to polynomial evaluation
(define (polyfn poly)
  (lambda (x) (horner-eval x poly)))

; Given an *exact* number n, returns the sequence of terms in simple continued
; fraction expansion of n
(define (simple-c-frac n)
  (let ((int-part (floor n)))
    (cond ((= int-part n)
           (list n))
          (else
           (cons int-part
                 (simple-c-frac
                  (/ (- n int-part))))))))

; Given a sequence corresponding to a simple continued fraction, returns the
; rational number it equals
(define (eval-c-frac li)
  (if (null? (cdr li))
      (car li)
      (+ (car li)
         (/ 1 (eval-c-frac (cdr li))))))

; Computes cont. fraction expansion of an algebraic number given only a
; fixed number of decimals. Arguments: No. of terms, an initial approximation
; of the number (must be an exact type), and the polynomial that it satisfies
; (given as a list of coefficients in degree-ascending order).
(define (cfrac terms approx f)
  (define (iter li i x y prx pry)
    (define (next-quot x y pry)
      (let ((complete-quot
             (- (abs (/ ((polyfn (deriv f)) (/ x y))
                        (* (sq y)
                           ((polyfn f) (/ x y)))))
                (/ pry y))))
        (if (= (denominator complete-quot) 1)
            (- complete-quot 1)
            (floor complete-quot))))
    (let ((quot (next-quot x y pry)))
      (if (< i 0)
          li
          (iter (append li (list quot))
                (- i 1)
                (+ (* quot x) prx)
                (+ (* quot y) pry)
                x y))))
  (define (first2 li)
    (list (eval-c-frac (take li 1))
          (eval-c-frac (take li 2))))
  (let ((init (first2 (simple-c-frac approx))))
    (let ((prev (cadr init))
          (pprev (car init)))
      (append (list pprev)
              (list (/ (- prev pprev)))
              (iter null
                    (- terms 3)
                    (numerator prev)
                    (denominator prev)
                    (numerator pprev)
                    (denominator pprev))))))
\$\endgroup\$
1
  • \$\begingroup\$ For future visitors to this page, I've expanded on this code here \$\endgroup\$ – Chris Brooks Feb 7 '17 at 0:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.