Commit 95f5e98f authored by Marco Antoniotti's avatar Marco Antoniotti

Initial checkin.

parents
*.*fasl
docs/
Copyright (c) 2018 Marco Antoniotti
All rights reserved.
Permission is hereby granted, without written agreement and without
license or royalty fees, to use, copy, modify, and distribute this
software and its documentation for any purpose, provided that the
above copyright notice and the following two paragraphs appear in all
copies of this software.
IN NO EVENT SHALL THE AUTHOR(S) BE LIABLE TO ANY PARTY FOR DIRECT,
INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF
THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR(S),
HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
THE AUTHOR(S) UNIVERSITY, COMPANY AND/OR AFFILIATION SPECIFICALLY
DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND
THE AUTHOR(S) HAVE NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT,
UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
Common Lisp Extension: Math
===========================
Copyright (c) 2018 Marco Antoniotti
See file COPYING for licensing information
DESCRIPTION
-----------
This library provides a layer of generic functions and some numerical
hooks that can be used to further "extend" Common Lisp mathematical
capabilities.
;;;; -*- Mode: Lisp -*-
;;;; common-math-package.lisp
;;;;
;;;; See file COPYING in the main folder for licensing information.
(defpackage "IT.UNIMIB.DISCO.MA.CL.EXT.MATH" (:use "CL")
(:nicknames
"COMMON-MATH"
"MATH"
"CL.MATH"
"CL.EXTENSIONS.MATH"
"CL.EXT.MATH")
;; IEEE names.
(:export
"NAN")
(:export
"LONG-FLOAT-POSITIVE-INFINITY"
"LONG-FLOAT-NEGATIVE-INFINITY"
"DOUBLE-FLOAT-POSITIVE-INFINITY"
"DOUBLE-FLOAT-NEGATIVE-INFINITY"
"SINGLE-FLOAT-POSITIVE-INFINITY"
"SINGLE-FLOAT-NEGATIVE-INFINITY"
"SHORT-FLOAT-POSITIVE-INFINITY"
"SHORT-FLOAT-NEGATIVE-INFINITY"
)
(:export "+POSITIVE-INFINITY+" "+NEGATIVE-INFINITY+")
(:export
"IS-NAN"
"NAN-P"
"IS-INFINITY"
"INFINITY-P"
)
(:export
"ROUNDING-MODES"
)
(:export
"UNIT-ROUNDOFF"
"+UNIT-ROUNDOFF-RATIONAL+"
"+UNIT-ROUNDOFF-SHORT-RATIONAL+"
"+UNIT-ROUNDOFF-SINGLE-RATIONAL+"
"+UNIT-ROUNDOFF-DOUBLE-RATIONAL+"
"+UNIT-ROUNDOFF-LONG-RATIONAL+"
"+UNIT-ROUNDOFF-FLOAT+"
"+UNIT-ROUNDOFF-SHORT-FLOAT+"
"+UNIT-ROUNDOFF-SINGLE-FLOAT+"
"+UNIT-ROUNDOFF-DOUBLE-FLOAT+"
"+UNIT-ROUNDOFF-LONG-FLOAT+"
)
;; Generic names (shadowing CL symbols).
(:shadow "=" "+" "-" "*" "/" ">" "<" ">=" "<=" "/=")
(:shadow "PLUSP" "MINUSP" "ZEROP")
(:shadow "EXPT" "GCD" "LCM")
(:export "=" "+" "-" "*" "/" ">" "<" ">=" "<=" "/=")
(:export "PLUSP" "MINUSP" "ZEROP")
(:export "EXPT" "GCD" "LCM")
(:export "INCREASING" "DECREASING")
(:export "MINIMUM" "MAXIMUM" "MINMAX")
(:export "OUTER-PRODUCT")
(:export "INNER-PRODUCT" "DOT")
(:export "SUM")
(:export "PRODUCT")
(:export ; Fixed arity operators.
"+." ".+."
"*." ".*."
"-." ".-."
"/." "./."
"=." ".=."
".<."
".>."
".<=."
".>=."
"/=."
"./=."
)
(:export "GCD." ".GCD." "LCM." ".LCM.")
(:export "/*" "./*.") ; Inner product.
(:export "*/" ".*/.") ; Outer product.
(:export "/+" "./+.") ; Inner sum.
(:export "+/" ".+/.") ; Outer sum.
;; Conditions.
(:export "UNDEFINED-OPERATION")
)
;;;; end of file -- common-math-package.lisp
;;;; -*- Mode: Lisp -*-
;;;; common-math.asd
;;;;
;;;; See file COPYING in the main folder for licensing information.
(asdf:defsystem "COMMON-MATH"
:author "Marco Antoniotti"
:license "BSD"
:components ((:file "common-math-package")
(:file "prologue"
:depends-on ("common-math-package"))
(:module "impl-dependent"
:depends-on ("prologue" "common-math-package")
:components (
#+lispworks
(:module "lispworks"
:components ((:file "lispworks")))
#+cmucl
(:module "cmucl"
:components ((:file "cmucl")))
#+sbcl
(:module "sbcl"
:components ((:file "sbcl")))
#+ccl
(:module "ccl"
:components ((:file "ccl")))
#+allegro
(:module "acl"
:components ((:file "acl")))
#-(or lispworks cmucl sbcl ccl allegro)
(:file "no-impl-dependent-code")
))
(:file "ieee-base"
:depends-on ("common-math-package" "impl-dependent"))
(:file "numerics-base"
:depends-on ("common-math-package"))
(:file "common-math"
:depends-on ("common-math-package"
"ieee-base"
"numerics-base"))
)
)
;;;; end of file -- common-math.asd
<
;;;; -*- 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))))
)))
)