head     1.1;
access   ;
symbols  ;
locks    ; strict;
comment  @@;


1.1
date     90.02.06.17.26.40;  author ram;  state Exp;
branches ;
next     ;


desc
@@



1.1
log
@Initial revision
@
text
@;;; -*- Mode: Lisp; Package: Lisp; Log: code.log -*-
;;;
;;; **********************************************************************
;;; This code was written as part of the Spice Lisp project at
;;; Carnegie-Mellon University, and has been placed in the public domain.
;;; If you want to use this code or any part of Spice Lisp, please contact
;;; Scott Fahlman (FAHLMAN@@CMUC). 
;;; **********************************************************************
;;;
;;; Routines to handle arithmetic functions not supported in the assembler
;;; routines, especially ratio arithmetic.  Some of the routines have been
;;; swiped from spnum.lisp.

(in-package "LISP")

(export '(numerator denominator gcd lcm logcount float-radix
		    float-sign float-digits float-precision integer-decode-float
		    decode-float scale-float rationalize))


(in-package "SYSTEM")
(export 'primep)

(in-package "EXTENSIONS")
(export '*ignore-floating-point-underflow*)

(in-package "LISP")

(defvar *ignore-floating-point-underflow* NIL
  "The variable *Ignore-floating-point-underflow* specifies whether the
  system should generate an error when floating point underflow occurs.
  The default value is NIL which means generate an error.  A non-NIL value
  means return 0.0 (short or long depending on the operands) when an
  underflow occurs without generating an error.")

;;; Macro to call the gcd miscop.

(defmacro %sp-gcd (x y)
  `(%primitive gcd ,x ,y))

(eval-when (compile)
  (defmacro %sp-make-ratio (x y)
    `(%primitive make-ratio ,x ,y))
)

;;;; Stuff used by the bignum printer:

(defun copy-xnum (x)
  "Might copy a bignum, but doesn't since bignum-fixnum-divide-inplace isn't
  currently destructive."
  x)

;;; Special division for printing et al.
;;;
;;; Both divisor and dividend must be positive, and the divisor must
;;; be of integer-length <= 19 (the closer the better).
;;; 
(defun bignum-fixnum-divide-inplace (x y)
  (truncate x y))

;;;; Logcount:

(defun logcount (x)
  "If X is negative, then the number of 0 bits is returned, otherwise
  number of 1 bits is returned."
  (%primitive logcount x))


;;;; Float stuff:

(defun float-sign (float1 &optional (float2 (float 1 float1)))
  "Returns a floating-point number that has the same sign as
   float1 and, if float2 is given, has the same absolute value
   as float2."
  (if (eq (minusp float1) (minusp float2))
      float2
      (- float2)))

(defun float-digits (f)
  "Returns a non-negative number of radix-b digits used in the
   representation of it's argument.  See Common Lisp: The Language
   by Guy Steele for more details."
  (typecase f
    (short-float 20)
;    (single-float 24)
    (double-float 53)
    (long-float 53)
    (t (error "Float-digits: ~A not a float" f))))

(defun float-precision (f)
  "Returns a non-negative number of significant radix-b digits
   in it's argument."
  (if (zerop f)
      0
      (float-digits f)))

(defun float-radix (f)
  "Returns (as an integer) the radix b of its floating-point
   argument."
  (declare (ignore f))
  2)

(defun integer-decode-float (x)
  "Returns three values:
   1) an integer-scaled version of  the significand.
   2) an exponent to which b must be raised to produce
       the appropriate power for division.
   3) -1 or 1 (i.e. the sign of the argument.)"
  (let ((precision (float-precision x)))
    (multiple-value-bind (f ex s) (decode-float x)
      (values (truncate (scale-float f precision))
	      (- ex precision)
	      s))))

(defun decode-float (f)
  "Returns three values:
   1) a floating-point number representing the
       significand.
   2) an integer representing the exponent.
   3) -1 or 1 (i.e. the sign of the argument.)"
  (%primitive decode-float f))

(defun scale-float (f ex)
  "Returns the value (* f (expt (float b f) e))"
  (%primitive scale-float f ex))

;;; Absolute value of various random numbers.

;;; %sp-abs-ratio accepts a single ratio and returns its absolute value.

(defun %sp-abs-ratio (x)
  (let ((n (numerator x))
	(d (denominator x)))
    (if (minusp n)
	(%sp-make-ratio (- n) d)
	x)))

