#|

A Fast Square Root Program

Copyright (C) 2004 David Greve <TheBeaNerd@yahoo.com>

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

A copy of the GNU General Public License is available at the GNU
website: http://gnu.org

Mon Oct 25 2004,  11:29:37 PM

|#
;;
;; Code for rapidly computing highly accurate rational approximations
;; to square roots of rational numbers.
;;
;; If by convention we say that:
;;
;; s = isqrt(C)
;;
;; d = C - s^2
;;
;; Then the square root of C can be expressed as as the infinite
;; continued fraction:
;;
;;        d
;; s + -----------------
;;             d
;;     2s + ------------
;;                  d
;;          2s + -------
;;               
;;               2s + ..
;;
;; If we designate the tail of this continued fraction using
;;
;; e(n) = nth error term (between 0 and 1)
;;
;; Then a pretty good first order approximation for a square root can
;; be computed as follows:
;;
;;                   d  
;; sqrt(C) ~= s + --------
;;                2s + e(1)
;;
;; We say that the nth term of the continued fraction representation
;; has the form:
;;
;;             d
;; e(n) =  ----------
;;         2s + e(n+1)
;;
;; Where we often drop the subscript on the error term for notational
;; convenience.
;;
;; A generalized expression for a partial evaluation of a continued
;; fraction can be represented as:
;;
;;               A + Be
;; sqrt(C) = s + -------
;;               C + De
;;
;; It is easy to see that the first order approximation given above is
;; an instance of this expression when A = d, B = 0, C = 2s, D = 1.
;; The generalized representation is useful for representing the
;; result of evaluating some number of sucessive terms in the
;; continued fraction representation.
;;
;; Assuming that the above representation is the result of evaluating
;; n terms of the continued fraction for C, then the n+1 term would be
;; computed by substituting the next term into the error expression in
;; the representation.
;;
;;     A + B(d/(2s + e))
;; s + ---------------
;;     C + D(d/(2s + e))
;;
;;     A(2s + e) + Bd
;; s + --------------
;;     C(2s + e) + Dd
;;
;;     (A2s + Bd) + Ae
;; s + -----------------
;;     (C2s + Dd) + Ce
;;
;; Which, we observe, is once again in the general representational
;; form.
;;
;; If the evaluation of the first n terms produced
;;
;; A + Be
;; ------
;; C + De
;;
;; and the evaluation of the next m terms were to produce
;;
;; W + Xe
;; ------
;; Y + Ze
;;
;; Then the evaluation of the first n+m terms would be:
;;
;;       W + Xe 
;; A + B(------)
;;       Y + Ze      A(Y + Ze) + B(W + Xe)
;; ------------  =   ---------------------
;;       W + Xe      C(Y + Ze) + D(W + Xe)
;; C + D(------)
;;       Y + Ze
;;
;; (AY + BW) + (AZ + BX)e
;; ----------------------
;; (CY + DW) + (CZ + DX)e
;;
;; Because the continued fraction representation of the square root is
;; uniform, the evaluation of any n sucessive terms will always produce
;; the same result.
;;
;; We can take advantage of this fact to refine our first order
;; approximation by substituting A,B,C, and D in for W,X,Y and Z in
;; the above expression.  The will double the number of terms we have
;; evaluated.  Of course, this procedure can be repeated again and
;; again, with each iteration of the algorithm doubling the number of
;; significant digits in our representation.


(defun symbol-listp (list)
  (declare (type t list))
  (if (consp list)
      (and (symbolp (car list))
	   (symbol-listp (cdr list)))
    (null list)))

(defun wf-binding (binding)
  (declare (type t binding))
  (and (consp binding)
       (symbol-listp (car binding))
       (consp (cdr binding))
       (null (cddr binding))))

