Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

[e90b2f]: src / lsp / numlib.lsp Maximize Restore History

Download this file

numlib.lsp    333 lines (296 with data), 11.7 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
;;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: SYSTEM -*-
;;;;
;;;; Copyright (c) 1984, Taiichi Yuasa and Masami Hagiya.
;;;; Copyright (c) 1990, Giuseppe Attardi.
;;;; Copyright (c) 2001, Juan Jose Garcia Ripoll.
;;;;
;;;; This program is free software; you can redistribute it and/or
;;;; modify it under the terms of the GNU Library General Public
;;;; License as published by the Free Software Foundation; either
;;;; version 2 of the License, or (at your option) any later version.
;;;;
;;;; See file '../Copyright' for full details.
;;;; number routines
(in-package "SYSTEM")
#-ecl-min
(ffi:clines "#include <math.h>")
#.
(flet ((binary-search (f min max)
(do ((new (/ (+ min max) 2) (/ (+ min max) 2)))
((>= min max)
max)
(if (funcall f new)
(if (= new max)
(return max)
(setq max new))
(if (= new min)
(return max)
(setq min new)))))
(epsilon+ (x)
(/= (float 1 x) (+ (float 1 x) x)))
(epsilon- (x)
(/= (float 1 x) (- (float 1 x) x))))
#+ecl-min
(si::trap-fpe 'last nil)
`(eval-when (compile load eval)
(defconstant short-float-epsilon
,(binary-search #'epsilon+ (coerce 0 'short-float) (coerce 1 'short-float))
"The smallest postive short-float E that satisfies
(not (= (float 1 E) (+ (float 1 E) E)))")
(defconstant single-float-epsilon
,(binary-search #'epsilon+ (coerce 0 'single-float) (coerce 1 'single-float))
"The smallest postive single-float E that satisfies
(not (= (float 1 E) (+ (float 1 E) E)))")
(defconstant double-float-epsilon
,(binary-search #'epsilon+ (coerce 0 'double-float) (coerce 1 'double-float))
"The smallest postive double-float E that satisfies
(not (= (float 1 E) (+ (float 1 E) E)))")
(defconstant long-float-epsilon
,(binary-search #'epsilon+ (coerce 0 'long-float) (coerce 1 'long-float))
"The smallest postive long-float E that satisfies
(not (= (float 1 E) (+ (float 1 E) E)))")
(defconstant short-float-negative-epsilon
,(binary-search #'epsilon- (coerce 0 'short-float) (coerce 1 'short-float))
"The smallest positive short-float E that satisfies
(not (= (float 1 E) (- (float 1 E) E)))")
(defconstant single-float-negative-epsilon
,(binary-search #'epsilon- (coerce 0 'single-float) (coerce 1 'single-float))
"The smallest positive single-float E that satisfies
(not (= (float 1 E) (- (float 1 E) E)))")
(defconstant double-float-negative-epsilon
,(binary-search #'epsilon- (coerce 0 'double-float) (coerce 1 'double-float))
"The smallest positive double-float E that satisfies
(not (= (float 1 E) (- (float 1 E) E)))")
(defconstant long-float-negative-epsilon
,(binary-search #'epsilon- (coerce 0 'long-float) (coerce 1 'long-float))
"The smallest positive long-float E that satisfies
(not (= (float 1 E) (- (float 1 E) E)))")
))
#+IEEE-FLOATING-POINT
(locally
(declare (notinline -))
(let* ((bits (si::trap-fpe 'last nil)))
(let ((a (/ (coerce 1 'short-float) (coerce 0.0 'short-float))))
(defconstant short-float-positive-infinity a)
(defconstant short-float-negative-infinity (- a)))
(let ((a (/ (coerce 1 'single-float) (coerce 0.0 'single-float))))
(defconstant single-float-positive-infinity a)
(defconstant single-float-negative-infinity (- a)))
(let ((a (/ (coerce 1 'double-float) (coerce 0.0 'double-float))))
(defconstant double-float-positive-infinity a)
(defconstant double-float-negative-infinity (- a)))
(let ((a (/ (coerce 1 'long-float) (coerce 0.0 'long-float))))
(defconstant long-float-positive-infinity a)
(defconstant long-float-negative-infinity (- a)))
(si::trap-fpe bits t)))
(defconstant imag-one #C(0.0 1.0))
(defun isqrt (i)
"Args: (integer)
Returns the integer square root of INTEGER."
(unless (and (integerp i) (>= i 0))
(error 'type-error :datum i :expected-type 'unsigned-byte))
(if (zerop i)
0
(let ((n (integer-length i)))
(do ((x (ash 1 (ceiling n 2)))
(y))
(nil)
(setq y (floor i x))
(when (<= x y)
(return x))
(setq x (floor (+ x y) 2))))))
(defun phase (x)
"Args: (number)
Returns the angle part (in radians) of the polar representation of NUMBER.
Returns zero for non-complex numbers."
(if (zerop x)
(if (eq x 0) 0.0 (float 0 (realpart x)))
(atan (imagpart x) (realpart x))))
(defun signum (x)
"Args: (number)
Returns a number that represents the sign of NUMBER. Returns NUMBER If it is
zero. Otherwise, returns the value of (/ NUMBER (ABS NUMBER))"
(if (zerop x) x (/ x (abs x))))
(defun cis (x)
"Args: (radians)
Returns a complex number whose realpart and imagpart are the values of (COS
RADIANS) and (SIN RADIANS) respectively."
(exp (* imag-one x)))
#-ecl-min
(eval-when (:compile-toplevel)
(defmacro c-num-op (name arg)
#+long-float
`(ffi::c-inline (,arg) (:long-double) :long-double
,(format nil "~al(#0)" name)
:one-liner t)
#-long-float
`(ffi::c-inline (,arg) (:double) :double
,(format nil "~a(#0)" name)
:one-liner t)))
(defun asin (x)
"Args: (number)
Returns the arc sine of NUMBER."
(if #+ecl-min t #-ecl-min (complexp x)
(complex-asin x)
#-ecl-min
(let* ((x (float x))
(xr (float x 1l0)))
(declare (long-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (c-num-op "asin" xr) x)
(complex-asin x)))))
;; Ported from CMUCL
(defun complex-asin (z)
(declare (number z)
(si::c-local))
(let ((sqrt-1-z (sqrt (- 1 z)))
(sqrt-1+z (sqrt (+ 1 z))))
(complex (atan (realpart z) (realpart (* sqrt-1-z sqrt-1+z)))
(asinh (imagpart (* (conjugate sqrt-1-z)
sqrt-1+z))))))
(defun acos (x)
"Args: (number)
Returns the arc cosine of NUMBER."
(if #+ecl-min t #-ecl-min (complexp x)
(complex-acos x)
#-ecl-min
(let* ((x (float x))
(xr (float x 1l0)))
(declare (long-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (c-num-op "acos" xr) (float x))
(complex-acos x)))))
;; Ported from CMUCL
(defun complex-acos (z)
(declare (number z)
(si::c-local))
(let ((sqrt-1+z (sqrt (+ 1 z)))
(sqrt-1-z (sqrt (- 1 z))))
(complex (* 2 (atan (realpart sqrt-1-z) (realpart sqrt-1+z)))
(asinh (imagpart (* (conjugate sqrt-1+z)
sqrt-1-z))))))
#+(and (not ecl-min) win32 (not mingw32))
(progn
(ffi:clines "double asinh(double x) { return log(x+sqrt(1.0+x*x)); }")
(ffi:clines "double acosh(double x) { return log(x+sqrt((x-1)*(x+1))); }")
(ffi:clines "double atanh(double x) { return log((1+x)/(1-x))/2; }"))
#+(and long-float (not ecl-min) win32 (not mingw32))
(progn
(ffi:clines "double asinhl(long double x) { return logl(x+sqrtl(1.0+x*x)); }")
(ffi:clines "double acoshl(long double x) { return logl(x+sqrtl((x-1)*(x+1))); }")
(ffi:clines "double atanhl(long double x) { return logl((1+x)/(1-x))/2; }"))
;; Ported from CMUCL
(defun asinh (x)
"Args: (number)
Returns the hyperbolic arc sine of NUMBER."
;(log (+ x (sqrt (+ 1.0 (* x x)))))
(if #+(or ecl-min) t #-(or ecl-min) (complexp x)
(let* ((iz (complex (- (imagpart x)) (realpart x)))
(result (complex-asin iz)))
(complex (imagpart result)
(- (realpart result))))
#-(or ecl-min)
(float (c-num-op "asinh" x) (float x))))
;; Ported from CMUCL
(defun acosh (x)
"Args: (number)
Returns the hyperbolic arc cosine of NUMBER."
;(log (+ x (sqrt (* (1- x) (1+ x)))))
(if #+(or ecl-min) t #-(or ecl-min) (complexp x)
(complex-acosh x)
#-(or ecl-min)
(let* ((x (float x))
(xr (float x 1d0)))
(declare (double-float xr))
(if (<= 1.0 xr)
(float (c-num-op "acosh" xr) (float x))
(complex-acosh x)))))
(defun complex-acosh (z)
(declare (number z) (si::c-local))
(let ((sqrt-z-1 (sqrt (- z 1)))
(sqrt-z+1 (sqrt (+ z 1))))
(complex (asinh (realpart (* (conjugate sqrt-z-1)
sqrt-z+1)))
(* 2 (atan (imagpart sqrt-z-1) (realpart sqrt-z+1))))))
(defun atanh (x)
"Args: (number)
Returns the hyperbolic arc tangent of NUMBER."
;(/ (- (log (1+ x)) (log (- 1 x))) 2)
(if #+(or ecl-min) t #-(or ecl-min) (complexp x)
(complex-atanh x)
#-(or ecl-min)
(let* ((x (float x))
(xr (float x 1d0)))
(declare (double-float xr))
(if (and (<= -1.0 xr) (<= xr 1.0))
(float (c-num-op "atanh" xr) (float x))
(complex-atanh x)))))
(defun complex-atanh (z)
(declare (number z) (si::c-local))
(/ (- (log (1+ z)) (log (- 1 z))) 2))
(defun ffloor (x &optional (y 1))
"Args: (number &optional (divisor 1))
Same as FLOOR, but returns a float as the first value."
(multiple-value-bind (i r) (floor x y)
(values (if (floatp r) (float i r) (float i)) r)))
(defun fceiling (x &optional (y 1.0f0))
"Args: (number &optional (divisor 1))
Same as CEILING, but returns a float as the first value."
(multiple-value-bind (i r) (ceiling x y)
(values (if (floatp r) (float i r) (float i)) r)))
(defun ftruncate (x &optional (y 1.0f0))
"Args: (number &optional (divisor 1))
Same as TRUNCATE, but returns a float as the first value."
(multiple-value-bind (i r) (truncate x y)
(values (if (floatp r) (float i r) (float i)) r)))
(defun fround (x &optional (y 1.0f0))
"Args: (number &optional (divisor 1))
Same as ROUND, but returns a float as the first value."
(multiple-value-bind (i r) (round x y)
(values (if (floatp r) (float i r) (float i)) r)))
(defun logtest (x y)
"Args: (integer1 integer2)
Equivalent to (NOT (ZEROP (LOGAND INTEGER1 INTEGER2)))."
(not (zerop (logand x y))))
(defun byte (size position)
"Args: (size position)
Returns a byte specifier of integers. The value specifies the SIZE-bits byte
starting the least-significant-bit but POSITION bits of integers. In ECL, a
byte specifier is represented by a dotted pair (SIZE . POSITION)."
(cons size position))
(defun byte-size (bytespec)
"Args: (byte)
Returns the size part (in ECL, the car part) of the byte specifier BYTE."
(car bytespec))
(defun byte-position (bytespec)
"Args: (byte)
Returns the position part (in ECL, the cdr part) of the byte specifier BYTE."
(cdr bytespec))
(defun ldb (bytespec integer)
"Args: (bytespec integer)
Extracts a byte from INTEGER at the specified byte position, right-justifies
the byte, and returns the result as an integer."
(logandc2 (ash integer (- (byte-position bytespec)))
(- (ash 1 (byte-size bytespec)))))
(defun ldb-test (bytespec integer)
"Args: (bytespec integer)
Returns T if at least one bit of the specified byte is 1; NIL otherwise."
(not (zerop (ldb bytespec integer))))
(defun mask-field (bytespec integer)
"Args: (bytespec integer)
Extracts the specified byte from INTEGER and returns the result as an integer."
(ash (ldb bytespec integer) (byte-position bytespec)))
(defun dpb (newbyte bytespec integer)
"Args: (newbyte bytespec integer)
Replaces the specified byte of INTEGER with NEWBYTE (an integer) and returns
the result."
(logxor integer
(mask-field bytespec integer)
(ash (logandc2 newbyte
(- (ash 1 (byte-size bytespec))))
(byte-position bytespec))))
(defun deposit-field (newbyte bytespec integer)
"Args: (integer1 bytespec integer2)
Returns an integer represented by the bit sequence obtained by replacing the
specified bits of INTEGER2 with the specified bits of INTEGER1."
(dpb (ash newbyte (- (byte-position bytespec))) bytespec integer))