(defun %sp-abs-complex (x)
  (let ((rx (realpart x))
	(ix (imagpart x)))
    (sqrt (+ (expt rx 2) (expt ix 2)))))

(defun %sp-negate-ratio (x)
  (%sp-make-ratio (- (numerator x)) (denominator x)))

(defun %sp-negate-complex (x)
  (%primitive make-complex (- (realpart x)) (- (imagpart x))))

;;; + Operator on various numbers.

(defun %sp-integer+ratio (x y)
  (let* ((dy (denominator y))
	 (n (+ (* x dy) (numerator y))))
    (%sp-make-ratio n dy)))

(defun %sp-ratio+ratio (x y)
  (let* ((nx (numerator x))
	 (dx (denominator x))
	 (ny (numerator y))
	 (dy (denominator y))
	 (g1 (%sp-gcd dx dy)))
    (if (= g1 1)
	(%sp-make-ratio (+ (* nx dy) (* dx ny)) (* dx dy))
	(let* ((t1 (+ (* nx (truncate dy g1)) (* (truncate dx g1) ny)))
	       (g2 (%sp-gcd t1 g1))
	       (t2 (truncate dx g1)))
	  (cond ((= t1 0) 0)
		((= g2 1)
		 (%sp-make-ratio t1 (* t2 dy)))
		(T (let* ((nn (truncate t1 g2))
			  (t3 (truncate dy g2))
			  (nd (if (= t2 1) t3 (* t2 t3))))
		     (if (= nd 1) nn (%sp-make-ratio nn nd)))))))))

(defun %sp-complex+number (x y)
  (%primitive make-complex
	      (+ (realpart x) y)
	      (imagpart x)))

(defun %sp-number+complex (x y)
  (%primitive make-complex
	      (+ (realpart y) x)
	      (imagpart y)))

(defun %sp-complex+complex (x y)
  (%primitive make-complex
	      (+ (realpart x) (realpart y))
	      (+ (imagpart x) (imagpart y))))

(defun %sp-1+ratio (x)
  (let ((dx (denominator x)))
    (%sp-make-ratio (+ (numerator x) dx) dx)))

(defun %sp-1+complex (x)
  (%primitive make-complex
	      (+ (%primitive realpart x) 1)
	      (%primitive imagpart x)))

;;; - Operator on various numbers.

(defun %sp-integer-ratio (x y)
  (let* ((dy (denominator y))
	 (n (- (* x dy) (numerator y))))
    (%sp-make-ratio n dy)))

(defun %sp-ratio-integer (x y)
  (let* ((dx (denominator x))
	 (n (- (numerator x) (* y dx))))
    (%sp-make-ratio n dx)))

(defun %sp-ratio-ratio (x y)
  (let* ((nx (numerator x))
	 (dx (denominator x))
	 (ny (numerator y))
	 (dy (denominator y))
	 (g1 (%sp-gcd dx dy)))
    (if (= g1 1)
	(%sp-make-ratio (- (* nx dy) (* dx ny)) (* dx dy))
	(let* ((t1 (- (* nx (truncate dy g1)) (* (truncate dx g1) ny)))
	       (g2 (%sp-gcd t1 g1))
	       (t2 (truncate dx g1)))
	  (cond ((= t1 0) 0)
		((= g2 1)
		 (%sp-make-ratio t1 (* t2 dy)))
		(T (let* ((nn (truncate t1 g2))
			  (t3 (truncate dy g2))
			  (nd (if (= t2 1) t3 (* t2 t3))))
		     (if (= nd 1) nn (%sp-make-ratio nn nd)))))))))

(defun %sp-complex-number (x y)
  (%primitive make-complex
	      (- (realpart x) y)
	      (imagpart x)))

(defun %sp-number-complex (x y)
  (%primitive make-complex
	      (- x (realpart y))
	      (- (imagpart y))))

(defun %sp-complex-complex (x y)
  (%primitive make-complex
	      (- (realpart x) (realpart y))
	      (- (imagpart x) (imagpart y))))

(defun %sp-1-ratio (x)
  (let ((dx (denominator x)))
    (%sp-make-ratio (- (numerator x) dx) dx)))

(defun %sp-1-complex (x)
  (%primitive make-complex
	      (- (realpart x) 1)
	      (imagpart x)))

;;; * operator on various random numbers.

