;; -*- Mode:Common-Lisp;Package:user; Base:10 -*-
;;; Common Lisp AFFINE rational number package based on IEEE 754 standard
;;; for binary floating point arithmetic.
;;; (c) copyright 1991, Richard J. Fateman
;;; Part of our intention here is to show how little
;;; a change would be needed to replace the CL rationals by this.
;;; This package includes operations and representations for rationals and
;;; Infinities (=1/0, -1/0)
;;; NotaNumber (0/0)
;;; signed zeros 0, -0 (= 0/-1)
;;; the option of trapping on divide by zero or invalid is provided.
;;; this "signed zero" means that the denominator of an erat might be
;;; negative, but then the denominator is in fact -1 and
;;; the numerator is 0.
;; version 4.0 ACL lisp should have "CLOS" built in, making this unnecessary.
(require "pcl")(use-package "pcl")
;; the user shouldn't use these two global variables
;; except via programs rat-dbyz-catch or rat-invalid-catch
;; or similar programs. So we made their names long.
(defvar *rat-dbyz-trap-enable* nil)
(defvar *rat-invalid-trap-enable* nil)
;; these two global variables can be queried by the user to
;; see if an un-trapped NaN was produced in a sequence of
;; operations.
(defvar *rat-dbyz-flag* nil)
(defvar *rat-invalid-flag* nil)
;; define a class holding the extended rationals
;; and how to print them.
(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.
#+ignore ; see redefinition below
(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
;; non-trapping versions ...
;(defmethod reciprocal((x erat)) (create-erat-simple (denom x)(numer x)))
;(defmethod reciprocal((x rational)) (create-erat-simple (denom x)(numer x)))
;; traps in *rat-dbyz-trap-enable* is non-nil
(defmethod reciprocal((x erat))(recip2 x))
(defmethod reciprocal((x rational))(recip2 x))
(defmethod reciprocal(x) `(Power ,x -1))
(defun recip2(x)
(cond((and *rat-dbyz-trap-enable* (eql (numer x) 0))
(throw 'rat-dbyz 'rat-dbyz))
(t (setf *rat-dbyz-flag* t)
(create-erat-simple (denom x)(numer x)))))
#+ignore ;; if we never trapped or looked at the flags, this would do..
(defun recip2 (x)(create-erat-simple (denom x)(numer x)))
;; 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
#+ignore ;; see redefinition below
(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))))))
#+ignore ;; see redefinition below
(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....
;; IEEE 754 defines invalid, underflow, overflow, divide-by-zero, inexact.
;; for this model, only divide-by-zero and invalid can even happen.
;; in medieval towns, rat-catchers roamed the streets...
;; rat-dbyz-catch is an example of a program which can be
;; used to interact with the exception handling. In this
;; case, (rat-dbyz-catch e val) evaluates the expression
;; e and in case there is a divide by zero, returns the value
;; represented by val. If for some reason you wish a divide-by-zero in (f x)
;; to return "1", use (rat-dbyz-catch (f x) 1) instead of (f x).
(defmacro rat-dbyz-catch (e val)
`(let* ((*rat-dbyz-trap-enable* t)
(res (catch 'rat-dbyz ,e)))
(if (eq res 'rat-dbyz) ,val res)))
;; the only place divide by zero can occur is reciprocal
;; this shows how to catch "invalid"
(defmacro rat-invalid-catch (e val)
`(let* ((*rat-invalid-trap-enable* t)
(res (catch 'rat-invalid ,e)))
(if (eq res 'rat-invalid) ,val res)))
(defun greaterp(x y) ;; sample greaterp version with invalid
(cond((eql (denom x) 0)
(if (eql (denom y) 0)
(> (numer x) 0 (numer y)) ; 1>0> -1 passes. others are invalid
(check-invalid-trap nil)))
((eq y nan-rat)
(check-invalid-trap nil))
(t (> (* (numer x)(denom y))(* (numer y)(denom x))))))
;; this program is called when we are using or producing
;; a NaN. If we are set up to trap on invalid,
;; it throws an to appropriate error catch program.
;; otherwise it sets the global invalid flag and returns x.
;; x is usually nil, but may be NaN.
(defun check-invalid-trap(x)
(cond(*rat-invalid-trap-enable*
(throw 'rat-invalid 'rat-invalid))
(t (setf *rat-invalid-flag* t)
x)))
(defun lessp (x y) ;; checks for invalid
(cond((eql (denom x) 0)
(if (eql (denom y) 0)
(< (numer x) 0 (numer y)) ; -1 <0< 1 passes. others are invalid
(check-invalid-trap nil)))
((eq y nan-rat) (check-invalid-trap nil))
(t (< (* (numer x)(denom y))(* (numer y)(denom x))))))
(defun create-erat-simple (x y) ;; sample create-erat-simple with invalid
;; this assumes integers x and y, where x/y is in lowest terms etc.
(case y
(0 (cond ((eql x 0) (check-invalid-trap 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))))
;; trapping on invalid can happen on any NaN producing operation. This
;; is always done by create-erat-simple. The other place an invalid
;; signal can occur is when an unordered operand is used in comparisons.
;; samples of greaterp and create-erat-simple are shown above.
;; IEEE754 also specifies the setting of global flags on invalid, dbyz
;; and other exceptions so that even if they don't stop execution,
;; they can be tested after a series of operations. Here the flags
;; are called *rat-dbyz-flag* and *rat-invalid-flag* .
;; bottom line: what haven't we done here?
;; THINGS WE MIGHT DO.
;; signalling vs. quiet nans;
;; depending on their purpose, signalling nans could be non-numbers
;; like "uninitialized array element" or a real interval or another
;; number format like arbitrary precision floating point.
;; the other 23 comparisons or a "compare and branch".
;; remainder, difference, copysign, finite, isnan, unordered, class.
;;
;; THINGS WE WOULD NOT DO (irrelevant in the rationals)
;; extra traps, namely inexact, over/underflow.
;; directed rounding modes
;; square-root (not closed in the rationals)
;; binary <-> decimal conversion (done by CL)
;; single/double/extended formats and their combinations
;; denormalized numbers and their consequences
;; scalb, logb, nextafter.
;; next: intervals.