;; Common Lisp interval arithmetic and plotting ;;; (c) 1992 by Cindy Ruan and Richard Fateman ;;; Real Intervals are represented as pairs of floats in a class. ;;; left and right endpoints are assumed to be numbers. This is ;;; checked when creating intervals with make-int. ;;; this package assumes left <= right. No exterior intervals for now. ;;; an interval can have endpoints that are +/- infinity or NaN (defclass interval() ((left-value :initarg :left :reader left-value :initform 0) (right-value :initarg :right :reader right-value :initform 0))) (defmethod left ((x interval)) (left-value x)) (defmethod right ((x interval)) (right-value x)) ;; for anything else, (left x) == x etc (defmethod left(x) x) (defmethod right(x) x) (defmethod print-object ((x interval) s) (format s "[~s ~s]" (left x) (right x))) (defmethod make-int((l number) (r number)) (make-instance 'interval :left l :right r)) (defvar onemse (- 1.0 lisp:single-float-epsilon)) (defvar onepse(+ 1.0 lisp:single-float-epsilon)) (defvar onemde (- 1.0d0 lisp:double-float-epsilon)) (defvar onepde(+ 1.0d0 lisp:double-float-epsilon)) ;; to make a floating point number into an interval, bump it ;; up and down by 1 unit in last place. Could be parameterized by epsilon.. (defmethod make-int1((l float)) (typecase l (single-float (if (> l 0.0) (make-int (* onemse l)(* onepse l )) (make-int (* onepse l )(* onemse l)))) (double-float (if (> l 0.0d0) (make-int (* onemde l)(* onepde l )) (make-int (* onepde l)(* onemde l)))))) (defmethod make-int1((l number)) (make-int l l)) ;non-floats are exact (defmethod intervalp ((x interval)) ;; true if X is an interval whose leftvalue <= rightvalue (<= (left X) (right X))) (defmethod intervalp (y) nil) ;;; Interval Computation ;; The following procedures assume that the interval ;; [a, b] is a valid interval ;; i.e. [a, b] = {t | a <= t <= b for all real numbers t} ;; basic operations. We use Allegro Common Lisp conventions here. ;; other CL supporting IEEE arithmetic must have similar conventions. (defvar infinity #.EXCL::*INFINITY-DOUBLE*) (defvar neginfinity #.EXCL::*NEGATIVE-INFINITY-DOUBLE*) (defun efpnp (x) ;; returns nil if x is not exceptional. otherwise returns ;; a SYMBOL with names as indicated below. (excl::exceptional-floating-point-number-p x)) (defun infinityp(x) (let ((ty (efpnp x))) (case ty ((EXCL::*INFINITY-SINGLE* EXCL::*INFINITY-DOUBLE*) t) (t nil)))) (defun nanp(x) (let ((ty (efpnp x))) (case ty ((EXCL::*NAN-SINGLE* EXCL::*NAN-DOUBLE*) t) (t nil)))) (defmethod intadd((A interval) (B interval)) (let ((a1 (left A)) (a2 (right A)) (b1 (left B)) (b2 (right B))) (make-int (+ a1 b1) (+ a2 b2)))) (defmethod intminus((A interval) (B interval)) (let ((a1 (left A)) (a2 (right A)) (b1 (left B)) (b2 (right B))) (make-int (- a1 b2) (- a2 b1)))) (defmethod intmult((A interval) (B interval)) (let* ((a1 (left A)) (a2 (right A)) (b1 (left B)) (b2 (right B)) (v1 (* a1 b1)) (v2 (* a1 b2)) (v3 (* a2 b1)) (v4 (* a2 b2))) (make-int (min v1 v2 v3 v4) (max v1 v2 v3 v4)))) (defmethod intdiv((A interval)(B interval)) (intmult A (inverse B))) (defmethod inverse((B interval)) ; 1/B (let ((b1 (left B)) (b2 (right B))) (cond ((and (>= 0 b1) (<= 0 b2)) ;; 0 is in [b1, b2] (make-int neginfinity infinity)) (t (make-int (/ 1 b1) (/ 1 b2)))))) ;; sin and cos (defvar m1to1 (make-int -1 1)) ;; the interval -1 to 1 (defvar twopi (* 2 pi)) (defvar piby2 (/ pi 2)) ;double precision pi / 2 (defvar intpiby2 (make-int piby2 piby2)) (defmethod intsin((z interval)) (intcos (intminus z intpiby2))) ; sin(x) = cos(x - pi/2) (defmethod intcos ((z interval)) (let ((low (left z)) (hi (right z))) (cond ;; if the interval contains infinity or nan, return [-1 1] ((or (efpnp low) (efpnp hi)) m1to1) ((> low hi) m1to1) ;external interval ((eql low hi) (make-int1 (sin low))) ; convert to an interval (t (let* ((u (cos low)) (v (cos hi)) (m (ceiling low pi)) (n (floor hi pi)) (minval (min u v)) (maxval (max u v))) (if (>= (- n m) 1) (return-from intcos m1to1)) (if (eql m n) ;; 1 or -1 exists between cos(low) and cos(hi) (case (mod m 2) (0 (setq maxval 1)) (1 (setq minval -1)))) (make-int minval maxval)))))) (defmethod inttan ((z interval)) (intdiv (intsin z) (intcos z))) (defmethod intcot ((z interval)) (intdiv (intcos z) (intsin z))) ; need more work..... (defmethod intcsc ((z interval)) ; need more work..... (inverse (intsin z))) (defmethod intsec ((z interval)) (inverse (intcos z))) ;; interval version of exp (defmethod intexp ((A interval)) (make-int (exp (left A)) (exp (right A)))) ;; interval version of log (defmethod intlog ((A interval) &optional (base (exp 1))) ;; assumes base is a positive number. (let ((a1 (left A)) (a2 (right A))) (cond ((< a1 0) (make-int neginfinity infinity)) ((= a1 0) (make-int neginfinity (log a2 base))) (t (make-int (log a1 base) (log a2 base)))))) ;; interval version of expt ;; not really correct right now... (defmethod intexpt ((A interval) (B interval)) (intexp (intmult B (intlog A)))) (defmethod intabs ((A interval)) (let ((a1 (left A)) (a2 (right A))) (cond ((> 0 a2) ;A < 0 (make-int (abs a2) (abs a1))) ((> 0 a1) ;0 is in A (make-int 0 (max (abs a1) a2))) (t A)))) ;; interval plot (defmacro intplot (f (xvar xmin xmax) &key (points 100)) `(block () (with-open-file (xtmp "/tmp/xtmp" :direction :output :if-exists :supersede) (multiple-value-bind (xlist ylist ymin ymax) (intplot1 ',f ',xvar ,xmin ,xmax ,points) (drawboxes xlist ylist ymin ymax xtmp)) (run-shell-command "xgraph /tmp/xtmp" :wait nil) 'XgraphOutput))) (defun drawboxes (xlist ylist ymin ymax xtmp &aux (neginf (- ymin (* 2 (abs ymin)))) ;in case of infinity (posinf (+ ymax (* 2 (abs ymax))))) (do* ((xl xlist (cdr xl)) (yl ylist (cdr yl)) (xint (car xl) (car xl)) (yint (car yl) (car yl))) ((null xl) nil) (let ((xl (coerce (left xint) 'single-float)) (xh (coerce (right xint) 'single-float)) (yl (coerce (left yint) 'single-float)) (yh (coerce (right yint) 'single-float))) (cond ((and (infinityp yl) (infinityp yh)) ;; y goes form +inf to -inf ;; don't draw top and bottom of the box (format xtmp "move ~12,4G ~12,4G ~12,4G ~12,4G move ~12,4G ~12,4G ~12,4G ~12,4G~%" xl posinf xl neginf xh posinf xh neginf)) ((or (nanp yl)(nanp yh)));don't draw anything ((infinityp yl) ;; y goes to -inf (format xtmp "move ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G~%" xl neginf xl yh xh yh xh neginf)) ((infinityp yh) ;; y goes to +inf (format xtmp "move ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G~%" xl posinf xl yl xh yl xh posinf)) (t ;; draw a complete box (format xtmp "move ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G ~12,4G~%" xl yl xl yh xh yh xh yl xl yl)))))) (defun intplot1 (f xvar xmin xmax points) ;; Let y=f(xvar) where xmin <= xvar <= xmax ;; Returns a list of x-intervals, a list of y-intervals, ;; min and max (except +infinity and -infinity) of the range (let (step xlist ylist xint yint ylo yhi (ymax neginfinity) ;; initialize ymax and ymin (ymin infinity) (alist (list (cons xvar xmin)))) ;;check: that xmin, xmax are both finite numbers and that they differ (if (and (numberp xmin)(numberp xmax)(not(>= xmin xmax))) nil (error "Bad interval to intplot ~s, ~s" xmin xmax)) ;;see if rescaling is necessary so that the plot axes etc are ;; useful (if (> (abs (/ (max xmin xmax) (- xmax xmin))) 100) (return-from intplot1 (rescale-intplot1 f xvar xmin xmax points))) (setq step (/ (- xmax xmin) points)) (do ((xlow xmin xhigh) (xhigh (+ xmin step) (+ xhigh step))) ((>= xhigh xmax) ;; set the right value of the last interval to be exactly xmax (setq xint (make-int xlow xmax)) (setq yint (ieval f (list (cons xvar xint)))) (setq ylo (left yint)) (setq yhi (right yint)) (if (and (not (infinityp yhi)) (> yhi ymax)) (setq ymax yhi)) (if (and (not (infinityp ylo)) (< ylo ymin)) (setq ymin ylo)) (setq xlist (cons xint xlist)) (setq ylist (cons yint ylist)) (values xlist ylist ymin ymax)) (setq xint (make-int xlow xhigh)) (setf (cdar alist) xint) ; update alist binding for xvar (setq yint (ieval f alist)) (setq ylo (left yint)) (setq yhi (right yint)) (if (and (not (infinityp yhi)) (> yhi ymax)) (setq ymax yhi)) (if (and (not (infinityp ylo)) (< ylo ymin)) (setq ymin ylo)) (setq xlist (cons xint xlist)) (setq ylist (cons yint ylist))))) (defun rescale-intplot1 (f xvar xmin xmax points) (let ((newf (subst `(+ ,xvar ,xmin) xvar f)) (newxmin 0.0) (newxmax (- xmax xmin))) (format t "~% To make the plot more meaningful~% we rescale ~s ~%to~% ~s" `(intplot ,f (,xvar ,xmin ,xmax)) `(intplot ,newf (,xvar ,newxmin ,newxmax))) (intplot1 newf xvar newxmin newxmax points))) ;; ieval... (setf (get '+ 'intervalmethod) #'intadd) (setf (get '- 'intervalmethod) #'intminus) (setf (get `* 'intervalmethod) #'intmult) (setf (get '/ 'intervalmethod) 'intdiv) (setf (get 'sin 'intervalmethod) 'intsin) (setf (get 'cos 'intervalmethod) 'intcos) (setf (get 'exp 'intervalmethod) 'intexp) (setf (get 'log 'intervalmethod) 'intlog) (setf (get 'expt 'intervalmethod) 'intexpt) (setf (get '^ 'intervalmethod) 'intexpt) (setf (get 'abs 'intervalmethod) 'intabs) (setf (get 'tan 'intervalmethod) 'inttan) (setf (get 'cot 'intervalmethod) 'intcot) (setf (get 'csc 'intervalmethod) 'intcsc) (setf (get 'sec 'intervalmethod) 'intsec) (defun intervalmathp(op) (get op 'intervalmethod)) ;; default case for r being atoms (defmethod ieval (r alist);r is not an interval (cond ((numberp r) (make-int1 r)) ;is r a number? ((let ((s (assoc r alist))) ;is r on alist? (cdr s))) (t r))) ;if r is not on alist return r itself ;; r is an interval (defmethod ieval ((r interval) alist) r) ;; r is a function expression (defmethod ieval ((r cons) alist) (let* ((op (car r)) (intmeth (get op 'intervalmethod))) (if intmeth (apply intmeth (ievalargs (cdr r) alist)) (error "Unknown interval math function.")))) ;; could do this as a local function.. (defun ievalargs (l alist) (mapcar #'(lambda(z) (ieval z alist)) l)) (defun foo (xmin xmax points) (declare (fixnum points)) (let ((width (- xmax xmin)) xhigh xint (result (make-array points))) (do ((count 1 (1+ count)) (xlow xmin xhigh)) ((> count points) result) (declare (fixnum count)) (setq xhigh (+ xmin (/ (* width count) points))) (setq xint (make-int xlow xhigh)) (format t "~%~s" xint) (setf (svref result (1- count)) xint) )))