(defun %sp-integer*ratio (x y)
  (if (zerop x) 0
      (let* ((ny (numerator y))
	     (dy (denominator y))
	     (gcd (%sp-gcd x dy)))
	(if (= gcd 1)
	    (%sp-make-ratio (* x ny) dy)
	    (let ((nn (* (truncate x gcd) ny))
		  (nd (truncate dy gcd)))
	      (if (= nd 1) nn
		  (%sp-make-ratio nn nd)))))))

(defun %sp-ratio*ratio (x y)
  (let* ((nx (numerator x))
	 (dx (denominator x))
	 (ny (numerator y))
	 (dy (denominator y))
	 (g1 (%sp-gcd nx dy))
	 (g2 (%sp-gcd dx ny))
	 (nn (* (if (= g1 1) nx (truncate nx g1))
		(if (= g2 1) ny (truncate ny g2))))
	 (nd (* (if (= g2 1) dx (truncate dx g2))
		(if (= g1 1) dy (truncate dy g1)))))
    (if (= nd 1) nn (%sp-make-ratio nn nd))))

(defun %sp-number*complex (x y)
  (%primitive make-complex
	      (* x (realpart y))
	      (* x (imagpart y))))

(defun %sp-complex*number (x y)
  (%primitive make-complex
	      (* (realpart x) y)
	      (* (imagpart x) y)))

(defun %sp-complex*complex (x y)
  (let* ((rx (realpart x))
	 (ix (imagpart x))
	 (ry (realpart y))
	 (iy (imagpart y)))
    (%primitive make-complex
		(- (* rx ry) (* ix iy))
		(+ (* rx iy) (* ix ry)))))

;;; / operator on various random numbers.

(defun %sp-bignum/fixnum (x y r)
  (let* ((g (%sp-gcd y r)))
    (if (eq g 1)
	(if (minusp y)
	    (%sp-make-ratio (- x) (- y))
	    (%sp-make-ratio x y))
	(let ((nn (truncate x g))
	      (nd (truncate y g)))
	  (when (minusp nd)
	    (setq nn (- nn))
	    (setq nd (- nd)))
	  (%sp-make-ratio nn nd)))))

(defun %sp-fixnum/bignum (x y)
  (if (zerop x)
      0
      (let ((g (%sp-gcd x y))
	    nn nd)
	(if (eq g 1)
	    (setq nn x
		  nd y)
	    (setq nn (truncate x g)
		  nd (truncate y g)))
	(when (minusp nd)
	  (setq nn (- nn))
	  (setq nd (- nd)))
	(%sp-make-ratio nn nd))))

(defun %sp-bignum/bignum (x y r)
  (let ((g (%sp-gcd y r))
	nn nd)
    (if (eq g 1)
	(setq nn x nd y)
	(setq nn (truncate x g) nd (truncate y g)))
    (when (minusp nd)
      (setq nn (- nn))
      (setq nd (- nd)))
    (%sp-make-ratio nn nd)))

;;; %sp-fixnum/ratio (x y) divides the fixnum x by the ratio y.

(defun %sp-integer/ratio (x y)
  (if (zerop x)
      0
      (let* ((ny (numerator y))
	     (dy (denominator y))
	     (gcd (%sp-gcd x ny))
	     (nn (* (if (eq gcd 1) x (truncate x gcd)) dy))
	     (nd (if (eq gcd 1) ny (truncate ny gcd))))
	(when (minusp nd)
	  (setq nn (- nn))
	  (setq nd (- nd)))
	(if (eq nd 1) nn (%sp-make-ratio nn nd)))))


(defun %sp-ratio/integer (x y)
  (let* ((nx (numerator x))
	 (gcd (%sp-gcd nx y))
	 (nn (if (eq gcd 1) nx (truncate nx gcd)))
	 (nd (* (if (eq gcd 1) y (truncate y gcd)) (denominator x))))
    (when (minusp nd)
      (setq nn (- nn))
      (setq nd (- nd)))
    (%sp-make-ratio nn nd)))