(defmacro met (binding &rest args)
  (declare (type (satisfies wf-binding) binding))
  `(multiple-value-bind 
       ,(car binding)
       ,(cadr binding)
     ,@args))

;; The following code doubles the prescision of a given sqrt
;; representation a number of times dictated by i.  
;;
;; Caution: the size of the values grow exponentially with i

(defun rsqrt-rec (i a b c d)
  (declare (type (integer 0 *) i)
	   (type integer a b c d))
  (if (zerop i) (values a b c d)
    (let ((w a)
	  (x b)
	  (y c)
	  (z d))
      (let ((a (+ (* a y) (* b w)))
	    (b (+ (* a z) (* b x)))
	    (c (+ (* c y) (* d w)))
	    (d (+ (* c z) (* d x))))
	(rsqrt-rec (1- i) a b c d)))))

(defun div (n d)
  (declare (type (integer 0 *) n)
	   (type (integer 1 *) d))
  (let ((m (mod n d)))
    (/ (- n m) d)))

(defun fraction-string-rec (i n d)
  (declare (type (integer 0 *) i n)
	   (type (integer 1 *) d)
	   (type (satisfies stringp) r))
  (if (zerop i) ""
    (let ((q (div n d)))
      (let ((qs (write-to-string q :base 10)))
	(let ((n (* 10 (- n (* q d)))))
	  (concatenate 'string qs (fraction-string-rec (1- i) n d)))))))

(defun fraction-string (i r)
  (declare (type (integer 0 *) i)
	   (type rational r))
  (let ((n (numerator   r))
	(d (denominator r)))
    (let ((q (div n d)))
      (let ((qs (concatenate 'string (write-to-string q :base 10) ".")))
	(let ((n (* 10 (- n (* q d)))))
	  (concatenate 'string qs (fraction-string-rec i n d)))))))

;; The following function computes the rational approximation to the
;; square root (integer) with a precision dictated by i.  

(defun rsqrti (i c)
  (declare (type (integer 0 *) i c))
  (let ((s (isqrt c)))
    (let ((d (- c (* s s))))
      (met ((n x d y) (rsqrt-rec i d 0 (* 2 s) 1))
	   (declare (ignore x y))
	   (+ s (/ n d))))))

;; The following function computes the rational approximation to the
;; square root (rational) with a precision dictated by i.  

(defun rsqrtr (i r)
  (declare (type (integer 0 *) i)
	   (type rational r))
  (let ((n (numerator r))
	(d (denominator r)))
    (let ((sn (rsqrti i n))
	  (sd (rsqrti i d)))
      (/ sn sd))))

(defun newton-rec (i s c)
  (if (zerop i) s
    (let ((s (* (/ 1 2) (+ s (/ c s)))))
      (newton-rec (1- i) s c))))

(defun newton (i c)
  (let ((s (isqrt c)))
    (let ((d (- c (* s s))))
      (newton-rec i (+ s (/ d (* 2 s))) c))))

;; Comparing this algorithm with Newton's method, startig at the same
;; point.

(defun prints (digits i c)
  (declare (type (integer 0 *) digits i)
	   (type rational c))
  (progn
    (format t "S ~2D : ~A~%" i (fraction-string digits (rsqrtr i c)))
    (format t "N ~2D : ~A~%" i (fraction-string digits (newton i c)))
    (format t "~%")
    (if (zerop i) nil
      (prints digits (1- i) c))))

;; It would appear that Newton is ahead by an iteration overall.
;; Phase skewing the results by one iteration, the rsqrtr method
;; appears to gain a few digits per iteration on Newton.

#|
(prints 100 6 1973)

S  6 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087
N  6 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087

S  5 : 44.41846462902561876438107965740906053959497442704659903610246205761940066180_6805283703947239542091268
N  5 : 44.4184646290256187643810796574090605395949744270465990361024620576194006618043686917147360058911830087

S  4 : 44.4184646290256187643810796574090605395_83327294135496111246850857337544379108155113103260155775597834
N  4 : 44.41846462902561876438107965740906053959497442704659903610246205761940066_1610650540744782742034584371

S  3 : 44.418464629025618764_752431232557204024098591484297090342230548605659661647016608600579504981095676558
N  3 : 44.4184646290256187643810796574090605_52241670431822943946992852924285398121711157344287940969052547226

S  2 : 44.41846462_8146721950566598272783899534327268025085249602382778806857857085449370627285274658896791048
N  2 : 44.4184646290256187_67435523789036838423352929675436605397104803285266972754092945518662856464686180991

S  1 : 44.4184_52114124148567022233646060917619843207813905667651972754144711476673949363834982650044981364863
N  1 : 44.4184646_35970603967534128700667457382729830926300611642131212353775669201609339752087257843205655945

S  0 : 44.4_04545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454
N  0 : 44.4_04545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454545454

|#