;;; -*- Package: USER; Mode: LISP; Syntax: Common-lisp -*-
(in-package "USER")

;;; *****************
;;; Reload this file.
;;; *****************

(defparameter *workfile* "~/lisp/fractions")

(defmacro L (&optional (filename *workfile*))
   `(load ,filename))

(defmacro CL (&optional (filename *workfile*))
    `(progn (compile-file ,filename) 
            (load ,filename)))

;;; ***********************************************************
;;; Continued fractions of rational and floating-point numbers.
;;; ***********************************************************

;; Given a finite list of partial quotients, (a0, a1, ... an), computes 
;; the rational number [a0, a1, ... an].
;; 
;; This function isn't used, but we've got to have it on principle.

(defun eval-finite-cf (partial-quotients)
  (+ (first partial-quotients)
     (if (null (rest partial-quotients))
         0
       (/ (eval-finite-cf (rest partial-quotients))))))


;; The rule for computing P[n+1] from P[n-1] and P[n].
;; n-1 and n are the names of those variables.  This macro
;; updates them to variables for n and n-1, respectively.

(defmacro PQ-advance (n-1 n a)
  `(psetf ,n-1 ,n                       ; parallel assignment
          ,n   (+ ,n-1 (* ,n ,a))))


;; Finds the continued fraction for number, and returns it as a list of partial quotients.
;; There are several optional keyword arguments:
;;    If print? is non-nil (the default), we'll print all the convergents as we go.
;;    epsilon is the accuracy of the number (defaults to its precision).

(defun continued-fraction (number &optional (print? t)
                                  (epsilon (if (rationalp number) 0 (precision number))))
  (assert (> number 0))
  (assert (>= epsilon 0))
  (when print? (format t "Finding continued fraction for ~A  (plus or minus ~A)." 
                       number epsilon))
  (loop ;; The square brackets and dashes below are just part of the variable names, and 
        ;; have no meaning to Common Lisp.  Think of n as -1 at initialization time.
        ;; After we compute a_0 (the first partial quotient), we will update these
        ;; variables to have appropriate values for n=0.

        with epsilon-lo = epsilon
        with epsilon-hi = epsilon

        with P[n-1] = 0     with Q[n-1] = 1   
        with P[n]   = 1     with Q[n]   = 0   

        ;; Get integer part and remainder.  The integer part is
        ;; our next partial quotient, so we must make sure that
        ;; it's uniquely determined within our current precision.

        for (partial-quotient leftover) = (multiple-value-list (floor number))
        for valid? = (= (floor (- number epsilon-lo))
                        (floor (+ number epsilon-hi)))


        when valid?
        collect partial-quotient

        when print?
        do (PQ-advance P[n-1] P[n] partial-quotient)   
           (PQ-advance Q[n-1] Q[n] partial-quotient)
           (format t "~&pq ~3A  conv ~18@A = ~F ~:[~;<--- not guaranteed; between about ~D and ~D~]"
                   partial-quotient  (/ P[n] Q[n]) (float (/ P[n] Q[n])) 
                   (not valid?) (floor (- number epsilon-lo)) (floor (+ number epsilon-hi)))

        until (or (zerop leftover) (not valid?))

        do (psetf epsilon-hi (/ epsilon-lo leftover (- leftover epsilon-lo))
                  epsilon-lo (/ epsilon-hi leftover (+ leftover epsilon-hi))
                  number     (/ leftover))))


;;; *******************************
;;; Continued fractions of sqrt(d).
;;; *******************************

;; Is d a perfect square?

(defun square? (d)
  (= d (expt (isqrt d) 2)))     ;; (isqrt d) gives (floor (sqrt d))


;; Attempts to find the repeating continued fraction expansion of
;; (sqrt d), from its floating-point representation.  
;;
;; Returns NIL if the representation isn't sufficiently precise to do this.
;;
;; Otherwise returns a list (a0 a1 a2 ... am), with am = twice a0,
;; representing the continued fraction [a0, a1, ... am, a1, ... am, a1, ...].

(defun approx-sqrt-cf (d)
  (assert (not (square? d)))
  (loop with cf = (continued-fraction (sqrt d) nil)
        for a in (rest cf)
        for end? = (= a (* 2 (first cf)))
        
        collect a into repeating-part
        until end?
        
        finally (return (when end? (cons (first cf) repeating-part)))))


;; Just like approx-sqrt-cf, but doesn't use floating-point 
;; approximations, and so always comes up with an answer.
;;
;; We represent the successive values of x by formal 
;; expressions of the form (a + (b sqrt(d))) / c.

(defun true-sqrt-cf (d)
  (assert (not (square? d)))
  (loop with (a b c) = (list 0 1 1)
        with first-partial-quotient = (isqrt d)
        with gcd

        for partial-quotient = (truncate (+ a (isqrt (* b b d))) c)

        collect partial-quotient
        
        do (decf a (* c partial-quotient))   ;; subtract partial-quotient from x
           (psetf a (* -1 a c)               ;; take reciprocal via parallel assignment
                  b (* b c)
                  c (- (* b b d) (* a a)))
           (setf gcd (gcd a b c)             ;; put into lowest terms
                 a (/ a gcd)
                 b (/ b gcd)
                 c (/ c gcd))

        until (= partial-quotient (* 2 first-partial-quotient))))  ;; about to repeat


;; Print the output cf of (approx-sqrt-cf d) or (true-sqrt-cf d),
;; in a nice format.

(defun print-sqrt-cf (cf d)
  (if cf
      (format t "~&sqrt(~A) ~10T= [~A~16T/~{~3D~} / ...]" 
              d (first cf) (rest cf))
    (format t "~&sqrt(~3D) = ???   (insufficient floating-point precision)" d))
  (values))       ; return nothing      


;;; ****************
;;; Pell's equation.
;;; ****************

;; Fundamental solution to the Pell or negative Pell equation,
;; depending on the second argument.  Returns a list (x y), 
;; or NIL if there is no solution.
;;
(defun Pell (d &optional (negative? nil))
  (loop with repeating-cf = (true-sqrt-cf d)
        with m = (1- (length repeating-cf))                 ;; length of repeating subsequence
        with target-n = (if negative?
                            (when (oddp m) (- m 1))         ;; NIL if no solution exists
                          (if (evenp m) (- m 1) (- (* 2 m) 1)))
        with cf = (append repeating-cf (rest repeating-cf)) ;; ensures we have terms thru target-m

        with P[n-1] = 0     with Q[n-1] = 1
        with P[n]   = 1     with Q[n]   = 0

        for partial-quotient in cf
        for n from 0
        
        do (PQ-advance P[n-1] P[n] partial-quotient)  ; compute P[n], Q[n]
           (PQ-advance Q[n-1] Q[n] partial-quotient)

        when (eql n target-n)
        do (assert (= (- (* P[n] P[n]) (* d Q[n] Q[n]))
                      (if negative? -1 1)))  ; check solution
           (return (list P[n] Q[n]))))
        

(defun print-Pell (solution d &optional (negative? nil))
  (cond ((not negative?)
         (format t "~&Pell with d = ~3D: fundamental soln x = ~A, ~57Ty = ~A"
                 d (first solution) (second solution)))
        (solution
         (format t "~&Negative Pell with d = ~3D: fundamental soln x = ~A, ~60Ty = ~A"
                 d (first solution) (second solution)))
        (t 
         (format t "~&Negative Pell with d = ~3D: (no solutions)" d))))


;;; *************************
;;; Floating-point precision.
;;; *************************

;; Finds the exponent (an integer) in the internal floating-point 
;; representation.

(defun floating-expt (float)
    (second (multiple-value-list (decode-float float))))


;; Precision (+ or -) of a floating-point number.
;; 
;; The system constant single-float-epsilon gives the smallest
;; positive number whose sum with 1.0 exceeds 1.0.  The precision of
;; 1.0 is exactly this -- that is, the floating-point number 1.0 is
;; used to represent any real between 1.0 - epsilon and 1.0 + epsilon.
;;
;; (One might argue that we could use epsilon/2, but this depends on 
;; whether routines like sqrt are very careful.)
;; 
;; The precision of a general number is found by bitshifting the
;; precision of 1.0.  The function scale-float accomplishes this
;; (taking radix into account -- we might be running on a ternary
;; instead of a binary machine).  
;; 
;; Note: The expression (float 1 float) gives a representation of 1.0 
;; in the same floating-point format as the argument, float.

(defun precision (float)
    (let ((expt (- (floating-expt float) 
                   (floating-expt (float 1 float)))))   
      (scale-float single-float-epsilon expt)))