(defun %sp-ratio/ratio (x y)
  (let* ((nx (numerator x))
	 (dx (denominator x))
	 (ny (numerator y))
	 (dy (denominator y))
	 (g1 (%sp-gcd nx ny))
	 (g2 (%sp-gcd dx dy))
	 (nn (* (if (= g1 1) nx (truncate nx g1))
		(if (= g2 1) dy (truncate dy g2))))
	 (nd (* (if (= g2 1) dx (truncate dx g2))
		(if (= g1 1) ny (truncate ny g1)))))
    (cond ((eq nd 1) nn)
	  ((eq nd -1) (- nn))
	  ((minusp nd)
	   (%sp-make-ratio (- nn) (- nd)))
	  (T (%sp-make-ratio nn nd)))))

(defun %sp-number/complex (x y)
  (let* ((ry (realpart y))
	 (iy (imagpart y))
	 (dn (+ (* ry ry) (* iy iy))))
    (%primitive make-complex
		(/ (* x ry) dn)
		(/ (- (* x iy)) dn))))

(defun %sp-complex/number (x y)
  (%primitive make-complex
	      (/ (realpart x) y)
	      (/ (imagpart x) y)))

(defun %sp-complex/complex (x y)
  (let* ((rx (realpart x))
	 (ix (imagpart x))
	 (ry (realpart y))
	 (iy (imagpart y))
	 (dn (+ (* ry ry) (* iy iy))))
    (%primitive make-complex
		(/ (+ (* rx ry) (* ix iy)) dn)
		(/ (- (* ix ry) (* rx iy)) dn))))

;;; Truncate operator

(defun %sp-integer-truncate-ratio (x y)
  (let* ((dy (denominator y)))
    (multiple-value-bind (q r) (truncate (* x dy) (numerator y))
      (if (zerop r)
	  (values q 0)
	  (let* ((g (%sp-gcd r dy)))
	    (values q (if (= g 1) (%sp-make-ratio r dy)
			  (%sp-make-ratio (truncate r g)
					  (truncate dy g)))))))))

(defun %sp-ratio-truncate-integer (x y)
  (let* ((dx (denominator x)))
    (multiple-value-bind (q r) (truncate (numerator x) (* dx y))
      (let* ((g (%sp-gcd r dx)))
	(values q (if (= g 1) (%sp-make-ratio r dx)
		      (%sp-make-ratio (truncate r g)
				      (truncate dx g))))))))

(defun %sp-ratio-truncate-ratio (x y)
  (let* ((nx (numerator x))
	 (dx (denominator x))
	 (ny (numerator y))
	 (dy (denominator y)))
    (multiple-value-bind (q r) (truncate (* nx dy) (* dx ny))
      (if (zerop r)
	  (values q 0)
	  (let* ((nd (* dx dy))
		 (g (%sp-gcd r (* dx dy))))
	    (values q (if (= g 1) (%sp-make-ratio r nd)
			  (%sp-make-ratio (truncate r g)
					  (truncate nd g)))))))))

(defun numerator (x)
  "Returns the numerator of a rational."
  (%primitive numerator x))

(defun denominator (x)
  "Returns the denominator of a rational."
  (%primitive denominator x))

(defun gcd (&rest args)
  "Returns the greatest common divisor of the arguments, which must be
  integers.  Gcd with no arguments is defined to be 0."
  (if (null args) 0
      (do* ((gcd (abs (car args)))
	    (args (cdr args) (cdr args)))
	   ((null args) gcd)
	(cond ((= gcd 1)
	       (dolist (x args)
		 (if (not (integerp x))
		     (error "Gcd - argument not an integer: ~A" x)))
	       (return 1)))
	(if (not (integerp (car args)))
	    (error "Gcd - argument not an integer: ~A" (car args))
	    (setq gcd (%sp-gcd gcd (car args)))))))

;;; Lcm.
;;; Lcm2 is defined this way so that operations won't unnecessarily bignumify.

(defun lcm2 (n m)
  "Least common multiple of two nonzero integers.  No type checking."
  (* (/ (max n m) (%sp-gcd n m)) (min n m)))


(defun lcm (&rest args)
  "Returns the least common multiple of one or more integers."
  (do ((args args (cdr args))
       (lcm 1 (lcm2 lcm (car args))))
      ((null args) lcm)
    (cond ((not (integerp (car args)))
	   (error "Lcm: argument not an integer, ~A"
		   (car args)))
	  ((zerop (car args))				;Result is zero.
	   (dolist (arg (cdr args))
	     (if (not (integerp arg))
		 (error "Lcm: argument not an integer, ~A" arg)))
	   (return 0)))))
;;; Rational and friends.

