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
8 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
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
24 ;;;; No actual modifications to the transcendentals are needed if
25 ;;;; single-float is the same as stat-float. The fearure
26 ;;;; :stat-float-is-double-float is used to indicate this.
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
41 ;;; ======== From the CLHS ==========
43 ;; Type SHORT-FLOAT, SINGLE-FLOAT, DOUBLE-FLOAT, LONG-FLOAT
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
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
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
93 ;; Three internal representations can be arranged in either of the
96 ;; Three types are provided: short-float, single-float, and
97 ;; double-float. An object can simultaneously be of type double-float and
100 ;; Three types are provided: single-float, double-float, and
101 ;; long-float. An object can simultaneously be of types single-float and
105 (in-package :cl-user
)
107 (defpackage :lisp-stat-float
109 (:export
+stat-float-typing
+ +stat-cfloat-typing
+ +stat-float-template
+
110 machine-epsilon base-float makedouble
112 make-base-trans-fun-2 make-base-trans-fun
114 BASE-LOG BASE-EXP BASE-EXPT BASE-SQRT BASE-SIN BASE-COS
115 BASE-TAN BASE-ASIN BASE-ACOS BASE-ATAN BASE-SINH
116 BASE-COSH BASE-TANH BASE-ASINH BASE-ACOSH BASE-ATANH
117 BASE-ABS BASE-PHASE BASE-FFLOOR BASE-FCEILING BASE-FTRUNCATE
118 BASE-FROUND BASE-SIGNUM BASE-CIS
))
120 (in-package #:lisp-stat-float
)
122 ;; This should technically be conditionalized to the Lisp
123 ;; implementation, i.e.
124 #+sbcl
(pushnew :stat-float-is-double-float
*features
*)
125 #+cmu
(pushnew :stat-float-is-double-float
*features
*)
126 #+clisp
(pushnew :stat-float-is-double-float
*features
*)
127 ;; etc... the above need to be verified!
129 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
130 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
132 ;;;; Constants determining default floating point type for statistical
133 ;;;; operations. This type generally corresponds to a C double.
135 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
138 ;; Why do we do this instead of defparameter? mostly to allow for doc
139 ;; defs as well as to minimize dereferencing mistakes.
140 (defmacro define-constant
(name value
&optional doc
)
141 `(defparameter ,name
(if (boundp ',name
) (symbol-value ',name
) ,value
)
142 ,@(when doc
(list doc
))))
144 (define-constant +stat-float-typing
+ 'long-float
)
145 (define-constant +stat-cfloat-typing
+ '(complex long-float
))
146 (define-constant +stat-float-template
+ 0.d0
)
148 ;(defparameter +stat-float-typing+ 'long-float)
149 ;(defparameter +stat-cfloat-typing+ '(complex long-float))
150 ;(defparameter +stat-float-template+ 0.d0)
152 (deftype stat-float
() +stat-float-typing
+)
153 (deftype stat-cfloat
() +stat-cfloat-typing
+)
155 (defparameter machine-epsilon
156 (do ((epsilon (float 1.0 +stat-float-template
+)
158 ((= (+ 1.0 (/ epsilon
2.0)) 1.0) epsilon
)))
160 (defmacro declare-double
(&rest vars
) `(declare (long-float ,@vars
)))
162 (setf *read-default-float-format
* +stat-float-typing
+)
164 ;;; FIX-BASE-DOC adds note about modification to the end of the
165 ;;; documentation string argument.
166 (defmacro fix-base-doc
(doc)
168 "~a~%Modified to coerce arguments(s) to stat-float or stat-cfloat."
171 ;;; FIX-BASE-FUN-DOC fixes documentation of SYM and installs in BASE-SYM
172 (defun fix-base-fun-doc (sym base-sym
)
173 (let ((doc (documentation sym
'function
)))
174 (if doc
(setf (documentation base-sym
'function
) (fix-base-doc doc
)))))
176 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
177 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
179 ;;;; Functions and Macros for modifying functions to coerce to standard
180 ;;;; floating point type.
182 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
183 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
185 ;;; FFIX - coerces its arguments to standard real or complex floating
186 ;;; point number if needed.
187 #-stat-float-is-double-float
190 (coerce ,x
+stat-cfloat-typing
+)
191 (float ,x
+stat-float-template
+)))
193 #+stat-float-is-double-float
194 (defmacro ffix
(x) x
)
196 ;;; MAKEDOUBLE coerces its argument to the standard floating point
197 ;;; type specified by +stat-float-template+
199 (defun makedouble (x) (float x
+stat-float-template
+))
201 ;; Why doesn't this optimization work?
202 ;; #+stat-float-is-double-float
203 ;; (eval-when (compile)
204 ;; (proclaim '(function makedouble (t) long-float)))
207 ;;; MAKE-BASE-TRANS-FUN Macro for re-defining one argument transcendental
210 #-
:stat-float-is-double-float
211 (defmacro make-base-trans-fun
(sym)
212 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
213 (doc (documentation sym
'function
))
214 (doc-list (if doc
(list (fix-base-doc doc
)))))
215 `(defun ,base-sym
(x)
217 (declare (inline ,sym coerce float
))
220 #+:stat-float-is-double-float
221 (defmacro make-base-trans-fun
(sym)
222 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
223 (doc (documentation sym
'function
)))
224 `(progn (setf (symbol-function ',base-sym
) (symbol-function ',sym
))
225 (if ,doc
(setf (documentation ',base-sym
'function
) ,doc
)))))
228 ;;; MAKE-BASE-TRANS-FUN-2 Macro for re-defining transcendental functions
229 ;;; with an optional second argument
231 #-
:stat-float-is-double-float
232 (defmacro make-base-trans-fun-2
(sym)
233 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
234 (doc (documentation sym
'function
))
235 (doc-list (if doc
(list (fix-base-doc doc
)))))
236 `(defun ,base-sym
(x &optional y
)
238 (declare (inline ,sym coerce float
))
239 (if y
(,sym
(ffix x
) (ffix y
)) (,sym
(ffix x
))))))
241 #+:stat-float-is-double-float
242 (defmacro make-base-trans-fun-2
(sym)
243 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
244 (doc (documentation sym
'function
)))
245 `(progn (setf (symbol-function ',base-sym
) (symbol-function ',sym
))
246 (if ,doc
(setf (documentation ',base-sym
'function
) ,doc
)))))
248 ;;; Modified base functions to coerce to standard floating point type
253 #-
:stat-float-is-double-float
255 (defun base-float (x &optional
(y +stat-float-template
+)) (float x y
))
256 (fix-base-fun-doc 'float
'base-float
))
258 #+:stat-float-is-double-float
259 (make-base-trans-fun float
)
262 #-
:stat-float-is-double-float
264 (defun base-expt (x y
)
265 (declare (inline expt coerce float integerp
))
266 (if (integerp y
) (expt x y
) (expt (ffix x
) (ffix y
))))
268 (fix-base-fun-doc 'expt
'base-expt
))
270 #+:stat-float-is-double-float
271 (make-base-trans-fun expt
)
274 (make-base-trans-fun-2 log
)
275 (make-base-trans-fun exp
)
276 (make-base-trans-fun sqrt
)
277 (make-base-trans-fun sin
)
278 (make-base-trans-fun cos
)
279 (make-base-trans-fun tan
)
280 (make-base-trans-fun asin
)
281 (make-base-trans-fun acos
)
282 (make-base-trans-fun-2 atan
)
283 (make-base-trans-fun sinh
)
284 (make-base-trans-fun cosh
)
285 (make-base-trans-fun tanh
)
286 (make-base-trans-fun asinh
)
287 (make-base-trans-fun acosh
)
288 (make-base-trans-fun atanh
)
289 (make-base-trans-fun abs
)
290 (make-base-trans-fun phase
)
291 (make-base-trans-fun-2 ffloor
)
292 (make-base-trans-fun-2 fceiling
)
293 (make-base-trans-fun-2 ftruncate
)
294 (make-base-trans-fun-2 fround
)
295 (make-base-trans-fun signum
)
296 (make-base-trans-fun cis
)