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 (in-package #:lisp-stat-float
)
13 ;;; Common Lisp allows for four different floating point types that need
14 ;;; not be distinct. For statistical work, the type I prefer to use is
15 ;;; the one that is closest to a C double. This type is named stat-float.
16 ;;; By setting the variable *read-default-float-format* to this type, you
17 ;;; insure that data entered as floating point data is read in with this
18 ;;; type. The problem arises with data read as integers that is passed to
19 ;;; a transcendental, like sqrt. Floating point contagion rules say these
20 ;;; integers are to be converted to type single-float. Unless single-float
21 ;;; is equivalent to C double, as it is in Mac CL and KCL, this is not
22 ;;; what I want. Hence this file redefines the transcendentals to first
23 ;;; coerce their arguments to stat-float before applying the built-in
26 ;;; No actual modifications to the transcendentals are needed if
27 ;;; single-float is the same as stat-float. The feature
28 ;;; :stat-float-is-double-float is used to indicate this.
31 ;;; In (A)KCL the type short-float corresponds to C float and the types
32 ;;; single-float, double-float and long-float correspond to C double.
33 ;;; But in the implementation of the transcendentals (A)KCL coerces
34 ;;; rationals to short-float, not single-float. CLtL1 is a little vague
35 ;;; on this (it talks about "single precision") but CLtL2 clarifies that
36 ;;; rationals should produce single-float results. So (A)KCL is wrong, at
37 ;;; least relative to the clarification in CLtL2. I therefore decided
38 ;;; to fix (A)KCL in files c/num_sfun.c and lsp/numlib.lsp. If these
39 ;;; fixes are applied, the feature :stat-float-is-double-float should be
43 ;;; ======== From the CLHS ==========
45 ;; Type SHORT-FLOAT, SINGLE-FLOAT, DOUBLE-FLOAT, LONG-FLOAT
49 ;; short-float: short-float, float, real, number, t
52 ;; single-float: single-float, float, real, number, t
54 ;; double-float: double-float, float, real, number, t
56 ;; long-float: long-float, float, real, number, t
58 ;; Description: For the four defined subtypes of type float, it is true
59 ;; that intermediate between the type short-float and the type long-float
60 ;; are the type single-float and the type double-float. The precise
61 ;; definition of these categories is implementation-defined. The
62 ;; precision (measured in ``bits'', computed as p log 2b) and the
63 ;; exponent size (also measured in ``bits,'' computed as log 2(n+1),
64 ;; where n is the maximum exponent value) is recommended to be at least
65 ;; as great as the values in the next figure. Each of the defined
66 ;; subtypes of type float might or might not have a minus zero.
69 ;; Format Minimum Precision Minimum Exponent Size
71 ;; Short 13 bits 5 bits
72 ;; Single 24 bits 8 bits
73 ;; Double 50 bits 8 bits
74 ;; Long 50 bits 8 bits
76 ;; Figure 12-12. Recommended Minimum Floating-Point Precision and Exponent Size
78 ;; There can be fewer than four internal representations for floats. If
79 ;; there are fewer distinct representations, the following rules apply:
81 ;; If there is only one, it is the type single-float. In this
82 ;; representation, an object is simultaneously of types single-float,
83 ;; double-float, short-float, and long-float.
85 ;; Two internal representations can be arranged in either of the
88 ;; Two types are provided: single-float and short-float. An object is
89 ;; simultaneously of types single-float, double-float, and long-float.
91 ;; Two types are provided: single-float and double-float. An object is
92 ;; simultaneously of types single-float and short-float, or double-float
95 ;; Three internal representations can be arranged in either of the
98 ;; Three types are provided: short-float, single-float, and
99 ;; double-float. An object can simultaneously be of type double-float and
102 ;; Three types are provided: single-float, double-float, and
103 ;; long-float. An object can simultaneously be of types single-float and
106 ;; This should technically be conditionalized to the Lisp
107 ;; implementation, i.e.
108 #+sbcl
(pushnew :stat-float-is-double-float
*features
*)
109 #+cmu
(pushnew :stat-float-is-double-float
*features
*)
110 #+clisp
(pushnew :stat-float-is-double-float
*features
*)
111 ;; etc... the above need to be verified!
113 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
114 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
116 ;;;; Constants determining default floating point type for statistical
117 ;;;; operations. This type generally corresponds to a C double.
119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
120 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
122 ;; Why do we do this instead of defparameter? mostly to allow for doc
123 ;; defs as well as to minimize dereferencing mistakes.
124 (defmacro define-constant
(name value
&optional doc
)
125 `(defparameter ,name
(if (boundp ',name
) (symbol-value ',name
) ,value
)
126 ,@(when doc
(list doc
))))
128 (define-constant +stat-float-typing
+ 'long-float
)
129 (define-constant +stat-cfloat-typing
+ '(complex long-float
))
130 (define-constant +stat-float-template
+ 0.d0
)
132 ;(defparameter +stat-float-typing+ 'long-float)
133 ;(defparameter +stat-cfloat-typing+ '(complex long-float))
134 ;(defparameter +stat-float-template+ 0.d0)
136 (deftype stat-float
() +stat-float-typing
+)
137 (deftype stat-cfloat
() +stat-cfloat-typing
+)
139 (defparameter machine-epsilon
140 (do ((epsilon (float 1.0 +stat-float-template
+)
142 ((= (+ 1.0 (/ epsilon
2.0)) 1.0) epsilon
)))
144 (defmacro declare-double
(&rest vars
) `(declare (long-float ,@vars
)))
146 (setf *read-default-float-format
* +stat-float-typing
+)
148 ;;; FIX-BASE-DOC adds note about modification to the end of the
149 ;;; documentation string argument.
150 (defmacro fix-base-doc
(doc)
152 "~a~%Modified to coerce arguments(s) to stat-float or stat-cfloat."
155 ;;; FIX-BASE-FUN-DOC fixes documentation of SYM and installs in BASE-SYM
156 (defun fix-base-fun-doc (sym base-sym
)
157 (let ((doc (documentation sym
'function
)))
158 (if doc
(setf (documentation base-sym
'function
) (fix-base-doc doc
)))))
160 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
161 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
163 ;;;; Functions and Macros for modifying functions to coerce to standard
164 ;;;; floating point type.
166 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
167 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
169 ;;; FFIX - coerces its arguments to standard real or complex floating
170 ;;; point number if needed.
171 #-stat-float-is-double-float
174 (coerce ,x
+stat-cfloat-typing
+)
175 (float ,x
+stat-float-template
+)))
177 #+stat-float-is-double-float
178 (defmacro ffix
(x) x
)
180 ;;; MAKEDOUBLE coerces its argument to the standard floating point
181 ;;; type specified by +stat-float-template+
183 (defun makedouble (x) (float x
+stat-float-template
+))
185 ;; Why doesn't this optimization work?
186 ;; #+stat-float-is-double-float
187 ;; (eval-when (compile)
188 ;; (proclaim '(function makedouble (t) long-float)))
191 ;;; MAKE-BASE-TRANS-FUN Macro for re-defining one argument transcendental
194 #-
:stat-float-is-double-float
195 (defmacro make-base-trans-fun
(sym)
196 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
197 (doc (documentation sym
'function
))
198 (doc-list (if doc
(list (fix-base-doc doc
)))))
199 `(defun ,base-sym
(x)
201 (declare (inline ,sym coerce float
))
204 #+:stat-float-is-double-float
205 (defmacro make-base-trans-fun
(sym)
206 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
207 (doc (documentation sym
'function
)))
208 `(progn (setf (symbol-function ',base-sym
) (symbol-function ',sym
))
209 (if ,doc
(setf (documentation ',base-sym
'function
) ,doc
)))))
212 ;;; MAKE-BASE-TRANS-FUN-2 Macro for re-defining transcendental functions
213 ;;; with an optional second argument
215 #-
:stat-float-is-double-float
216 (defmacro make-base-trans-fun-2
(sym)
217 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
218 (doc (documentation sym
'function
))
219 (doc-list (if doc
(list (fix-base-doc doc
)))))
220 `(defun ,base-sym
(x &optional y
)
222 (declare (inline ,sym coerce float
))
223 (if y
(,sym
(ffix x
) (ffix y
)) (,sym
(ffix x
))))))
225 #+:stat-float-is-double-float
226 (defmacro make-base-trans-fun-2
(sym)
227 (let* ((base-sym (intern (string-upcase (format nil
"BASE-~s" sym
))))
228 (doc (documentation sym
'function
)))
229 `(progn (setf (symbol-function ',base-sym
) (symbol-function ',sym
))
230 (if ,doc
(setf (documentation ',base-sym
'function
) ,doc
)))))
232 ;;; Modified base functions to coerce to standard floating point type
237 #-
:stat-float-is-double-float
239 (defun base-float (x &optional
(y +stat-float-template
+)) (float x y
))
240 (fix-base-fun-doc 'float
'base-float
))
242 #+:stat-float-is-double-float
243 (make-base-trans-fun float
)
246 #-
:stat-float-is-double-float
248 (defun base-expt (x y
)
249 (declare (inline expt coerce float integerp
))
250 (if (integerp y
) (expt x y
) (expt (ffix x
) (ffix y
))))
252 (fix-base-fun-doc 'expt
'base-expt
))
254 #+:stat-float-is-double-float
255 (make-base-trans-fun expt
)
258 (make-base-trans-fun-2 log
)
259 (make-base-trans-fun exp
)
260 (make-base-trans-fun sqrt
)
261 (make-base-trans-fun sin
)
262 (make-base-trans-fun cos
)
263 (make-base-trans-fun tan
)
264 (make-base-trans-fun asin
)
265 (make-base-trans-fun acos
)
266 (make-base-trans-fun-2 atan
)
267 (make-base-trans-fun sinh
)
268 (make-base-trans-fun cosh
)
269 (make-base-trans-fun tanh
)
270 (make-base-trans-fun asinh
)
271 (make-base-trans-fun acosh
)
272 (make-base-trans-fun atanh
)
273 (make-base-trans-fun abs
)
274 (make-base-trans-fun phase
)
275 (make-base-trans-fun-2 ffloor
)
276 (make-base-trans-fun-2 fceiling
)
277 (make-base-trans-fun-2 ftruncate
)
278 (make-base-trans-fun-2 fround
)
279 (make-base-trans-fun signum
)
280 (make-base-trans-fun cis
)