/[cmucl]/src/code/float.lisp
ViewVC logotype

Contents of /src/code/float.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.49 - (show annotations)
Sat Sep 3 05:19:03 2011 UTC (2 years, 7 months ago) by rtoy
Branch: MAIN
CVS Tags: GIT-CONVERSION, snapshot-2011-09, HEAD
Changes since 1.48: +36 -31 lines
Fix rounding for large numbers.

Bug was pointed by Christophe in private email.  Fix is based on his
suggested solution.  Some examples that should work now:

(round 100000000002.9d0) -> 100000000003

(round (+ most-positive-fixnum 1.5w0)) -> 536870912
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written as part of the CMU Common Lisp project at
5 ;;; Carnegie Mellon University, and has been placed in the public domain.
6 ;;;
7 (ext:file-comment
8 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/float.lisp,v 1.49 2011/09/03 05:19:03 rtoy Exp $")
9 ;;;
10 ;;; **********************************************************************
11 ;;;
12 ;;; This file contains the definitions of float specific number support
13 ;;; (other than irrational stuff, which is in irrat.) There is code in here
14 ;;; that assumes there are only two float formats: IEEE single and double.
15 ;;;
16 ;;; Author: Rob MacLachlan
17 ;;; Long-float support by Douglas Crosher, 1998.
18 ;;;
19 (in-package "KERNEL")
20 (intl:textdomain "cmucl")
21
22 (export '(%unary-truncate %unary-round %unary-ftruncate
23 %unary-ftruncate/single-float %unary-ftruncate/double-float))
24
25 #-x87
26 (export '(%unary-fround/single-float %unary-fround/double-float))
27
28 (in-package "LISP")
29 (export '(least-positive-normalized-short-float
30 least-positive-normalized-single-float
31 least-positive-normalized-double-float
32 least-positive-normalized-long-float
33 least-negative-normalized-short-float
34 least-negative-normalized-single-float
35 least-negative-normalized-double-float
36 least-negative-normalized-long-float
37 least-positive-single-float
38 least-positive-short-float
39 least-negative-single-float
40 least-negative-short-float
41 least-positive-double-float
42 least-positive-long-float
43 least-negative-double-float
44 least-negative-long-float
45 most-positive-single-float
46 most-positive-short-float
47 most-negative-single-float
48 most-negative-short-float
49 most-positive-double-float
50 most-positive-long-float
51 most-negative-double-float
52 most-negative-long-float))
53
54 (in-package "EXTENSIONS")
55 (export '(single-float-positive-infinity short-float-positive-infinity
56 double-float-positive-infinity long-float-positive-infinity
57 single-float-negative-infinity short-float-negative-infinity
58 double-float-negative-infinity long-float-negative-infinity
59 set-floating-point-modes float-denormalized-p float-nan-p
60 float-trapping-nan-p float-infinity-p))
61
62 (in-package "KERNEL")
63
64
65 ;;;; Utilities:
66
67 ;;; SINGLE-FROM-BITS, DOUBLE-FROM-BITS -- Internal
68 ;;;
69 ;;; These functions let us create floats from bits with the significand
70 ;;; uniformly represented as an integer. This is less efficient for double
71 ;;; floats, but is more convenient when making special values, etc.
72 ;;;
73 (defun single-from-bits (sign exp sig)
74 (declare (type bit sign) (type (unsigned-byte 24) sig)
75 (type (unsigned-byte 8) exp))
76 (make-single-float
77 (dpb exp vm:single-float-exponent-byte
78 (dpb sig vm:single-float-significand-byte
79 (if (zerop sign) 0 -1)))))
80 ;;;
81 (defun double-from-bits (sign exp sig)
82 (declare (type bit sign) (type (unsigned-byte 53) sig)
83 (type (unsigned-byte 11) exp))
84 (make-double-float (dpb exp vm:double-float-exponent-byte
85 (dpb (ash sig -32) vm:double-float-significand-byte
86 (if (zerop sign) 0 -1)))
87 (ldb (byte 32 0) sig)))
88
89 #+double-double
90 (defun double-double-from-bits (sign exp sig)
91 (declare (type bit sign) (type (unsigned-byte #.vm:double-double-float-digits) sig)
92 (type (unsigned-byte 11) exp)
93 (ignore exp))
94 (let ((lo (ldb (byte vm:double-float-digits 0) sig))
95 (hi (ldb (byte vm:double-float-digits vm:double-float-digits) sig))
96 (fsign (- 1 (* 2 sign))))
97 ;; We can't return a double-double because exp might be very
98 ;; small, and the low part of the double-double would get scaled
99 ;; to 0. So we return both parts and expect the caller to figure
100 ;; out what to do.
101 (values (scale-float (* fsign (float hi 1d0)) #.(- vm:double-float-digits))
102 (scale-float (* fsign (float lo 1d0)) #.(- vm:double-float-digits)))))
103 ;;;
104 #+(and long-float x86)
105 (defun long-from-bits (sign exp sig)
106 (declare (type bit sign) (type (unsigned-byte 64) sig)
107 (type (unsigned-byte 15) exp))
108 (make-long-float (logior (ash sign 15) exp)
109 (ldb (byte 32 32) sig)
110 (ldb (byte 32 0) sig)))
111
112
113 ;;;; Float parameters:
114
115 (defconstant least-positive-single-float (single-from-bits 0 0 1))
116 (defconstant least-positive-short-float least-positive-single-float)
117 (defconstant least-negative-single-float (single-from-bits 1 0 1))
118 (defconstant least-negative-short-float least-negative-single-float)
119 (defconstant least-positive-double-float (double-from-bits 0 0 1))
120 #-long-float
121 (defconstant least-positive-long-float least-positive-double-float)
122 #+(and long-float x86)
123 (defconstant least-positive-long-float (long-from-bits 0 0 1))
124 (defconstant least-negative-double-float (double-from-bits 1 0 1))
125 #-long-float
126 (defconstant least-negative-long-float least-negative-double-float)
127 #+(and long-float x86)
128 (defconstant least-negative-long-float (long-from-bits 1 0 1))
129
130 (defconstant least-positive-normalized-single-float
131 (single-from-bits 0 vm:single-float-normal-exponent-min 0))
132 (defconstant least-positive-normalized-short-float
133 least-positive-normalized-single-float)
134 (defconstant least-negative-normalized-single-float
135 (single-from-bits 1 vm:single-float-normal-exponent-min 0))
136 (defconstant least-negative-normalized-short-float
137 least-negative-normalized-single-float)
138 (defconstant least-positive-normalized-double-float
139 (double-from-bits 0 vm:double-float-normal-exponent-min 0))
140 #-long-float
141 (defconstant least-positive-normalized-long-float
142 least-positive-normalized-double-float)
143 #+(and long-float x86)
144 (defconstant least-positive-normalized-long-float
145 (long-from-bits 0 vm:long-float-normal-exponent-min
146 (ash vm:long-float-hidden-bit 32)))
147 (defconstant least-negative-normalized-double-float
148 (double-from-bits 1 vm:double-float-normal-exponent-min 0))
149 #-long-float
150 (defconstant least-negative-normalized-long-float
151 least-negative-normalized-double-float)
152 #+(and long-float x86)
153 (defconstant least-negative-normalized-long-float
154 (long-from-bits 1 vm:long-float-normal-exponent-min
155 (ash vm:long-float-hidden-bit 32)))
156
157 (defconstant most-positive-single-float
158 (single-from-bits 0 vm:single-float-normal-exponent-max
159 (ldb vm:single-float-significand-byte -1)))
160 (defconstant most-positive-short-float most-positive-single-float)
161 (defconstant most-negative-single-float
162 (single-from-bits 1 vm:single-float-normal-exponent-max
163 (ldb vm:single-float-significand-byte -1)))
164 (defconstant most-negative-short-float most-negative-single-float)
165 (defconstant most-positive-double-float
166 (double-from-bits 0 vm:double-float-normal-exponent-max
167 (ldb (byte vm:double-float-digits 0) -1)))
168 #-long-float
169 (defconstant most-positive-long-float most-positive-double-float)
170 #+(and long-float x86)
171 (defconstant most-positive-long-float
172 (long-from-bits 0 vm:long-float-normal-exponent-max
173 (ldb (byte vm:long-float-digits 0) -1)))
174 (defconstant most-negative-double-float
175 (double-from-bits 1 vm:double-float-normal-exponent-max
176 (ldb (byte vm:double-float-digits 0) -1)))
177 #-long-float
178 (defconstant most-negative-long-float most-negative-double-float)
179 #+(and long-float x86)
180 (defconstant most-negative-long-float
181 (long-from-bits 1 vm:long-float-normal-exponent-max
182 (ldb (byte vm:long-float-digits 0) -1)))
183
184 (defconstant single-float-positive-infinity
185 (single-from-bits 0 (1+ vm:single-float-normal-exponent-max) 0))
186 (defconstant short-float-positive-infinity single-float-positive-infinity)
187 (defconstant single-float-negative-infinity
188 (single-from-bits 1 (1+ vm:single-float-normal-exponent-max) 0))
189 (defconstant short-float-negative-infinity single-float-negative-infinity)
190 (defconstant double-float-positive-infinity
191 (double-from-bits 0 (1+ vm:double-float-normal-exponent-max) 0))
192 #-long-float
193 (defconstant long-float-positive-infinity double-float-positive-infinity)
194 #+(and long-float x86)
195 (defconstant long-float-positive-infinity
196 (long-from-bits 0 (1+ vm:long-float-normal-exponent-max)
197 (ash vm:long-float-hidden-bit 32)))
198 (defconstant double-float-negative-infinity
199 (double-from-bits 1 (1+ vm:double-float-normal-exponent-max) 0))
200 #-long-float
201 (defconstant long-float-negative-infinity double-float-negative-infinity)
202 #+(and long-float x86)
203 (defconstant long-float-negative-infinity
204 (long-from-bits 1 (1+ vm:long-float-normal-exponent-max)
205 (ash vm:long-float-hidden-bit 32)))
206
207 (defconstant single-float-epsilon
208 (single-from-bits 0 (- vm:single-float-bias (1- vm:single-float-digits)) 1))
209 (defconstant short-float-epsilon single-float-epsilon)
210 (defconstant single-float-negative-epsilon
211 (single-from-bits 0 (- vm:single-float-bias vm:single-float-digits) 1))
212 (defconstant short-float-negative-epsilon single-float-negative-epsilon)
213 #-(and long-float x86)
214 (defconstant double-float-epsilon
215 (double-from-bits 0 (- vm:double-float-bias (1- vm:double-float-digits)) 1))
216 #+(and long-float x86)
217 (defconstant double-float-epsilon
218 (double-from-bits 0 (- vm:double-float-bias (1- vm:double-float-digits))
219 (expt 2 42)))
220 #-long-float
221 (defconstant long-float-epsilon double-float-epsilon)
222 #+(and long-float x86)
223 (defconstant long-float-epsilon
224 (long-from-bits 0 (- vm:long-float-bias (1- vm:long-float-digits))
225 (+ 1 (ash vm:long-float-hidden-bit 32))))
226 #-(and long-float x86)
227 (defconstant double-float-negative-epsilon
228 (double-from-bits 0 (- vm:double-float-bias vm:double-float-digits) 1))
229 #+(and long-float x86)
230 (defconstant double-float-negative-epsilon
231 (double-from-bits 0 (- vm:double-float-bias vm:double-float-digits)
232 (expt 2 42)))
233 #-long-float
234 (defconstant long-float-negative-epsilon double-float-negative-epsilon)
235 #+(and long-float x86)
236 (defconstant long-float-negative-epsilon
237 (long-from-bits 0 (- vm:long-float-bias vm:long-float-digits)
238 (+ 1 (ash vm:long-float-hidden-bit 32))))
239
240
241 ;;;; Float predicates and environment query:
242
243 (declaim (maybe-inline float-denormalized-p float-infinity-p float-nan-p
244 float-trapping-nan-p))
245
246 ;;; FLOAT-DENORMALIZED-P -- Public
247 ;;;
248 (defun float-denormalized-p (x)
249 "Return true if the float X is denormalized."
250 (number-dispatch ((x float))
251 ((single-float)
252 (and (zerop (ldb vm:single-float-exponent-byte (single-float-bits x)))
253 (not (zerop x))))
254 ((double-float)
255 (and (zerop (ldb vm:double-float-exponent-byte
256 (double-float-high-bits x)))
257 (not (zerop x))))
258 #+(and long-float x86)
259 ((long-float)
260 (and (zerop (ldb vm:long-float-exponent-byte (long-float-exp-bits x)))
261 (not (zerop x))))))
262
263 (macrolet ((frob (name doc single double #+(and long-float x86) long
264 #+double-double double-double)
265 `(defun ,name (x)
266 ,doc
267 (number-dispatch ((x float))
268 ((single-float)
269 (let ((bits (single-float-bits x)))
270 (and (> (ldb vm:single-float-exponent-byte bits)
271 vm:single-float-normal-exponent-max)
272 ,single)))
273 ((double-float)
274 (let ((hi (double-float-high-bits x))
275 (lo (double-float-low-bits x)))
276 (declare (ignorable lo))
277 (and (> (ldb vm:double-float-exponent-byte hi)
278 vm:double-float-normal-exponent-max)
279 ,double)))
280 #+(and long-float x86)
281 ((long-float)
282 (let ((exp (long-float-exp-bits x))
283 (hi (long-float-high-bits x))
284 (lo (long-float-low-bits x)))
285 (declare (ignorable lo))
286 (and (> (ldb vm:long-float-exponent-byte exp)
287 vm:long-float-normal-exponent-max)
288 ,long)))
289 #+double-double
290 ((double-double-float)
291 ,double-double)))))
292
293 (frob float-infinity-p "Return true if the float X is an infinity (+ or -)."
294 (zerop (ldb vm:single-float-significand-byte bits))
295 (and (zerop (ldb vm:double-float-significand-byte hi))
296 (zerop lo))
297 #+(and long-float x86)
298 (and (zerop (ldb vm:long-float-significand-byte hi))
299 (zerop lo))
300 #+double-double
301 (float-infinity-p (double-double-hi x)))
302
303 (frob float-nan-p "Return true if the float X is a NaN (Not a Number)."
304 (not (zerop (ldb vm:single-float-significand-byte bits)))
305 (or (not (zerop (ldb vm:double-float-significand-byte hi)))
306 (not (zerop lo)))
307 #+(and long-float x86)
308 (or (not (zerop (ldb vm:long-float-significand-byte hi)))
309 (not (zerop lo)))
310 #+double-double
311 (float-nan-p (double-double-hi x)))
312
313 (frob float-trapping-nan-p
314 "Return true if the float X is a trapping NaN (Not a Number)."
315 (zerop (logand (ldb vm:single-float-significand-byte bits)
316 vm:single-float-trapping-nan-bit))
317 (zerop (logand (ldb vm:double-float-significand-byte hi)
318 vm:double-float-trapping-nan-bit))
319 #+(and long-float x86)
320 (zerop (logand (ldb vm:long-float-significand-byte hi)
321 vm:long-float-trapping-nan-bit))
322 #+double-double
323 (float-trapping-nan-p (double-double-hi x))))
324
325
326 ;;; FLOAT-PRECISION -- Public
327 ;;;
328 ;;; If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
329 ;;; actual exponent (and hence how denormalized it is), otherwise we just
330 ;;; return the number of digits or 0.
331 ;;;
332 (declaim (maybe-inline float-precision))
333 (defun float-precision (f)
334 "Returns a non-negative number of significant digits in it's float argument.
335 Will be less than FLOAT-DIGITS if denormalized or zero."
336 (macrolet ((frob (digits bias decode)
337 `(cond ((zerop f) 0)
338 ((float-denormalized-p f)
339 (multiple-value-bind (ignore exp)
340 (,decode f)
341 (declare (ignore ignore))
342 (truly-the fixnum
343 (+ ,digits (1- ,digits) ,bias exp))))
344 (t
345 ,digits))))
346 (number-dispatch ((f float))
347 ((single-float)
348 (frob vm:single-float-digits vm:single-float-bias
349 integer-decode-single-denorm))
350 ((double-float)
351 (frob vm:double-float-digits vm:double-float-bias
352 integer-decode-double-denorm))
353 #+long-float
354 ((long-float)
355 (frob vm:long-float-digits vm:long-float-bias
356 integer-decode-long-denorm))
357 #+double-double
358 ((double-double-float)
359 ;; What exactly is the precision for a double-double? We make
360 ;; it the sum of the precisions of the two components.
361 (let ((hi (double-double-hi f)))
362 (cond ((zerop hi)
363 ;; ANSI CL says the precision is 0
364 0)
365 ((< (abs hi) (scale-float least-positive-normalized-double-float 53))
366 ;; For every power of 2 below this, we lose a bit of
367 ;; precision. More or less.
368 (let ((cutoff
369 (nth-value
370 1
371 (decode-float
372 (scale-float least-positive-normalized-double-float 53)))))
373 (multiple-value-bind (f exp)
374 (decode-float hi)
375 (declare (ignore f))
376 (- 106 (- cutoff exp)))))
377 (t
378 ;; Normally we have 106 bits of precision (twice the
379 ;; double-float precision)
380 106)))))))
381
382
383 #+nil
384 (defun float-sign (float1 &optional (float2 (float 1 float1)))
385 "Returns a floating-point number that has the same sign as
386 float1 and, if float2 is given, has the same absolute value
387 as float2."
388 (declare (float float1 float2))
389 (* (if (etypecase float1
390 (single-float (minusp (single-float-bits float1)))
391 (double-float (minusp (double-float-high-bits float1)))
392 #+long-float
393 (long-float (minusp (long-float-exp-bits float1)))
394 #+double-double
395 (double-double-float (minusp (double-double-hi float1))))
396 (float -1 float1)
397 (float 1 float1))
398 (abs float2)))
399
400 (defun float-sign (float1 &optional float2)
401 "Returns a floating-point number that has the same sign as
402 float1 and, if float2 is given, has the same absolute value
403 as float2."
404 (declare (float float1)
405 (type (or null float) float2))
406 (let ((f1-sign (if (etypecase float1
407 (single-float (minusp (single-float-bits float1)))
408 (double-float (minusp (double-float-high-bits float1)))
409 #+long-float
410 (long-float (minusp (long-float-exp-bits float1)))
411 #+double-double
412 (double-double-float (minusp (float-sign (double-double-hi float1)))))
413 (float -1 float1)
414 (float 1 float1))))
415 ;; Multiplication of double-double-float's doesn't preserve the
416 ;; sign of signed-zeroes, so we split the case of float2 this way
417 ;; so we can get the right sign.
418 (if float2
419 (if (minusp f1-sign)
420 (- (abs float2))
421 (abs float2))
422 f1-sign)))
423
424 (defun float-format-digits (format)
425 (ecase format
426 ((short-float single-float) vm:single-float-digits)
427 ((double-float #-long-float long-float) vm:double-float-digits)
428 #+long-float
429 (long-float vm:long-float-digits)
430 #+double-double
431 (double-double-float vm:double-double-float-digits)))
432
433 (declaim (inline float-digits float-radix))
434
435 (defun float-digits (f)
436 "Returns a non-negative number of radix-b digits used in the
437 representation of it's argument. See Common Lisp: The Language
438 by Guy Steele for more details."
439 (number-dispatch ((f float))
440 ((single-float) vm:single-float-digits)
441 ((double-float) vm:double-float-digits)
442 #+long-float
443 ((long-float) vm:long-float-digits)
444 #+double-double
445 ((double-double-float) vm:double-double-float-digits)))
446
447 (defun float-radix (f)
448 "Returns (as an integer) the radix b of its floating-point
449 argument."
450 (number-dispatch ((f float))
451 ((float) 2)))
452
453
454
455 ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT:
456
457 (declaim (maybe-inline integer-decode-single-float
458 integer-decode-double-float))
459
460 ;;; INTEGER-DECODE-SINGLE-DENORM -- Internal
461 ;;;
462 ;;; Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
463 ;;;
464 (defun integer-decode-single-denorm (x)
465 (declare (type single-float x))
466 (let* ((bits (single-float-bits (abs x)))
467 (sig (ash (ldb vm:single-float-significand-byte bits) 1))
468 (extra-bias 0))
469 (declare (type (unsigned-byte 24) sig)
470 (type (integer 0 23) extra-bias))
471 (loop
472 (unless (zerop (logand sig vm:single-float-hidden-bit))
473 (return))
474 (setq sig (ash sig 1))
475 (incf extra-bias))
476 (values sig
477 (- (- vm:single-float-bias) vm:single-float-digits extra-bias)
478 (if (minusp (float-sign x)) -1 1))))
479
480
481 ;;; INTEGER-DECODE-SINGLE-FLOAT -- Internal
482 ;;;
483 ;;; Handle the single-float case of INTEGER-DECODE-FLOAT. If an infinity or
484 ;;; NAN, error. If a denorm, call i-d-s-DENORM to handle it.
485 ;;;
486 (defun integer-decode-single-float (x)
487 (declare (single-float x))
488 (let* ((bits (single-float-bits (abs x)))
489 (exp (ldb vm:single-float-exponent-byte bits))
490 (sig (ldb vm:single-float-significand-byte bits))
491 (sign (if (minusp (float-sign x)) -1 1))
492 (biased (- exp vm:single-float-bias vm:single-float-digits)))
493 (declare (fixnum biased))
494 (unless (<= exp vm:single-float-normal-exponent-max)
495 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
496 (cond ((and (zerop exp) (zerop sig))
497 (values 0 biased sign))
498 ((< exp vm:single-float-normal-exponent-min)
499 (integer-decode-single-denorm x))
500 (t
501 (values (logior sig vm:single-float-hidden-bit) biased sign)))))
502
503
504 ;;; INTEGER-DECODE-DOUBLE-DENORM -- Internal
505 ;;;
506 ;;; Like INTEGER-DECODE-SINGLE-DENORM, only doubly so.
507 ;;;
508 (defun integer-decode-double-denorm (x)
509 (declare (type double-float x))
510 (let* ((high-bits (double-float-high-bits (abs x)))
511 (sig-high (ldb vm:double-float-significand-byte high-bits))
512 (low-bits (double-float-low-bits x))
513 (sign (if (minusp (float-sign x)) -1 1))
514 (biased (- (- vm:double-float-bias) vm:double-float-digits)))
515 (if (zerop sig-high)
516 (let ((sig low-bits)
517 (extra-bias (- vm:double-float-digits 33))
518 (bit (ash 1 31)))
519 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
520 (loop
521 (unless (zerop (logand sig bit)) (return))
522 (setq sig (ash sig 1))
523 (incf extra-bias))
524 (values (ash sig (- vm:double-float-digits 32))
525 (truly-the fixnum (- biased extra-bias))
526 sign))
527 (let ((sig (ash sig-high 1))
528 (extra-bias 0))
529 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
530 (loop
531 (unless (zerop (logand sig vm:double-float-hidden-bit))
532 (return))
533 (setq sig (ash sig 1))
534 (incf extra-bias))
535 (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
536 (truly-the fixnum (- biased extra-bias))
537 sign)))))
538
539
540 ;;; INTEGER-DECODE-DOUBLE-FLOAT -- Internal
541 ;;;
542 ;;; Like INTEGER-DECODE-SINGLE-FLOAT, only doubly so.
543 ;;;
544 (defun integer-decode-double-float (x)
545 (declare (double-float x))
546 (let* ((abs (abs x))
547 (hi (double-float-high-bits abs))
548 (lo (double-float-low-bits abs))
549 (exp (ldb vm:double-float-exponent-byte hi))
550 (sig (ldb vm:double-float-significand-byte hi))
551 (sign (if (minusp (float-sign x)) -1 1))
552 (biased (- exp vm:double-float-bias vm:double-float-digits)))
553 (declare (fixnum biased))
554 (unless (<= exp vm:double-float-normal-exponent-max)
555 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
556 (cond ((and (zerop exp) (zerop sig) (zerop lo))
557 (values 0 biased sign))
558 ((< exp vm:double-float-normal-exponent-min)
559 (integer-decode-double-denorm x))
560 (t
561 (values
562 (logior (ash (logior (ldb vm:double-float-significand-byte hi)
563 vm:double-float-hidden-bit)
564 32)
565 lo)
566 biased sign)))))
567
568
569 ;;; INTEGER-DECODE-LONG-DENORM -- Internal
570 ;;;
571 #+(and long-float x86)
572 (defun integer-decode-long-denorm (x)
573 (declare (type long-float x))
574 (let* ((high-bits (long-float-high-bits (abs x)))
575 (sig-high (ldb vm:long-float-significand-byte high-bits))
576 (low-bits (long-float-low-bits x))
577 (sign (if (minusp (float-sign x)) -1 1))
578 (biased (- (- vm:long-float-bias) vm:long-float-digits)))
579 (if (zerop sig-high)
580 (let ((sig low-bits)
581 (extra-bias (- vm:long-float-digits 33))
582 (bit (ash 1 31)))
583 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
584 (loop
585 (unless (zerop (logand sig bit)) (return))
586 (setq sig (ash sig 1))
587 (incf extra-bias))
588 (values (ash sig (- vm:long-float-digits 32))
589 (truly-the fixnum (- biased extra-bias))
590 sign))
591 (let ((sig (ash sig-high 1))
592 (extra-bias 0))
593 (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
594 (loop
595 (unless (zerop (logand sig vm:long-float-hidden-bit))
596 (return))
597 (setq sig (ash sig 1))
598 (incf extra-bias))
599 (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
600 (truly-the fixnum (- biased extra-bias))
601 sign)))))
602
603 #+double-double
604 (defun integer-decode-double-double-float (x)
605 (declare (type double-double-float x))
606 (multiple-value-bind (hi-int hi-exp sign)
607 (integer-decode-float (double-double-hi x))
608 (if (zerop (double-double-lo x))
609 (values (ash hi-int 53) (- hi-exp 53) sign)
610 (multiple-value-bind (lo-int lo-exp lo-sign)
611 (integer-decode-float (double-double-lo x))
612 ;; We have x = 2^e1*i1 + 2^e2*i2
613 ;; = 2^e2*(2^(e1-e2)*i1 + i2)
614 ;;
615 ;; NOTE: The hi and lo parts could actually have different
616 ;; signs, so we need to add the two parts together with the
617 ;; right sign!
618 (values (+ (* (* sign lo-sign) lo-int)
619 (ash hi-int (- hi-exp lo-exp)))
620 lo-exp
621 sign)))))
622
623 ;;; INTEGER-DECODE-LONG-FLOAT -- Internal
624 ;;;
625 #+(and long-float x86)
626 (defun integer-decode-long-float (x)
627 (declare (long-float x))
628 (let* ((hi (long-float-high-bits x))
629 (lo (long-float-low-bits x))
630 (exp-bits (long-float-exp-bits x))
631 (exp (ldb vm:long-float-exponent-byte exp-bits))
632 (sign (if (minusp exp-bits) -1 1))
633 (biased (- exp vm:long-float-bias vm:long-float-digits)))
634 (declare (fixnum biased))
635 (unless (<= exp vm:long-float-normal-exponent-max)
636 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
637 (cond ((and (zerop exp) (zerop hi) (zerop lo))
638 (values 0 biased sign))
639 ((< exp vm:long-float-normal-exponent-min)
640 (integer-decode-long-denorm x))
641 (t
642 (values (logior (ash hi 32) lo) biased sign)))))
643
644
645 ;;; INTEGER-DECODE-FLOAT -- Public
646 ;;;
647 ;;; Dispatch to the correct type-specific i-d-f function.
648 ;;;
649 (defun integer-decode-float (x)
650 "Returns three values:
651 1) an integer representation of the significand.
652 2) the exponent for the power of 2 that the significand must be multiplied
653 by to get the actual value. This differs from the DECODE-FLOAT exponent
654 by FLOAT-DIGITS, since the significand has been scaled to have all its
655 digits before the radix point.
656 3) -1 or 1 (i.e. the sign of the argument.)"
657 (number-dispatch ((x float))
658 ((single-float)
659 (integer-decode-single-float x))
660 ((double-float)
661 (integer-decode-double-float x))
662 #+long-float
663 ((long-float)
664 (integer-decode-long-float x))
665 #+double-double
666 ((double-double-float)
667 (integer-decode-double-double-float x))))
668
669
670 (declaim (maybe-inline decode-single-float decode-double-float))
671
672 ;;; DECODE-SINGLE-DENORM -- Internal
673 ;;;
674 ;;; Handle the denormalized case of DECODE-SINGLE-FLOAT. We call
675 ;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
676 ;;;
677 (defun decode-single-denorm (x)
678 (declare (type single-float x))
679 (multiple-value-bind (sig exp sign)
680 (integer-decode-single-denorm x)
681 (values (make-single-float
682 (dpb sig vm:single-float-significand-byte
683 (dpb vm:single-float-bias vm:single-float-exponent-byte 0)))
684 (truly-the fixnum (+ exp vm:single-float-digits))
685 (float sign x))))
686
687
688 ;;; DECODE-SINGLE-FLOAT -- Internal
689 ;;;
690 ;;; Handle the single-float case of DECODE-FLOAT. If an infinity or NAN,
691 ;;; error. If a denorm, call d-s-DENORM to handle it.
692 ;;;
693 (defun decode-single-float (x)
694 (declare (single-float x))
695 (let* ((bits (single-float-bits (abs x)))
696 (exp (ldb vm:single-float-exponent-byte bits))
697 (sign (float-sign x))
698 (biased (truly-the single-float-exponent
699 (- exp vm:single-float-bias))))
700 (unless (<= exp vm:single-float-normal-exponent-max)
701 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
702 (cond ((zerop x)
703 (values 0.0f0 biased sign))
704 ((< exp vm:single-float-normal-exponent-min)
705 (decode-single-denorm x))
706 (t
707 (values (make-single-float
708 (dpb vm:single-float-bias
709 vm:single-float-exponent-byte
710 bits))
711 biased sign)))))
712
713
714 ;;; DECODE-DOUBLE-DENORM -- Internal
715 ;;;
716 ;;; Like DECODE-SINGLE-DENORM, only doubly so.
717 ;;;
718 (defun decode-double-denorm (x)
719 (declare (double-float x))
720 (multiple-value-bind (sig exp sign)
721 (integer-decode-double-denorm x)
722 (values (make-double-float
723 (dpb (logand (ash sig -32) (lognot vm:double-float-hidden-bit))
724 vm:double-float-significand-byte
725 (dpb vm:double-float-bias vm:double-float-exponent-byte 0))
726 (ldb (byte 32 0) sig))
727 (truly-the fixnum (+ exp vm:double-float-digits))
728 (float sign x))))
729
730
731 ;;; DECODE-DOUBLE-FLOAT -- Public
732 ;;;
733 ;;; Like DECODE-SINGLE-FLOAT, only doubly so.
734 ;;;
735 (defun decode-double-float (x)
736 (declare (double-float x))
737 (let* ((abs (abs x))
738 (hi (double-float-high-bits abs))
739 (lo (double-float-low-bits abs))
740 (exp (ldb vm:double-float-exponent-byte hi))
741 (sign (float-sign x))
742 (biased (truly-the double-float-exponent
743 (- exp vm:double-float-bias))))
744 (unless (<= exp vm:double-float-normal-exponent-max)
745 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
746 (cond ((zerop x)
747 (values 0.0d0 biased sign))
748 ((< exp vm:double-float-normal-exponent-min)
749 (decode-double-denorm x))
750 (t
751 (values (make-double-float
752 (dpb vm:double-float-bias vm:double-float-exponent-byte hi)
753 lo)
754 biased sign)))))
755
756
757 ;;; DECODE-LONG-DENORM -- Internal
758 ;;;
759 #+(and long-float x86)
760 (defun decode-long-denorm (x)
761 (declare (long-float x))
762 (multiple-value-bind (sig exp sign)
763 (integer-decode-long-denorm x)
764 (values (make-long-float vm:long-float-bias (ash sig -32)
765 (ldb (byte 32 0) sig))
766 (truly-the fixnum (+ exp vm:long-float-digits))
767 (float sign x))))
768
769
770 ;;; DECODE-LONG-FLOAT -- Public
771 ;;;
772 #+(and long-float x86)
773 (defun decode-long-float (x)
774 (declare (long-float x))
775 (let* ((hi (long-float-high-bits x))
776 (lo (long-float-low-bits x))
777 (exp-bits (long-float-exp-bits x))
778 (exp (ldb vm:long-float-exponent-byte exp-bits))
779 (sign (if (minusp exp-bits) -1l0 1l0))
780 (biased (truly-the long-float-exponent (- exp vm:long-float-bias))))
781 (unless (<= exp vm:long-float-normal-exponent-max)
782 (error (intl:gettext "Can't decode NAN or infinity: ~S.") x))
783 (cond ((zerop x)
784 (values 0.0l0 biased sign))
785 ((< exp vm:long-float-normal-exponent-min)
786 (decode-long-denorm x))
787 (t
788 (values (make-long-float
789 (dpb vm:long-float-bias vm:long-float-exponent-byte
790 exp-bits)
791 hi
792 lo)
793 biased sign)))))
794
795 ;;; DECODE-DOUBLE-DOUBLE-FLOAT -- Public
796 #+double-double
797 (defun decode-double-double-float (x)
798 (declare (type double-double-float x))
799 (multiple-value-bind (hi-frac hi-exp sign)
800 (decode-float (double-double-hi x))
801 (values (make-double-double-float hi-frac
802 (scale-float (double-double-lo x) (- hi-exp)))
803 hi-exp
804 (coerce sign 'double-double-float))))
805
806 ;;; DECODE-FLOAT -- Public
807 ;;;
808 ;;; Dispatch to the appropriate type-specific function.
809 ;;;
810 (defun decode-float (f)
811 "Returns three values:
812 1) a floating-point number representing the significand. This is always
813 between 0.5 (inclusive) and 1.0 (exclusive).
814 2) an integer representing the exponent.
815 3) -1.0 or 1.0 (i.e. the sign of the argument.)"
816 (number-dispatch ((f float))
817 ((single-float)
818 (decode-single-float f))
819 ((double-float)
820 (decode-double-float f))
821 #+long-float
822 ((long-float)
823 (decode-long-float f))
824 #+double-double
825 ((double-double-float)
826 (decode-double-double-float f))))
827
828
829 ;;;; SCALE-FLOAT:
830
831 (declaim (maybe-inline scale-single-float scale-double-float))
832
833 ;;; SCALE-FLOAT-MAYBE-UNDERFLOW -- Internal
834 ;;;
835 ;;; Handle float scaling where the X is denormalized or the result is
836 ;;; denormalized or underflows to 0.
837 ;;;
838 (defun scale-float-maybe-underflow (x exp)
839 (multiple-value-bind (sig old-exp)
840 (integer-decode-float x)
841 (let* ((digits (float-digits x))
842 (new-exp (+ exp old-exp digits
843 (etypecase x
844 (single-float vm:single-float-bias)
845 (double-float vm:double-float-bias))))
846 (sign (if (minusp (float-sign x)) 1 0)))
847 (cond
848 ((< new-exp
849 (etypecase x
850 (single-float vm:single-float-normal-exponent-min)
851 (double-float vm:double-float-normal-exponent-min)))
852 (when (vm:current-float-trap :inexact)
853 (error 'floating-point-inexact :operation 'scale-float
854 :operands (list x exp)))
855 (when (vm:current-float-trap :underflow)
856 (error 'floating-point-underflow :operation 'scale-float
857 :operands (list x exp)))
858 (let ((shift (1- new-exp)))
859 ;; Is it necessary to have this IF here? Is there any case
860 ;; where (ash sig shift) won't return 0 when
861 ;; shift < -(digits-1)?
862 (if (< shift (- (1- digits)))
863 (etypecase x
864 (single-float (single-from-bits sign 0 0))
865 (double-float (double-from-bits sign 0 0)))
866 (etypecase x
867 (single-float (single-from-bits sign 0 (ash sig shift)))
868 (double-float (double-from-bits sign 0 (ash sig shift)))))))
869 (t
870 (etypecase x
871 (single-float (single-from-bits sign new-exp sig))
872 (double-float (double-from-bits sign new-exp sig))))))))
873
874
875 ;;; SCALE-FLOAT-MAYBE-OVERFLOW -- Internal
876 ;;;
877 ;;; Called when scaling a float overflows, or the oringinal float was a NaN
878 ;;; or infinity. If overflow errors are trapped, then error, otherwise return
879 ;;; the appropriate infinity. If a NaN, signal or not as appropriate.
880 ;;;
881 (defun scale-float-maybe-overflow (x exp)
882 (cond
883 ((float-infinity-p x)
884 ;; Infinity is infinity, no matter how small...
885 x)
886 ((float-nan-p x)
887 (when (and (float-trapping-nan-p x)
888 (vm:current-float-trap :invalid))
889 (error 'floating-point-invalid-operation :operation 'scale-float
890 :operands (list x exp)))
891 x)
892 (t
893 (when (vm:current-float-trap :overflow)
894 (error 'floating-point-overflow :operation 'scale-float
895 :operands (list x exp)))
896 (when (vm:current-float-trap :inexact)
897 (error 'floating-point-inexact :operation 'scale-float
898 :operands (list x exp)))
899 (* (float-sign x)
900 (etypecase x
901 (single-float single-float-positive-infinity)
902 (double-float double-float-positive-infinity))))))
903
904
905 ;;; SCALE-SINGLE-FLOAT, SCALE-DOUBLE-FLOAT -- Internal
906 ;;;
907 ;;; Scale a single or double float, calling the correct over/underflow
908 ;;; functions.
909 ;;;
910 (defun scale-single-float (x exp)
911 (declare (single-float x) (fixnum exp))
912 (let* ((bits (single-float-bits x))
913 (old-exp (ldb vm:single-float-exponent-byte bits))
914 (new-exp (+ old-exp exp)))
915 (cond
916 ((zerop x) x)
917 ((or (< old-exp vm:single-float-normal-exponent-min)
918 (< new-exp vm:single-float-normal-exponent-min))
919 (scale-float-maybe-underflow x exp))
920 ((or (> old-exp vm:single-float-normal-exponent-max)
921 (> new-exp vm:single-float-normal-exponent-max))
922 (scale-float-maybe-overflow x exp))
923 (t
924 (make-single-float (dpb new-exp vm:single-float-exponent-byte bits))))))
925 ;;;
926 (declaim (maybe-inline scale-double-float))
927 (defun scale-double-float (x exp)
928 (declare (double-float x) (fixnum exp))
929 (let* ((hi (double-float-high-bits x))
930 (lo (double-float-low-bits x))
931 (old-exp (ldb vm:double-float-exponent-byte hi))
932 (new-exp (+ old-exp exp)))
933 (cond
934 ((zerop x) x)
935 ((or (< old-exp vm:double-float-normal-exponent-min)
936 (< new-exp vm:double-float-normal-exponent-min))
937 (scale-float-maybe-underflow x exp))
938 ((or (> old-exp vm:double-float-normal-exponent-max)
939 (> new-exp vm:double-float-normal-exponent-max))
940 (scale-float-maybe-overflow x exp))
941 (t
942 (make-double-float (dpb new-exp vm:double-float-exponent-byte hi)
943 lo)))))
944
945 #+(and x86 long-float)
946 (defun scale-long-float (x exp)
947 (declare (long-float x) (fixnum exp))
948 (scale-float x exp))
949
950 #+double-double
951 (defun scale-double-double-float (x exp)
952 (declare (type double-double-float x) (fixnum exp))
953 (let ((hi (double-double-hi x))
954 (lo (double-double-lo x)))
955 (make-double-double-float (scale-double-float hi exp)
956 (scale-double-float lo exp))))
957
958 ;;; SCALE-FLOAT -- Public
959 ;;;
960 ;;; Dispatch to the correct type-specific scale-float function.
961 ;;;
962 (defun scale-float (f ex)
963 "Returns the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
964 of precision or overflow."
965 (number-dispatch ((f float))
966 ((single-float)
967 (scale-single-float f ex))
968 ((double-float)
969 (scale-double-float f ex))
970 #+long-float
971 ((long-float)
972 (scale-long-float f ex))
973 #+double-double
974 ((double-double-float)
975 (scale-double-double-float f ex))))
976
977
978 ;;;; Converting to/from floats:
979
980 (defun float (number &optional (other () otherp))
981 "Converts any REAL to a float. If OTHER is not provided, it returns a
982 SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the
983 result is the same float format as OTHER."
984 (if otherp
985 (number-dispatch ((number real) (other float))
986 (((foreach rational single-float double-float #+long-float long-float
987 #+double-double double-double-float)
988 (foreach single-float double-float #+long-float long-float
989 #+double-double double-double-float))
990 (coerce number '(dispatch-type other))))
991 (if (floatp number)
992 number
993 (coerce number 'single-float))))
994
995
996 (macrolet ((frob (name type)
997 `(defun ,name (x)
998 (number-dispatch ((x real))
999 (((foreach single-float double-float
1000 #+long-float long-float
1001 #+double-double double-double-float
1002 fixnum))
1003 (coerce x ',type))
1004 ((bignum)
1005 (bignum-to-float x ',type))
1006 ((ratio)
1007 (float-ratio x ',type))))))
1008 (frob %single-float single-float)
1009 (frob %double-float double-float)
1010 #+long-float
1011 (frob %long-float long-float)
1012 #+(and nil double-double)
1013 (frob %double-double-float double-double-float))
1014
1015
1016 ;; Note: There might be bootstrap issues when adding double-double
1017 ;; floats. Use the above definition if that's the case.
1018 #+double-double
1019 (defun %double-double-float (n)
1020 (typecase n
1021 (fixnum
1022 (%make-double-double-float (float n 1d0) 0d0))
1023 (single-float
1024 (%make-double-double-float (float n 1d0) 0d0))
1025 (double-float
1026 (%make-double-double-float (float n 1d0) 0d0))
1027 (double-double-float
1028 n)
1029 (bignum
1030 (bignum:bignum-to-float n 'double-double-float))
1031 (ratio
1032 (kernel::float-ratio n 'double-double-float))))
1033
1034 ;; Convert rational to double-double-float. The core algorithm is
1035 ;; from Richard Fateman. It was slightly modified to handle rationals
1036 ;; instead of just bignums.
1037 (defun rational-to-dd (rat)
1038 (declare (rational rat))
1039 (let* ((p (coerce rat 'double-float))
1040 (ans (%make-double-double-float p 0d0))
1041 (remainder rat))
1042 (declare (double-float p)
1043 (rational remainder)
1044 (type double-double-float ans))
1045 ;; Once or twice? I think once is enough.
1046 (dotimes (k 1 ans)
1047 (decf remainder (rational p))
1048 (setf p (coerce remainder 'double-float))
1049 (setf ans (+ ans p)))))
1050
1051 (defun float-ratio-float (x format)
1052 (let* ((signed-num (numerator x))
1053 (plusp (plusp signed-num))
1054 (num (if plusp signed-num (- signed-num)))
1055 (den (denominator x))
1056 (digits (float-format-digits format))
1057 (scale 0))
1058 (declare (fixnum digits scale))
1059 ;;
1060 ;; Strip any trailing zeros from the denominator and move it into the scale
1061 ;; factor (to minimize the size of the operands.)
1062 (let ((den-twos (1- (integer-length (logxor den (1- den))))))
1063 (declare (fixnum den-twos))
1064 (decf scale den-twos)
1065 (setq den (ash den (- den-twos))))
1066 ;;
1067 ;; Guess how much we need to scale by from the magnitudes of the numerator
1068 ;; and denominator. We want one extra bit for a guard bit.
1069 (let* ((num-len (integer-length num))
1070 (den-len (integer-length den))
1071 (delta (- den-len num-len))
1072 (shift (1+ (the fixnum (+ delta digits))))
1073 (shifted-num (ash num shift)))
1074 (declare (fixnum delta shift))
1075 (decf scale delta)
1076 (labels ((float-and-scale (bits)
1077 (let* ((bits (ash bits -1))
1078 (len (integer-length bits)))
1079 (cond ((> len digits)
1080 (assert (= len (the fixnum (1+ digits))))
1081 (multiple-value-bind (f0)
1082 (floatit (ash bits -1))
1083 #+nil
1084 (progn
1085 (format t "1: f0, f1 = ~A ~A~%" f0 f1)
1086 (format t " scale = ~A~%" (1+ scale)))
1087
1088 (scale-float f0 (1+ scale))))
1089 (t
1090 (multiple-value-bind (f0)
1091 (floatit bits)
1092 #+nil
1093 (progn
1094 (format t "2: f0, f1 = ~A ~A~%" f0 f1)
1095 (format t " scale = ~A~%" scale)
1096 (format t "scale-float f0 = ~A~%" (scale-float f0 scale))
1097 (when f1
1098 (format t "scale-float f1 = ~A~%"
1099 (scale-float f1 (- scale 53)))))
1100
1101 (scale-float f0 scale))))))
1102 (floatit (bits)
1103 (let ((sign (if plusp 0 1)))
1104 (case format
1105 (single-float
1106 (single-from-bits sign vm:single-float-bias bits))
1107 (double-float
1108 (double-from-bits sign vm:double-float-bias bits))
1109 #+long-float
1110 (long-float
1111 (long-from-bits sign vm:long-float-bias bits))))))
1112 (loop
1113 (multiple-value-bind (fraction-and-guard rem)
1114 (truncate shifted-num den)
1115 (let ((extra (- (integer-length fraction-and-guard) digits)))
1116 (declare (fixnum extra))
1117 (cond ((/= extra 1)
1118 (assert (> extra 1)))
1119 ((oddp fraction-and-guard)
1120 (return
1121 (if (zerop rem)
1122 (float-and-scale
1123 (if (zerop (logand fraction-and-guard 2))
1124 fraction-and-guard
1125 (1+ fraction-and-guard)))
1126 (float-and-scale (1+ fraction-and-guard)))))
1127 (t
1128 (return (float-and-scale fraction-and-guard)))))
1129 (setq shifted-num (ash shifted-num -1))
1130 (incf scale)))))))
1131
1132 ;;; FLOAT-RATIO -- Internal
1133 ;;;
1134 ;;; Convert a ratio to a float. We avoid any rounding error by doing an
1135 ;;; integer division. Accuracy is important to preserve read/print
1136 ;;; consistency, since this is ultimately how the reader reads a float. We
1137 ;;; scale the numerator by a power of two until the division results in the
1138 ;;; desired number of fraction bits, then do round-to-nearest.
1139 ;;;
1140 #+double-double
1141 (defun float-ratio (x format)
1142 (if (eq format 'double-double-float)
1143 (rational-to-dd x)
1144 (float-ratio-float x format)))
1145
1146 #-double-double
1147 (defun float-ratio (x format)
1148 (float-ratio-float x format))
1149
1150 #|
1151 These might be useful if we ever have a machine w/o float/integer conversion
1152 hardware. For now, we'll use special ops that uninterruptibly frob the
1153 rounding modes & do ieee round-to-integer.
1154
1155 ;;; %UNARY-TRUNCATE-SINGLE-FLOAT/FIXNUM -- Interface
1156 ;;;
1157 ;;; The compiler compiles a call to this when we are doing %UNARY-TRUNCATE
1158 ;;; and the result is known to be a fixnum. We can avoid some generic
1159 ;;; arithmetic in this case.
1160 ;;;
1161 (defun %unary-truncate-single-float/fixnum (x)
1162 (declare (single-float x) (values fixnum))
1163 (locally (declare (optimize (speed 3) (safety 0)))
1164 (let* ((bits (single-float-bits x))
1165 (exp (ldb vm:single-float-exponent-byte bits))
1166 (frac (logior (ldb vm:single-float-significand-byte bits)
1167 vm:single-float-hidden-bit))
1168 (shift (- exp vm:single-float-digits vm:single-float-bias)))
1169 (when (> exp vm:single-float-normal-exponent-max)
1170 (error 'floating-point-invalid-operation :operator 'truncate
1171 :operands (list x)))
1172 (if (<= shift (- vm:single-float-digits))
1173 0
1174 (let ((res (ash frac shift)))
1175 (declare (type (unsigned-byte 31) res))
1176 (if (minusp bits)
1177 (- res)
1178 res))))))
1179
1180
1181 ;;; %UNARY-TRUNCATE-DOUBLE-FLOAT/FIXNUM -- Interface
1182 ;;;
1183 ;;; Double-float version of this operation (see above single op).
1184 ;;;
1185 (defun %unary-truncate-double-float/fixnum (x)
1186 (declare (double-float x) (values fixnum))
1187 (locally (declare (optimize (speed 3) (safety 0)))
1188 (let* ((hi-bits (double-float-high-bits x))
1189 (exp (ldb vm:double-float-exponent-byte hi-bits))
1190 (frac (logior (ldb vm:double-float-significand-byte hi-bits)
1191 vm:double-float-hidden-bit))
1192 (shift (- exp (- vm:double-float-digits vm:word-bits)
1193 vm:double-float-bias)))
1194 (when (> exp vm:double-float-normal-exponent-max)
1195 (error 'floating-point-invalid-operation :operator 'truncate
1196 :operands (list x)))
1197 (if (<= shift (- vm:word-bits vm:double-float-digits))
1198 0
1199 (let* ((res-hi (ash frac shift))
1200 (res (if (plusp shift)
1201 (logior res-hi
1202 (the fixnum
1203 (ash (double-float-low-bits x)
1204 (- shift vm:word-bits))))
1205 res-hi)))
1206 (declare (type (unsigned-byte 31) res-hi res))
1207 (if (minusp hi-bits)
1208 (- res)
1209 res))))))
1210 |#
1211
1212
1213 ;;; %UNARY-TRUNCATE -- Interface
1214 ;;;
1215 ;;; This function is called when we are doing a truncate without any funky
1216 ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do
1217 ;;; *not* return the second value of truncate, so it must be computed by the
1218 ;;; caller if needed.
1219 ;;;
1220 ;;; In the float case, we pick off small arguments so that compiler can use
1221 ;;; special-case operations. We use an exclusive test, since (due to round-off
1222 ;;; error), (float most-positive-fixnum) may be greater than
1223 ;;; most-positive-fixnum.
1224 ;;;
1225 (defun %unary-truncate (number)
1226 (number-dispatch ((number real))
1227 ((integer) number)
1228 ((ratio) (values (truncate (numerator number) (denominator number))))
1229 (((foreach single-float double-float #+long-float long-float))
1230 (if (< (float most-negative-fixnum number)
1231 number
1232 (float most-positive-fixnum number))
1233 (truly-the fixnum (%unary-truncate number))
1234 (multiple-value-bind (bits exp)
1235 (integer-decode-float number)
1236 (let ((res (ash bits exp)))
1237 (if (minusp number)
1238 (- res)
1239 res)))))
1240 #+double-double
1241 ((double-double-float)
1242 (multiple-value-bind (bits exp)
1243 (integer-decode-double-double-float number)
1244 (let ((res (ash bits exp)))
1245 (if (minusp number)
1246 (- res)
1247 res))))))
1248
1249
1250 ;;; %UNARY-ROUND -- Interface
1251 ;;;
1252 ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we
1253 ;;; can't use the round primitive, then we do our own round-to-nearest on the
1254 ;;; result of i-d-f. [Note that this rounding will really only happen with
1255 ;;; double floats, since the whole single-float fraction will fit in a fixnum,
1256 ;;; so all single-floats larger than most-positive-fixnum can be precisely
1257 ;;; represented by an integer.]
1258 ;;;
1259 (defun %unary-round (number)
1260 (flet ((round-integer (bits exp sign)
1261 (let* ((shifted (ash bits exp))
1262 (roundup-p
1263 ;; Round if the are fraction bits (exp is
1264 ;; negative).
1265 (when (minusp exp)
1266 (let ((fraction (ldb (byte (- exp) 0) bits))
1267 (half (ash 1 (- -1 exp))))
1268 ;; If the fraction is less than half, then no
1269 ;; rounding. Otherwise, round up if the
1270 ;; fraction is greater than half or the
1271 ;; integer part is odd (for round-to-even).
1272 (cond ((> fraction half) t)
1273 ((< fraction half) nil)
1274 ((oddp shifted)
1275 t)))))
1276 (rounded (if roundup-p
1277 (1+ shifted)
1278 shifted)))
1279
1280 (if (minusp sign)
1281 (- rounded)
1282 rounded))))
1283 (number-dispatch ((number real))
1284 ((integer) number)
1285 ((ratio) (values (round (numerator number) (denominator number))))
1286 (((foreach single-float double-float #+long-float long-float))
1287 (if (< (float most-negative-fixnum number)
1288 number
1289 (float most-positive-fixnum number))
1290 (truly-the fixnum (%unary-round number))
1291 (multiple-value-bind (bits exp sign)
1292 (integer-decode-float number)
1293 (round-integer bits exp sign))))
1294 #+double-double
1295 ((double-double-float)
1296 (multiple-value-bind (bits exp sign)
1297 (integer-decode-float number)
1298 (round-integer bits exp sign))))))
1299
1300 (declaim (maybe-inline %unary-ftruncate/single-float
1301 %unary-ftruncate/double-float))
1302 ;; %UNARY-FTRUNCATE/SINGLE-FLOAT
1303 ;;
1304 ;; Basically the same as ftruncate, but specialized to handle only
1305 ;; single-floats and to only return the first value.
1306 (defun %unary-ftruncate/single-float (x)
1307 (declare (single-float x)
1308 (optimize (speed 3) (safety 0)))
1309 ;; We diddle the bits directly to truncate the number to reduce
1310 ;; consing.
1311 #+x87
1312 (let* ((bits (kernel:single-float-bits x))
1313 (exp (ldb vm:single-float-exponent-byte bits))
1314 (biased (truly-the kernel:single-float-exponent
1315 (- exp vm:single-float-bias))))
1316 (declare (type (signed-byte 32) bits))
1317 ;; At this point, we have the number represented as
1318 ;;
1319 ;; x = sign * 2^exp * frac
1320 ;;
1321 ;; where 0.5 <= frac < 1.0
1322 ;;
1323 ;; So if the exp <= 0, |x| < 1 and we know the result is 0.
1324 ;;
1325 ;; If exp >= (float-digits x), we know that all of the bits
1326 ;; in the fraction are integer bits, so the result is x.
1327 ;;
1328 ;; For everything else, some of the bits in frac are
1329 ;; fractional. Figure which ones they are and set them to
1330 ;; zero.
1331 ;;
1332 (cond ((> exp vm:single-float-normal-exponent-max)
1333 ;; Infinity or NaN. Convert NaN to quiet NaN.
1334 (let ((signif (ldb vm:single-float-significand-byte bits)))
1335 (if (zerop signif)
1336 ;; Infinity
1337 x
1338 ;; NaN. Convert to quiet Nan
1339 (make-single-float (logior bits vm:single-float-trapping-nan-bit)))))
1340 ((<= biased 0)
1341 ;; Number is less than 1. IEEE 754 says it should have the
1342 ;; same sign. Make it so.
1343 (* x 0f0))
1344 ((>= biased (float-digits x))
1345 ;; Number is an integer
1346 x)
1347 (t
1348 ;; Somewhere in between. Zap the bits that
1349 ;; represent the fraction part.
1350 (let ((frac-bits (- (float-digits x) biased)))
1351 (setf bits (logandc2 bits (- (ash 1 frac-bits) 1)))
1352 (kernel:make-single-float bits)))))
1353 ;; We have a fast fround, so use that and fix up the result to
1354 ;; truncate.
1355 #-x87
1356 (let ((r (%unary-fround/single-float x)))
1357 (if (> (abs r) (abs x))
1358 (if (> r 0)
1359 (- r 1)
1360 (+ r 1))
1361 r)))
1362
1363 ;; %UNARY-FTRUNCATE/DOUBLE-FLOAT
1364 ;;
1365 ;; Like %UNARY-FTRUNCATE/SINGLE-FLOAT, except for double-float.
1366
1367 (defun %unary-ftruncate/double-float (x)
1368 (declare (double-float x)
1369 (optimize (speed 3) (safety 0)))
1370 #+x87
1371 (let* ((hi (kernel:double-float-high-bits x))
1372 (lo (kernel:double-float-low-bits x))
1373 (exp (ldb vm:double-float-exponent-byte hi))
1374 (biased (truly-the kernel:double-float-exponent
1375 (- exp vm:double-float-bias))))
1376 (declare (type (signed-byte 32) hi)
1377 (type (unsigned-byte 32) lo))
1378 (cond ((> exp vm:double-float-normal-exponent-max)
1379 ;; Infinity or NaN. Convert NaN to quiet NaN.
1380 (if (and (zerop (ldb vm:double-float-significand-byte hi))
1381 (zerop lo))
1382 ;; Infinity
1383 x
1384 ;; NaN. Convert to quiet NaN
1385 (make-double-float (logior hi vm:double-float-trapping-nan-bit)
1386 lo)))
1387 ((<= biased 0)
1388 ;; Number is less than 1. IEEE 754 says it should have the
1389 ;; same sign. Make it so.
1390 (* x 0d0))
1391 ((>= biased (float-digits x))
1392 ;; Number is an integer
1393 x)
1394 (t
1395 ;; Somewhere in between. Zap the bits that
1396 ;; represent the fraction part.
1397 (let ((frac-bits (- (float-digits x) biased)))
1398 (cond ((< frac-bits 32)
1399 (setf lo (logandc2 lo (- (ash 1 frac-bits) 1))))
1400 (t
1401 (setf lo 0)
1402 (setf hi (logandc2 hi (- (ash 1 (- frac-bits 32)) 1)))))
1403 (kernel:make-double-float hi lo)))))
1404 ;; We have a fast fround, so use that and fix up the result to
1405 ;; truncate.
1406 #-x87
1407 (let ((r (%unary-fround/double-float x)))
1408 (if (> (abs r) (abs x))
1409 (if (> r 0)
1410 (- r 1)
1411 (+ r 1))
1412 r)))
1413
1414 #+double-double
1415 (defun %unary-ftruncate/double-double-float (x)
1416 ;; We should do something like what we do for double-double-float,
1417 ;; but I'm lazy right now.
1418 (coerce (truncate x) 'double-double-float))
1419
1420 ;;; %UNARY-FTRUNCATE -- Interface
1421 ;;;
1422 ;;; This function is called when we are doing an ftruncate without any
1423 ;;; funky divisor. Note that we do *not* return the second value of
1424 ;;; truncate, so it must be computed by the caller if needed.
1425 ;;;
1426 ;;; In the float case, we pick off small arguments so that compiler
1427 ;;; can use special-case operations. We use an exclusive test, since
1428 ;;; (due to round-off error), (float most-positive-signed-byte-32) may
1429 ;;; be greater than most-positive-signed-byte-32. We chose this
1430 ;;; because %unary-truncate can be optimized when the result is known
1431 ;;; to be a (signed-byte 32).
1432 ;;;
1433 (defun %unary-ftruncate (number)
1434 (number-dispatch ((number real))
1435 ((integer)
1436 (float number))
1437 ((ratio)
1438 (float (truncate (numerator number) (denominator number))))
1439 ((single-float)
1440 (%unary-ftruncate/single-float number))
1441 ((double-float)
1442 (%unary-ftruncate/double-float number))
1443 #+double-double
1444 ((double-double-float)
1445 (%unary-ftruncate/double-double-float number))))
1446
1447
1448 ;; %UNARY-FROUND not implemented for x87. I think there are potential
1449 ;; roundoff problems due to the 80-bit FPU registers.
1450 #-x87
1451 (progn
1452 (declaim (inline %unary-fround/single-float %unary-fround/double-float))
1453
1454 ;; %UNARY-FROUND/DOUBLE-FLOAT
1455 ;;
1456 ;; Basically the same as fround, but specialized to handle only
1457 ;; double-floats and to return the first value. This is the algorithm
1458 ;; given by Anton Ertl in comp.arch.arithmetic, Oct 26, 2002.
1459 (defun %unary-fround/double-float (x)
1460 (declare (double-float x))
1461 (let ((const (scale-float 1d0 53)))
1462 (cond ((> x 0)
1463 (+ (- x const) const))
1464 ((< x 0)
1465 (- (+ x const) const))
1466 (t
1467 x))))
1468
1469
1470 ;; %UNARY-FROUND/SINLGE-FLOAT
1471 ;;
1472 ;; Same as %UNARY-FROUND/DOUBLE-FLOAT, except specialized for single-float.
1473 (defun %unary-fround/single-float (x)
1474 (declare (single-float x))
1475 (let ((const (scale-float 1f0 24)))
1476 (cond ((> x 0)
1477 (+ (- x const) const))
1478 ((< x 0)
1479 (- (+ x const) const))
1480 (t
1481 x))))
1482
1483 ;;; %UNARY-FROUND -- Interface
1484 ;;;
1485 ;;; This function is called when we are doing an fround without any
1486 ;;; funky divisor. Note that we do *not* return the second value of
1487 ;;; round, so it must be computed by the caller if needed.
1488 ;;;
1489 (defun %unary-fround (number)
1490 (number-dispatch ((number real))
1491 ((integer)
1492 (float number))
1493 ((ratio)
1494 (float (round (numerator number) (denominator number))))
1495 ((single-float)
1496 (%unary-fround/single-float number))
1497 ((double-float)
1498 (%unary-fround/double-float number))))
1499
1500 ) ; not x87
1501
1502 (defun rational (x)
1503 "RATIONAL produces a rational number for any real numeric argument. This is
1504 more efficient than RATIONALIZE, but it assumes that floating-point is
1505 completely accurate, giving a result that isn't as pretty."
1506 (number-dispatch ((x real))
1507 (((foreach single-float double-float #+long-float long-float
1508 #+double-double double-double-float))
1509 (multiple-value-bind (bits exp)
1510 (integer-decode-float x)
1511 (if (eql bits 0)
1512 0
1513 (let* ((int (if (minusp x) (- bits) bits))
1514 (digits (float-digits x))
1515 (ex (+ exp digits)))
1516 (if (minusp ex)
1517 (integer-/-integer int (ash 1 (+ digits (- ex))))
1518 (integer-/-integer (ash int ex) (ash 1 digits)))))))
1519 ((rational) x)))
1520
1521
1522 #+nil
1523 (defun rationalize (x)
1524 "Converts any REAL to a RATIONAL. Floats are converted to a simple rational
1525 representation exploiting the assumption that floats are only accurate to
1526 their precision. RATIONALIZE (and also RATIONAL) preserve the invariant:
1527 (= x (float (rationalize x) x))"
1528 (number-dispatch ((x real))
1529 (((foreach single-float double-float #+long-float long-float))
1530 ;; Thanks to Kim Fateman, who stole this function rationalize-float
1531 ;; from macsyma's rational. Macsyma'a rationalize was written
1532 ;; by the legendary Gosper (rwg). Gosper is now working for Symbolics.
1533 ;; Guy Steele said about Gosper, "He has been called the
1534 ;; only living 17th century mathematician and is also the best
1535 ;; pdp-10 hacker I know." So, if you can understand or debug this
1536 ;; code you win big.
1537 (cond ((minusp x) (- (rationalize (- x))))
1538 ((zerop x) 0)
1539 (t
1540 (let ((eps (etypecase x
1541 (single-float single-float-epsilon)
1542 (double-float double-float-epsilon)
1543 #+long-float
1544 (long-float long-float-epsilon)))
1545 (y ())
1546 (a ()))
1547 (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
1548 (num (setq a (truncate x))
1549 (+ (* (setq a (truncate y)) num) onum))
1550 (den 1 (+ (* a den) oden))
1551 (onum 1 num)
1552 (oden 0 den))
1553 ((or (= xx (float a x))
1554 (and (not (zerop den))
1555 (not (> (abs (/ (- x (/ (float num x)
1556 (float den x)))
1557 x))
1558 eps))))
1559 (integer-/-integer num den))
1560 (declare ((dispatch-type x) xx)))))))
1561 ((rational) x)))
1562
1563 ;;; RATIONALIZE -- Public
1564 ;;;
1565 ;;; The algorithm here is the method described in CLISP. Bruno Haible has
1566 ;;; graciously given permission to use this algorithm. He says, "You can use
1567 ;;; it, if you present the following explanation of the algorithm."
1568 ;;;
1569 ;;; Algorithm (recursively presented):
1570 ;;; If x is a rational number, return x.
1571 ;;; If x = 0.0, return 0.
1572 ;;; If x < 0.0, return (- (rationalize (- x))).
1573 ;;; If x > 0.0:
1574 ;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
1575 ;;; exponent, sign).
1576 ;;; If m = 0 or e >= 0: return x = m*2^e.
1577 ;;; Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
1578 ;;; with smallest possible numerator and denominator.
1579 ;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
1580 ;;; But in this case the result will be x itself anyway, regardless of
1581 ;;; the choice of a. Therefore we can simply ignore this case.
1582 ;;; Note 2: At first, we need to consider the closed interval [a,b].
1583 ;;; but since a and b have the denominator 2^(|e|+1) whereas x itself
1584 ;;; has a denominator <= 2^|e|, we can restrict the seach to the open
1585 ;;; interval (a,b).
1586 ;;; So, for given a and b (0 < a < b) we are searching a rational number
1587 ;;; y with a <= y <= b.
1588 ;;; Recursive algorithm fraction_between(a,b):
1589 ;;; c := (ceiling a)
1590 ;;; if c < b
1591 ;;; then return c ; because a <= c < b, c integer
1592 ;;; else
1593 ;;; ; a is not integer (otherwise we would have had c = a < b)
1594 ;;; k := c-1 ; k = floor(a), k < a < b <= k+1
1595 ;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
1596 ;;; ; note 1 <= 1/(b-k) < 1/(a-k)
1597 ;;;
1598 ;;; You can see that we are actually computing a continued fraction expansion.
1599 ;;;
1600 ;;; Algorithm (iterative):
1601 ;;; If x is rational, return x.
1602 ;;; Call (integer-decode-float x). It returns a m,e,s (mantissa,
1603 ;;; exponent, sign).
1604 ;;; If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
1605 ;;; Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
1606 ;;; (positive and already in lowest terms because the denominator is a
1607 ;;; power of two and the numerator is odd).
1608 ;;; Start a continued fraction expansion
1609 ;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
1610 ;;; Loop
1611 ;;; c := (ceiling a)
1612 ;;; if c >= b
1613 ;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
1614 ;;; goto Loop
1615 ;;; finally partial_quotient(c).
1616 ;;; Here partial_quotient(c) denotes the iteration
1617 ;;; i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
1618 ;;; At the end, return s * (p[i]/q[i]).
1619 ;;; This rational number is already in lowest terms because
1620 ;;; p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
1621 ;;;
1622 (defun rationalize (x)
1623 "Converts any REAL to a RATIONAL. Floats are converted to a simple rational
1624 representation exploiting the assumption that floats are only accurate to
1625 their precision. RATIONALIZE (and also RATIONAL) preserve the invariant:
1626 (= x (float (rationalize x) x))"
1627 (number-dispatch ((x real))
1628 (((foreach single-float double-float #+long-float long-float
1629 #+double-double double-double-float))
1630 ;; This is a fairly straigtforward implementation of the iterative
1631 ;; algorithm above.
1632 (multiple-value-bind (frac expo sign)
1633 (integer-decode-float x)
1634 (cond ((or (zerop frac) (>= expo 0))
1635 (if (minusp sign)
1636 (- (ash frac expo))
1637 (ash frac expo)))
1638 (t
1639 ;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e),
1640 ;; so build the fraction up immediately, without having to do
1641 ;; a gcd.
1642 (let ((a (build-ratio (- (* 2 frac) 1) (ash 1 (- 1 expo))))
1643 (b (build-ratio (+ (* 2 frac) 1) (ash 1 (- 1 expo))))
1644 (p0 0)
1645 (q0 1)
1646 (p1 1)
1647 (q1 0))
1648 (do ((c (ceiling a) (ceiling a)))
1649 ((< c b)
1650 (let ((top (+ (* c p1) p0))
1651 (bot (+ (* c q1) q0)))
1652 (build-ratio (if (minusp sign)
1653 (- top)
1654 top)
1655 bot)))
1656 (let* ((k (- c 1))
1657 (p2 (+ (* k p1) p0))
1658 (q2 (+ (* k q1) q0)))
1659 (psetf a (/ (- b k))
1660 b (/ (- a k)))
1661 (setf p0 p1
1662 q0 q1
1663 p1 p2
1664 q1 q2))))))))
1665 ((rational) x)))
1666

  ViewVC Help
Powered by ViewVC 1.1.5