Scheme Lisp Continued Fraction для квадратных корней

Вероятно, лучший способ сделать это — умножить матрицу, поскольку квадратные корни положительных чисел представлены тривиально. Если a является целым квадратным корнем из N и b = N-a^2, то цепная дробь равна a+b/(2a+b/(2a+b...)). Это может быть представлено произведением бесконечной матрицы ((a b)(1 0)) на бесконечное произведение ((2a b)(1 0)) с любой точностью. Когда у вас будет столько терминов, сколько вы хотите, просто возьмите рациональное в качестве первого столбца.

Вот пример, который должно быть легко расширить.

#lang racket/base
(require math/number-theory)

(define (times left right)
  (define (transpose m)
    (apply map list m))
  (define r (transpose right))
  (for/list ((row left))
    (for/list ((col r))
      (apply + (map * row col)))))

(define iterations (make-parameter 50))

(define (square-root N)
  (define a (integer-root N 2))
  (define b (- N (* a a)))
  (if (zero? b)
      a
      (let* ((first-matrix (list (list a b) (list 1 0)))
             (2a (* 2 a))
             (rest-matrix (list (list 2a b) (list 1 0))))
        (let ((result-matrix
               (for/fold ((accum first-matrix))
                         ((i (in-range (iterations))))
                 (times accum rest-matrix))))
          (apply / (map car result-matrix))))))
 ;;;;;;
Welcome to DrRacket, version 6.12 [3m].
Language: racket/base [custom]; memory limit: 16384 MB.
> (square-root 2)
1 4866752642924153522/11749380235262596085
> (- (sqrt 2) (square-root 2))
0.0
> (square-root 23)
;;;4 and a huge fraction that breaks formatting here
> (- (sqrt 23) (square-root 23))
0.0
> 

Кстати, доступен пакет с непрерывной дробью, но он немного больше сложнее, чем это и построено на последовательностях.

EDIT Поскольку это было помечено как r5rs, я должен был написать это именно так. Я не уверен на 100% в соответствии с Racket/R5RS, но это работает.

#lang r5rs

(define (transpose m)
  (apply map list m))

(define (row*col r c)
  (apply + (map * r c)))

(define (make-row* r)
  (lambda(c) (row*col r c)))

(define (times left right)
  (let ((r (transpose right)))
    (let loop ((rows left))
      (if (null? rows)
          '()
          (cons (map (make-row* (car rows))
                     r)
                (loop (cdr rows)))))))

(define (matrix-expt m e)
  (cond ((= 1 e) m)
        ((even? e) (matrix-expt (times m m)
                                (/ e 2)))
        (else (times m (matrix-expt (times m m)
                                    (/ (- e 1) 2))))))

(define (integer-root N exp)
  (letrec ((recur (lambda(i)
                    (let ((next (+ 1 i)))
                      (if (< N (expt next exp))
                          i
                          (recur next))))))
    (recur 1)))

(define (square-root N iters)
  (let* ((a (integer-root N 2))
         (b (- N (* a a))))
    (if (zero? b)
        a
        (let* ((first-matrix (list (list a b)
                                   (list 1 0)))
               (2a (* 2 a))
               (rest-matrix (list (list 2a b)
                                  (list 1 0)))
               (result-matrix
                (times first-matrix
                       (matrix-expt rest-matrix
                                    iters))))
          (apply / (map car result-matrix))))))