;;; -*- Mode:LISP; Package:NEW-MATH; Base:10; Readtable:CL -*-
;;;
 
;;; Greatest Common Denominator
;;; Steins algorithm, see Knuth Volume 2, page 321
;;;
;;; based on folowing four rules for positive integers U and V:
;;;
;;;	1) If U and V are both even  then gcd(U,V) = 2 gcd(U/2,V/2)
;;;	2) If U is even and V is odd then gcd(U,V) = gcd(U/2,V)
;;;     3) Euclids algorithm ------------ gcd(U,V) = gcd(U-V,V)
;;;     4) if U and V are odd then U-V is even and abs(U-V) < max(U,V)
;-----------------------------------------------------------------------------

(defun gcd-fixnum (u v)
  (when (not (plusp u))
    (if (minusp u)
	(setq u (- u))
      (return-from gcd-fixnum 1)))
  (when (not (plusp v))
    (if (minusp v)
	(setq v (- v))
      (return-from gcd-fixnum 1)))
  (let ((k 0)
	tt)
    (tagbody
	B1
	   (when (not (hw:32logbitp 0 (hw:32logior u v)))
	     (setq u (ash u -1))
	     (setq v (ash v -1))
	     (setq k (1+ k))
	     (go B1))
	B2
	   (when (hw:32logbitp 0 u)
	       (setq tt (- v))
	       (go B4))
	   (setq tt u)
	B3
	   (setq tt (ash tt -1))
	B4
	   (when (not (hw:32logbitp 0 tt))
	     (go B3))
	B5
	   (if (plusp tt)
	       (setq u tt)
	     (setq v (- tt)))
	B6
	   (setq tt (- u v))
	   (when (not (zerop tt))
	     (go B3))
	   (return-from gcd-fixnum
	     (hw:ldb (hw:32logical-shift-up u k) vinc:%%fixnum-field 0)))))

;-----------------------------------------------------------------------------
;   GCD Bignum
;-----------------------------------------------------------------------------
;;;; @@@ This will cons up the wazoo,,, if you get inspired then write a better one!
(defun gcd-bignum (u v)
  (unless (plusp u)
    (if (minusp u)
	(setq u (- u))
      (return-from gcd-bignum 1)))
  (unless (plusp v)
    (if (minusp v)
	(setq v (- v))
      (return-from gcd-bignum 1)))
  (let ((k 0)
	tt)
    (tagbody
	B1
	   (when (and (li:evenp u) (li:evenp v))
	     (setq u (ash u -1))
	     (setq v (ash v -1))
	     (setq k (1+ k))
	     (go B1))
	B2
	   (when (li:oddp u)
	       (setq tt (- v))
	       (go B4))
	   (setq tt u)
	B3
	   (setq tt (ash tt -1))
	B4
	   (when (li:evenp tt)
	     (go B3))
	B5
	   (if (plusp tt)
	       (setq u tt)
	     (setq v (- tt)))
	B6
	   (setq tt (- u v))
	   (when (not (zerop tt))
	     (go B3))
	   (return-from gcd-bignum
	     (* u (li::^ 2 k))))))
;;; Changed to avoid ASH bug ||| Jim 10/10/88
;;; Please fix the ASH bug.
;	     (ash u k)))))

;-----------------------------------------------------------------------------
;   GCD
;-----------------------------------------------------------------------------
(defun gcd (&rest numbers)
  (if (null numbers)
      0
    (do ((value (li:abs (cons:car numbers)))
	 (rest (cons:cdr numbers) (cons:cdr rest)))
	((null rest) value)
      (setq value (gcd-generic value (li:abs (cons:car rest)))))))      

;-----------------------------------------------------------------------------
;   Least Common Multiple
;-----------------------------------------------------------------------------
(defun lcm (number &rest numbers)
  (do ((value (li:abs number))
       (rest numbers (cons:cdr rest)))
      ((null rest) value)
    (let ((next-number (li:abs (cons:car rest))))
      (setq value (if (or (zerop value) (zerop next-number))
		      (return 0)
		    (truncate-generic (* value next-number)
				      (gcd-generic value next-number)))))))

