;; A basis for rational arithmetic with infinities ;; constructed by overloading generic arithmetic. ;; Richard Fateman, November, 2005 ;; ;; this is the AFFINE model ;; There are 2 infinities 1/0 and -1/0 and 2 zeros: 0, 0/-1, as well as 0/0. ;; It is kind of unfair to equate these rational numbers to their ;; corresponding names/concepts in floating point. The float infinity ;; is used for numbers that are "too large" or overflow. The rational ;; infinity is a limit. The NaN (not a number) floating point name is ;; a misfit in the rational model, because in a symbolic system like ;; Lisp, MOST objects are not numbers. Calling the 0/0 undefined is ;; not fair either, because it is defined to some extent: (at least we ;; know it is not the same as some particular number, say 43). It is ;; indeterminate, perhaps. But as a consideration for analogical ;; reasoning, we will use the name NaN, but let's let it stand for ;; something else. iNdetermiNAte. I know that's NNA, but to ;; a dyslexic, ohw serac? (defpackage :affine ;uses generic arithmetic (:use :cl ) (:shadowing-import-from :ga "+" "-" "/" "*" "expt" ;binary arith "=" "/=" ">" "<" "<=" ">=" ;binary comparisons "sin" "cos" "tan" ;... more trig "atan" "asin" "acos" ;... more inverse trig "sinh" "cosh" "atanh" ;... more hyperbolic "expt" "log" "exp" "sqrt" ;... more exponential, powers ;; "1-" "1+" "abs" ;; maybe need these too ;; "numerator" "denominator" "two-arg-/" ;; "two-arg-/=" ) ;(:shadow "numerator" "denominator") (:export "aff" ) ) (eval-when (load eval)(require "ga")) ;generic arithmetic (in-package :affine) ;package is projective rationals (defstruct aff n d) ;an affine rat. n is 0 or 1 ;; the constructor for aff is just (/ a b) (defconstant NaN (make-aff :n 0 :d 0)) ; 0/0 (defconstant inf (make-aff :n 1 :d 0)) ; 1/0 (defconstant minf (make-aff :n -1 :d 0)) ; -1/0 (defconstant minz (make-aff :n 0 :d -1)) ; 0/-1 (defmethod print-object ((a aff) stream) (format stream "~a/~a" (aff-n a) (aff-d a))) ;; Non-rational operations on extended rationals are not very useful, though ;; they do provide some help with locating branch cuts. ;; Algebraic and transcendental functions are of course ;; not so great for regular rationals either --- unless you ;; have an approximation in mind. ;; Here are some, anyway. (defmethod ga::sin ((s aff)) NaN) (defmethod ga::cos ((s aff)) NaN) (defmethod ga::tan ((s aff)) NaN) (defmethod ga::exp ((s aff)) NaN) (defmethod ga::log ((s aff)) NaN) (defmethod ga::two-arg-expt((s aff)(n t)) (if (zerop n) 1 s)) ;NaN^0=1 (defmethod ga::two-arg-expt((s t)(n aff)) (cond ((= s 1) 1) ;1^inf, 1^nan ((= s 0) 0) ;0^inf, 0^nan ;; maybe |x|<1 x^inf? (t n))) ;x^inf, x^nan (defmethod ga::atan ((s aff)) NaN) ;;; these all were copied from projective, which is simpler. ;;; but they are not accurate for affine model. ;;; Here is the table of operations for + #| 1 2 3 4 5 6 ------------------------------------- + 0 -0 1/0 -1/0 0/0 x --------------------------------------- 0 0 0 1/0 -1/0 || x -0 0 -0 1/0 -1/0 || x 1/0 1/0 1/0 1/0 0/0 || 1/0 -1/0 -1/0 -1/0 0/0 -1/0 || -1/0 0/0 - - - - - - - - - - - - - - - 0/0 all the same.. y y y 1/0 -1/0 0/0 x+y |# (defun minzp(x)(and (= (aff-n x) 0)(= (aff-d x) -1))) ;minus zero predicate (defun pinfp(x)(and (= (aff-n x) 1)(= (aff-d x) 0))) ;pos inf predicate (defun minfp(x)(and (= (aff-n x) -1)(= (aff-d x) 0))) ;minus inf predicate (defun NaNp(x)(and (= (aff-n x) 0)(= (aff-d x) 0))) ;minus inf predicate (defmethod ga::two-arg-+ ((r aff) (s aff)) (cond ((or (NaNp r)(NaNp s)) r) ;col 5, row 5 ((minzp r) s) ;col 2 ((minzp s) r) ;row 2 ((pinfp r)(if (minfp s) NaN r)) ;col 3 ((pinfp s)(if (minfp r) NaN s)) ;row 3 ((minfp r)(if (pinfp s) NaN r)) ; row 4 ((minfp s)(if (pinfp r) NaN s)) ; col 4 (t (error "add ~s and ~s?" r s)))) (defmethod ga::two-arg-+ ((r aff) (s t)) (if (minzp r)(if (= s 0) 0 s) r)) ; (defmethod ga::two-arg-+ ((s t) (r aff)) (if (minzp r)(if (= s 0) 0 s) r)) ;; checked 1/2/06 #| Here is the table of operations for * 1 2 3 4 5 6 ------------------------------ * 0 -0 1/0 -1/0 0/0 x --------------------------------------- 0 0 -0 0/0 0/0 || 0*s ; s=sign(x) -0 -0 0 0/0 0/0 || -0*s 1/0 0/0 0/0 1/0 -1/0 || s/0 -1/0 0/0 0/0 -1/0 1/0 || -s/0 0/0 - - - - - - - - - - - - - - - 0/0 all the same.. y 0*t -0*t t/0 -t/0 0/0 x*y t=sign y |# (defmethod ga::two-arg-* ((r aff) (s aff)) (let ((n(* (aff-n r)(aff-n s))) (d (* (aff-d r)(aff-d s)))) (if (= d 1) n ;; from (-1) *( -1) (make-aff :n n :d d)))) (defmethod ga::two-arg-* ((r aff) (s t)) (make-aff :n (* (aff-n r)(signum s)) :d (aff-d r))) (defmethod ga::two-arg-* ((s t)(r aff)) (ga::two-arg-* r s)) ;same as above (defmethod ga::two-arg-/ ((r rational) (s (eql 0))) ;;norm/norm. check for norm/0 . This is how we construct INF and UND (make-aff :n (signum r) :d 0)) ;; this does not do division by 0.0 ;;More could be done; check with IEEE float rules affine ;; The names for nan are not standardized. These are the names for allegro cl. #+allegro (defmethod ga::two-arg-/ ((r number) (s (eql 0.0d0))) ;; e.g. floats. ;;norm/norm. check for norm/0 (cond ((zerop r) #.excl::*nan-double*) ;0/0 ((excl::nan-p r) #.excl::*nan-double*) ; nan/0 (t #.excl::*infinity-double*))) (defmethod ga::two-arg-/ ((r t) (s aff)) (* r (invert-aff s))) (defmethod invert-aff((x aff)) (let ((n (aff-n x))) (if (= (aff-d x) -1) minf ; 1/(0/-1) is minus inf (cond ((= n 1) 0) ; 1/ (1/0) = 0 ((= n -1) minz) ; 1/ (-1/0) = -0 (t NaN)) ; 1/(0/0) is 0/0 ))) (defmethod invert-aff((x (eql 0)))inf) (defmethod invert-aff((x t))(/ 1 x)) (defmethod negative-aff((x t))(- x)) (defmethod negative-aff((x aff))(if (= -1 (aff-d x)) 0 ; -0 -> 0 (make-aff :n (- (aff-n x)) ; 0/0 -> 0/0, -1/0 -> 1/0, 1/0 -> -1/0 :d (aff-d x)))) (defmethod ga::two-arg-- ((r t) (s aff)) (+ r (negative-aff s))) ;; These are new: we can't shadow these in :ga. ;;(defmethod aff-numerator((s aff))(aff-n s)) ;;(defmethod aff-numerator((s t))(cl::numerator s)) ;;(defmethod aff-denominator((s aff))(aff-d s)) ;;(defmethod aff-denominator((s t))(cl::denominator s)) ;;(defun numerator(x)(if (aff-p x)(aff-n x) (cl::numerator x))) ;;(defun denominator(x)(if (aff-p x)(aff-d x) (cl::denominator x))) ;; this was copied from projective. It is wrong here. #+ignore (defmacro defcomparison (op) (let ((two-arg (intern (concatenate 'string "two-arg-" (symbol-name op)) :ga ))) `(progn ;; very few compares work. Just notequal. See below (format t "~%defining ~s" ',two-arg) (defmethod ,two-arg ((arg1 aff) (arg2 number)) nil) (defmethod ,two-arg ((arg1 number) (arg2 aff)) nil) (defmethod ,two-arg ((arg1 aff) (arg2 aff)) nil) (compile ',two-arg) ',op))) ;;(defcomparison >) ;;(defcomparison <) ;;(defcomparison <=) ;;(defcomparison >=) ;; ;; No regular or extended number is equal to 1/0 or 0/0, but 0 = -minz (defmethod ga::two-arg-/= ((arg1 aff) (arg2 number)) (not (ga::two-arg-= arg1 arg2))) (defmethod ga::two-arg-/= ((arg1 number) (arg2 aff)) (not (ga::two-arg-= arg1 arg2))) (defmethod ga::two-arg-= ((arg1 (eql 0)) (arg2 aff)) (=(aff-d arg2) -1)) ; 0=-0 ;; what about 0.0? (defmethod ga::two-arg-= ((arg1 number) (arg2 aff)) nil) (defmethod ga::two-arg-= ((arg1 aff) (arg2 number))(ga::two-arg-= arg2 arg1)) (defmethod ga::two-arg-= ((arg1 aff) (arg2 aff)) (and (= (aff-n arg1)(aff-n arg2)) ;1=1 or -1=-1 (= (aff-d arg1)(aff-d arg2)) ;-1=-1 or 0=0 (not (= (aff-d arg1) (aff-n arg1) 0)))) ; neither is 0/0. ;; we can do this: ;; a/b < c/d if ad) (defcomparison <) (defmethod ga::two-arg-<= ((a t)(b t)) (or (= a b)(< a b))) (defmethod ga::two-arg->= ((a t)(b t)) (or (= a b)(> a b))) (defmethod ga::two-arg-/= ((arg1 t) (arg2 t)) (not (ga::two-arg-= arg1 arg2))) ;; IEEE 754 has 26 distinct comparisons, >, <, =, unordered (X settings of flags for exceptions) ;; still to do. min, max via two-arg-min?