1 | ;;; -*- Mode: Lisp; Package: CCL -*- |
---|
2 | ;;; |
---|
3 | ;;; Copyright (C) 2009 Clozure Associates |
---|
4 | ;;; Copyright (C) 1994-2001 Digitool, Inc |
---|
5 | ;;; This file is part of Clozure CL. |
---|
6 | ;;; |
---|
7 | ;;; Clozure CL is licensed under the terms of the Lisp Lesser GNU Public |
---|
8 | ;;; License , known as the LLGPL and distributed with Clozure CL as the |
---|
9 | ;;; file "LICENSE". The LLGPL consists of a preamble and the LGPL, |
---|
10 | ;;; which is distributed with Clozure CL as the file "LGPL". Where these |
---|
11 | ;;; conflict, the preamble takes precedence. |
---|
12 | ;;; |
---|
13 | ;;; Clozure CL is referenced in the preamble as the "LIBRARY." |
---|
14 | ;;; |
---|
15 | ;;; The LLGPL is also available online at |
---|
16 | ;;; http://opensource.franz.com/preamble.html |
---|
17 | ;;; |
---|
18 | ;;; level-0;l0-float.lisp |
---|
19 | |
---|
20 | (in-package "CCL") |
---|
21 | |
---|
22 | (eval-when (:compile-toplevel :execute) |
---|
23 | (require "NUMBER-MACROS") |
---|
24 | (require :number-case-macro) |
---|
25 | (defconstant single-float-pi (coerce pi 'single-float)) |
---|
26 | (defconstant double-float-half-pi (/ pi 2)) |
---|
27 | (defconstant single-float-half-pi (coerce (/ pi 2) 'single-float)) |
---|
28 | (defconstant single-float-log2 0.6931472) ; (log 2) |
---|
29 | (defconstant double-float-log2 0.6931471805599453d0) ; (log 2.0d0) |
---|
30 | (defconstant double-float-log2^23 15.942385152878742d0) ; (log (expt 2 23)) |
---|
31 | ) |
---|
32 | |
---|
33 | ;;; used by float reader |
---|
34 | (defun make-float-from-fixnums (hi lo exp sign &optional result) |
---|
35 | ;;(require-null-or-double-float-sym result) |
---|
36 | ;; maybe nuke all these require-types? |
---|
37 | ;;(setq hi (require-type hi 'fixnum)) |
---|
38 | ;;(setq lo (require-type lo 'fixnum)) |
---|
39 | ;;(setq exp (require-type exp 'fixnum)) |
---|
40 | ;;(setq sign (require-type sign 'fixnum)) |
---|
41 | (let ((the-float (or result (%make-dfloat)))) |
---|
42 | (%make-float-from-fixnums the-float hi lo exp sign) |
---|
43 | the-float)) |
---|
44 | |
---|
45 | |
---|
46 | #+32-bit-target |
---|
47 | (defun make-short-float-from-fixnums (significand biased-exp sign &optional result) |
---|
48 | (%make-short-float-from-fixnums (or result (%make-sfloat)) significand biased-exp sign)) |
---|
49 | |
---|
50 | #+64-bit-target |
---|
51 | (defun make-short-float-from-fixnums (significand biased-exp sign) |
---|
52 | (declare (fixnum significand biased-exp sign)) |
---|
53 | (host-single-float-from-unsigned-byte-32 |
---|
54 | (logior |
---|
55 | (the fixnum (if (< sign 0) (ash 1 31) 0)) |
---|
56 | (the fixnum (ash biased-exp IEEE-single-float-exponent-offset)) |
---|
57 | (the fixnum (logand significand |
---|
58 | (1- (ash 1 IEEE-single-float-hidden-bit))))))) |
---|
59 | |
---|
60 | |
---|
61 | (defun float-sign (n1 &optional n2) ; second arg silly |
---|
62 | "Return a floating-point number that has the same sign as |
---|
63 | FLOAT1 and, if FLOAT2 is given, has the same absolute value |
---|
64 | as FLOAT2." |
---|
65 | (if (and n2 (not (typep n2 'float))) |
---|
66 | (setq n2 (require-type n2 'float))) |
---|
67 | (number-case n1 |
---|
68 | (double-float |
---|
69 | (if (%double-float-sign n1) |
---|
70 | (if n2 |
---|
71 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2)) |
---|
72 | -1.0d0) |
---|
73 | (if n2 |
---|
74 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2) |
---|
75 | 1.0d0))) |
---|
76 | (short-float |
---|
77 | (if (%short-float-sign n1) |
---|
78 | (if n2 |
---|
79 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2)) |
---|
80 | -1.0s0) |
---|
81 | (if n2 |
---|
82 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2) |
---|
83 | 1.0s0))))) |
---|
84 | |
---|
85 | |
---|
86 | |
---|
87 | (defun %double-float-minusp (n) |
---|
88 | (and (%double-float-sign n)(not (%double-float-zerop n)))) |
---|
89 | |
---|
90 | (defun %short-float-minusp (n) |
---|
91 | (and (%short-float-sign n) (not (%short-float-zerop n)))) |
---|
92 | |
---|
93 | (defun %double-float-abs (n) |
---|
94 | (if (not (%double-float-sign n)) |
---|
95 | n |
---|
96 | (%%double-float-abs! n (%make-dfloat)))) |
---|
97 | |
---|
98 | #+32-bit-target |
---|
99 | (defun %short-float-abs (n) |
---|
100 | (if (not (%short-float-sign n)) |
---|
101 | n |
---|
102 | (%%short-float-abs! n (%make-sfloat)))) |
---|
103 | |
---|
104 | (defun fixnum-decode-float (n) |
---|
105 | (etypecase n |
---|
106 | (double-float (%integer-decode-double-float n)))) |
---|
107 | |
---|
108 | (defun nan-or-infinity-p (n) |
---|
109 | (etypecase n |
---|
110 | (double-float (eq 2047 (%double-float-exp n))) |
---|
111 | (short-float (eq 255 (%short-float-exp n))))) |
---|
112 | |
---|
113 | ; not sure this is right |
---|
114 | (defun infinity-p (n) |
---|
115 | (etypecase n |
---|
116 | (double-float (multiple-value-bind (hi lo exp)(fixnum-decode-float n) |
---|
117 | (and (eq 2047 exp) |
---|
118 | (eq #x1000000 hi) |
---|
119 | (eq 0 lo)))) |
---|
120 | (short-float |
---|
121 | #+32-bit-target |
---|
122 | (multiple-value-bind (high low)(%sfloat-hwords n) |
---|
123 | (let* ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F)))) |
---|
124 | (exp (%ilsr 7 (%ilogand2 high #x7F80)))) |
---|
125 | (and (eq exp 255) |
---|
126 | (eq 0 mantissa)))) |
---|
127 | #+64-bit-target |
---|
128 | (let* ((bits (single-float-bits n)) |
---|
129 | (exp (ldb (byte IEEE-single-float-exponent-width |
---|
130 | IEEE-single-float-exponent-offset) |
---|
131 | bits)) |
---|
132 | (mantissa (ldb (byte IEEE-single-float-mantissa-width |
---|
133 | IEEE-single-float-mantissa-offset) |
---|
134 | bits))) |
---|
135 | (declare (fixnum bits exp mantissa)) |
---|
136 | (and (= exp 255) |
---|
137 | (zerop mantissa)))))) |
---|
138 | |
---|
139 | #+32-bit-target |
---|
140 | (defun fixnum-decode-short-float (float) |
---|
141 | (multiple-value-bind (high low)(%sfloat-hwords float) |
---|
142 | (let* ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F)))) |
---|
143 | (exp (%ilsr 7 (%ilogand2 high #x7F80)))) |
---|
144 | (if (and (neq exp 0)(neq exp 255))(setq mantissa (%ilogior mantissa #x800000))) |
---|
145 | (values mantissa exp (%ilsr 15 high))))) |
---|
146 | |
---|
147 | #+64-bit-target |
---|
148 | (defun fixnum-decode-short-float (float) |
---|
149 | (let* ((bits (single-float-bits float))) |
---|
150 | (declare (fixnum bits)) |
---|
151 | (let* ((mantissa (ldb (byte IEEE-single-float-mantissa-width |
---|
152 | IEEE-single-float-mantissa-offset) |
---|
153 | bits)) |
---|
154 | (exp (ldb (byte IEEE-single-float-exponent-width |
---|
155 | IEEE-single-float-exponent-offset) |
---|
156 | bits)) |
---|
157 | (sign (lsh bits -31))) |
---|
158 | (declare (fixnum mantissa exp sign)) |
---|
159 | (unless (or (= exp 0) (= exp 255)) |
---|
160 | (setq mantissa (logior mantissa (ash 1 IEEE-single-float-hidden-bit)))) |
---|
161 | (values mantissa exp sign)))) |
---|
162 | |
---|
163 | |
---|
164 | |
---|
165 | #+32-bit-target |
---|
166 | (defun integer-decode-double-float (n) |
---|
167 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
168 | ; is only 53 bits and positive so should be easy |
---|
169 | ;(values (logior (ash hi 28) lo) exp sign))) |
---|
170 | ; if denormalized, may fit in a fixnum |
---|
171 | (setq exp (- exp (if (< hi #x1000000) |
---|
172 | (+ IEEE-double-float-mantissa-width IEEE-double-float-bias) |
---|
173 | (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))) |
---|
174 | (if (< hi (ash 1 (1- target::fixnumshift))) ; aka 2 |
---|
175 | (values (logior (ash hi 28) lo) exp sign) |
---|
176 | ; might fit in 1 word? |
---|
177 | (let ((big (%alloc-misc 2 target::subtag-bignum))) |
---|
178 | (make-big-53 hi lo big) |
---|
179 | (if (< hi #x1000000) (%normalize-bignum big)) |
---|
180 | (values big exp sign))))) |
---|
181 | |
---|
182 | #+64-bit-target |
---|
183 | (defun integer-decode-double-float (n) |
---|
184 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
185 | (setq exp (- exp (if (< hi #x1000000) |
---|
186 | (+ IEEE-double-float-mantissa-width IEEE-double-float-bias) |
---|
187 | (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))) |
---|
188 | (values (logior (ash hi 28) lo) exp sign))) |
---|
189 | |
---|
190 | |
---|
191 | ;;; actually only called when magnitude bigger than a fixnum |
---|
192 | #+32-bit-target |
---|
193 | (defun %truncate-double-float (n) |
---|
194 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
195 | (if (< exp (1+ IEEE-double-float-bias)) ; this is false in practice |
---|
196 | 0 |
---|
197 | (progn |
---|
198 | (setq exp (- exp (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))) |
---|
199 | (if (eq sign 1) ; positive |
---|
200 | (logior (ash hi (+ 28 exp))(ash lo exp)) |
---|
201 | (if (<= exp 0) ; exp positive - negate before shift - else after |
---|
202 | (let ((poo (logior (ash hi (+ 28 exp))(ash lo exp)))) |
---|
203 | (- poo)) |
---|
204 | (let ((poo (logior (ash hi 28) lo))) |
---|
205 | (ash (- poo) exp)))))))) |
---|
206 | |
---|
207 | #+64-bit-target |
---|
208 | (defun %truncate-double-float (n) |
---|
209 | (multiple-value-bind (mantissa exp sign) (integer-decode-float n) |
---|
210 | (* sign (ash mantissa exp)))) |
---|
211 | |
---|
212 | |
---|
213 | |
---|
214 | ; actually only called when bigger than a fixnum |
---|
215 | (defun %truncate-short-float (n) |
---|
216 | (multiple-value-bind (mantissa exp sign)(fixnum-decode-short-float n) |
---|
217 | (if (< exp (1+ IEEE-single-float-bias)) ; is magnitude less than 1 - false in practice |
---|
218 | 0 |
---|
219 | (progn |
---|
220 | (setq exp (- exp (+ IEEE-single-float-mantissa-width (1+ IEEE-single-float-bias)))) |
---|
221 | (ash (if (eq sign 0) mantissa (- mantissa)) exp))))) |
---|
222 | |
---|
223 | (defun decode-float (n) |
---|
224 | "Return three values: |
---|
225 | 1) a floating-point number representing the significand. This is always |
---|
226 | between 0.5 (inclusive) and 1.0 (exclusive). |
---|
227 | 2) an integer representing the exponent. |
---|
228 | 3) -1.0 or 1.0 (i.e. the sign of the argument.)" |
---|
229 | (number-case n |
---|
230 | (double-float |
---|
231 | (let* ((old-exp (%double-float-exp n)) |
---|
232 | (sign (if (%double-float-sign n) -1.0d0 1.0d0))) |
---|
233 | (if (eq 0 old-exp) |
---|
234 | (if (%double-float-zerop n) |
---|
235 | (values 0.0d0 0 sign) |
---|
236 | (let* ((val (%make-dfloat)) |
---|
237 | (zeros (dfloat-significand-zeros n))) |
---|
238 | (%%double-float-abs! n val) |
---|
239 | (%%scale-dfloat! val (+ 2 IEEE-double-float-bias zeros) val) ; get it normalized |
---|
240 | (set-%double-float-exp val IEEE-double-float-bias) ; then bash exponent |
---|
241 | (values val (- old-exp zeros IEEE-double-float-bias) sign))) |
---|
242 | (if (> old-exp IEEE-double-float-normal-exponent-max) |
---|
243 | (error "Can't decode NAN or infinity ~s" n) |
---|
244 | (let ((val (%make-dfloat))) |
---|
245 | (%%double-float-abs! n val) |
---|
246 | (set-%double-float-exp val IEEE-double-float-bias) |
---|
247 | (values val (- old-exp IEEE-double-float-bias) sign)))))) |
---|
248 | (short-float |
---|
249 | (let* ((old-exp (%short-float-exp n)) |
---|
250 | (sign (if (%short-float-sign n) -1.0s0 1.0s0))) |
---|
251 | (if (eq 0 old-exp) |
---|
252 | (if (%short-float-zerop n) |
---|
253 | (values 0.0s0 0 sign) |
---|
254 | #+32-bit-target |
---|
255 | (let* ((val (%make-sfloat)) |
---|
256 | (zeros (sfloat-significand-zeros n))) |
---|
257 | (%%short-float-abs! n val) |
---|
258 | (%%scale-sfloat! val (+ 2 IEEE-single-float-bias zeros) val) ; get it normalized |
---|
259 | (set-%short-float-exp val IEEE-single-float-bias) ; then bash exponent |
---|
260 | (values val (- old-exp zeros IEEE-single-float-bias) sign)) |
---|
261 | #+64-bit-target |
---|
262 | (let* ((zeros (sfloat-significand-zeros n)) |
---|
263 | (val (%%scale-sfloat (%short-float-abs n) |
---|
264 | (+ 2 IEEE-single-float-bias zeros)))) |
---|
265 | (values (set-%short-float-exp val IEEE-single-float-bias) |
---|
266 | (- old-exp zeros IEEE-single-float-bias) sign))) |
---|
267 | (if (> old-exp IEEE-single-float-normal-exponent-max) |
---|
268 | (error "Can't decode NAN or infinity ~s" n) |
---|
269 | #+32-bit-target |
---|
270 | (let ((val (%make-sfloat))) |
---|
271 | (%%short-float-abs! n val) |
---|
272 | (set-%short-float-exp val IEEE-single-float-bias) |
---|
273 | (values val (- old-exp IEEE-single-float-bias) sign)) |
---|
274 | #+64-bit-target |
---|
275 | (values (set-%short-float-exp (%short-float-abs n) |
---|
276 | IEEE-single-float-bias) |
---|
277 | (- old-exp IEEE-single-float-bias) sign))))))) |
---|
278 | |
---|
279 | ; (* float (expt 2 int)) |
---|
280 | (defun scale-float (float int) |
---|
281 | "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss |
---|
282 | of precision or overflow." |
---|
283 | (unless (fixnump int)(setq int (require-type int 'fixnum))) |
---|
284 | (number-case float |
---|
285 | (double-float |
---|
286 | (let* ((float-exp (%double-float-exp float)) |
---|
287 | (new-exp (+ float-exp int))) |
---|
288 | (if (eq 0 float-exp) ; already denormalized? |
---|
289 | (if (%double-float-zerop float) |
---|
290 | float |
---|
291 | (let ((result (%make-dfloat))) |
---|
292 | (%%scale-dfloat! float (+ (1+ IEEE-double-float-bias) int) result))) |
---|
293 | (if (<= new-exp 0) ; maybe going denormalized |
---|
294 | (if (<= new-exp (- IEEE-double-float-digits)) |
---|
295 | 0.0d0 ; should this be underflow? - should just be normal and result is fn of current fpu-mode |
---|
296 | ;(error "Can't scale ~s by ~s." float int) ; should signal something |
---|
297 | (let ((result (%make-dfloat))) |
---|
298 | (%copy-double-float float result) |
---|
299 | (set-%double-float-exp result 1) ; scale by float-exp -1 |
---|
300 | (%%scale-dfloat! result (+ IEEE-double-float-bias (+ float-exp int)) result) |
---|
301 | result)) |
---|
302 | (if (> new-exp IEEE-double-float-normal-exponent-max) |
---|
303 | (error (make-condition 'floating-point-overflow |
---|
304 | :operation 'scale-float |
---|
305 | :operands (list float int))) |
---|
306 | (let ((new-float (%make-dfloat))) |
---|
307 | (%copy-double-float float new-float) |
---|
308 | (set-%double-float-exp new-float new-exp) |
---|
309 | new-float)))))) |
---|
310 | (short-float |
---|
311 | (let* ((float-exp (%short-float-exp float)) |
---|
312 | (new-exp (+ float-exp int))) |
---|
313 | (if (eq 0 float-exp) ; already denormalized? |
---|
314 | (if (%short-float-zerop float) |
---|
315 | float |
---|
316 | #+32-bit-target |
---|
317 | (let ((result (%make-sfloat))) |
---|
318 | (%%scale-sfloat! float (+ (1+ IEEE-single-float-bias) int) result)) |
---|
319 | #+64-bit-target |
---|
320 | (%%scale-sfloat float (+ (1+ IEEE-single-float-bias) int))) |
---|
321 | (if (<= new-exp 0) ; maybe going denormalized |
---|
322 | (if (<= new-exp (- IEEE-single-float-digits)) |
---|
323 | ;; should this be underflow? - should just be normal and |
---|
324 | ;; result is fn of current fpu-mode (error "Can't scale |
---|
325 | ;; ~s by ~s." float int) ; should signal something |
---|
326 | 0.0s0 |
---|
327 | #+32-bit-target |
---|
328 | (let ((result (%make-sfloat))) |
---|
329 | (%copy-short-float float result) |
---|
330 | (set-%short-float-exp result 1) ; scale by float-exp -1 |
---|
331 | (%%scale-sfloat! result (+ IEEE-single-float-bias (+ float-exp int)) result) |
---|
332 | result) |
---|
333 | #+64-bit-target |
---|
334 | (%%scale-sfloat (set-%short-float-exp float 1) |
---|
335 | (+ IEEE-single-float-bias (+ float-exp int)))) |
---|
336 | (if (> new-exp IEEE-single-float-normal-exponent-max) |
---|
337 | (error (make-condition 'floating-point-overflow |
---|
338 | :operation 'scale-float |
---|
339 | :operands (list float int))) |
---|
340 | #+32-bit-target |
---|
341 | (let ((new-float (%make-sfloat))) |
---|
342 | (%copy-short-float float new-float) |
---|
343 | (set-%short-float-exp new-float new-exp) |
---|
344 | new-float) |
---|
345 | #+64-bit-target |
---|
346 | (set-%short-float-exp float new-exp)))))))) |
---|
347 | |
---|
348 | (defun %copy-float (f) |
---|
349 | ;Returns a freshly consed float. float can also be a macptr. |
---|
350 | (cond ((double-float-p f) (%copy-double-float f (%make-dfloat))) |
---|
351 | ((macptrp f) |
---|
352 | (let ((float (%make-dfloat))) |
---|
353 | (%copy-ptr-to-ivector f 0 float (* 4 target::double-float.value-cell) 8) |
---|
354 | float)) |
---|
355 | (t (error "Illegal arg ~s to %copy-float" f)))) |
---|
356 | |
---|
357 | (defun float-precision (float) ; not used - not in cltl2 index ? |
---|
358 | "Return a non-negative number of significant digits in its float argument. |
---|
359 | Will be less than FLOAT-DIGITS if denormalized or zero." |
---|
360 | (number-case float |
---|
361 | (double-float |
---|
362 | (if (eq 0 (%double-float-exp float)) |
---|
363 | (if (not (%double-float-zerop float)) |
---|
364 | ; denormalized |
---|
365 | (- IEEE-double-float-mantissa-width (dfloat-significand-zeros float)) |
---|
366 | 0) |
---|
367 | IEEE-double-float-digits)) |
---|
368 | (short-float |
---|
369 | (if (eq 0 (%short-float-exp float)) |
---|
370 | (if (not (%short-float-zerop float)) |
---|
371 | ; denormalized |
---|
372 | (- IEEE-single-float-mantissa-width (sfloat-significand-zeros float)) |
---|
373 | 0) |
---|
374 | IEEE-single-float-digits)))) |
---|
375 | |
---|
376 | |
---|
377 | (defun %double-float (number &optional result) |
---|
378 | ;(require-null-or-double-float-sym result) |
---|
379 | ; use number-case when macro is common |
---|
380 | (number-case number |
---|
381 | (double-float |
---|
382 | (if result |
---|
383 | (%copy-double-float number result) |
---|
384 | number)) |
---|
385 | (short-float |
---|
386 | (%short-float->double-float number (or result (%make-dfloat)))) |
---|
387 | (fixnum |
---|
388 | (%fixnum-dfloat number (or result (%make-dfloat)))) |
---|
389 | (bignum (%bignum-dfloat number result)) |
---|
390 | (ratio |
---|
391 | (if (not result)(setq result (%make-dfloat))) |
---|
392 | (let* ((num (%numerator number)) |
---|
393 | (den (%denominator number))) |
---|
394 | ; dont error if result is floatable when either top or bottom is not. |
---|
395 | ; maybe do usual first, catching error |
---|
396 | (if (not (or (bignump num)(bignump den))) |
---|
397 | (with-stack-double-floats ((fnum num) |
---|
398 | (fden den)) |
---|
399 | (%double-float/-2! fnum fden result)) |
---|
400 | (let* ((numlen (integer-length num)) |
---|
401 | (denlen (integer-length den)) |
---|
402 | (exp (- numlen denlen)) |
---|
403 | (minusp (minusp num))) |
---|
404 | (if (and (<= numlen IEEE-double-float-bias) |
---|
405 | (<= denlen IEEE-double-float-bias) |
---|
406 | #|(not (minusp exp))|# |
---|
407 | (<= (abs exp) IEEE-double-float-mantissa-width)) |
---|
408 | (with-stack-double-floats ((fnum num) |
---|
409 | (fden den)) |
---|
410 | |
---|
411 | (%double-float/-2! fnum fden result)) |
---|
412 | (if (> exp IEEE-double-float-mantissa-width) |
---|
413 | (progn (%double-float (round num den) result)) |
---|
414 | (if (>= exp 0) |
---|
415 | ; exp between 0 and 53 and nums big |
---|
416 | (let* ((shift (- IEEE-double-float-digits exp)) |
---|
417 | (num (if minusp (- num) num)) |
---|
418 | (int (round (ash num shift) den)) ; gaak |
---|
419 | (intlen (integer-length int)) |
---|
420 | (new-exp (+ intlen (- IEEE-double-float-bias shift)))) |
---|
421 | |
---|
422 | (when (> intlen IEEE-double-float-digits) |
---|
423 | (setq shift (1- shift)) |
---|
424 | (setq int (round (ash num shift) den)) |
---|
425 | (setq intlen (integer-length int)) |
---|
426 | (setq new-exp (+ intlen (- IEEE-double-float-bias shift)))) |
---|
427 | (when (> new-exp 2046) |
---|
428 | (error (make-condition 'floating-point-overflow |
---|
429 | :operation 'double-float |
---|
430 | :operands (list number)))) |
---|
431 | (make-float-from-fixnums (ldb (byte 25 (- intlen 25)) int) |
---|
432 | (ldb (byte 28 (max (- intlen 53) 0)) int) |
---|
433 | new-exp ;(+ intlen (- IEEE-double-float-bias 53)) |
---|
434 | (if minusp -1 1) |
---|
435 | result)) |
---|
436 | ; den > num - exp negative |
---|
437 | (progn |
---|
438 | (float-rat-neg-exp num den (if minusp -1 1) result))))))))))) |
---|
439 | |
---|
440 | |
---|
441 | #+32-bit-target |
---|
442 | (defun %short-float-ratio (number &optional result) |
---|
443 | (if (not result)(setq result (%make-sfloat))) |
---|
444 | (let* ((num (%numerator number)) |
---|
445 | (den (%denominator number))) |
---|
446 | ;; dont error if result is floatable when either top or bottom is |
---|
447 | ;; not. maybe do usual first, catching error |
---|
448 | (if (not (or (bignump num)(bignump den))) |
---|
449 | (target::with-stack-short-floats ((fnum num) |
---|
450 | (fden den)) |
---|
451 | (%short-float/-2! fnum fden result)) |
---|
452 | (let* ((numlen (integer-length num)) |
---|
453 | (denlen (integer-length den)) |
---|
454 | (exp (- numlen denlen)) |
---|
455 | (minusp (minusp num))) |
---|
456 | (if (and (<= numlen IEEE-single-float-bias) |
---|
457 | (<= denlen IEEE-single-float-bias) |
---|
458 | #|(not (minusp exp))|# |
---|
459 | (<= (abs exp) IEEE-single-float-mantissa-width)) |
---|
460 | (target::with-stack-short-floats ((fnum num) |
---|
461 | (fden den)) |
---|
462 | (%short-float/-2! fnum fden result)) |
---|
463 | (if (> exp IEEE-single-float-mantissa-width) |
---|
464 | (progn (%short-float (round num den) result)) |
---|
465 | (if (>= exp 0) |
---|
466 | ; exp between 0 and 23 and nums big |
---|
467 | (let* ((shift (- IEEE-single-float-digits exp)) |
---|
468 | (num (if minusp (- num) num)) |
---|
469 | (int (round (ash num shift) den)) ; gaak |
---|
470 | (intlen (integer-length int)) |
---|
471 | (new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
472 | (when (> intlen IEEE-single-float-digits) |
---|
473 | (setq shift (1- shift)) |
---|
474 | (setq int (round (ash num shift) den)) |
---|
475 | (setq intlen (integer-length int)) |
---|
476 | (setq new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
477 | (when (> new-exp IEEE-single-float-normal-exponent-max) |
---|
478 | (error (make-condition 'floating-point-overflow |
---|
479 | :operation 'short-float |
---|
480 | :operands (list number)))) |
---|
481 | (make-short-float-from-fixnums |
---|
482 | (ldb (byte IEEE-single-float-digits (- intlen IEEE-single-float-digits)) int) |
---|
483 | new-exp |
---|
484 | (if minusp -1 1) |
---|
485 | result)) |
---|
486 | ; den > num - exp negative |
---|
487 | (progn |
---|
488 | (float-rat-neg-exp num den (if minusp -1 1) result t))))))))) |
---|
489 | |
---|
490 | #+64-bit-target |
---|
491 | (defun %short-float-ratio (number) |
---|
492 | (let* ((num (%numerator number)) |
---|
493 | (den (%denominator number))) |
---|
494 | ;; dont error if result is floatable when either top or bottom is |
---|
495 | ;; not. maybe do usual first, catching error |
---|
496 | (if (not (or (bignump num)(bignump den))) |
---|
497 | (/ (the short-float (%short-float num)) |
---|
498 | (the short-float (%short-float den))) |
---|
499 | (let* ((numlen (integer-length num)) |
---|
500 | (denlen (integer-length den)) |
---|
501 | (exp (- numlen denlen)) |
---|
502 | (minusp (minusp num))) |
---|
503 | (if (and (<= numlen IEEE-single-float-bias) |
---|
504 | (<= denlen IEEE-single-float-bias) |
---|
505 | #|(not (minusp exp))|# |
---|
506 | (<= (abs exp) IEEE-single-float-mantissa-width)) |
---|
507 | (/ (the short-float (%short-float num)) |
---|
508 | (the short-float (%short-float den))) |
---|
509 | (if (> exp IEEE-single-float-mantissa-width) |
---|
510 | (progn (%short-float (round num den))) |
---|
511 | (if (>= exp 0) |
---|
512 | ; exp between 0 and 23 and nums big |
---|
513 | (let* ((shift (- IEEE-single-float-digits exp)) |
---|
514 | (num (if minusp (- num) num)) |
---|
515 | (int (round (ash num shift) den)) ; gaak |
---|
516 | (intlen (integer-length int)) |
---|
517 | (new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
518 | (when (> intlen IEEE-single-float-digits) |
---|
519 | (setq shift (1- shift)) |
---|
520 | (setq int (round (ash num shift) den)) |
---|
521 | (setq intlen (integer-length int)) |
---|
522 | (setq new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
523 | (when (> new-exp IEEE-single-float-normal-exponent-max) |
---|
524 | (error (make-condition 'floating-point-overflow |
---|
525 | :operation 'short-float |
---|
526 | :operands (list number)))) |
---|
527 | (make-short-float-from-fixnums |
---|
528 | (ldb (byte IEEE-single-float-digits (- intlen IEEE-single-float-digits)) int) |
---|
529 | new-exp |
---|
530 | (if minusp 1 0))) |
---|
531 | ; den > num - exp negative |
---|
532 | (progn |
---|
533 | (float-rat-neg-exp num den (if minusp -1 1) nil t))))))))) |
---|
534 | |
---|
535 | |
---|
536 | #+32-bit-target |
---|
537 | (defun %short-float (number &optional result) |
---|
538 | (number-case number |
---|
539 | (short-float |
---|
540 | (if result (%copy-short-float number result) number)) |
---|
541 | (double-float |
---|
542 | (%double-float->short-float number (or result (%make-sfloat)))) |
---|
543 | (fixnum |
---|
544 | (%fixnum-sfloat number (or result (%make-sfloat)))) |
---|
545 | (bignum |
---|
546 | (%bignum-sfloat number (or result (%make-sfloat)))) |
---|
547 | (ratio |
---|
548 | (%short-float-ratio number result)))) |
---|
549 | |
---|
550 | #+64-bit-target |
---|
551 | (defun %short-float (number) |
---|
552 | (number-case number |
---|
553 | (short-float number) |
---|
554 | (double-float (%double-float->short-float number)) |
---|
555 | (fixnum (%fixnum-sfloat number)) |
---|
556 | (bignum (%bignum-sfloat number)) |
---|
557 | (ratio (%short-float-ratio number)))) |
---|
558 | |
---|
559 | |
---|
560 | (defun float-rat-neg-exp (integer divisor sign &optional result short) |
---|
561 | (if (minusp sign)(setq integer (- integer))) |
---|
562 | (let* ((integer-length (integer-length integer)) |
---|
563 | ;; make sure we will have enough bits in the quotient |
---|
564 | ;; (and a couple extra for rounding) |
---|
565 | (shift-factor (+ (- (integer-length divisor) integer-length) (if short 28 60))) ; fix |
---|
566 | (scaled-integer integer)) |
---|
567 | (if (plusp shift-factor) |
---|
568 | (setq scaled-integer (ash integer shift-factor)) |
---|
569 | (setq divisor (ash divisor (- shift-factor))) ; assume div > num |
---|
570 | ) |
---|
571 | ;(pprint (list shift-factor scaled-integer divisor)) |
---|
572 | (multiple-value-bind (quotient remainder)(floor scaled-integer divisor) |
---|
573 | (unless (zerop remainder) ; whats this - tells us there's junk below |
---|
574 | (setq quotient (logior quotient 1))) |
---|
575 | ;; why do it return 2 values? |
---|
576 | (values (float-and-scale-and-round sign quotient (- shift-factor) short result))))) |
---|
577 | |
---|
578 | |
---|
579 | |
---|
580 | ;;; when is (negate-bignum (bignum-ashift-right big)) ; can't negate |
---|
581 | ;;; in place cause may get bigger cheaper than (negate-bignum big) - 6 |
---|
582 | ;;; 0r 8 digits ; 8 longs so win if digits > 7 or negate it on the |
---|
583 | ;;; stack |
---|
584 | |
---|
585 | (defun %bignum-dfloat (big &optional result) |
---|
586 | (let* ((minusp (bignum-minusp big))) |
---|
587 | (flet |
---|
588 | ((doit (new-big) |
---|
589 | (let* ((int-len (bignum-integer-length new-big))) |
---|
590 | (when (>= int-len (- 2047 IEEE-double-float-bias)) ; args? |
---|
591 | (error (make-condition 'floating-point-overflow |
---|
592 | :operation 'float :operands (list big)))) |
---|
593 | (if (> int-len 53) |
---|
594 | (let* ((hi (ldb (byte 25 (- int-len 25)) new-big)) |
---|
595 | (lo (ldb (byte 28 (- int-len 53)) new-big))) |
---|
596 | ;(print (list new-big hi lo)) |
---|
597 | (when (and (logbitp (- int-len 54) new-big) ; round bit |
---|
598 | (or (%ilogbitp 0 lo) ; oddp |
---|
599 | ;; or more bits below round |
---|
600 | (%i< (one-bignum-factor-of-two new-big) (- int-len 54)))) |
---|
601 | (if (eq lo #xfffffff) |
---|
602 | (setq hi (1+ hi) lo 0) |
---|
603 | (setq lo (1+ lo))) |
---|
604 | (when (%ilogbitp 25 hi) ; got bigger |
---|
605 | (setq int-len (1+ int-len)) |
---|
606 | (let ((bit (%ilogbitp 0 hi))) |
---|
607 | (setq hi (%ilsr 1 hi)) |
---|
608 | (setq lo (%ilsr 1 lo)) |
---|
609 | (if bit (setq lo (%ilogior #x8000000 lo)))))) |
---|
610 | (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len)(if minusp -1 1) result)) |
---|
611 | (let* ((hi (ldb (byte 25 (- int-len 25)) new-big)) |
---|
612 | (lobits (min (- int-len 25) 28)) |
---|
613 | (lo (ldb (byte lobits (- int-len (+ lobits 25))) new-big))) |
---|
614 | (if (< lobits 28) (setq lo (ash lo (- 28 lobits)))) |
---|
615 | (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len) (if minusp -1 1) result)))))) |
---|
616 | (declare (dynamic-extent #'doit)) |
---|
617 | (with-one-negated-bignum-buffer big doit)))) |
---|
618 | |
---|
619 | #+32-bit-target |
---|
620 | (defun %bignum-sfloat (big &optional result) |
---|
621 | (let* ((minusp (bignum-minusp big))) |
---|
622 | (flet |
---|
623 | ((doit (new-big) |
---|
624 | (let* ((int-len (bignum-integer-length new-big))) |
---|
625 | (when (>= int-len (- 255 IEEE-single-float-bias)) ; args? |
---|
626 | (error (make-condition 'floating-point-overflow |
---|
627 | :operation 'float :operands (list big 1.0s0)))) |
---|
628 | (if t ;(> int-len IEEE-single-float-digits) ; always true |
---|
629 | (let* ((lo (ldb (byte IEEE-single-float-digits (- int-len IEEE-single-float-digits)) new-big))) |
---|
630 | (when (and (logbitp (- int-len 25) new-big) ; round bit |
---|
631 | (or (%ilogbitp 0 lo) ; oddp |
---|
632 | ; or more bits below round |
---|
633 | (%i< (one-bignum-factor-of-two new-big) (- int-len 25)))) |
---|
634 | (setq lo (1+ lo)) |
---|
635 | (when (%ilogbitp 24 lo) ; got bigger |
---|
636 | (setq int-len (1+ int-len)) |
---|
637 | (setq lo (%ilsr 1 lo)))) |
---|
638 | (make-short-float-from-fixnums lo (+ IEEE-single-float-bias int-len)(if minusp -1 1) result)) |
---|
639 | )))) |
---|
640 | (declare (dynamic-extent #'doit)) |
---|
641 | (with-one-negated-bignum-buffer big doit)))) |
---|
642 | |
---|
643 | |
---|
644 | #+64-bit-target |
---|
645 | (defun %bignum-sfloat (big) |
---|
646 | (let* ((minusp (bignum-minusp big))) |
---|
647 | (flet |
---|
648 | ((doit (new-big) |
---|
649 | (let* ((int-len (bignum-integer-length new-big))) |
---|
650 | (when (>= int-len (- 255 IEEE-single-float-bias)) ; args? |
---|
651 | (error (make-condition 'floating-point-overflow |
---|
652 | :operation 'float :operands (list big 1.0s0)))) |
---|
653 | (if t ;(> int-len IEEE-single-float-digits) ; always true |
---|
654 | (let* ((lo (ldb (byte IEEE-single-float-digits (- int-len IEEE-single-float-digits)) new-big))) |
---|
655 | (when (and (logbitp (- int-len 25) new-big) ; round bit |
---|
656 | (or (%ilogbitp 0 lo) ; oddp |
---|
657 | ; or more bits below round |
---|
658 | (%i< (one-bignum-factor-of-two new-big) (- int-len 25)))) |
---|
659 | (setq lo (1+ lo)) |
---|
660 | (when (%ilogbitp 24 lo) ; got bigger |
---|
661 | (setq int-len (1+ int-len)) |
---|
662 | (setq lo (%ilsr 1 lo)))) |
---|
663 | (make-short-float-from-fixnums lo (+ IEEE-single-float-bias int-len)(if minusp -1 1))) |
---|
664 | )))) |
---|
665 | (declare (dynamic-extent #'doit)) |
---|
666 | (with-one-negated-bignum-buffer big doit)))) |
---|
667 | |
---|
668 | |
---|
669 | |
---|
670 | |
---|
671 | (defun %fixnum-dfloat (fix &optional result) |
---|
672 | (if (eq 0 fix) |
---|
673 | (if result (%copy-double-float 0.0d0 result) 0.0d0) |
---|
674 | (progn |
---|
675 | (when (not result)(setq result (%make-dfloat))) |
---|
676 | ; it better return result |
---|
677 | (%int-to-dfloat fix result)))) |
---|
678 | |
---|
679 | |
---|
680 | #+32-bit-target |
---|
681 | (defun %fixnum-sfloat (fix &optional result) |
---|
682 | (if (eq 0 fix) |
---|
683 | (if result (%copy-short-float 0.0s0 result) 0.0s0) |
---|
684 | (%int-to-sfloat! fix (or result (%make-sfloat))))) |
---|
685 | |
---|
686 | #+64-bit-target |
---|
687 | (defun %fixnum-sfloat (fix) |
---|
688 | (if (eq 0 fix) |
---|
689 | 0.0s0 |
---|
690 | (%int-to-sfloat fix))) |
---|
691 | |
---|
692 | ;;; Transcendental functions. |
---|
693 | (defun sin (x) |
---|
694 | "Return the sine of NUMBER." |
---|
695 | (if (complexp x) |
---|
696 | (let* ((r (realpart x)) |
---|
697 | (i (imagpart x))) |
---|
698 | (complex (* (sin r) (cosh i)) |
---|
699 | (* (cos r) (sinh i)))) |
---|
700 | (if (typep x 'double-float) |
---|
701 | (%double-float-sin! x (%make-dfloat)) |
---|
702 | #+32-bit-target |
---|
703 | (target::with-stack-short-floats ((sx x)) |
---|
704 | (%single-float-sin! sx (%make-sfloat))) |
---|
705 | #+64-bit-target |
---|
706 | (%single-float-sin (%short-float x))))) |
---|
707 | |
---|
708 | (defun cos (x) |
---|
709 | "Return the cosine of NUMBER." |
---|
710 | (if (complexp x) |
---|
711 | (let* ((r (realpart x)) |
---|
712 | (i (imagpart x))) |
---|
713 | (complex (* (cos r) (cosh i)) |
---|
714 | (- (* (sin r) (sinh i))))) |
---|
715 | (if (typep x 'double-float) |
---|
716 | (%double-float-cos! x (%make-dfloat)) |
---|
717 | #+32-bit-target |
---|
718 | (target::with-stack-short-floats ((sx x)) |
---|
719 | (%single-float-cos! sx (%make-sfloat))) |
---|
720 | #+64-bit-target |
---|
721 | (%single-float-cos (%short-float x))))) |
---|
722 | |
---|
723 | (defun tan (x) |
---|
724 | "Return the tangent of NUMBER." |
---|
725 | (if (complexp x) |
---|
726 | (/ (sin x) (cos x)) |
---|
727 | (if (typep x 'double-float) |
---|
728 | (%double-float-tan! x (%make-dfloat)) |
---|
729 | #+32-bit-target |
---|
730 | (target::with-stack-short-floats ((sx x)) |
---|
731 | (%single-float-tan! sx (%make-sfloat))) |
---|
732 | #+64-bit-target |
---|
733 | (%single-float-tan (%short-float x)) |
---|
734 | ))) |
---|
735 | |
---|
736 | |
---|
737 | ;;; Multiply by i (with additional optional scale factor) |
---|
738 | ;;; Does the "right thing" with minus zeroes (see CLTL2) |
---|
739 | (defun i* (number &optional (scale 1)) |
---|
740 | (complex (* (- scale) (imagpart number)) |
---|
741 | (* scale (realpart number)))) |
---|
742 | |
---|
743 | ;;; complex atanh |
---|
744 | (defun %complex-atanh (z) |
---|
745 | (let* ((rx (realpart z)) |
---|
746 | (ix (imagpart z)) |
---|
747 | (sign (typecase rx |
---|
748 | (double-float (%double-float-sign rx)) |
---|
749 | (short-float (%short-float-sign rx)) |
---|
750 | (t (minusp rx)))) |
---|
751 | (x rx) |
---|
752 | (y ix) |
---|
753 | (y1 (abs y)) |
---|
754 | ra ia) |
---|
755 | ;; following code requires non-negative x |
---|
756 | (when sign |
---|
757 | (setf x (- x)) |
---|
758 | (setf y (- y))) |
---|
759 | (cond ((> (max x y1) 1.8014399e+16) |
---|
760 | ;; large value escape |
---|
761 | (setq ra (if (> x y1) |
---|
762 | (let ((r (/ y x))) |
---|
763 | (/ (/ x) (1+ (* r r)))) |
---|
764 | (let ((r (/ x y))) |
---|
765 | (/ (/ r y) (1+ (* r r)))))) |
---|
766 | (setq ia (typecase y |
---|
767 | (double-float (float-sign y double-float-half-pi)) |
---|
768 | (single-float (float-sign y single-float-half-pi)) |
---|
769 | (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi))))) |
---|
770 | ((= x 1) |
---|
771 | (setq ra (if (< y1 1e-9) |
---|
772 | (/ (log-e (/ 2 y1)) 2) |
---|
773 | (/ (log1+ (/ 4 (* y y))) 4))) |
---|
774 | (setq ia (/ (atan (/ 2 y) -1) 2))) |
---|
775 | (t |
---|
776 | (let ((r2 (+ (* x x) (* y y)))) |
---|
777 | (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4)) |
---|
778 | (setq ia (/ (atan (* 2 y) (- 1 r2)) 2))))) |
---|
779 | ;; fixup signs, with special case for real arguments |
---|
780 | (cond (sign |
---|
781 | (setq ra (- ra)) |
---|
782 | (when (typep z 'complex) |
---|
783 | (setq ia (- ia)))) |
---|
784 | (t |
---|
785 | (unless (typep z 'complex) |
---|
786 | (setq ia (- ia))))) |
---|
787 | (complex ra ia))) |
---|
788 | |
---|
789 | (defun atan (y &optional (x nil x-p)) |
---|
790 | "Return the arc tangent of Y if X is omitted or Y/X if X is supplied." |
---|
791 | (cond (x-p |
---|
792 | (cond ((or (typep x 'double-float) |
---|
793 | (typep y 'double-float)) |
---|
794 | (with-stack-double-floats ((dy y) |
---|
795 | (dx x)) |
---|
796 | (%df-atan2 dy dx))) |
---|
797 | (t |
---|
798 | (when (and (rationalp x) (rationalp y)) |
---|
799 | ;; rescale arguments so that the maximum absolute value is 1 |
---|
800 | (let ((x1 (abs x)) (y1 (abs y))) |
---|
801 | (cond ((> y1 x1) |
---|
802 | (setf x (/ x y1)) |
---|
803 | (setf y (signum y))) |
---|
804 | ((not (zerop x)) |
---|
805 | (setf y (/ y x1)) |
---|
806 | (setf x (signum x)))))) |
---|
807 | #+32-bit-target |
---|
808 | (target::with-stack-short-floats ((sy y) |
---|
809 | (sx x)) |
---|
810 | (%sf-atan2! sy sx)) |
---|
811 | #+64-bit-target |
---|
812 | (%sf-atan2 (%short-float y) (%short-float x))))) |
---|
813 | ((typep y 'double-float) |
---|
814 | (%double-float-atan! y (%make-dfloat))) |
---|
815 | ((typep y 'single-float) |
---|
816 | #+32-bit-target |
---|
817 | (%single-float-atan! y (%make-sfloat)) |
---|
818 | #+64-bit-target |
---|
819 | (%single-float-atan y)) |
---|
820 | ((typep y 'rational) |
---|
821 | (cond ((<= (abs y) most-positive-short-float) |
---|
822 | #+32-bit-target |
---|
823 | (target::with-stack-short-floats ((sy y)) |
---|
824 | (%single-float-atan! sy (%make-sfloat))) |
---|
825 | #+64-bit-target |
---|
826 | (%single-float-atan (%short-float y))) |
---|
827 | ((minusp y) |
---|
828 | #.(- single-float-half-pi)) |
---|
829 | (t |
---|
830 | single-float-half-pi))) |
---|
831 | (t |
---|
832 | (let ((r (realpart y)) |
---|
833 | (i (imagpart y))) |
---|
834 | (if (zerop i) |
---|
835 | (complex (atan r) i) |
---|
836 | (i* (%complex-atanh (complex (- i) r)) -1)))))) |
---|
837 | |
---|
838 | |
---|
839 | |
---|
840 | (defun log (x &optional (b nil b-p)) |
---|
841 | "Return the logarithm of NUMBER in the base BASE, which defaults to e." |
---|
842 | (if b-p |
---|
843 | (if (zerop b) |
---|
844 | (if (zerop x) |
---|
845 | (report-bad-arg x '(not (satisfies zerop) )) |
---|
846 | ;; ** CORRECT THE CONTAGION for complex args ** |
---|
847 | (+ 0 (* x b))) |
---|
848 | ;; do the float/rational contagion before the division |
---|
849 | ;; but take care with negative zeroes |
---|
850 | (let ((x1 (realpart x)) |
---|
851 | (b1 (realpart b))) |
---|
852 | (if (and (typep x1 'float) |
---|
853 | (typep b1 'float)) |
---|
854 | (/ (log-e (* x (float 1.0 b1))) |
---|
855 | (log-e (* b (float 1.0 x1)))) |
---|
856 | (let ((r (/ (cond ((typep x 'rational) |
---|
857 | (%rational-log x 1.0d0)) |
---|
858 | ((typep x1 'rational) |
---|
859 | (%rational-complex-log x 1.0d0)) |
---|
860 | (t |
---|
861 | (log-e (* x 1.0d0)))) |
---|
862 | (cond ((typep b 'rational) |
---|
863 | (%rational-log b 1.0d0)) |
---|
864 | ((typep b1 'rational) |
---|
865 | (%rational-complex-log b 1.0d0)) |
---|
866 | (t |
---|
867 | (log-e (* b 1.0d0))))))) |
---|
868 | (cond ((or (typep x1 'double-float) |
---|
869 | (typep b1 'double-float)) |
---|
870 | r) |
---|
871 | ((complexp r) |
---|
872 | (complex (%short-float (realpart r)) |
---|
873 | (%short-float (imagpart r)))) |
---|
874 | (t |
---|
875 | (%short-float r))))))) |
---|
876 | (log-e x))) |
---|
877 | |
---|
878 | |
---|
879 | |
---|
880 | (defun log-e (x) |
---|
881 | (cond |
---|
882 | ((typep x 'double-float) |
---|
883 | (if (%double-float-sign x) |
---|
884 | (with-stack-double-floats ((dx x)) |
---|
885 | (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi)) |
---|
886 | (%double-float-log! x (%make-dfloat)))) |
---|
887 | ((typep x 'short-float) |
---|
888 | #+32-bit-target |
---|
889 | (if (%short-float-sign x) |
---|
890 | (target::with-stack-short-floats ((sx x)) |
---|
891 | (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat)) |
---|
892 | single-float-pi)) |
---|
893 | (%single-float-log! x (%make-sfloat))) |
---|
894 | #+64-bit-target |
---|
895 | (if (%short-float-sign x) |
---|
896 | (complex (%single-float-log (%short-float-abs (%short-float x))) |
---|
897 | single-float-pi) |
---|
898 | (%single-float-log (%short-float x)))) |
---|
899 | ((typep x 'complex) |
---|
900 | (if (typep (realpart x) 'rational) |
---|
901 | (%rational-complex-log x 1.0s0) |
---|
902 | ;; take care that intermediate results do not overflow/underflow: |
---|
903 | ;; pre-scale argument by an appropriate power of two; |
---|
904 | ;; we only need to scale for very large and very small values - |
---|
905 | ;; hence the various 'magic' numbers (values may not be too |
---|
906 | ;; critical but do depend on the sqrt update to fix abs's operation) |
---|
907 | (let ((m (max (abs (realpart x)) |
---|
908 | (abs (imagpart x)))) |
---|
909 | (log-s 0) |
---|
910 | (s 1)) |
---|
911 | (if (typep m 'short-float) |
---|
912 | (let ((expon (- (%short-float-exp m) IEEE-single-float-bias))) |
---|
913 | (cond ((> expon 126) |
---|
914 | (setq log-s double-float-log2^23) |
---|
915 | (setq s #.(ash 1 23))) |
---|
916 | ((< expon -124) |
---|
917 | (setq log-s #.(- double-float-log2^23)) |
---|
918 | (setq s #.(/ 1.0s0 (ash 1 23)))))) |
---|
919 | (let ((expon (- (%double-float-exp m) IEEE-double-float-bias))) |
---|
920 | (cond ((> expon 1022) |
---|
921 | (setq log-s double-float-log2^23) |
---|
922 | (setq s #.(ash 1 23))) |
---|
923 | ((< expon -1020) |
---|
924 | (setq log-s #.(- double-float-log2^23)) |
---|
925 | (setq s #.(/ 1.0d0 (ash 1 23))))))) |
---|
926 | (if (eql s 1) |
---|
927 | (complex (log-abs x) (phase x)) |
---|
928 | (let ((temp (log-abs (/ x s)))) |
---|
929 | (complex (float (+ log-s temp) temp) (phase x))))))) |
---|
930 | (t |
---|
931 | (%rational-log x 1.0s0)))) |
---|
932 | |
---|
933 | ;;; helpers for rational log |
---|
934 | (defun %rational-log (x prototype) |
---|
935 | (cond ((minusp x) |
---|
936 | (complex (%rational-log (- x) prototype) |
---|
937 | (if (typep prototype 'short-float) |
---|
938 | single-float-pi |
---|
939 | pi))) |
---|
940 | ((bignump x) |
---|
941 | ;(let* ((base1 3) |
---|
942 | ; (guess (floor (1- (integer-length x)) |
---|
943 | ; (log base1 2))) |
---|
944 | ; (guess1 (* guess (log-e base1)))) |
---|
945 | ; (+ guess1 (log-e (/ x (expt base1 guess))))) |
---|
946 | ; Using base1 = 2 is *much* faster. Was there a reason for 3? |
---|
947 | (let* ((guess (1- (integer-length x))) |
---|
948 | (guess1 (* guess double-float-log2))) |
---|
949 | (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype))) |
---|
950 | ((and (ratiop x) |
---|
951 | ;; Rational arguments near +1 can be specified with great precision: don't lose it |
---|
952 | (cond ((< 0.5 x 3) |
---|
953 | (log1+ (float (- x 1) prototype))) |
---|
954 | ( |
---|
955 | ;; Escape out small values as well as large |
---|
956 | (or (> x most-positive-short-float) |
---|
957 | (< x least-positive-normalized-short-float)) |
---|
958 | ;; avoid loss of precision due to subtracting logs of numerator and denominator |
---|
959 | (let* ((n (%numerator x)) |
---|
960 | (d (%denominator x)) |
---|
961 | (sn (1- (integer-length n))) |
---|
962 | (sd (1- (integer-length d)))) |
---|
963 | (float (+ (* (- sn sd) double-float-log2) |
---|
964 | (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0)) |
---|
965 | (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0)))) |
---|
966 | prototype)))))) |
---|
967 | ((typep prototype 'short-float) |
---|
968 | #+32-bit-target |
---|
969 | (target::with-stack-short-floats ((sx x)) |
---|
970 | (%single-float-log! sx (%make-sfloat))) |
---|
971 | #+64-bit-target |
---|
972 | (%single-float-log (%short-float x))) |
---|
973 | (t |
---|
974 | (with-stack-double-floats ((dx x)) |
---|
975 | (%double-float-log! dx (%make-dfloat)))))) |
---|
976 | |
---|
977 | ;;; (log1+ x) = (log (1+ x)) |
---|
978 | ;;; but is much more precise near x = 0 |
---|
979 | (defun log1+ (x) |
---|
980 | ;;(cond ((typep x 'complex) |
---|
981 | ;; (let ((r (realpart x)) |
---|
982 | ;; (i (imagpart x))) |
---|
983 | ;; (if (and (< (abs r) 0.5) |
---|
984 | ;; (< (abs i) 3)) |
---|
985 | ;; (let* ((n (+ (* r (+ 2 r)) (* i i))) |
---|
986 | ;; (d (1+ (sqrt (1+ n))))) |
---|
987 | ;; (complex (log1+ (/ n d)) (atan i (1+ r)))) |
---|
988 | ;; (log (1+ x))))) |
---|
989 | ;; (t |
---|
990 | (if (and (typep x 'ratio) |
---|
991 | (< -0.5 x 2)) |
---|
992 | (setq x (%short-float x))) |
---|
993 | (let ((y (1+ x))) |
---|
994 | (if (eql y x) |
---|
995 | (log-e y) |
---|
996 | (let ((e (1- y))) |
---|
997 | (if (zerop e) |
---|
998 | (* x 1.0) |
---|
999 | (- (log-e y) (/ (- e x) y))))))) |
---|
1000 | |
---|
1001 | ;;; helper for complex log |
---|
1002 | ;;; uses more careful approach when (abs x) is near 1 |
---|
1003 | (defun log-abs (x) |
---|
1004 | (let ((a (abs x))) |
---|
1005 | (if (< 0.5 a 3) |
---|
1006 | (let* ((r (realpart x)) |
---|
1007 | (i (imagpart x)) |
---|
1008 | (n (if (> (abs r) (abs i)) |
---|
1009 | (+ (* (1+ r) (1- r)) (* i i)) |
---|
1010 | (+ (* r r) (* (1+ i) (1- i)))))) |
---|
1011 | (log1+ (/ n (1+ a)))) |
---|
1012 | (log-e a)))) |
---|
1013 | |
---|
1014 | (defun %rational-complex-log (x prototype &aux ra ia) |
---|
1015 | (let* ((rx (realpart x)) |
---|
1016 | (ix (imagpart x)) |
---|
1017 | (x (abs rx)) |
---|
1018 | (y (abs ix))) |
---|
1019 | (if (> y x) |
---|
1020 | (let ((r (float (/ rx y) 1.0d0))) |
---|
1021 | (setq ra (+ (%rational-log y 1.0d0) |
---|
1022 | (/ (log1+ (* r r)) 2))) |
---|
1023 | (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r))) |
---|
1024 | (let ((r (float (/ ix x) 1.0d0))) |
---|
1025 | (setq ra (+ (%rational-log x 1.0d0) |
---|
1026 | (/ (log1+ (* r r)) 2))) |
---|
1027 | (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0))))) |
---|
1028 | (if (typep prototype 'short-float) |
---|
1029 | (complex (%short-float ra) (%short-float ia)) |
---|
1030 | (complex ra ia)))) |
---|
1031 | |
---|
1032 | (defun exp (x) |
---|
1033 | "Return e raised to the power NUMBER." |
---|
1034 | (typecase x |
---|
1035 | (complex (* (exp (realpart x)) (cis (imagpart x)))) |
---|
1036 | (double-float (%double-float-exp! x (%make-dfloat))) |
---|
1037 | (t |
---|
1038 | #+32-bit-target |
---|
1039 | (target::with-stack-short-floats ((sx x)) |
---|
1040 | (%single-float-exp! sx (%make-sfloat))) |
---|
1041 | #+64-bit-target |
---|
1042 | (%single-float-exp (%short-float x))))) |
---|
1043 | |
---|
1044 | |
---|
1045 | (defun positive-realpart-p (n) |
---|
1046 | (> (realpart n) 0)) |
---|
1047 | |
---|
1048 | (defun expt (b e) |
---|
1049 | "Return BASE raised to the POWER." |
---|
1050 | (cond ((zerop e) (1+ (* b e))) |
---|
1051 | ((integerp e) |
---|
1052 | (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e))) |
---|
1053 | ((zerop b) |
---|
1054 | (if (plusp (realpart e)) b (report-bad-arg e '(satisfies positive-realpart-p)))) |
---|
1055 | ((and (realp b) (plusp b) (realp e)) |
---|
1056 | (if (or (typep b 'double-float) |
---|
1057 | (typep e 'double-float)) |
---|
1058 | (with-stack-double-floats ((b1 b) |
---|
1059 | (e1 e)) |
---|
1060 | (%double-float-expt! b1 e1 (%make-dfloat))) |
---|
1061 | #+32-bit-target |
---|
1062 | (target::with-stack-short-floats ((b1 b) |
---|
1063 | (e1 e)) |
---|
1064 | (%single-float-expt! b1 e1 (%make-sfloat))) |
---|
1065 | #+64-bit-target |
---|
1066 | (%single-float-expt (%short-float b) (%short-float e)) |
---|
1067 | )) |
---|
1068 | ((typep (realpart e) 'double-float) |
---|
1069 | ;; Avoid intermediate single-float result from LOG |
---|
1070 | (let ((promoted-base (* 1d0 b))) |
---|
1071 | (exp (* e (log promoted-base))))) |
---|
1072 | (t (exp (* e (log b)))))) |
---|
1073 | |
---|
1074 | |
---|
1075 | |
---|
1076 | (defun sqrt (x &aux a b) |
---|
1077 | "Return the square root of NUMBER." |
---|
1078 | (cond ((zerop x) x) |
---|
1079 | ((complexp x) |
---|
1080 | (let ((rx (realpart x)) |
---|
1081 | (ix (imagpart x))) |
---|
1082 | (cond ((rationalp rx) |
---|
1083 | (if (zerop rx) |
---|
1084 | (let ((s (sqrt (/ (abs ix) 2)))) |
---|
1085 | (complex s (if (minusp ix) (- s) s))) |
---|
1086 | (let* ((s (+ (* rx rx) (* ix ix))) |
---|
1087 | (d (if (ratiop s) |
---|
1088 | (/ (isqrt (%numerator s)) |
---|
1089 | (isqrt (%denominator s))) |
---|
1090 | (isqrt s)))) |
---|
1091 | (unless (eql s (* d d)) |
---|
1092 | (setf d (%double-float-hypot (%double-float rx) |
---|
1093 | (%double-float ix)))) |
---|
1094 | (cond ((minusp rx) |
---|
1095 | (setq b (sqrt (/ (- d rx) 2))) |
---|
1096 | (when (minusp ix) |
---|
1097 | (setq b (- b))) |
---|
1098 | (setq a (/ ix (+ b b)))) |
---|
1099 | (t |
---|
1100 | (setq a (sqrt (/ (+ rx d) 2))) |
---|
1101 | (setq b (/ ix (+ a a))))) |
---|
1102 | (if (rationalp a) |
---|
1103 | (complex a b) |
---|
1104 | (complex (%short-float a) (%short-float b)))))) |
---|
1105 | ((minusp rx) |
---|
1106 | (if (zerop ix) |
---|
1107 | (complex 0 (float-sign ix (sqrt (- rx)))) |
---|
1108 | (let ((shift (cond ((< rx -1) -3) |
---|
1109 | ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25) |
---|
1110 | (t -1)))) |
---|
1111 | (setq rx (scale-float rx shift)) |
---|
1112 | (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx)))) |
---|
1113 | (setq b (scale-float s (ash (- -1 shift) -1))) |
---|
1114 | (when (minusp ix) |
---|
1115 | (setq b (- b))) |
---|
1116 | (setq a (/ ix (scale-float b 1))) |
---|
1117 | (complex a b))))) |
---|
1118 | (t |
---|
1119 | (if (zerop ix) |
---|
1120 | (complex (sqrt rx) ix) |
---|
1121 | (let ((shift (cond ((> rx 1) -3) |
---|
1122 | ((and (< rx 5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25) |
---|
1123 | (t -1)))) |
---|
1124 | (setq rx (scale-float rx shift)) |
---|
1125 | (let ((s (fsqrt (+ rx (abs (complex rx (scale-float ix shift))))))) |
---|
1126 | (setq a (scale-float s (ash (- -1 shift) -1))) |
---|
1127 | (setq b (/ ix (scale-float a 1))) |
---|
1128 | (complex a b)))))))) |
---|
1129 | ((minusp x) (complex 0 (sqrt (- x)))) |
---|
1130 | ((floatp x) |
---|
1131 | (fsqrt x)) |
---|
1132 | ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a) |
---|
1133 | ((and (ratiop x) |
---|
1134 | (let ((n (numerator x)) |
---|
1135 | d) |
---|
1136 | (and (eql n (* (setq a (isqrt n)) a)) |
---|
1137 | (eql (setq d (denominator x)) |
---|
1138 | (* (setq b (isqrt d)) b))))) |
---|
1139 | (/ a b)) |
---|
1140 | (t |
---|
1141 | (float (fsqrt (float x 0.0d0)) 1.0s0)))) |
---|
1142 | |
---|
1143 | |
---|
1144 | |
---|
1145 | (defun asin (x) |
---|
1146 | "Return the arc sine of NUMBER." |
---|
1147 | (number-case x |
---|
1148 | (complex |
---|
1149 | (let ((sqrt-1-x (sqrt (- 1 x))) |
---|
1150 | (sqrt-1+x (sqrt (+ 1 x)))) |
---|
1151 | (complex (atan (realpart x) |
---|
1152 | (realpart (* sqrt-1-x sqrt-1+x))) |
---|
1153 | (asinh (imagpart (* (conjugate sqrt-1-x) |
---|
1154 | sqrt-1+x)))))) |
---|
1155 | (double-float |
---|
1156 | (locally (declare (type double-float x)) |
---|
1157 | (if (and (<= -1.0d0 x) |
---|
1158 | (<= x 1.0d0)) |
---|
1159 | (%double-float-asin! x (%make-dfloat)) |
---|
1160 | (let* ((temp (+ (complex -0.0d0 x) |
---|
1161 | (sqrt (- 1.0d0 (the double-float (* x x))))))) |
---|
1162 | (complex (phase temp) (- (log (abs temp)))))))) |
---|
1163 | ((short-float rational) |
---|
1164 | #+32-bit-target |
---|
1165 | (let* ((x1 (%make-sfloat))) |
---|
1166 | (declare (dynamic-extent x1)) |
---|
1167 | (if (and (realp x) |
---|
1168 | (<= -1.0s0 (setq x (%short-float x x1))) |
---|
1169 | (<= x 1.0s0)) |
---|
1170 | (%single-float-asin! x1 (%make-sfloat)) |
---|
1171 | (progn |
---|
1172 | (setq x (+ (complex (- (imagpart x)) (realpart x)) |
---|
1173 | (sqrt (- 1 (* x x))))) |
---|
1174 | (complex (phase x) (- (log (abs x))))))) |
---|
1175 | #+64-bit-target |
---|
1176 | (if (and (realp x) |
---|
1177 | (<= -1.0s0 (setq x (%short-float x))) |
---|
1178 | (<= x 1.0s0)) |
---|
1179 | (%single-float-asin x) |
---|
1180 | (progn |
---|
1181 | (setq x (+ (complex (- (imagpart x)) (realpart x)) |
---|
1182 | (sqrt (- 1 (* x x))))) |
---|
1183 | (complex (phase x) (- (log (abs x)))))) |
---|
1184 | ))) |
---|
1185 | |
---|
1186 | |
---|
1187 | (defun acos (x) |
---|
1188 | "Return the arc cosine of NUMBER." |
---|
1189 | (number-case x |
---|
1190 | (complex |
---|
1191 | (if (typep (realpart x) 'double-float) |
---|
1192 | (- double-float-half-pi (asin x)) |
---|
1193 | (- single-float-half-pi (asin x)))) |
---|
1194 | (double-float |
---|
1195 | (locally (declare (type double-float x)) |
---|
1196 | (if (and (<= -1.0d0 x) |
---|
1197 | (<= x 1.0d0)) |
---|
1198 | (%double-float-acos! x (%make-dfloat)) |
---|
1199 | (- double-float-half-pi (asin x))))) |
---|
1200 | ((short-float rational) |
---|
1201 | #+32-bit-target |
---|
1202 | (target::with-stack-short-floats ((sx x)) |
---|
1203 | (locally |
---|
1204 | (declare (type short-float sx)) |
---|
1205 | (if (and (<= -1.0s0 sx) |
---|
1206 | (<= sx 1.0s0)) |
---|
1207 | (%single-float-acos! sx (%make-sfloat)) |
---|
1208 | (- single-float-half-pi (asin sx))))) |
---|
1209 | #+64-bit-target |
---|
1210 | (let* ((sx (%short-float x))) |
---|
1211 | (declare (type short-float sx)) |
---|
1212 | (if (and (<= -1.0s0 sx) |
---|
1213 | (<= sx 1.0s0)) |
---|
1214 | (%single-float-acos sx) |
---|
1215 | (- single-float-half-pi (asin sx)))) |
---|
1216 | ))) |
---|
1217 | |
---|
1218 | |
---|
1219 | (defun fsqrt (x) |
---|
1220 | (etypecase x |
---|
1221 | (double-float (%double-float-sqrt! x (%make-dfloat))) |
---|
1222 | (single-float |
---|
1223 | #+32-bit-target |
---|
1224 | (%single-float-sqrt! x (%make-sfloat)) |
---|
1225 | #+64-bit-target |
---|
1226 | (%single-float-sqrt x)))) |
---|
1227 | |
---|
1228 | |
---|
1229 | |
---|
1230 | (defun %df-atan2 (y x) |
---|
1231 | (if (zerop x) |
---|
1232 | (if (zerop y) |
---|
1233 | (if (plusp (float-sign x)) |
---|
1234 | (if (eql y -0.0d0) |
---|
1235 | -0.0d0 |
---|
1236 | 0.0d0) |
---|
1237 | (float-sign y pi)) |
---|
1238 | (float-sign y double-float-half-pi)) |
---|
1239 | (%double-float-atan2! y x (%make-dfloat)))) |
---|
1240 | |
---|
1241 | #+32-bit-target |
---|
1242 | (defun %sf-atan2! (y x) |
---|
1243 | (if (zerop x) |
---|
1244 | (if (zerop y) |
---|
1245 | (if (plusp (float-sign x)) |
---|
1246 | ;; Don't return Y (which may be stack-consed) here. |
---|
1247 | ;; We know that (ZEROP Y) is true, so: |
---|
1248 | (if (eql y -0.0s0) |
---|
1249 | -0.0s0 |
---|
1250 | 0.0s0) |
---|
1251 | (float-sign y single-float-pi)) |
---|
1252 | (float-sign y single-float-half-pi)) |
---|
1253 | (%single-float-atan2! y x (%make-sfloat)))) |
---|
1254 | |
---|
1255 | #+64-bit-target |
---|
1256 | (defun %sf-atan2 (y x) |
---|
1257 | (if (zerop x) |
---|
1258 | (if (zerop y) |
---|
1259 | (if (plusp (float-sign x)) |
---|
1260 | y |
---|
1261 | (float-sign y single-float-pi)) |
---|
1262 | (float-sign y single-float-half-pi)) |
---|
1263 | (%single-float-atan2 y x))) |
---|
1264 | |
---|
1265 | #+64-bit-target |
---|
1266 | (defun %short-float-exp (n) |
---|
1267 | (let* ((bits (single-float-bits n))) |
---|
1268 | (declare (type (unsigned-byte 32) bits)) |
---|
1269 | (ldb (byte IEEE-single-float-exponent-width IEEE-single-float-exponent-offset) bits))) |
---|
1270 | |
---|
1271 | |
---|
1272 | #+64-bit-target |
---|
1273 | (defun set-%short-float-exp (float exp) |
---|
1274 | (host-single-float-from-unsigned-byte-32 |
---|
1275 | (dpb exp |
---|
1276 | (byte IEEE-single-float-exponent-width |
---|
1277 | IEEE-single-float-exponent-offset) |
---|
1278 | (the (unsigned-byte 32) (single-float-bits float))))) |
---|
1279 | |
---|
1280 | #+64-bit-target |
---|
1281 | (defun %%scale-sfloat (float int) |
---|
1282 | (* (the single-float float) |
---|
1283 | (the single-float (host-single-float-from-unsigned-byte-32 |
---|
1284 | (dpb int |
---|
1285 | (byte IEEE-single-float-exponent-width |
---|
1286 | IEEE-single-float-exponent-offset) |
---|
1287 | 0))))) |
---|
1288 | |
---|
1289 | #+64-bit-target |
---|
1290 | (defun %double-float-exp (n) |
---|
1291 | (let* ((highword (double-float-bits n))) |
---|
1292 | (declare (fixnum highword)) |
---|
1293 | (logand (1- (ash 1 IEEE-double-float-exponent-width)) |
---|
1294 | (ash highword (- (- IEEE-double-float-exponent-offset 32)))))) |
---|
1295 | |
---|
1296 | #+64-bit-target |
---|
1297 | (defun set-%double-float-exp (float exp) |
---|
1298 | (let* ((highword (double-float-bits float))) |
---|
1299 | (declare (fixnum highword)) |
---|
1300 | (setf (uvref float target::double-float.val-high-cell) |
---|
1301 | (dpb exp |
---|
1302 | (byte IEEE-double-float-exponent-width |
---|
1303 | (- IEEE-double-float-exponent-offset 32)) |
---|
1304 | highword)) |
---|
1305 | exp)) |
---|
1306 | |
---|
1307 | #+64-bit-target |
---|
1308 | (defun %integer-decode-double-float (f) |
---|
1309 | (multiple-value-bind (hiword loword) (double-float-bits f) |
---|
1310 | (declare (type (unsigned-byte 32) hiword loword)) |
---|
1311 | (let* ((exp (ldb (byte IEEE-double-float-exponent-width |
---|
1312 | (- IEEE-double-float-exponent-offset 32)) |
---|
1313 | hiword)) |
---|
1314 | (mantissa (logior |
---|
1315 | (the fixnum |
---|
1316 | (dpb (ldb (byte (- IEEE-double-float-mantissa-width 32) |
---|
1317 | IEEE-double-float-mantissa-offset) |
---|
1318 | hiword) |
---|
1319 | (byte (- IEEE-double-float-mantissa-width 32) |
---|
1320 | 32) |
---|
1321 | loword)) |
---|
1322 | (if (zerop exp) |
---|
1323 | 0 |
---|
1324 | (ash 1 IEEE-double-float-hidden-bit)))) |
---|
1325 | (sign (if (logbitp 31 hiword) -1 1))) |
---|
1326 | (declare (fixnum exp mantissa sign)) |
---|
1327 | (values (ldb (byte 25 28) mantissa) |
---|
1328 | (ldb (byte 28 0) mantissa) |
---|
1329 | exp |
---|
1330 | sign)))) |
---|
1331 | |
---|
1332 | ;;; end of l0-float.lisp |
---|