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 | ) |
---|
26 | |
---|
27 | ;;; used by float reader |
---|
28 | (defun make-float-from-fixnums (hi lo exp sign &optional result) |
---|
29 | ;;(require-null-or-double-float-sym result) |
---|
30 | ;; maybe nuke all these require-types? |
---|
31 | ;;(setq hi (require-type hi 'fixnum)) |
---|
32 | ;;(setq lo (require-type lo 'fixnum)) |
---|
33 | ;;(setq exp (require-type exp 'fixnum)) |
---|
34 | ;;(setq sign (require-type sign 'fixnum)) |
---|
35 | (let ((the-float (or result (%make-dfloat)))) |
---|
36 | (%make-float-from-fixnums the-float hi lo exp sign) |
---|
37 | the-float)) |
---|
38 | |
---|
39 | |
---|
40 | #+32-bit-target |
---|
41 | (defun make-short-float-from-fixnums (significand biased-exp sign &optional result) |
---|
42 | (%make-short-float-from-fixnums (or result (%make-sfloat)) significand biased-exp sign)) |
---|
43 | |
---|
44 | #+64-bit-target |
---|
45 | (defun make-short-float-from-fixnums (significand biased-exp sign) |
---|
46 | (declare (fixnum significand biased-exp sign)) |
---|
47 | (host-single-float-from-unsigned-byte-32 |
---|
48 | (logior |
---|
49 | (the fixnum (if (< sign 0) (ash 1 31) 0)) |
---|
50 | (the fixnum (ash biased-exp IEEE-single-float-exponent-offset)) |
---|
51 | (the fixnum (logand significand |
---|
52 | (1- (ash 1 IEEE-single-float-hidden-bit))))))) |
---|
53 | |
---|
54 | |
---|
55 | (defun float-sign (n1 &optional n2) ; second arg silly |
---|
56 | "Return a floating-point number that has the same sign as |
---|
57 | FLOAT1 and, if FLOAT2 is given, has the same absolute value |
---|
58 | as FLOAT2." |
---|
59 | (if (and n2 (not (typep n2 'float))) |
---|
60 | (setq n2 (require-type n2 'float))) |
---|
61 | (number-case n1 |
---|
62 | (double-float |
---|
63 | (if (%double-float-sign n1) |
---|
64 | (if n2 |
---|
65 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2)) |
---|
66 | -1.0d0) |
---|
67 | (if n2 |
---|
68 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2) |
---|
69 | 1.0d0))) |
---|
70 | (short-float |
---|
71 | (if (%short-float-sign n1) |
---|
72 | (if n2 |
---|
73 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) n2 (- n2)) |
---|
74 | -1.0s0) |
---|
75 | (if n2 |
---|
76 | (if (if (typep n2 'double-float) (%double-float-minusp n2) (%short-float-minusp n2)) (- n2) n2) |
---|
77 | 1.0s0))))) |
---|
78 | |
---|
79 | |
---|
80 | |
---|
81 | (defun %double-float-minusp (n) |
---|
82 | (and (%double-float-sign n)(not (%double-float-zerop n)))) |
---|
83 | |
---|
84 | (defun %short-float-minusp (n) |
---|
85 | (and (%short-float-sign n) (not (%short-float-zerop n)))) |
---|
86 | |
---|
87 | (defun %double-float-abs (n) |
---|
88 | (if (not (%double-float-sign n)) |
---|
89 | n |
---|
90 | (%%double-float-abs! n (%make-dfloat)))) |
---|
91 | |
---|
92 | #+32-bit-target |
---|
93 | (defun %short-float-abs (n) |
---|
94 | (if (not (%short-float-sign n)) |
---|
95 | n |
---|
96 | (%%short-float-abs! n (%make-sfloat)))) |
---|
97 | |
---|
98 | (defun fixnum-decode-float (n) |
---|
99 | (etypecase n |
---|
100 | (double-float (%integer-decode-double-float n)))) |
---|
101 | |
---|
102 | (defun nan-or-infinity-p (n) |
---|
103 | (etypecase n |
---|
104 | (double-float (eq 2047 (%double-float-exp n))) |
---|
105 | (short-float (eq 255 (%short-float-exp n))))) |
---|
106 | |
---|
107 | ; not sure this is right |
---|
108 | (defun infinity-p (n) |
---|
109 | (etypecase n |
---|
110 | (double-float (multiple-value-bind (hi lo exp)(fixnum-decode-float n) |
---|
111 | (and (eq 2047 exp) |
---|
112 | (eq #x1000000 hi) |
---|
113 | (eq 0 lo)))) |
---|
114 | (short-float |
---|
115 | #+32-bit-target |
---|
116 | (multiple-value-bind (high low)(%sfloat-hwords n) |
---|
117 | (let* ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F)))) |
---|
118 | (exp (%ilsr 7 (%ilogand2 high #x7F80)))) |
---|
119 | (and (eq exp 255) |
---|
120 | (eq 0 mantissa)))) |
---|
121 | #+64-bit-target |
---|
122 | (let* ((bits (single-float-bits n)) |
---|
123 | (exp (ldb (byte IEEE-single-float-exponent-width |
---|
124 | IEEE-single-float-exponent-offset) |
---|
125 | bits)) |
---|
126 | (mantissa (ldb (byte IEEE-single-float-mantissa-width |
---|
127 | IEEE-single-float-mantissa-offset) |
---|
128 | bits))) |
---|
129 | (declare (fixnum bits exp mantissa)) |
---|
130 | (and (= exp 255) |
---|
131 | (zerop mantissa)))))) |
---|
132 | |
---|
133 | #+32-bit-target |
---|
134 | (defun fixnum-decode-short-float (float) |
---|
135 | (multiple-value-bind (high low)(%sfloat-hwords float) |
---|
136 | (let* ((mantissa (%ilogior2 low (%ilsl 16 (%ilogand2 high #x007F)))) |
---|
137 | (exp (%ilsr 7 (%ilogand2 high #x7F80)))) |
---|
138 | (if (and (neq exp 0)(neq exp 255))(setq mantissa (%ilogior mantissa #x800000))) |
---|
139 | (values mantissa exp (%ilsr 15 high))))) |
---|
140 | |
---|
141 | #+64-bit-target |
---|
142 | (defun fixnum-decode-short-float (float) |
---|
143 | (let* ((bits (single-float-bits float))) |
---|
144 | (declare (fixnum bits)) |
---|
145 | (let* ((mantissa (ldb (byte IEEE-single-float-mantissa-width |
---|
146 | IEEE-single-float-mantissa-offset) |
---|
147 | bits)) |
---|
148 | (exp (ldb (byte IEEE-single-float-exponent-width |
---|
149 | IEEE-single-float-exponent-offset) |
---|
150 | bits)) |
---|
151 | (sign (lsh bits -31))) |
---|
152 | (declare (fixnum mantissa exp sign)) |
---|
153 | (unless (or (= exp 0) (= exp 255)) |
---|
154 | (setq mantissa (logior mantissa (ash 1 IEEE-single-float-hidden-bit)))) |
---|
155 | (values mantissa exp sign)))) |
---|
156 | |
---|
157 | |
---|
158 | |
---|
159 | #+32-bit-target |
---|
160 | (defun integer-decode-double-float (n) |
---|
161 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
162 | ; is only 53 bits and positive so should be easy |
---|
163 | ;(values (logior (ash hi 28) lo) exp sign))) |
---|
164 | ; if denormalized, may fit in a fixnum |
---|
165 | (setq exp (- exp (if (< hi #x1000000) |
---|
166 | (+ IEEE-double-float-mantissa-width IEEE-double-float-bias) |
---|
167 | (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))) |
---|
168 | (if (< hi (ash 1 (1- target::fixnumshift))) ; aka 2 |
---|
169 | (values (logior (ash hi 28) lo) exp sign) |
---|
170 | ; might fit in 1 word? |
---|
171 | (let ((big (%alloc-misc 2 target::subtag-bignum))) |
---|
172 | (make-big-53 hi lo big) |
---|
173 | (if (< hi #x1000000) (%normalize-bignum big)) |
---|
174 | (values big exp sign))))) |
---|
175 | |
---|
176 | #+64-bit-target |
---|
177 | (defun integer-decode-double-float (n) |
---|
178 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
179 | (setq exp (- exp (if (< hi #x1000000) |
---|
180 | (+ IEEE-double-float-mantissa-width IEEE-double-float-bias) |
---|
181 | (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias))))) |
---|
182 | (values (logior (ash hi 28) lo) exp sign))) |
---|
183 | |
---|
184 | |
---|
185 | ;;; actually only called when magnitude bigger than a fixnum |
---|
186 | #+32-bit-target |
---|
187 | (defun %truncate-double-float (n) |
---|
188 | (multiple-value-bind (hi lo exp sign)(%integer-decode-double-float n) |
---|
189 | (if (< exp (1+ IEEE-double-float-bias)) ; this is false in practice |
---|
190 | 0 |
---|
191 | (progn |
---|
192 | (setq exp (- exp (+ IEEE-double-float-mantissa-width (1+ IEEE-double-float-bias)))) |
---|
193 | (if (eq sign 1) ; positive |
---|
194 | (logior (ash hi (+ 28 exp))(ash lo exp)) |
---|
195 | (if (<= exp 0) ; exp positive - negate before shift - else after |
---|
196 | (let ((poo (logior (ash hi (+ 28 exp))(ash lo exp)))) |
---|
197 | (- poo)) |
---|
198 | (let ((poo (logior (ash hi 28) lo))) |
---|
199 | (ash (- poo) exp)))))))) |
---|
200 | |
---|
201 | #+64-bit-target |
---|
202 | (defun %truncate-double-float (n) |
---|
203 | (multiple-value-bind (mantissa exp sign) (integer-decode-float n) |
---|
204 | (* sign (ash mantissa exp)))) |
---|
205 | |
---|
206 | |
---|
207 | |
---|
208 | ; actually only called when bigger than a fixnum |
---|
209 | (defun %truncate-short-float (n) |
---|
210 | (multiple-value-bind (mantissa exp sign)(fixnum-decode-short-float n) |
---|
211 | (if (< exp (1+ IEEE-single-float-bias)) ; is magnitude less than 1 - false in practice |
---|
212 | 0 |
---|
213 | (progn |
---|
214 | (setq exp (- exp (+ IEEE-single-float-mantissa-width (1+ IEEE-single-float-bias)))) |
---|
215 | (ash (if (eq sign 0) mantissa (- mantissa)) exp))))) |
---|
216 | |
---|
217 | (defun decode-float (n) |
---|
218 | "Return three values: |
---|
219 | 1) a floating-point number representing the significand. This is always |
---|
220 | between 0.5 (inclusive) and 1.0 (exclusive). |
---|
221 | 2) an integer representing the exponent. |
---|
222 | 3) -1.0 or 1.0 (i.e. the sign of the argument.)" |
---|
223 | (number-case n |
---|
224 | (double-float |
---|
225 | (let* ((old-exp (%double-float-exp n)) |
---|
226 | (sign (if (%double-float-sign n) -1.0d0 1.0d0))) |
---|
227 | (if (eq 0 old-exp) |
---|
228 | (if (%double-float-zerop n) |
---|
229 | (values 0.0d0 0 sign) |
---|
230 | (let* ((val (%make-dfloat)) |
---|
231 | (zeros (dfloat-significand-zeros n))) |
---|
232 | (%%double-float-abs! n val) |
---|
233 | (%%scale-dfloat! val (+ 2 IEEE-double-float-bias zeros) val) ; get it normalized |
---|
234 | (set-%double-float-exp val IEEE-double-float-bias) ; then bash exponent |
---|
235 | (values val (- old-exp zeros IEEE-double-float-bias) sign))) |
---|
236 | (if (> old-exp IEEE-double-float-normal-exponent-max) |
---|
237 | (error "Can't decode NAN or infinity ~s" n) |
---|
238 | (let ((val (%make-dfloat))) |
---|
239 | (%%double-float-abs! n val) |
---|
240 | (set-%double-float-exp val IEEE-double-float-bias) |
---|
241 | (values val (- old-exp IEEE-double-float-bias) sign)))))) |
---|
242 | (short-float |
---|
243 | (let* ((old-exp (%short-float-exp n)) |
---|
244 | (sign (if (%short-float-sign n) -1.0s0 1.0s0))) |
---|
245 | (if (eq 0 old-exp) |
---|
246 | (if (%short-float-zerop n) |
---|
247 | (values 0.0s0 0 sign) |
---|
248 | #+32-bit-target |
---|
249 | (let* ((val (%make-sfloat)) |
---|
250 | (zeros (sfloat-significand-zeros n))) |
---|
251 | (%%short-float-abs! n val) |
---|
252 | (%%scale-sfloat! val (+ 2 IEEE-single-float-bias zeros) val) ; get it normalized |
---|
253 | (set-%short-float-exp val IEEE-single-float-bias) ; then bash exponent |
---|
254 | (values val (- old-exp zeros IEEE-single-float-bias) sign)) |
---|
255 | #+64-bit-target |
---|
256 | (let* ((zeros (sfloat-significand-zeros n)) |
---|
257 | (val (%%scale-sfloat (%short-float-abs n) |
---|
258 | (+ 2 IEEE-single-float-bias zeros)))) |
---|
259 | (values (set-%short-float-exp val IEEE-single-float-bias) |
---|
260 | (- old-exp zeros IEEE-single-float-bias) sign))) |
---|
261 | (if (> old-exp IEEE-single-float-normal-exponent-max) |
---|
262 | (error "Can't decode NAN or infinity ~s" n) |
---|
263 | #+32-bit-target |
---|
264 | (let ((val (%make-sfloat))) |
---|
265 | (%%short-float-abs! n val) |
---|
266 | (set-%short-float-exp val IEEE-single-float-bias) |
---|
267 | (values val (- old-exp IEEE-single-float-bias) sign)) |
---|
268 | #+64-bit-target |
---|
269 | (values (set-%short-float-exp (%short-float-abs n) |
---|
270 | IEEE-single-float-bias) |
---|
271 | (- old-exp IEEE-single-float-bias) sign))))))) |
---|
272 | |
---|
273 | ; (* float (expt 2 int)) |
---|
274 | (defun scale-float (float int) |
---|
275 | "Return the value (* f (expt (float 2 f) ex)), but with no unnecessary loss |
---|
276 | of precision or overflow." |
---|
277 | (unless (fixnump int)(setq int (require-type int 'fixnum))) |
---|
278 | (number-case float |
---|
279 | (double-float |
---|
280 | (let* ((float-exp (%double-float-exp float)) |
---|
281 | (new-exp (+ float-exp int))) |
---|
282 | (if (eq 0 float-exp) ; already denormalized? |
---|
283 | (if (%double-float-zerop float) |
---|
284 | float |
---|
285 | (let ((result (%make-dfloat))) |
---|
286 | (%%scale-dfloat! float (+ (1+ IEEE-double-float-bias) int) result))) |
---|
287 | (if (<= new-exp 0) ; maybe going denormalized |
---|
288 | (if (<= new-exp (- IEEE-double-float-digits)) |
---|
289 | 0.0d0 ; should this be underflow? - should just be normal and result is fn of current fpu-mode |
---|
290 | ;(error "Can't scale ~s by ~s." float int) ; should signal something |
---|
291 | (let ((result (%make-dfloat))) |
---|
292 | (%copy-double-float float result) |
---|
293 | (set-%double-float-exp result 1) ; scale by float-exp -1 |
---|
294 | (%%scale-dfloat! result (+ IEEE-double-float-bias (+ float-exp int)) result) |
---|
295 | result)) |
---|
296 | (if (> new-exp IEEE-double-float-normal-exponent-max) |
---|
297 | (error (make-condition 'floating-point-overflow |
---|
298 | :operation 'scale-float |
---|
299 | :operands (list float int))) |
---|
300 | (let ((new-float (%make-dfloat))) |
---|
301 | (%copy-double-float float new-float) |
---|
302 | (set-%double-float-exp new-float new-exp) |
---|
303 | new-float)))))) |
---|
304 | (short-float |
---|
305 | (let* ((float-exp (%short-float-exp float)) |
---|
306 | (new-exp (+ float-exp int))) |
---|
307 | (if (eq 0 float-exp) ; already denormalized? |
---|
308 | (if (%short-float-zerop float) |
---|
309 | float |
---|
310 | #+32-bit-target |
---|
311 | (let ((result (%make-sfloat))) |
---|
312 | (%%scale-sfloat! float (+ (1+ IEEE-single-float-bias) int) result)) |
---|
313 | #+64-bit-target |
---|
314 | (%%scale-sfloat float (+ (1+ IEEE-single-float-bias) int))) |
---|
315 | (if (<= new-exp 0) ; maybe going denormalized |
---|
316 | (if (<= new-exp (- IEEE-single-float-digits)) |
---|
317 | ;; should this be underflow? - should just be normal and |
---|
318 | ;; result is fn of current fpu-mode (error "Can't scale |
---|
319 | ;; ~s by ~s." float int) ; should signal something |
---|
320 | 0.0s0 |
---|
321 | #+32-bit-target |
---|
322 | (let ((result (%make-sfloat))) |
---|
323 | (%copy-short-float float result) |
---|
324 | (set-%short-float-exp result 1) ; scale by float-exp -1 |
---|
325 | (%%scale-sfloat! result (+ IEEE-single-float-bias (+ float-exp int)) result) |
---|
326 | result) |
---|
327 | #+64-bit-target |
---|
328 | (%%scale-sfloat (set-%short-float-exp float 1) |
---|
329 | (+ IEEE-single-float-bias (+ float-exp int)))) |
---|
330 | (if (> new-exp IEEE-single-float-normal-exponent-max) |
---|
331 | (error (make-condition 'floating-point-overflow |
---|
332 | :operation 'scale-float |
---|
333 | :operands (list float int))) |
---|
334 | #+32-bit-target |
---|
335 | (let ((new-float (%make-sfloat))) |
---|
336 | (%copy-short-float float new-float) |
---|
337 | (set-%short-float-exp new-float new-exp) |
---|
338 | new-float) |
---|
339 | #+64-bit-target |
---|
340 | (set-%short-float-exp float new-exp)))))))) |
---|
341 | |
---|
342 | (defun %copy-float (f) |
---|
343 | ;Returns a freshly consed float. float can also be a macptr. |
---|
344 | (cond ((double-float-p f) (%copy-double-float f (%make-dfloat))) |
---|
345 | ((macptrp f) |
---|
346 | (let ((float (%make-dfloat))) |
---|
347 | (%copy-ptr-to-ivector f 0 float (* 4 target::double-float.value-cell) 8) |
---|
348 | float)) |
---|
349 | (t (error "Illegal arg ~s to %copy-float" f)))) |
---|
350 | |
---|
351 | (defun float-precision (float) ; not used - not in cltl2 index ? |
---|
352 | "Return a non-negative number of significant digits in its float argument. |
---|
353 | Will be less than FLOAT-DIGITS if denormalized or zero." |
---|
354 | (number-case float |
---|
355 | (double-float |
---|
356 | (if (eq 0 (%double-float-exp float)) |
---|
357 | (if (not (%double-float-zerop float)) |
---|
358 | ; denormalized |
---|
359 | (- IEEE-double-float-mantissa-width (dfloat-significand-zeros float)) |
---|
360 | 0) |
---|
361 | IEEE-double-float-digits)) |
---|
362 | (short-float |
---|
363 | (if (eq 0 (%short-float-exp float)) |
---|
364 | (if (not (%short-float-zerop float)) |
---|
365 | ; denormalized |
---|
366 | (- IEEE-single-float-mantissa-width (sfloat-significand-zeros float)) |
---|
367 | 0) |
---|
368 | IEEE-single-float-digits)))) |
---|
369 | |
---|
370 | |
---|
371 | (defun %double-float (number &optional result) |
---|
372 | ;(require-null-or-double-float-sym result) |
---|
373 | ; use number-case when macro is common |
---|
374 | (number-case number |
---|
375 | (double-float |
---|
376 | (if result |
---|
377 | (%copy-double-float number result) |
---|
378 | number)) |
---|
379 | (short-float |
---|
380 | (%short-float->double-float number (or result (%make-dfloat)))) |
---|
381 | (fixnum |
---|
382 | (%fixnum-dfloat number (or result (%make-dfloat)))) |
---|
383 | (bignum (%bignum-dfloat number result)) |
---|
384 | (ratio |
---|
385 | (if (not result)(setq result (%make-dfloat))) |
---|
386 | (let* ((num (%numerator number)) |
---|
387 | (den (%denominator number))) |
---|
388 | ; dont error if result is floatable when either top or bottom is not. |
---|
389 | ; maybe do usual first, catching error |
---|
390 | (if (not (or (bignump num)(bignump den))) |
---|
391 | (with-stack-double-floats ((fnum num) |
---|
392 | (fden den)) |
---|
393 | (%double-float/-2! fnum fden result)) |
---|
394 | (let* ((numlen (integer-length num)) |
---|
395 | (denlen (integer-length den)) |
---|
396 | (exp (- numlen denlen)) |
---|
397 | (minusp (minusp num))) |
---|
398 | (if (and (<= numlen IEEE-double-float-bias) |
---|
399 | (<= denlen IEEE-double-float-bias) |
---|
400 | #|(not (minusp exp))|# |
---|
401 | (<= (abs exp) IEEE-double-float-mantissa-width)) |
---|
402 | (with-stack-double-floats ((fnum num) |
---|
403 | (fden den)) |
---|
404 | |
---|
405 | (%double-float/-2! fnum fden result)) |
---|
406 | (if (> exp IEEE-double-float-mantissa-width) |
---|
407 | (progn (%double-float (round num den) result)) |
---|
408 | (if (>= exp 0) |
---|
409 | ; exp between 0 and 53 and nums big |
---|
410 | (let* ((shift (- IEEE-double-float-digits exp)) |
---|
411 | (num (if minusp (- num) num)) |
---|
412 | (int (round (ash num shift) den)) ; gaak |
---|
413 | (intlen (integer-length int)) |
---|
414 | (new-exp (+ intlen (- IEEE-double-float-bias shift)))) |
---|
415 | |
---|
416 | (when (> intlen IEEE-double-float-digits) |
---|
417 | (setq shift (1- shift)) |
---|
418 | (setq int (round (ash num shift) den)) |
---|
419 | (setq intlen (integer-length int)) |
---|
420 | (setq new-exp (+ intlen (- IEEE-double-float-bias shift)))) |
---|
421 | (when (> new-exp 2046) |
---|
422 | (error (make-condition 'floating-point-overflow |
---|
423 | :operation 'double-float |
---|
424 | :operands (list number)))) |
---|
425 | (make-float-from-fixnums (ldb (byte 25 (- intlen 25)) int) |
---|
426 | (ldb (byte 28 (max (- intlen 53) 0)) int) |
---|
427 | new-exp ;(+ intlen (- IEEE-double-float-bias 53)) |
---|
428 | (if minusp -1 1) |
---|
429 | result)) |
---|
430 | ; den > num - exp negative |
---|
431 | (progn |
---|
432 | (float-rat-neg-exp num den (if minusp -1 1) result))))))))))) |
---|
433 | |
---|
434 | |
---|
435 | #+32-bit-target |
---|
436 | (defun %short-float-ratio (number &optional result) |
---|
437 | (if (not result)(setq result (%make-sfloat))) |
---|
438 | (let* ((num (%numerator number)) |
---|
439 | (den (%denominator number))) |
---|
440 | ;; dont error if result is floatable when either top or bottom is |
---|
441 | ;; not. maybe do usual first, catching error |
---|
442 | (if (not (or (bignump num)(bignump den))) |
---|
443 | (target::with-stack-short-floats ((fnum num) |
---|
444 | (fden den)) |
---|
445 | (%short-float/-2! fnum fden result)) |
---|
446 | (let* ((numlen (integer-length num)) |
---|
447 | (denlen (integer-length den)) |
---|
448 | (exp (- numlen denlen)) |
---|
449 | (minusp (minusp num))) |
---|
450 | (if (and (<= numlen IEEE-single-float-bias) |
---|
451 | (<= denlen IEEE-single-float-bias) |
---|
452 | #|(not (minusp exp))|# |
---|
453 | (<= (abs exp) IEEE-single-float-mantissa-width)) |
---|
454 | (target::with-stack-short-floats ((fnum num) |
---|
455 | (fden den)) |
---|
456 | (%short-float/-2! fnum fden result)) |
---|
457 | (if (> exp IEEE-single-float-mantissa-width) |
---|
458 | (progn (%short-float (round num den) result)) |
---|
459 | (if (>= exp 0) |
---|
460 | ; exp between 0 and 23 and nums big |
---|
461 | (let* ((shift (- IEEE-single-float-digits exp)) |
---|
462 | (num (if minusp (- num) num)) |
---|
463 | (int (round (ash num shift) den)) ; gaak |
---|
464 | (intlen (integer-length int)) |
---|
465 | (new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
466 | (when (> intlen IEEE-single-float-digits) |
---|
467 | (setq shift (1- shift)) |
---|
468 | (setq int (round (ash num shift) den)) |
---|
469 | (setq intlen (integer-length int)) |
---|
470 | (setq new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
471 | (when (> new-exp IEEE-single-float-normal-exponent-max) |
---|
472 | (error (make-condition 'floating-point-overflow |
---|
473 | :operation 'short-float |
---|
474 | :operands (list number)))) |
---|
475 | (make-short-float-from-fixnums |
---|
476 | (ldb (byte IEEE-single-float-digits (- intlen IEEE-single-float-digits)) int) |
---|
477 | new-exp |
---|
478 | (if minusp -1 1) |
---|
479 | result)) |
---|
480 | ; den > num - exp negative |
---|
481 | (progn |
---|
482 | (float-rat-neg-exp num den (if minusp -1 1) result t))))))))) |
---|
483 | |
---|
484 | #+64-bit-target |
---|
485 | (defun %short-float-ratio (number) |
---|
486 | (let* ((num (%numerator number)) |
---|
487 | (den (%denominator number))) |
---|
488 | ;; dont error if result is floatable when either top or bottom is |
---|
489 | ;; not. maybe do usual first, catching error |
---|
490 | (if (not (or (bignump num)(bignump den))) |
---|
491 | (/ (the short-float (%short-float num)) |
---|
492 | (the short-float (%short-float den))) |
---|
493 | (let* ((numlen (integer-length num)) |
---|
494 | (denlen (integer-length den)) |
---|
495 | (exp (- numlen denlen)) |
---|
496 | (minusp (minusp num))) |
---|
497 | (if (and (<= numlen IEEE-single-float-bias) |
---|
498 | (<= denlen IEEE-single-float-bias) |
---|
499 | #|(not (minusp exp))|# |
---|
500 | (<= (abs exp) IEEE-single-float-mantissa-width)) |
---|
501 | (/ (the short-float (%short-float num)) |
---|
502 | (the short-float (%short-float den))) |
---|
503 | (if (> exp IEEE-single-float-mantissa-width) |
---|
504 | (progn (%short-float (round num den))) |
---|
505 | (if (>= exp 0) |
---|
506 | ; exp between 0 and 23 and nums big |
---|
507 | (let* ((shift (- IEEE-single-float-digits exp)) |
---|
508 | (num (if minusp (- num) num)) |
---|
509 | (int (round (ash num shift) den)) ; gaak |
---|
510 | (intlen (integer-length int)) |
---|
511 | (new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
512 | (when (> intlen IEEE-single-float-digits) |
---|
513 | (setq shift (1- shift)) |
---|
514 | (setq int (round (ash num shift) den)) |
---|
515 | (setq intlen (integer-length int)) |
---|
516 | (setq new-exp (+ intlen (- IEEE-single-float-bias shift)))) |
---|
517 | (when (> new-exp IEEE-single-float-normal-exponent-max) |
---|
518 | (error (make-condition 'floating-point-overflow |
---|
519 | :operation 'short-float |
---|
520 | :operands (list number)))) |
---|
521 | (make-short-float-from-fixnums |
---|
522 | (ldb (byte IEEE-single-float-digits (- intlen IEEE-single-float-digits)) int) |
---|
523 | new-exp |
---|
524 | (if minusp 1 0))) |
---|
525 | ; den > num - exp negative |
---|
526 | (progn |
---|
527 | (float-rat-neg-exp num den (if minusp -1 1) nil t))))))))) |
---|
528 | |
---|
529 | |
---|
530 | #+32-bit-target |
---|
531 | (defun %short-float (number &optional result) |
---|
532 | (number-case number |
---|
533 | (short-float |
---|
534 | (if result (%copy-short-float number result) number)) |
---|
535 | (double-float |
---|
536 | (%double-float->short-float number (or result (%make-sfloat)))) |
---|
537 | (fixnum |
---|
538 | (%fixnum-sfloat number (or result (%make-sfloat)))) |
---|
539 | (bignum |
---|
540 | (%bignum-sfloat number (or result (%make-sfloat)))) |
---|
541 | (ratio |
---|
542 | (%short-float-ratio number result)))) |
---|
543 | |
---|
544 | #+64-bit-target |
---|
545 | (defun %short-float (number) |
---|
546 | (number-case number |
---|
547 | (short-float number) |
---|
548 | (double-float (%double-float->short-float number)) |
---|
549 | (fixnum (%fixnum-sfloat number)) |
---|
550 | (bignum (%bignum-sfloat number)) |
---|
551 | (ratio (%short-float-ratio number)))) |
---|
552 | |
---|
553 | |
---|
554 | (defun float-rat-neg-exp (integer divisor sign &optional result short) |
---|
555 | (if (minusp sign)(setq integer (- integer))) |
---|
556 | (let* ((integer-length (integer-length integer)) |
---|
557 | ;; make sure we will have enough bits in the quotient |
---|
558 | ;; (and a couple extra for rounding) |
---|
559 | (shift-factor (+ (- (integer-length divisor) integer-length) (if short 28 60))) ; fix |
---|
560 | (scaled-integer integer)) |
---|
561 | (if (plusp shift-factor) |
---|
562 | (setq scaled-integer (ash integer shift-factor)) |
---|
563 | (setq divisor (ash divisor (- shift-factor))) ; assume div > num |
---|
564 | ) |
---|
565 | ;(pprint (list shift-factor scaled-integer divisor)) |
---|
566 | (multiple-value-bind (quotient remainder)(floor scaled-integer divisor) |
---|
567 | (unless (zerop remainder) ; whats this - tells us there's junk below |
---|
568 | (setq quotient (logior quotient 1))) |
---|
569 | ;; why do it return 2 values? |
---|
570 | (values (float-and-scale-and-round sign quotient (- shift-factor) short result))))) |
---|
571 | |
---|
572 | |
---|
573 | |
---|
574 | ;;; when is (negate-bignum (bignum-ashift-right big)) ; can't negate |
---|
575 | ;;; in place cause may get bigger cheaper than (negate-bignum big) - 6 |
---|
576 | ;;; 0r 8 digits ; 8 longs so win if digits > 7 or negate it on the |
---|
577 | ;;; stack |
---|
578 | |
---|
579 | (defun %bignum-dfloat (big &optional result) |
---|
580 | (let* ((minusp (bignum-minusp big))) |
---|
581 | (flet |
---|
582 | ((doit (new-big) |
---|
583 | (let* ((int-len (bignum-integer-length new-big))) |
---|
584 | (when (>= int-len (- 2047 IEEE-double-float-bias)) ; args? |
---|
585 | (error (make-condition 'floating-point-overflow |
---|
586 | :operation 'float :operands (list big)))) |
---|
587 | (if (> int-len 53) |
---|
588 | (let* ((hi (ldb (byte 25 (- int-len 25)) new-big)) |
---|
589 | (lo (ldb (byte 28 (- int-len 53)) new-big))) |
---|
590 | ;(print (list new-big hi lo)) |
---|
591 | (when (and (logbitp (- int-len 54) new-big) ; round bit |
---|
592 | (or (%ilogbitp 0 lo) ; oddp |
---|
593 | ;; or more bits below round |
---|
594 | (%i< (one-bignum-factor-of-two new-big) (- int-len 54)))) |
---|
595 | (if (eq lo #xfffffff) |
---|
596 | (setq hi (1+ hi) lo 0) |
---|
597 | (setq lo (1+ lo))) |
---|
598 | (when (%ilogbitp 25 hi) ; got bigger |
---|
599 | (setq int-len (1+ int-len)) |
---|
600 | (let ((bit (%ilogbitp 0 hi))) |
---|
601 | (setq hi (%ilsr 1 hi)) |
---|
602 | (setq lo (%ilsr 1 lo)) |
---|
603 | (if bit (setq lo (%ilogior #x8000000 lo)))))) |
---|
604 | (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len)(if minusp -1 1) result)) |
---|
605 | (let* ((hi (ldb (byte 25 (- int-len 25)) new-big)) |
---|
606 | (lobits (min (- int-len 25) 28)) |
---|
607 | (lo (ldb (byte lobits (- int-len (+ lobits 25))) new-big))) |
---|
608 | (if (< lobits 28) (setq lo (ash lo (- 28 lobits)))) |
---|
609 | (make-float-from-fixnums hi lo (+ IEEE-double-float-bias int-len) (if minusp -1 1) result)))))) |
---|
610 | (declare (dynamic-extent #'doit)) |
---|
611 | (with-one-negated-bignum-buffer big doit)))) |
---|
612 | |
---|
613 | #+32-bit-target |
---|
614 | (defun %bignum-sfloat (big &optional result) |
---|
615 | (let* ((minusp (bignum-minusp big))) |
---|
616 | (flet |
---|
617 | ((doit (new-big) |
---|
618 | (let* ((int-len (bignum-integer-length new-big))) |
---|
619 | (when (>= int-len (- 255 IEEE-single-float-bias)) ; args? |
---|
620 | (error (make-condition 'floating-point-overflow |
---|
621 | :operation 'float :operands (list big 1.0s0)))) |
---|
622 | (if t ;(> int-len IEEE-single-float-digits) ; always true |
---|
623 | (let* ((lo (ldb (byte IEEE-single-float-digits (- int-len IEEE-single-float-digits)) new-big))) |
---|
624 | (when (and (logbitp (- int-len 25) new-big) ; round bit |
---|
625 | (or (%ilogbitp 0 lo) ; oddp |
---|
626 | ; or more bits below round |
---|
627 | (%i< (one-bignum-factor-of-two new-big) (- int-len 25)))) |
---|
628 | (setq lo (1+ lo)) |
---|
629 | (when (%ilogbitp 24 lo) ; got bigger |
---|
630 | (setq int-len (1+ int-len)) |
---|
631 | (setq lo (%ilsr 1 lo)))) |
---|
632 | (make-short-float-from-fixnums lo (+ IEEE-single-float-bias int-len)(if minusp -1 1) result)) |
---|
633 | )))) |
---|
634 | (declare (dynamic-extent #'doit)) |
---|
635 | (with-one-negated-bignum-buffer big doit)))) |
---|
636 | |
---|
637 | |
---|
638 | #+64-bit-target |
---|
639 | (defun %bignum-sfloat (big) |
---|
640 | (let* ((minusp (bignum-minusp big))) |
---|
641 | (flet |
---|
642 | ((doit (new-big) |
---|
643 | (let* ((int-len (bignum-integer-length new-big))) |
---|
644 | (when (>= int-len (- 255 IEEE-single-float-bias)) ; args? |
---|
645 | (error (make-condition 'floating-point-overflow |
---|
646 | :operation 'float :operands (list big 1.0s0)))) |
---|
647 | (if t ;(> int-len IEEE-single-float-digits) ; always true |
---|
648 | (let* ((lo (ldb (byte IEEE-single-float-digits (- int-len IEEE-single-float-digits)) new-big))) |
---|
649 | (when (and (logbitp (- int-len 25) new-big) ; round bit |
---|
650 | (or (%ilogbitp 0 lo) ; oddp |
---|
651 | ; or more bits below round |
---|
652 | (%i< (one-bignum-factor-of-two new-big) (- int-len 25)))) |
---|
653 | (setq lo (1+ lo)) |
---|
654 | (when (%ilogbitp 24 lo) ; got bigger |
---|
655 | (setq int-len (1+ int-len)) |
---|
656 | (setq lo (%ilsr 1 lo)))) |
---|
657 | (make-short-float-from-fixnums lo (+ IEEE-single-float-bias int-len)(if minusp -1 1))) |
---|
658 | )))) |
---|
659 | (declare (dynamic-extent #'doit)) |
---|
660 | (with-one-negated-bignum-buffer big doit)))) |
---|
661 | |
---|
662 | |
---|
663 | |
---|
664 | |
---|
665 | (defun %fixnum-dfloat (fix &optional result) |
---|
666 | (if (eq 0 fix) |
---|
667 | (if result (%copy-double-float 0.0d0 result) 0.0d0) |
---|
668 | (progn |
---|
669 | (when (not result)(setq result (%make-dfloat))) |
---|
670 | ; it better return result |
---|
671 | (%int-to-dfloat fix result)))) |
---|
672 | |
---|
673 | |
---|
674 | #+32-bit-target |
---|
675 | (defun %fixnum-sfloat (fix &optional result) |
---|
676 | (if (eq 0 fix) |
---|
677 | (if result (%copy-short-float 0.0s0 result) 0.0s0) |
---|
678 | (%int-to-sfloat! fix (or result (%make-sfloat))))) |
---|
679 | |
---|
680 | #+64-bit-target |
---|
681 | (defun %fixnum-sfloat (fix) |
---|
682 | (if (eq 0 fix) |
---|
683 | 0.0s0 |
---|
684 | (%int-to-sfloat fix))) |
---|
685 | |
---|
686 | ;;; Transcendental functions. |
---|
687 | (defun sin (x) |
---|
688 | "Return the sine of NUMBER." |
---|
689 | (if (complexp x) |
---|
690 | (let* ((r (realpart x)) |
---|
691 | (i (imagpart x))) |
---|
692 | (complex (* (sin r) (cosh i)) |
---|
693 | (* (cos r) (sinh i)))) |
---|
694 | (if (typep x 'double-float) |
---|
695 | (%double-float-sin! x (%make-dfloat)) |
---|
696 | #+32-bit-target |
---|
697 | (target::with-stack-short-floats ((sx x)) |
---|
698 | (%single-float-sin! sx (%make-sfloat))) |
---|
699 | #+64-bit-target |
---|
700 | (%single-float-sin (%short-float x))))) |
---|
701 | |
---|
702 | (defun cos (x) |
---|
703 | "Return the cosine of NUMBER." |
---|
704 | (if (complexp x) |
---|
705 | (let* ((r (realpart x)) |
---|
706 | (i (imagpart x))) |
---|
707 | (complex (* (cos r) (cosh i)) |
---|
708 | (- (* (sin r) (sinh i))))) |
---|
709 | (if (typep x 'double-float) |
---|
710 | (%double-float-cos! x (%make-dfloat)) |
---|
711 | #+32-bit-target |
---|
712 | (target::with-stack-short-floats ((sx x)) |
---|
713 | (%single-float-cos! sx (%make-sfloat))) |
---|
714 | #+64-bit-target |
---|
715 | (%single-float-cos (%short-float x))))) |
---|
716 | |
---|
717 | (defun tan (x) |
---|
718 | "Return the tangent of NUMBER." |
---|
719 | (if (complexp x) |
---|
720 | (/ (sin x) (cos x)) |
---|
721 | (if (typep x 'double-float) |
---|
722 | (%double-float-tan! x (%make-dfloat)) |
---|
723 | #+32-bit-target |
---|
724 | (target::with-stack-short-floats ((sx x)) |
---|
725 | (%single-float-tan! sx (%make-sfloat))) |
---|
726 | #+64-bit-target |
---|
727 | (%single-float-tan (%short-float x)) |
---|
728 | ))) |
---|
729 | |
---|
730 | |
---|
731 | |
---|
732 | |
---|
733 | (defun atan (y &optional (x nil x-p)) |
---|
734 | "Return the arc tangent of Y if X is omitted or Y/X if X is supplied." |
---|
735 | (if x-p |
---|
736 | (if (or (typep x 'double-float) |
---|
737 | (typep y 'double-float)) |
---|
738 | (with-stack-double-floats ((dy y) |
---|
739 | (dx x)) |
---|
740 | (%df-atan2 dy dx)) |
---|
741 | #+32-bit-target |
---|
742 | (target::with-stack-short-floats ((sy y) |
---|
743 | (sx x)) |
---|
744 | (%sf-atan2! sy sx)) |
---|
745 | #+64-bit-target |
---|
746 | (%sf-atan2 (%short-float y) (%short-float x))) |
---|
747 | (if (typep y 'complex) |
---|
748 | (let* ((iy (* (sqrt -1) y))) |
---|
749 | (/ (- (log (+ 1 iy)) (log (- 1 iy))) |
---|
750 | #c(0 2))) |
---|
751 | (if (typep y 'double-float) |
---|
752 | (%double-float-atan! y (%make-dfloat)) |
---|
753 | #+32-bit-target |
---|
754 | (target::with-stack-short-floats ((sy y)) |
---|
755 | (%single-float-atan! sy (%make-sfloat))) |
---|
756 | #+64-bit-target |
---|
757 | (%single-float-atan (%short-float y)) |
---|
758 | )))) |
---|
759 | |
---|
760 | |
---|
761 | |
---|
762 | (defun log (x &optional (b nil b-p)) |
---|
763 | "Return the logarithm of NUMBER in the base BASE, which defaults to e." |
---|
764 | (if b-p |
---|
765 | (if (zerop b) |
---|
766 | (if (zerop x) |
---|
767 | (report-bad-arg x '(not (satisfies zerop) )) |
---|
768 | (if (floatp x) (float 0.0d0 x) 0)) |
---|
769 | (/ (log-e x) (log-e b))) |
---|
770 | (log-e x))) |
---|
771 | |
---|
772 | (defun log-e (x) |
---|
773 | (cond |
---|
774 | ((bignump x) |
---|
775 | (if (minusp x) |
---|
776 | (complex (log-e (- x)) pi) |
---|
777 | (let* ((base1 3) |
---|
778 | (guess (floor (1- (integer-length x)) |
---|
779 | (log base1 2))) |
---|
780 | (guess1 (* guess (log-e base1)))) |
---|
781 | (+ guess1 (log-e (/ x (expt base1 guess))))))) |
---|
782 | ((and (ratiop x) |
---|
783 | (or (> x most-positive-short-float) |
---|
784 | (< x most-negative-short-float))) |
---|
785 | (- (log-e (%numerator x)) (log-e (%denominator x)))) |
---|
786 | ((typep x 'complex) |
---|
787 | (complex (log-e (abs x)) (phase x))) |
---|
788 | ((typep x 'double-float) |
---|
789 | (with-stack-double-floats ((dx x)) |
---|
790 | (if (minusp x) |
---|
791 | (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi) |
---|
792 | (%double-float-log! dx (%make-dfloat))))) |
---|
793 | (t |
---|
794 | #+32-bit-target |
---|
795 | (target::with-stack-short-floats ((sx x)) |
---|
796 | (if (minusp x) |
---|
797 | (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat)) |
---|
798 | #.(coerce pi 'short-float)) |
---|
799 | (%single-float-log! sx (%make-sfloat)))) |
---|
800 | #+64-bit-target |
---|
801 | (if (minusp x) |
---|
802 | (complex (%single-float-log (%short-float-abs (%short-float x))) #.(coerce pi 'single-float)) |
---|
803 | (%single-float-log (%short-float x))) |
---|
804 | ))) |
---|
805 | |
---|
806 | |
---|
807 | |
---|
808 | (defun exp (x) |
---|
809 | "Return e raised to the power NUMBER." |
---|
810 | (typecase x |
---|
811 | (complex (* (exp (realpart x)) (cis (imagpart x)))) |
---|
812 | (double-float (%double-float-exp! x (%make-dfloat))) |
---|
813 | (t |
---|
814 | #+32-bit-target |
---|
815 | (target::with-stack-short-floats ((sx x)) |
---|
816 | (%single-float-exp! sx (%make-sfloat))) |
---|
817 | #+64-bit-target |
---|
818 | (%single-float-exp (%short-float x))))) |
---|
819 | |
---|
820 | |
---|
821 | |
---|
822 | (defun expt (b e) |
---|
823 | "Return BASE raised to the POWER." |
---|
824 | (cond ((zerop e) (1+ (* b e))) |
---|
825 | ((integerp e) |
---|
826 | (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e))) |
---|
827 | ((zerop b) |
---|
828 | (if (plusp (realpart e)) b (report-bad-arg e '(number (0) *)))) |
---|
829 | ((and (realp b) (plusp b) (realp e)) |
---|
830 | (if (or (typep b 'double-float) |
---|
831 | (typep e 'double-float)) |
---|
832 | (with-stack-double-floats ((b1 b) |
---|
833 | (e1 e)) |
---|
834 | (%double-float-expt! b1 e1 (%make-dfloat))) |
---|
835 | #+32-bit-target |
---|
836 | (target::with-stack-short-floats ((b1 b) |
---|
837 | (e1 e)) |
---|
838 | (%single-float-expt! b1 e1 (%make-sfloat))) |
---|
839 | #+64-bit-target |
---|
840 | (%single-float-expt (%short-float b) (%short-float e)) |
---|
841 | )) |
---|
842 | ((typep (realpart e) 'double-float) |
---|
843 | ;; Avoid intermediate single-float result from LOG |
---|
844 | (let ((promoted-base (* 1d0 b))) |
---|
845 | (exp (* e (log promoted-base))))) |
---|
846 | (t (exp (* e (log b)))))) |
---|
847 | |
---|
848 | |
---|
849 | |
---|
850 | (defun sqrt (x &aux a b) |
---|
851 | "Return the square root of NUMBER." |
---|
852 | (cond ((zerop x) x) |
---|
853 | ((complexp x) (* (sqrt (abs x)) (cis (/ (phase x) 2)))) |
---|
854 | ((minusp x) (complex 0 (sqrt (- x)))) |
---|
855 | ((floatp x) |
---|
856 | (fsqrt x)) |
---|
857 | ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a) |
---|
858 | ((and (ratiop x) |
---|
859 | (let ((n (numerator x)) |
---|
860 | d) |
---|
861 | (and (eql n (* (setq a (isqrt n)) a)) |
---|
862 | (eql (setq d (denominator x)) |
---|
863 | (* (setq b (isqrt d)) b))))) |
---|
864 | (/ a b)) |
---|
865 | (t |
---|
866 | #+32-bit-target |
---|
867 | (target::with-stack-short-floats ((f1)) |
---|
868 | (fsqrt (%short-float x f1))) |
---|
869 | #+64-bit-target |
---|
870 | (fsqrt (%short-float x))))) |
---|
871 | |
---|
872 | |
---|
873 | |
---|
874 | (defun asin (x) |
---|
875 | "Return the arc sine of NUMBER." |
---|
876 | (number-case x |
---|
877 | (complex |
---|
878 | (let ((sqrt-1-x (sqrt (- 1 x))) |
---|
879 | (sqrt-1+x (sqrt (+ 1 x)))) |
---|
880 | (complex (atan (/ (realpart x) |
---|
881 | (realpart (* sqrt-1-x sqrt-1+x)))) |
---|
882 | (asinh (imagpart (* (conjugate sqrt-1-x) |
---|
883 | sqrt-1+x)))))) |
---|
884 | (double-float |
---|
885 | (locally (declare (type double-float x)) |
---|
886 | (if (and (<= -1.0d0 x) |
---|
887 | (<= x 1.0d0)) |
---|
888 | (%double-float-asin! x (%make-dfloat)) |
---|
889 | (let* ((temp (+ (complex -0.0d0 x) |
---|
890 | (sqrt (- 1.0d0 (the double-float (* x x))))))) |
---|
891 | (complex (phase temp) (- (log (abs temp)))))))) |
---|
892 | ((short-float rational) |
---|
893 | #+32-bit-target |
---|
894 | (let* ((x1 (%make-sfloat))) |
---|
895 | (declare (dynamic-extent x1)) |
---|
896 | (if (and (realp x) |
---|
897 | (<= -1.0s0 (setq x (%short-float x x1))) |
---|
898 | (<= x 1.0s0)) |
---|
899 | (%single-float-asin! x1 (%make-sfloat)) |
---|
900 | (progn |
---|
901 | (setq x (+ (complex (- (imagpart x)) (realpart x)) |
---|
902 | (sqrt (- 1 (* x x))))) |
---|
903 | (complex (phase x) (- (log (abs x))))))) |
---|
904 | #+64-bit-target |
---|
905 | (if (and (realp x) |
---|
906 | (<= -1.0s0 (setq x (%short-float x))) |
---|
907 | (<= x 1.0s0)) |
---|
908 | (%single-float-asin x) |
---|
909 | (progn |
---|
910 | (setq x (+ (complex (- (imagpart x)) (realpart x)) |
---|
911 | (sqrt (- 1 (* x x))))) |
---|
912 | (complex (phase x) (- (log (abs x)))))) |
---|
913 | ))) |
---|
914 | |
---|
915 | |
---|
916 | (eval-when (:execute :compile-toplevel) |
---|
917 | (defconstant double-float-half-pi (asin 1.0d0)) |
---|
918 | (defconstant single-float-half-pi (asin 1.0f0)) |
---|
919 | ) |
---|
920 | |
---|
921 | |
---|
922 | |
---|
923 | (defun acos (x) |
---|
924 | "Return the arc cosine of NUMBER." |
---|
925 | (number-case x |
---|
926 | (complex |
---|
927 | (let ((sqrt-1+x (sqrt (+ 1 x))) |
---|
928 | (sqrt-1-x (sqrt (- 1 x)))) |
---|
929 | (complex (* 2 (atan (/ (realpart sqrt-1-x) |
---|
930 | (realpart sqrt-1+x)))) |
---|
931 | (asinh (imagpart (* (conjugate sqrt-1+x) |
---|
932 | sqrt-1-x)))))) |
---|
933 | |
---|
934 | (double-float |
---|
935 | (locally (declare (type double-float x)) |
---|
936 | (if (and (<= -1.0d0 x) |
---|
937 | (<= x 1.0d0)) |
---|
938 | (%double-float-acos! x (%make-dfloat)) |
---|
939 | (- double-float-half-pi (asin x))))) |
---|
940 | ((short-float rational) |
---|
941 | #+32-bit-target |
---|
942 | (target::with-stack-short-floats ((sx x)) |
---|
943 | (locally |
---|
944 | (declare (type short-float sx)) |
---|
945 | (if (and (<= -1.0s0 sx) |
---|
946 | (<= sx 1.0s0)) |
---|
947 | (%single-float-acos! sx (%make-sfloat)) |
---|
948 | (- single-float-half-pi (asin sx))))) |
---|
949 | #+64-bit-target |
---|
950 | (let* ((sx (%short-float x))) |
---|
951 | (declare (type short-float sx)) |
---|
952 | (if (and (<= -1.0s0 sx) |
---|
953 | (<= sx 1.0s0)) |
---|
954 | (%single-float-acos sx) |
---|
955 | (- single-float-half-pi (asin sx)))) |
---|
956 | ))) |
---|
957 | |
---|
958 | |
---|
959 | (defun fsqrt (x) |
---|
960 | (etypecase x |
---|
961 | (double-float (%double-float-sqrt! x (%make-dfloat))) |
---|
962 | (single-float |
---|
963 | #+32-bit-target |
---|
964 | (%single-float-sqrt! x (%make-sfloat)) |
---|
965 | #+64-bit-target |
---|
966 | (%single-float-sqrt x)))) |
---|
967 | |
---|
968 | |
---|
969 | |
---|
970 | (defun %df-atan2 (y x &optional result) |
---|
971 | (if (zerop x) |
---|
972 | (if (zerop y) |
---|
973 | (if (plusp (float-sign x)) |
---|
974 | y |
---|
975 | (float-sign y pi)) |
---|
976 | (float-sign y double-float-half-pi)) |
---|
977 | (%double-float-atan2! y x (or result (%make-dfloat))))) |
---|
978 | |
---|
979 | #+32-bit-target |
---|
980 | (defun %sf-atan2! (y x &optional result) |
---|
981 | (if (zerop x) |
---|
982 | (if (zerop y) |
---|
983 | (if (plusp (float-sign x)) |
---|
984 | y |
---|
985 | (float-sign y pi)) |
---|
986 | (float-sign y single-float-half-pi)) |
---|
987 | (%single-float-atan2! y x (or result (%make-sfloat))))) |
---|
988 | |
---|
989 | #+64-bit-target |
---|
990 | (defun %sf-atan2 (y x) |
---|
991 | (if (zerop x) |
---|
992 | (if (zerop y) |
---|
993 | (if (plusp (float-sign x)) |
---|
994 | y |
---|
995 | (float-sign y pi)) |
---|
996 | (float-sign y single-float-half-pi)) |
---|
997 | (%single-float-atan2 y x))) |
---|
998 | |
---|
999 | #+64-bit-target |
---|
1000 | (defun %short-float-exp (n) |
---|
1001 | (let* ((bits (single-float-bits n))) |
---|
1002 | (declare (type (unsigned-byte 32) bits)) |
---|
1003 | (ldb (byte IEEE-single-float-exponent-width IEEE-single-float-exponent-offset) bits))) |
---|
1004 | |
---|
1005 | |
---|
1006 | #+64-bit-target |
---|
1007 | (defun set-%short-float-exp (float exp) |
---|
1008 | (host-single-float-from-unsigned-byte-32 |
---|
1009 | (dpb exp |
---|
1010 | (byte IEEE-single-float-exponent-width |
---|
1011 | IEEE-single-float-exponent-offset) |
---|
1012 | (the (unsigned-byte 32) (single-float-bits float))))) |
---|
1013 | |
---|
1014 | #+64-bit-target |
---|
1015 | (defun %%scale-sfloat (float int) |
---|
1016 | (* (the single-float float) |
---|
1017 | (the single-float (host-single-float-from-unsigned-byte-32 |
---|
1018 | (dpb int |
---|
1019 | (byte IEEE-single-float-exponent-width |
---|
1020 | IEEE-single-float-exponent-offset) |
---|
1021 | 0))))) |
---|
1022 | |
---|
1023 | #+64-bit-target |
---|
1024 | (defun %double-float-exp (n) |
---|
1025 | (let* ((highword (double-float-bits n))) |
---|
1026 | (declare (fixnum highword)) |
---|
1027 | (logand (1- (ash 1 IEEE-double-float-exponent-width)) |
---|
1028 | (ash highword (- (- IEEE-double-float-exponent-offset 32)))))) |
---|
1029 | |
---|
1030 | #+64-bit-target |
---|
1031 | (defun set-%double-float-exp (float exp) |
---|
1032 | (let* ((highword (double-float-bits float))) |
---|
1033 | (declare (fixnum highword)) |
---|
1034 | (setf (uvref float target::double-float.val-high-cell) |
---|
1035 | (dpb exp |
---|
1036 | (byte IEEE-double-float-exponent-width |
---|
1037 | (- IEEE-double-float-exponent-offset 32)) |
---|
1038 | highword)) |
---|
1039 | exp)) |
---|
1040 | |
---|
1041 | #+64-bit-target |
---|
1042 | (defun %integer-decode-double-float (f) |
---|
1043 | (multiple-value-bind (hiword loword) (double-float-bits f) |
---|
1044 | (declare (type (unsigned-byte 32) hiword loword)) |
---|
1045 | (let* ((exp (ldb (byte IEEE-double-float-exponent-width |
---|
1046 | (- IEEE-double-float-exponent-offset 32)) |
---|
1047 | hiword)) |
---|
1048 | (mantissa (logior |
---|
1049 | (the fixnum |
---|
1050 | (dpb (ldb (byte (- IEEE-double-float-mantissa-width 32) |
---|
1051 | IEEE-double-float-mantissa-offset) |
---|
1052 | hiword) |
---|
1053 | (byte (- IEEE-double-float-mantissa-width 32) |
---|
1054 | 32) |
---|
1055 | loword)) |
---|
1056 | (if (zerop exp) |
---|
1057 | 0 |
---|
1058 | (ash 1 IEEE-double-float-hidden-bit)))) |
---|
1059 | (sign (if (logbitp 31 hiword) -1 1))) |
---|
1060 | (declare (fixnum exp mantissa sign)) |
---|
1061 | (values (ldb (byte 25 28) mantissa) |
---|
1062 | (ldb (byte 28 0) mantissa) |
---|
1063 | exp |
---|
1064 | sign)))) |
---|
1065 | |
---|
1066 | ;;; end of l0-float.lisp |
---|