;;; -*- Mode:Common-Lisp;Package:user; Base:10 -*- ;;; extended projective rational package ;;; (rationals + Infinity (=1/0) + NaN (0/0)) ;; newer lisp should have "CLOS" built in. (require "pcl") (use-package "pcl") ;;; Define an extended set for rational numbers. ;;; This is a model of the PROJECTIVE LINE. ;;; We add two symbols, nan-rat (0/0) and inf (1/0) ;;; 0/0, meaning "no particular real number". ;;; Projective model has one Infinity. Imagine the ;;; projective numbers on a circle with +/- Infinity opposite 0. ;; infinity ;; _-^-_ ;; / \ ;; ( ) The projective "circle" ;; \__.__/ ;; - 0 + ;; 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)) (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 numbers or possibly erats. If the 2nd arg is missing, ;; it is taken to be 1. (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 'number))(not (typep y 'number))) (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. (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 (if (eql x 0) nan-rat inf-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)(create-erat-simple r s)) (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 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 (defun greaterp(x y) (cond((eql (denom x) 0) nil) ;; 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) nil) ;; 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 erat= (x y) (and (not (eql (denom y) 0)) ; checking one denom is enough (eql (denom y)(denom x)) (eql (numer y)(numer x)))) ;; other items needed include expt ;; coercion to other formats ;; zerop plusp minusp notequal etc. ;; for "normal" rationals, these are already defined. ;; because 1/0 is NOT SIGNED infinity, we'd better 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....