;;;; -*- Mode: Lisp -*- ;;;; common-math.lisp ;;;; ;;;; See file COPYING in the main folder for licensing information. (in-package "CL.EXT.MATH") ;;; Conditions ;;; ========== ;;; undefined-operation (define-condition undefined-operation (error) ; undefined-function ((operator :reader undefined-operation-operator ; name :initarg :operator) (arguments :reader undefined-operation-arguments :initarg :arguments) ) (:report (lambda (uoc stream) (format stream "Undefined operation ~S called with ~S." (undefined-operation-operator uoc) (undefined-operation-arguments uoc)))) (:default-initargs :arguments () :operator nil)) ;;; Constants and Variables ;;; ======================= (defconstant +negative-infinity+ '+negative-infinity+) (defconstant +positive-infinity+ '+positive-infinity+) ;;; (defparameter *ignore-comparison-errors-p* nil) ;;; Generic Operations Interface ;;; ============================ ;;; Boolean Operations ;;; ------------------ (defgeneric .<. (x y) (:method ((x real) (y real)) (cl:< x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (.<. double-float-negative-infinity +negative-infinity+) nil) (:method ((y (eql +negative-infinity+)) (x real)) T) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. nil) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) T) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) nil) (:method ((x real) (y (eql +positive-infinity+))) T) (:method ((y (eql +positive-infinity+)) (x real)) nil) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. nil) ) (defgeneric <. (x) (:method ((n real)) (cl:< n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) (defgeneric .>. (x y) (:method ((x real) (y real)) (cl:> x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (.>. double-float-negative-infinity +negative-infinity+) nil) (:method ((y (eql +negative-infinity+)) (x real)) T) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. nil) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) T) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) nil) (:method ((x real) (y (eql +positive-infinity+))) T) (:method ((y (eql +positive-infinity+)) (x real)) nil) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. nil) ) (defgeneric >. (x) (:method ((n real)) (cl:> n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) (defgeneric .=. (x y) (:method ((x number) (y number)) (cl:= x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (.=. double-float-negative-infinity +negative-infinity+) nil) (:method ((y (eql +negative-infinity+)) (x real)) nil) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. T) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) nil) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) nil) (:method ((x real) (y (eql +positive-infinity+))) nil) (:method ((y (eql +positive-infinity+)) (x real)) nil) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. t) ) (defgeneric =. (x) (:method ((n number)) (cl:= n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) (defgeneric ./=. (x y) (:method ((x number) (y number)) (cl:/= x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (./=. double-float-negative-infinity +negative-infinity+) T) (:method ((y (eql +negative-infinity+)) (x real)) T) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. nil) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) t) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) t) (:method ((x real) (y (eql +positive-infinity+))) t) (:method ((y (eql +positive-infinity+)) (x real)) t) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. nil) ) (defgeneric /=. (x) (:method ((n number)) (cl:/= n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) (defgeneric .<=. (x y) (:method ((x real) (y real)) (cl:<= x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (.<=. double-float-negative-infinity +negative-infinity+) nil) (:method ((y (eql +negative-infinity+)) (x real)) T) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. t) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) nil) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) nil) (:method ((x real) (y (eql +positive-infinity+))) t) (:method ((y (eql +positive-infinity+)) (x real)) nil) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. t) ) (defgeneric <=. (x) (:method ((n real)) (cl:<= n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) (defgeneric .>=. (x y) (:method ((x real) (y real)) (cl:>= x y)) (:method ((x real) (y (eql +negative-infinity+))) ;; Fix behavior for (.<=. double-float-negative-infinity +negative-infinity+) nil) (:method ((y (eql +negative-infinity+)) (x real)) T) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+))) ;; This is per "standard mathematics" and also per IEEE 754 (more or less). ;; "infinity" is defined to be greater that any number except NaN and itself. t) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+))) nil) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+))) t) (:method ((x real) (y (eql +positive-infinity+))) nil) (:method ((y (eql +positive-infinity+)) (x real)) t) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+))) ;; See above. t) ) (defgeneric >=. (x) (:method ((n real)) (cl:>= n)) (:method ((n (eql +negative-infinity+))) t) (:method ((n (eql +positive-infinity+))) t) ) ;;; Standard Dyadic and Monadic Arithmetic Operations ;;; ------------------------------------------------- ;;; ;;; Note that as an implentation choice, every dyadic operation takes ;;; an optional argument that can be used as a "deposit" of the ;;; result. This becomes useful for, e.g., matrix operations. (defgeneric .+. (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) (let ((r (cl:+ x y))) (if (and (is-nan r) *error-on-nan-returning-operations*) (error 'floating-point-invalid-operation :operation '.+. :operands (list x y)) r))) (:method ((x real) (y (eql +negative-infinity+)) &optional r) (declare (ignore r)) +negative-infinity+) (:method ((y (eql +negative-infinity+)) (x real) &optional r) (declare (ignore r)) +negative-infinity+) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) &optional r) (declare (ignore r)) +negative-infinity+) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)) &optional r) (declare (ignore r)) (if *error-on-nan-returning-operations* (error 'floating-point-invalid-operation :operation '.+. :operands (list y x)) nan)) (:method ((x number) (y (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((y (eql +positive-infinity+)) (x number) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) &optional r) (if *error-on-nan-returning-operations* (error 'floating-point-invalid-operation :operation '.+. :operands (list y x)) nan)) ) (defgeneric +. (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) (cl:+ x)) (:method ((x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((x (eql +negative-infinity+)) &optional r) (declare (ignore r)) +negative-infinity+) ) (defgeneric .*. (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) (let ((r (cl:* x y))) (if (and (is-nan r) *error-on-nan-returning-operations*) (error 'floating-point-invalid-operation :operation '.+. :operands (list x y)) r))) (:method ((x real) (y (eql +negative-infinity+)) &optional r) (declare (ignore r)) (let ((s (cl:signum x))) (cond ((cl:plusp s) +negative-infinity+) ((cl:zerop s) (if *error-on-nan-returning-operations* (error 'floating-point-invalid-operation :operation '.+. :operands (list x y)) nan)) (t ; (minusp s) +positive-infinity+)))) (:method ((y (eql +negative-infinity+)) (x real) &optional r) (declare (ignore r)) (.*. x y)) (:method ((y (eql +negative-infinity+)) (x (eql +negative-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((y (eql +negative-infinity+)) (x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((x number) (y (eql +positive-infinity+)) &optional r) (declare (ignore r)) (let ((s (cl:signum x))) (cond ((cl:plusp s) +positive-infinity+) ((cl:zerop s) (if *error-on-nan-returning-operations* (error 'floating-point-invalid-operation :operation '.*. :operands (list x y)) nan)) (t ; (minusp s) +negative-infinity+)))) (:method ((y (eql +positive-infinity+)) (x number) &optional r) (declare (ignore r)) (.*. x y)) (:method ((y (eql +positive-infinity+)) (x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((y (eql +positive-infinity+)) (x (eql +negative-infinity+)) &optional r) +negative-infinity+) ) (defgeneric *. (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) (cl:* x)) (:method ((x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((x (eql +negative-infinity+)) &optional r) (declare (ignore r)) +negative-infinity+) ) (defgeneric .-. (x y &optional r) (:method ((x number) (y number) &optional r) (let ((r (cl:- x y))) (if (and (is-nan r) *error-on-nan-returning-operations*) (error 'floating-point-invalid-operation :operation '.-. :operands (list x y)) r))) (:method ((x t) (y t) &optional r) (.+. x (-. y) r))) (defgeneric -. (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) (cl:- x)) (:method ((x (eql +positive-infinity+)) &optional r) (declare (ignore r)) +negative-infinity+) (:method ((x (eql +negative-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) ) (defgeneric ./. (x y &optional r) (:method ((x number) (y number) &optional r) (declare (ignore r)) (let ((r (cl:/ x y))) (if (and (is-nan r) *error-on-nan-returning-operations*) (error 'floating-point-invalid-operation :operation './. :operands (list x y)) r))) (:method ((x t) (y t) &optional r) (.*. x (/. y) r)) ) (defgeneric /. (x &optional r) (:method ((x number) &optional r) (declare (ignore r)) (cl:/ x)) (:method ((x (eql +positive-infinity+)) &optional (r 'single-float)) ;; Let's abuse the R argument. (coerce 0 r)) (:method ((x (eql +negative-infinity+)) &optional (r 'single-float)) ;; Let's abuse the R argument. (coerce 0 r)) ) ;;; Redefined Common Lisp N-adic operators ;;; -------------------------------------- (defun + (&rest args) (if (null args) (cl:+) (let ((n-args (list-length args))) (declare (type fixnum n-args)) (cond ((cl:= n-args 1) (+. (first args))) ((cl:= n-args 2) (.+. (first args) (second args))) (t (.+. (first args) (apply #'+ (rest args)))) )))) (define-compiler-macro + (&rest args) (if (null args) '(cl:+) (let ((n-args (list-length args))) (cond ((cl:= n-args 1) `(+. ,(first args))) ((cl:= n-args 2) `(.+. ,(first args) ,(second args))) (t `(.+. ,(first args) (+ ,@(rest args)))) ))) ) (defun - (arg1 &rest args) (if (null args) (-. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.-. arg1 (first args)) (.-. arg1 (apply #'+ args)))))) (define-compiler-macro - (arg1 &rest args) (if (null args) `(-. ,arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) `(.-. ,arg1 ,(first args)) `(.-. ,arg1 (+ ,@(rest args))))) )) (defun * (&rest args) (if (null args) (cl:*) (let ((n-args (list-length args))) (cond ((cl:= n-args 1) (*. (first args))) ((cl:= n-args 2) (.*. (first args) (second args))) (t (.*. (first args) (apply #'* (rest args)))) )))) (define-compiler-macro * (&rest args) (if (null args) `(cl:*) (let ((n-args (list-length args))) (cond ((cl:= n-args 1) `(*. ,(first args))) ((cl:= n-args 2) `(.*. ,(first args) ,(second args))) (t `(.*. ,(first args) (* ,@(rest args)))) )))) ;;; TODO: add other compiler macros or move all of them out to a ;;; separate file. (defun / (arg1 &rest args) (if (null args) (/. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (./. arg1 (first args)) (./. arg1 (apply #'* args)))))) (defun = (arg1 &rest args) (if (null args) (=. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.=. arg1 (first args)) (.=. arg1 (apply #'= args)))))) (defun /= (arg1 &rest args) (if (null args) (/=. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (./=. arg1 (first args)) (./=. arg1 (apply #'/= args)))))) (defun < (arg1 &rest args) (if (null args) (<. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.<. arg1 (first args)) (.<. arg1 (apply #'< args)))))) (defun <= (arg1 &rest args) (if (null args) (<=. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.<=. arg1 (first args)) (.<=. arg1 (apply #'<= args)))))) (defun > (arg1 &rest args) (if (null args) (>. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.>. arg1 (first args)) (.>. arg1 (apply #'> args)))))) (defun >= (arg1 &rest args) (if (null args) (>=. arg1) (let ((n-args (list-length args))) (if (cl:= n-args 1) (.>=. arg1 (first args)) (.>=. arg1 (apply #'>= args)))))) ;;; Other Redefined Common Lisp Operators ;;; ------------------------------------- (defgeneric expt (base power &optional r) (:method ((base number) (power number) &optional r) (declare (ignore r)) (cl:expt base power))) (defgeneric zerop (n) (:method ((n number)) (cl:zerop n)) (:method ((x (eql +positive-infinity+))) nil) (:method ((x (eql +negative-infinity+))) nil) ) (defgeneric gcd. (n &optional r) (:method ((n integer) &optional r) (declare (ignore r)) (cl:gcd n)) (:method ((n (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((n (eql +negative-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) ) (defgeneric .gcd. (n m &optional r) (:method ((n integer) (m integer) &optional r) (declare (ignore r)) (cl:gcd n m)) (:method ((n (eql +positive-infinity+)) (m integer) &optional r) (declare (ignore r)) (abs m)) (:method ((n (eql +negative-infinity+)) (m integer) &optional r) (declare (ignore r)) (abs m)) (:method ((m integer) (n (eql +positive-infinity+)) &optional r) (declare (ignore r)) (.gcd. n m)) (:method ((m integer) (n (eql +negative-infinity+)) &optional r) (declare (ignore r)) (.gcd. n m)) ) (defun gcd (&rest args) (case (list-length args) (0 (cl:gcd)) (1 (gcd. (first args))) (2 (.gcd. (first args) (second args))) (t (apply #'gcd (.gcd. (first args) (second args)) (cddr args))))) (defgeneric lcm. (n &optional r) (:method ((n integer) &optional r) (declare (ignore r)) (cl:lcm n)) (:method ((n (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((n (eql +negative-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) ) (defgeneric .lcm. (n m &optional r) (:method ((n integer) (m integer) &optional r) (declare (ignore r)) (cl:lcm n m)) (:method ((n (eql +positive-infinity+)) (m integer) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((n (eql +negative-infinity+)) (m integer) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((m integer) (n (eql +positive-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) (:method ((m integer) (n (eql +negative-infinity+)) &optional r) (declare (ignore r)) +positive-infinity+) ) (defun lcm (&rest args) (case (list-length args) (0 (cl:lcm)) (1 (lcm. (first args))) (2 (.lcm. (first args) (second args))) (t (apply #'lcm (.lcm. (first args) (second args)) (cddr args))))) (defun increasing (x) (declare (type sequence x)) (loop for i from 0 below (1- (length x)) for a1 = (elt x 0) then a2 for a2 = (elt x (1+ i)) always (.<. a1 a2))) (defun decreasing (x) (declare (type sequence x)) (loop for i from 0 below (1- (length x)) for a1 = (elt x 0) then a2 for a2 = (elt x (1+ i)) always (.>. a1 a2))) (defgeneric plusp (x) (:method ((x number)) (cl:plusp x)) (:method ((x (eql +positive-infinity+))) t) (:method ((x (eql +negative-infinity+))) nil)) (defgeneric minusp (x) (:method ((x number)) (cl:minusp x)) (:method ((x (eql +positive-infinity+))) nil) (:method ((x (eql +negative-infinity+))) t)) (defgeneric minimum (x)) (defgeneric maximum (x)) (defgeneric minmax (x)) (defun outer-product (x y &optional result) (.*/. x y result)) (defun inner-product (x y &optional result) (./*. x y result)) (defgeneric */ (x y &rest more-args)) (defgeneric .*/. (x y &optional result)) (defgeneric /* (x y &rest more-args)) (defgeneric ./*. (x y &optional result)) (defun outer-sum (x y &optional result) (.+/. x y result)) (defun inner-sum (x y &optional result) (./+. x y result)) (defgeneric +/ (x y &rest more-args)) (defgeneric .+/. (x y &optional result)) (defgeneric /+ (x y &rest more-args)) (defgeneric ./+. (x y &optional result)) ;;;; end of file -- common-math.lisp