;;; Build-ratio takes two integer arguments and builds
;;; the rational number which is their quotient.
;;; When reduction is implemented, it will go in here.

(defun build-ratio (x y)
  "Build-ratio takes two integer arguments and builds
   the rational number which is their quotient."
  (multiple-value-bind (q r) (truncate x y)
    (if (zerop r) q
	(let ((gcd (%sp-gcd x y)))
	  (unless (= gcd 1)
	    (setq x (/ x gcd))
	    (setq y (/ y gcd)))
	  (if (minusp y)
	      (%primitive make-ratio (- x) (- y))
	      (%primitive make-ratio x y)))))))

;;; Rational produces a rational number for any numeric argument.
;;; Rational assumed that the floating point is completely accurate.

(defun rational (x)
  "Rational produces a rational number for any numeric argument.
   It assumes that floating-point is completely accurate."
  (typecase x
    (float (multiple-value-bind (f ex) (decode-float x)
	     (let* ((precision (float-precision f))
		    (f (truncate (scale-float f precision))))
	       (if (minusp ex)
		   (build-ratio f (ash 1 (+ precision (abs ex))))
		   (build-ratio (ash f ex) (ash 1 precision))))))
    (rational x)
    (otherwise (error "Argument not a non complex number, ~A."
		      x))))

;;; Rationalize does a rational, but it assumes that floats
;;; are only accurate to their precision, and generates a good
;;; rational aproximation of them.
(defun rationalize (x)
  "Rationalize calls the function rational, but it assumes that
   floats are only accurate to their precision, and generates
   a good rational approximation of them."
  (typecase x
    (rational x)
    (short-float (rationalize-float x short-float-epsilon))
    (long-float (rationalize-float x long-float-epsilon))
    (otherwise (error "Argument not a non complex number, ~A."
		      x))))

;;; Thanks to Kim Fateman, who stole this function rationalize-float
;;; from macsyma's rational. Macsyma'a rationalize was written
;;; by the legendary Gosper (rwg). Gosper is now working for Symbolics.
;;; Guy Steele said about Gosper, "He has been called the
;;; only living 17th century mathematician and is also the best
;;; pdp-10 hacker I know." So, if you can understand or debug this
;;; code you win big.

(defun rationalize-float (x &optional (eps long-float-epsilon))
  (cond ((minusp x) (- (rationalize (- x))))
	((zerop x) 0)
	(t (let ((y ())
		 (a ()))
	     (do ((xx x (setq y (/ 1.0 (- xx (float a x)))))
		  (num (setq a (truncate x))
		       (+ (* (setq a (truncate y)) num) onum))
		  (den 1 (+ (* a den) oden))
		  (onum 1 num)
		  (oden 0 den))
		 ((and (not (zerop den))
		       (not (> (abs (/ (- x (/ (float num x)
					       (float den x)))
				       x))
			       eps)))
		  (/ num den)))))))


;;; ISQRT:  Integer square root - isqrt(n)**2 <= n
;;; Upper and lower bounds on the result are estimated using integer-length.
;;; On each iteration, one of the bounds is replaced by their mean.
;;; The lower bound is returned when the bounds meet or differ by only 1.
;;; Initial bounds guarantee that lg(sqrt(n)) = lg(n)/2 iterations suffice.

(defun isqrt (n)
  "Returns the root of the nearest integer less than
   n which is a perfect square."
  (if (and (integerp n) (not (minusp n)))
      (do* ((lg (integer-length n))
	    (lo (ash 1 (ash (1- lg) -1)))
	    (hi (+ lo (ash lo (if (oddp lg) -1 0))))) ;tighten by 3/4 if possible.
	   ((<= (1- hi) lo) lo)
	(let ((mid (ash (+ lo hi) -1)))
	  (if (<= (* mid mid) n) (setq lo mid) (setq hi mid))))
      (error "Isqrt: ~S argument must be a nonnegative integer" n)))

(defun primep (x)
  "Returns T iff X is a positive prime integer."
  (if (<= x 5)
      (and (>= x 2) (/= x 4))
      (and (not (evenp x))
	   (not (= 0 (rem x 3)))
	   (do ((q 6)
		(r 1)
		(inc 2 (logxor inc 6)) ;; 2,4,2,4...
		(d 5 (+ d inc)))
	       ((or (= r 0) (> d q)) (/= r 0))
	     (multiple-value-setq (q r) (truncate x d))))))
@