;-----------------------------------------------------------------------------
;   Numerator function
;-----------------------------------------------------------------------------
(defun numerator (x)
  (prims:dispatch vinc::%%data-type x
		  ($$dtp-fixnum x)
		  ($$dtp-bignum x)
		  ($$dtp-rational (hw:vma-start-read-vma-boxed-md-boxed x)
				  (hw:read-md))
		  (t (li:error "You can't get the numerator of that!"))))

;-----------------------------------------------------------------------------
;  Denominator function
;-----------------------------------------------------------------------------
(defun denominator (x)
  (prims:dispatch vinc::%%data-type x
		  ($$dtp-fixnum 1)
		  ($$dtp-bignum 1)
		  ($$dtp-rational (hw:vma-start-read-cdr-vma-boxed-md-boxed x)
				  (hw:read-md))
		  (t (li:error "You can't get the denominator of that!"))))

;-----------------------------------------------------------------------------
;   add rational
;-----------------------------------------------------------------------------
(defun add-rational (x y)
  (let* ((denominator-x (denominator x))
	 (denominator-y (denominator y)))
    (make-canonical-rational
      (+ (multiply-generic (numerator x) denominator-y)
	 (multiply-generic (numerator y) denominator-x))
      (multiply-generic denominator-x denominator-y))))

;-----------------------------------------------------------------------------
;   subtract rational
;-----------------------------------------------------------------------------
(defun subtract-rational (x y)
  (let* ((denominator-x (denominator x))
	 (denominator-y (denominator y)))
    (make-canonical-rational
      (- (multiply-generic (numerator x) denominator-y)
	 (multiply-generic (numerator y) denominator-x))
      (multiply-generic denominator-x denominator-y))))

;-----------------------------------------------------------------------------
;   multiply rational
;-----------------------------------------------------------------------------
(defun multiply-rational (x y)
  (make-canonical-rational
    (multiply-generic (numerator x) (numerator y))
    (multiply-generic (denominator x) (denominator y))))

;-----------------------------------------------------------------------------
;   divide rational
;-----------------------------------------------------------------------------
(defun divide-rational (x y)
  (make-canonical-rational
    (multiply-generic (numerator x) (denominator y))
    (multiply-generic (denominator x) (numerator y))))

;-----------------------------------------------------------------------------
;   make canonical rational
;-----------------------------------------------------------------------------
(defun make-canonical-rational (numerator denominator)
  (when (zerop denominator)
    (li:error "Can't make rational with zero denominator"))
  (when (minusp denominator)
    (setq numerator   (- numerator))
    (setq denominator (- denominator)))
  (if (zerop numerator)
      0
    (let* ((gcd             (gcd-generic numerator denominator))
	   (new-numerator   (truncate numerator   gcd))
	   (new-denominator (truncate denominator gcd)))
      (values
	(if (= 1 new-denominator)
	    new-numerator
	  (hw:dpb-boxed vinc:$$dtp-rational vinc::%%data-type
			(cons:cons new-numerator new-denominator)))
	(multiple-value-bind (crud status)
	    (test-generic new-numerator)
	  (declare (ignore crud))
	  status)))))

;-----------------------------------------------------------------------------
;   negate rational
;-----------------------------------------------------------------------------
(defun negate-rational (x)
  (let ((new-numerator (- (numerator x))))
    (values (hw:dpb-boxed vinc:$$dtp-rational vinc::%%data-type
			  (cons:cons new-numerator
				     (denominator x)))
	    (multiple-value-bind (crud status)
		(test-generic new-numerator)
	      (declare (ignore crud))
	      status))))

;-----------------------------------------------------------------------------
;   test rational
;-----------------------------------------------------------------------------
(defun test-rational (x)
  (test-generic (numerator x)))

;-----------------------------------------------------------------------------
;   compare  rational
;-----------------------------------------------------------------------------
(defun compare-rational (x y)
  (test-generic (numerator (subtract-rational x y))))