;;; -*- Mode:Common-Lisp;Package:user; Base:10 -*- ;;; extended AFFINE rational package based on IEEE 754 standard ;;; for binary floating point arithmetic. ;;; This includes rationals ;;; Infinities (=1/0, -1/0) ;;; NotaNumber (0/0) ;;; + signed zeros 0, -0 (= 0/-1) ;;; this last case means that the denominator of an erat might be ;;; negative, but only in the case that the denominator is -1 and ;;; the numerator is 0. ;; newer lisp should have "CLOS" built in. (require "pcl") (use-package "pcl") ;; define a class holding the extended rationals ;; and how to print them. The intention here is to show how little ;; a change would be needed to replace the CL rationals by this (defclass erat() ((numerator :initarg :numerator :reader num :initform 0) (denominator :initarg :denominator :reader den :initform 0))) (defvar nan-rat (make-instance 'erat :numerator 0 :denominator 0)) (defvar inf-rat (make-instance 'erat :numerator 1 :denominator 0)) (defvar minf-rat (make-instance 'erat :numerator -1 :denominator 0)) (defvar mzero-rat (make-instance 'erat :numerator 0 :denominator -1)) (defmethod print-object ((x erat) s) (if (eql (den x) 1) (format s "~s" (num x)) ;; an integer n/1 (format s "~s./~s" (num x)(den x)))) ;; use ./ instead of / in printing ;; create-erat is the top-level function to use to make a new erat. ;; it assumes nothing much about the inputs x and y except that they ;; are either CL real numbers or possibly erats. If the 2nd arg is missing, ;; it is taken to be 1. (deftype realnum()'(and number (not complex))) (defun create-erat(x &optional(y 1));;create x/y as erat, or ratio (cond((eql y 0)(create-erat-simple x y)) ((typep x 'erat)(quotient x (create-erat y))) ((typep y 'erat)(quotient (create-erat x) y)) ((or (not (typep x 'realnum))(not (typep y 'realnum))) (error "~s/~s cannot be converted to erat" x y)) ;; we could actually return a CL rational or integer by just ;; saying (rationalize (/ x y)) here.. ;; we should (but don't) check for float-nans etc. ;; rationalize will give an error for inf or nan stuff. ;; if we had a portable way of checking for this, we could ;; do the coercion. (t(let ((h (rationalize (/ x y)))) ;; handles normal floats etc (create-erat-simple (numerator h)(denominator h)))))) ;; this function should be used whenever possible internally ;; to avoid another gcd computation. (defun create-erat-simple (x y) ;; this assumes integers x and y, where x/y is in lowest terms etc. (case y (0 (cond ((eql x 0) nan-rat) ((> x 0) inf-rat) ((< x 0) minf-rat))) (1 x) ;just return the integer x (t (make-instance 'erat :numerator x :denominator y)))) ;; these accessors make it irrelevant as to whether x is ratio or erat (defmethod numer ((x rational))(numerator x)) (defmethod numer ((x erat))(num x)) (defmethod denom ((x rational))(denominator x)) (defmethod denom ((x erat))(den x)) ;; define extended plus ;; here's a backstop for nonsensical arguments (defmethod plus (a b) `(Plus ,a ,b)) (defmethod plus (x (y erat)) (plus y x)); put erat first if anywhere ;; adds both normal integers or ratios. not used for erats... (defmethod plus ((x rational)(y rational))(create-erat (+ x y))) (defmethod plus ((x erat) y) (let* ((a (numer x))(b (denom x)) (c (numer y))(d (denom y)) ;a/b + c/d = r/s (r(+ (* a d)(* b c))) (s (* b d)) g) (cond((eql s 0) (cond((> r 0) inf-rat) ((< r 0) minf-rat) (t ;; a*d+b*c = 0. Look at the cases. We know that b*d=0. ;; go through the cases (cond ;; b=0 ==> a*d=0. ;; case 1: a=0 ==> 0/0 + ?? -> 0/0 ((eql a 0) nan-rat) ;; case 2: a != 0 ==> d=0. But then b*c = 0. ;; case 3: c=0 ==> ?? + 0/0 -> 0/0 ((eql c 0) nan-rat) ;; that also takes care of the case d=0. ;; case 4: b=d=0 ==> a,c = + or - 1. ((not (eql a b)) nan-rat) ;; 1/0 + -1/0 -> 0/0 ((> a 0) inf-rat) (t minf-rat))))) (t (setq g (gcd r s)) (create-erat-simple (/ r g) (/ s g)))))) ;; define extended times ;; here's a backstop for nonsensical arguments (defmethod times (a b) `(Times ,a ,b)) (defmethod times (x (y erat)) (times y x)); put erat first if anywhere (defmethod times ((x rational)(y rational))(* x y)) (defmethod times ((x erat) y) (let* ((a (numer x))(b (denom x)) (c (numer y))(d (denom y)) ;a/b * c/d = r/s (r (* a c)) (s (* b d)) g) (cond((eql s 0)(create-erat-simple (signum r) s)) (t (setq g (gcd r s)) (create-erat-simple (/ r g) (/ s g)))))) ;; quotient (defun quotient(x y)(times x (reciprocal y))) ;; reciprocal (defmethod reciprocal((x erat)) (create-erat-simple (denom x)(numer x))) (defmethod reciprocal((x rational)) (create-erat-simple (denom x)(numer x))) (defmethod reciprocal(x) `(Power ,x -1)) ;; comparisons (taking IEEE float rules as gospel) ;; .. 1/0 > 0/-1 true ;; ... 0/1 > 0/-1 false ;; ... any > 0/0 false ;; ... 0/0 > any false (defun greaterp(x y) (cond((eql (denom x) 0) (if (eql (denom y) 0) (> (numer x) 0 (numer y)))) ; 1>0> -1 passes ;; do we need to check if (denom y) is zero? ;; no, because if that is the case, (numer y) is one, ;; and the test below is then (> 0 (denom x)). ;; but (denom x) is either 0 or positive, so result is nil (t (> (* (numer x)(denom y))(* (numer y)(denom x)))))) (defun lessp (x y) (cond((eql (denom x) 0) (if (eql (denom y) 0) (< (numer x) 0 (numer y)))) ; -1 <0< 1 passes ;; do we need to check if (denom y) is zero? ;; no, because if that is the case, (numer y) is one, ;; and the test below is then (< 0 (denom x)). ;; but (denom x) is either 0 or positive, so result is nil (t (< (* (numer x)(denom y))(* (numer y)(denom x)))))) ;; 1/0 = any false ;; -1/0 = any false ;; 0/0 = any false ;; 0/1 = 0/-1 true (defun erat= (x y) (and (not (eql (denom y) 0)) ; checking one denom is enough (eql (numer y)(numer x)) ;; if nums are equal then check if (or (eql (denom y)(denom x)) ;; denoms are equal or case of (eql (minusp (denom x))(denom y)) ;; 0/1 = 0/-1 ))) ;; actually IEEE 754 has 26 functionally distinct comparisons composed ;; from the setting of 4 relation flags and an exception flag. We ;; have provided only 3 of them. ;; greater less equal unordered (exception if invalid ;; than than or unorder) ;; > T - - - T ;; < - T - - T ;; = - - T - - ;; There are some subtleties because not-equal is not the same as ;; less-than-or-greater-than. The latter gives an exception for ;; unordered. ;; other items needed include expt ;; coercion to and from other formats ;; difference, remainder, zerop plusp minusp notequal etc. ;; for "normal" rationals, these are already defined. ;; we can use nans for results of ;; (coerce inf-rat 'float) (coerce inf-rat 'single-float) ;; (coerce inf-rat 'double-float) (coerce inf-rat 'long-float) ;; coercion of nan-rat to a float type produces nans of that type. ;; in Franz Inc's allegro CL we have the forms ;; #.excl::*nan-double* #.excl::*infinity-double* ;; #.excl::*negative-infinity-double* ;; #.excl::*nan-single* #.excl::*infinity-single* ;; #.excl::*negative-infinity-single* ;; semantics for the float-nans is not specified in the proposed CL ;; standard. It may be, in fact, something we should specify.... ;; not implemented in this version at all are exception flags. ;; IEEE 754 defines invalid, underflow, overflow, divide-by-zero, inexact. ;; for this model, only divide-by-zero and invalid can even happen.