Initial CLOS-based statistical data structure.
[CommonLispStat.git] / lsfloat.lsp
blob4db30b0ab312af232b29417217619d358d027ba0
1 ;;; -*- mode: lisp -*-
2 ;;; Copyright (c) 2005--2007, by A.J. Rossini <blindglobe@gmail.com>
3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
4 ;;; Since 1991, ANSI was finally finished. Edited for ANSI Common Lisp.
6 ;;; lsfloat -- Floating point specs and transcendental functions
7 ;;;
8 ;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
9 ;;; unrestricted use.
10 ;;;
11 ;;; Common Lisp allows for four different floating point types that need
12 ;;; not be distinct. For statistical work, the type I prefer to use is
13 ;;; the one that is closest to a C double. This type is named stat-float.
14 ;;; By setting the variable *read-default-float-format* to this type, you
15 ;;; insure that data entered as floating point data is read in with this
16 ;;; type. The problem arises with data read as integers that is passed to
17 ;;; a transcendental, like sqrt. Floating point contagion rules say these
18 ;;; integers are to be converted to type single-float. Unless single-float
19 ;;; is equivalent to C double, as it is in Mac CL and KCL, this is not
20 ;;; what I want. Hence this file redefines the transcendentals to first
21 ;;; coerce their arguments to stat-float before applying the built-in
22 ;;; functions.
23 ;;;
24 ;;; No actual modifications to the transcendentals are needed if
25 ;;; single-float is the same as stat-float. The feature
26 ;;; :stat-float-is-double-float is used to indicate this.
27 ;;;
28 ;;; KCL NOTE:
29 ;;; In (A)KCL the type short-float corresponds to C float and the types
30 ;;; single-float, double-float and long-float correspond to C double.
31 ;;; But in the implementation of the transcendentals (A)KCL coerces
32 ;;; rationals to short-float, not single-float. CLtL1 is a little vague
33 ;;; on this (it talks about "single precision") but CLtL2 clarifies that
34 ;;; rationals should produce single-float results. So (A)KCL is wrong, at
35 ;;; least relative to the clarification in CLtL2. I therefore decided
36 ;;; to fix (A)KCL in files c/num_sfun.c and lsp/numlib.lsp. If these
37 ;;; fixes are applied, the feature :stat-float-is-double-float should be
38 ;;; defined.
41 ;;; ======== From the CLHS ==========
43 ;; Type SHORT-FLOAT, SINGLE-FLOAT, DOUBLE-FLOAT, LONG-FLOAT
45 ;; Supertypes:
47 ;; short-float: short-float, float, real, number, t
50 ;; single-float: single-float, float, real, number, t
52 ;; double-float: double-float, float, real, number, t
54 ;; long-float: long-float, float, real, number, t
56 ;; Description: For the four defined subtypes of type float, it is true
57 ;; that intermediate between the type short-float and the type long-float
58 ;; are the type single-float and the type double-float. The precise
59 ;; definition of these categories is implementation-defined. The
60 ;; precision (measured in ``bits'', computed as p log 2b) and the
61 ;; exponent size (also measured in ``bits,'' computed as log 2(n+1),
62 ;; where n is the maximum exponent value) is recommended to be at least
63 ;; as great as the values in the next figure. Each of the defined
64 ;; subtypes of type float might or might not have a minus zero.
67 ;; Format Minimum Precision Minimum Exponent Size
68 ;; ----------
69 ;; Short 13 bits 5 bits
70 ;; Single 24 bits 8 bits
71 ;; Double 50 bits 8 bits
72 ;; Long 50 bits 8 bits
74 ;; Figure 12-12. Recommended Minimum Floating-Point Precision and Exponent Size
76 ;; There can be fewer than four internal representations for floats. If
77 ;; there are fewer distinct representations, the following rules apply:
79 ;; If there is only one, it is the type single-float. In this
80 ;; representation, an object is simultaneously of types single-float,
81 ;; double-float, short-float, and long-float.
83 ;; Two internal representations can be arranged in either of the
84 ;; following ways:
86 ;; Two types are provided: single-float and short-float. An object is
87 ;; simultaneously of types single-float, double-float, and long-float.
89 ;; Two types are provided: single-float and double-float. An object is
90 ;; simultaneously of types single-float and short-float, or double-float
91 ;; and long-float.
93 ;; Three internal representations can be arranged in either of the
94 ;; following ways:
96 ;; Three types are provided: short-float, single-float, and
97 ;; double-float. An object can simultaneously be of type double-float and
98 ;; long-float.
100 ;; Three types are provided: single-float, double-float, and
101 ;; long-float. An object can simultaneously be of types single-float and
102 ;; short-float.
104 ;;; Package Setup
106 (in-package :cl-user)
108 (defpackage :lisp-stat-float
109 (:use :common-lisp)
110 (:export +stat-float-typing+ +stat-cfloat-typing+ +stat-float-template+
111 machine-epsilon base-float makedouble
113 make-base-trans-fun-2 make-base-trans-fun
115 BASE-LOG BASE-EXP BASE-EXPT BASE-SQRT BASE-SIN BASE-COS
116 BASE-TAN BASE-ASIN BASE-ACOS BASE-ATAN BASE-SINH
117 BASE-COSH BASE-TANH BASE-ASINH BASE-ACOSH BASE-ATANH
118 BASE-ABS BASE-PHASE BASE-FFLOOR BASE-FCEILING BASE-FTRUNCATE
119 BASE-FROUND BASE-SIGNUM BASE-CIS))
121 (in-package #:lisp-stat-float)
123 ;; This should technically be conditionalized to the Lisp
124 ;; implementation, i.e.
125 #+sbcl(pushnew :stat-float-is-double-float *features*)
126 #+cmu(pushnew :stat-float-is-double-float *features*)
127 #+clisp(pushnew :stat-float-is-double-float *features*)
128 ;; etc... the above need to be verified!
130 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
131 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
132 ;;;;
133 ;;;; Constants determining default floating point type for statistical
134 ;;;; operations. This type generally corresponds to a C double.
135 ;;;;
136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
137 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
139 ;; Why do we do this instead of defparameter? mostly to allow for doc
140 ;; defs as well as to minimize dereferencing mistakes.
141 (defmacro define-constant (name value &optional doc)
142 `(defparameter ,name (if (boundp ',name) (symbol-value ',name) ,value)
143 ,@(when doc (list doc))))
145 (define-constant +stat-float-typing+ 'long-float)
146 (define-constant +stat-cfloat-typing+ '(complex long-float))
147 (define-constant +stat-float-template+ 0.d0)
149 ;(defparameter +stat-float-typing+ 'long-float)
150 ;(defparameter +stat-cfloat-typing+ '(complex long-float))
151 ;(defparameter +stat-float-template+ 0.d0)
153 (deftype stat-float () +stat-float-typing+)
154 (deftype stat-cfloat () +stat-cfloat-typing+)
156 (defparameter machine-epsilon
157 (do ((epsilon (float 1.0 +stat-float-template+)
158 (/ epsilon 2.0)))
159 ((= (+ 1.0 (/ epsilon 2.0)) 1.0) epsilon)))
161 (defmacro declare-double (&rest vars) `(declare (long-float ,@vars)))
163 (setf *read-default-float-format* +stat-float-typing+)
165 ;;; FIX-BASE-DOC adds note about modification to the end of the
166 ;;; documentation string argument.
167 (defmacro fix-base-doc (doc)
168 `(format nil
169 "~a~%Modified to coerce arguments(s) to stat-float or stat-cfloat."
170 ,doc))
172 ;;; FIX-BASE-FUN-DOC fixes documentation of SYM and installs in BASE-SYM
173 (defun fix-base-fun-doc (sym base-sym)
174 (let ((doc (documentation sym 'function)))
175 (if doc (setf (documentation base-sym 'function) (fix-base-doc doc)))))
177 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
178 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
179 ;;;;
180 ;;;; Functions and Macros for modifying functions to coerce to standard
181 ;;;; floating point type.
182 ;;;;
183 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
184 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
186 ;;; FFIX - coerces its arguments to standard real or complex floating
187 ;;; point number if needed.
188 #-stat-float-is-double-float
189 (defmacro ffix (x)
190 `(if (complexp ,x)
191 (coerce ,x +stat-cfloat-typing+)
192 (float ,x +stat-float-template+)))
194 #+stat-float-is-double-float
195 (defmacro ffix (x) x)
197 ;;; MAKEDOUBLE coerces its argument to the standard floating point
198 ;;; type specified by +stat-float-template+
200 (defun makedouble (x) (float x +stat-float-template+))
202 ;; Why doesn't this optimization work?
203 ;; #+stat-float-is-double-float
204 ;; (eval-when (compile)
205 ;; (proclaim '(function makedouble (t) long-float)))
208 ;;; MAKE-BASE-TRANS-FUN Macro for re-defining one argument transcendental
209 ;;; functions
211 #-:stat-float-is-double-float
212 (defmacro make-base-trans-fun (sym)
213 (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
214 (doc (documentation sym 'function))
215 (doc-list (if doc (list (fix-base-doc doc)))))
216 `(defun ,base-sym (x)
217 ,@doc-list
218 (declare (inline ,sym coerce float))
219 (,sym (ffix x)))))
221 #+:stat-float-is-double-float
222 (defmacro make-base-trans-fun (sym)
223 (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
224 (doc (documentation sym 'function)))
225 `(progn (setf (symbol-function ',base-sym) (symbol-function ',sym))
226 (if ,doc (setf (documentation ',base-sym 'function) ,doc)))))
229 ;;; MAKE-BASE-TRANS-FUN-2 Macro for re-defining transcendental functions
230 ;;; with an optional second argument
232 #-:stat-float-is-double-float
233 (defmacro make-base-trans-fun-2 (sym)
234 (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
235 (doc (documentation sym 'function))
236 (doc-list (if doc (list (fix-base-doc doc)))))
237 `(defun ,base-sym (x &optional y)
238 ,@doc-list
239 (declare (inline ,sym coerce float))
240 (if y (,sym (ffix x) (ffix y)) (,sym (ffix x))))))
242 #+:stat-float-is-double-float
243 (defmacro make-base-trans-fun-2 (sym)
244 (let* ((base-sym (intern (string-upcase (format nil "BASE-~s" sym))))
245 (doc (documentation sym 'function)))
246 `(progn (setf (symbol-function ',base-sym) (symbol-function ',sym))
247 (if ,doc (setf (documentation ',base-sym 'function) ,doc)))))
249 ;;; Modified base functions to coerce to standard floating point type
253 ;;; BASE-FLOAT
254 #-:stat-float-is-double-float
255 (progn
256 (defun base-float (x &optional (y +stat-float-template+)) (float x y))
257 (fix-base-fun-doc 'float 'base-float))
259 #+:stat-float-is-double-float
260 (make-base-trans-fun float)
262 ;;; BASE-EXPT
263 #-:stat-float-is-double-float
264 (progn
265 (defun base-expt (x y)
266 (declare (inline expt coerce float integerp))
267 (if (integerp y) (expt x y) (expt (ffix x) (ffix y))))
269 (fix-base-fun-doc 'expt 'base-expt))
271 #+:stat-float-is-double-float
272 (make-base-trans-fun expt)
274 ;;; Others
275 (make-base-trans-fun-2 log)
276 (make-base-trans-fun exp)
277 (make-base-trans-fun sqrt)
278 (make-base-trans-fun sin)
279 (make-base-trans-fun cos)
280 (make-base-trans-fun tan)
281 (make-base-trans-fun asin)
282 (make-base-trans-fun acos)
283 (make-base-trans-fun-2 atan)
284 (make-base-trans-fun sinh)
285 (make-base-trans-fun cosh)
286 (make-base-trans-fun tanh)
287 (make-base-trans-fun asinh)
288 (make-base-trans-fun acosh)
289 (make-base-trans-fun atanh)
290 (make-base-trans-fun abs)
291 (make-base-trans-fun phase)
292 (make-base-trans-fun-2 ffloor)
293 (make-base-trans-fun-2 fceiling)
294 (make-base-trans-fun-2 ftruncate)
295 (make-base-trans-fun-2 fround)
296 (make-base-trans-fun signum)
297 (make-base-trans-fun cis)