;; Some continued fraction functions.
;;
;; I can't tell you how many times I have written these.  I hope this
;; can be the last.

(defmacro met (bind &rest args)
  `(multiple-value-bind ,(car bind) ,(cadr bind)
			,@args))

(defmacro mv (&rest args)
  `(values ,@args))

;; This is the matrix operation being performed repeatedly:
;;
;; | a b || 0 1 |
;; | c d || 1 x |

(defun fwd-rec (a b c d list)
  (if (consp list)
      (let ((x (car list)))
	(let ((a b)
	      (b (+ a (* x b)))
	      (c d)
	      (d (+ c (* x d))))
	  (fwd-rec a b c d (cdr list))))
    (mv a b c d list)))

(defun fwd (list)
  (fwd-rec 1 0 0 1 list))

(defun frac-rec (n d)
  (if (zerop d) n
    (let ((m (mod n d)))
      (let ((q (/ (- n m) d)))
	(cons q (frac-rec d m))))))

(defun normalize-rec (s qm2 qm1 list)
  (if (consp list)
      (cons qm2 (normalize-rec (- s) qm1 (car list) (cdr list)))
    (if (< s 0) (if (< (abs qm1) 2) (list* (+ qm2 qm1) list) 
		  (list* qm2 (1- qm1) 1 list))
      (list* qm2 qm1 list))))

(defun normalize (list)
  (if (consp list)
      (let ((q1   (car list))
	    (list (cdr list)))
	(if (consp list)
	    (let ((q2   (car list))
		  (list (cdr list)))
	      (normalize-rec -1 q1 q2 list))
	  (list* q1 list)))
    list))

(defun frac (n d)
  (let ((list (frac-rec n d)))
    (normalize list)))

(defun defrac-rec (p n)
  (if (zerop n) nil
    (progn
      (format t "~D ~A~%" n (frac p n))
      (defrac-rec p (1- n)))))

;; We use defrac to make pretty pictures.  It prints the continued
;; fracion for p/n for every value of n less than p.  Note that the
;; continued fractions for modular inverses are list reversals of each
;; other.

(defun defrac (p)
  (defrac-rec p (1- p)))

;; We usually want to minimize the energy in x/y .. this becomes
;; an optimization problem :
;;
;; y + zA => 0
;; x - zB => 0
;;
;; We minimize M(z) = (y + zA)(x - zB) = xy + zAx - zBy - ABz^2
;;
;; To do this, set M'(z) = 0 and solve ..
;;
;; M'(z) = Ax - By - 2ABz = 0
;;
;;         z = (Ax - By)/2AB

(defun unbias (A x B y)
  (let ((z (truncate (/ (- (* A x) (* B y)) (* 2 A B)))))
    (mv (- x (* z B))
	(+ y (* z A)))))
	

;; Given A,B,C find x,y,g such that:
;;
;;   Ax + By = C/g

(defun solve (A B C)
  (met ((x Bx y Ay g) (fwd (frac A B)))
       (let ((x (- x)))
	 ;;(format t "(equal (+ (* ~D ~D) (* ~D ~D)) 1))" Ay x Bx y)
	 (let ((x (* (/ C g) x))
	       (y (* (/ C g) y)))
	   (met ((x y) (unbias Ay x Bx y))
		;;(format t "(equal (* ~D (+ (* ~D ~D) (* ~D ~D))) ~D)" g Ay x Bx y C)
		(mv x y g))))))

;; (solve 3 11 -17) => (equal (* 1 (+ (* 3 -2) (* 11 -1))) -17)
;;
;; (solve (* 7 3) (* 11 3) (* 19 3)) => (equal (* 3 (+ (* 7 -2) (* 11 3))) 57)
;;