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 | ;;; |
---|
19 | ;;; level-0;l0-numbers.lisp |
---|
20 | |
---|
21 | (in-package "CCL") |
---|
22 | |
---|
23 | (eval-when (:compile-toplevel :execute) |
---|
24 | (require "ARCH") |
---|
25 | (require "LISPEQU") |
---|
26 | (require "NUMBER-MACROS") |
---|
27 | (require "NUMBER-CASE-MACRO") |
---|
28 | |
---|
29 | |
---|
30 | |
---|
31 | (defvar *dfloat-dops* '((* . %double-float*-2!)(/ . %double-float/-2!) |
---|
32 | (+ . %double-float+-2!)(- . %double-float--2!))) |
---|
33 | |
---|
34 | (defvar *sfloat-dops* '((* . %short-float*-2!)(/ . %short-float/-2!) |
---|
35 | (+ . %short-float+-2!)(- . %short-float--2!))) |
---|
36 | |
---|
37 | (defmacro dfloat-rat (op x y &optional (destructive-op (cdr (assq op *dfloat-dops*)))) |
---|
38 | (if destructive-op |
---|
39 | (let ((f2 (gensym))) |
---|
40 | `(let ((,f2 (%double-float ,y (%make-dfloat)))) |
---|
41 | (,destructive-op ,x ,f2 ,f2))) |
---|
42 | `(,op (the double-float ,x) (the double-float (%double-float ,y))))) |
---|
43 | |
---|
44 | (defmacro rat-dfloat (op x y &optional (destructive-op (cdr (assq op *dfloat-dops*)))) |
---|
45 | (if destructive-op |
---|
46 | (let ((f1 (gensym))) |
---|
47 | `(let ((,f1 (%double-float ,x (%make-dfloat)))) |
---|
48 | (,destructive-op ,f1 ,y ,f1))) |
---|
49 | `(,op (the double-float (%double-float ,x)) (the double-float ,y)))) |
---|
50 | |
---|
51 | (defmacro sfloat-rat (op x y &optional (destructive-op (cdr (assq op *sfloat-dops*)))) |
---|
52 | (let* ((use-destructive-op |
---|
53 | (target-word-size-case |
---|
54 | (32 destructive-op) |
---|
55 | (64 nil)))) |
---|
56 | (if use-destructive-op |
---|
57 | (let ((f2 (gensym))) |
---|
58 | `(let ((,f2 (%short-float ,y (%make-sfloat)))) |
---|
59 | (,destructive-op ,x ,f2 ,f2))) |
---|
60 | `(,op (the short-float ,x) (the short-float (%short-float ,y)))))) |
---|
61 | |
---|
62 | (defmacro rat-sfloat (op x y &optional (destructive-op (cdr (assq op *sfloat-dops*)))) |
---|
63 | (let* ((use-destructive-op |
---|
64 | (target-word-size-case |
---|
65 | (32 destructive-op) |
---|
66 | (64 nil)))) |
---|
67 | (if use-destructive-op |
---|
68 | (let ((f1 (gensym))) |
---|
69 | `(let ((,f1 (%short-float ,x (%make-sfloat)))) |
---|
70 | (,destructive-op ,f1 ,y ,f1))) |
---|
71 | `(,op (the short-float (%short-float ,x)) (the short-float ,y))))) |
---|
72 | |
---|
73 | |
---|
74 | |
---|
75 | |
---|
76 | |
---|
77 | (declaim (inline %make-complex %make-ratio)) |
---|
78 | (declaim (inline canonical-complex)) |
---|
79 | (declaim (inline build-ratio)) |
---|
80 | (declaim (inline maybe-truncate))) |
---|
81 | |
---|
82 | |
---|
83 | |
---|
84 | (defun %make-complex (realpart imagpart) |
---|
85 | (gvector :complex realpart imagpart)) |
---|
86 | |
---|
87 | (defun %make-ratio (numerator denominator) |
---|
88 | (gvector :ratio numerator denominator)) |
---|
89 | |
---|
90 | |
---|
91 | |
---|
92 | ; this is no longer used |
---|
93 | (defun %integer-signum (num) |
---|
94 | (if (fixnump num) |
---|
95 | (%fixnum-signum num) |
---|
96 | ; there is no such thing as bignum zero we hope |
---|
97 | (if (bignum-minusp num) -1 1))) |
---|
98 | |
---|
99 | |
---|
100 | ; Destructive primitives. |
---|
101 | (macrolet ((defdestructive-df-op (non-destructive-name destructive-name op) |
---|
102 | `(progn |
---|
103 | (defun ,non-destructive-name (x y) |
---|
104 | (,destructive-name x y (%make-dfloat))) |
---|
105 | (defun ,destructive-name (x y result) |
---|
106 | (declare (double-float x y result)) |
---|
107 | (%setf-double-float result (the double-float (,op x y))))))) |
---|
108 | (defdestructive-df-op %double-float+-2 %double-float+-2! +) |
---|
109 | (defdestructive-df-op %double-float--2 %double-float--2! -) |
---|
110 | (defdestructive-df-op %double-float*-2 %double-float*-2! *) |
---|
111 | (defdestructive-df-op %double-float/-2 %double-float/-2! /)) |
---|
112 | |
---|
113 | #-64-bit-target |
---|
114 | (macrolet ((defdestructive-sf-op (non-destructive-name destructive-name op) |
---|
115 | `(progn |
---|
116 | (defun ,non-destructive-name (x y) |
---|
117 | (,destructive-name x y (%make-sfloat))) |
---|
118 | (defun ,destructive-name (x y result) |
---|
119 | (declare (short-float x y result)) |
---|
120 | (%setf-short-float result (the short-float (,op x y))))))) |
---|
121 | (defdestructive-sf-op %short-float+-2 %short-float+-2! +) |
---|
122 | (defdestructive-sf-op %short-float--2 %short-float--2! -) |
---|
123 | (defdestructive-sf-op %short-float*-2 %short-float*-2! *) |
---|
124 | (defdestructive-sf-op %short-float/-2 %short-float/-2! /)) |
---|
125 | |
---|
126 | |
---|
127 | (defun %negate (x) |
---|
128 | (number-case x |
---|
129 | (fixnum (- (the fixnum x))) |
---|
130 | (double-float (%double-float-negate! x (%make-dfloat))) |
---|
131 | (short-float |
---|
132 | #+32-bit-target (%short-float-negate! x (%make-sfloat)) |
---|
133 | #+64-bit-target (%short-float-negate x)) |
---|
134 | (bignum (negate-bignum x)) |
---|
135 | (ratio (%make-ratio (%negate (%numerator x)) (%denominator x))) |
---|
136 | (complex (%make-complex (%negate (%realpart X))(%negate (%imagpart X))) ))) |
---|
137 | |
---|
138 | (defun %double-float-zerop (n) |
---|
139 | (zerop (the double-float n))) |
---|
140 | |
---|
141 | (defun %short-float-zerop (n) |
---|
142 | (zerop (the single-float n))) |
---|
143 | |
---|
144 | (defun zerop (number) |
---|
145 | "Is this number zero?" |
---|
146 | (number-case number |
---|
147 | (integer (eq number 0)) |
---|
148 | (short-float (%short-float-zerop number)) |
---|
149 | (double-float (%double-float-zerop number)) |
---|
150 | (ratio nil) |
---|
151 | (complex |
---|
152 | (number-case (%realpart number) |
---|
153 | (short-float (and (%short-float-zerop (%realpart number)) |
---|
154 | (%short-float-zerop (%imagpart number)))) |
---|
155 | (double-float (and (%double-float-zerop (%realpart number)) |
---|
156 | (%double-float-zerop (%imagpart number)))) |
---|
157 | (t (and (eql 0 (%realpart number))(eql 0 (%imagpart number)))))))) |
---|
158 | |
---|
159 | (defun %short-float-plusp (x) |
---|
160 | (> (the single-float x) 0.0f0)) |
---|
161 | |
---|
162 | (defun %double-float-plusp (x) |
---|
163 | (> (the double-float x) 0.0d0)) |
---|
164 | |
---|
165 | (defun plusp (number) |
---|
166 | "Is this real number strictly positive?" |
---|
167 | (number-case number |
---|
168 | (fixnum (%i> number 0)) |
---|
169 | (bignum (bignum-plusp number)) |
---|
170 | (short-float (%short-float-plusp number)) |
---|
171 | (double-float (%double-float-plusp number)) |
---|
172 | (ratio (plusp (%numerator number))))) |
---|
173 | |
---|
174 | |
---|
175 | (defun minusp (number) |
---|
176 | "Is this real number strictly negative?" |
---|
177 | (number-case number |
---|
178 | (fixnum (%i< number 0)) |
---|
179 | (bignum (bignum-minusp number)) |
---|
180 | (short-float (%short-float-minusp number)) |
---|
181 | (double-float (%double-float-minusp number)) |
---|
182 | (ratio (minusp (%numerator number))))) |
---|
183 | |
---|
184 | |
---|
185 | (defun oddp (n) |
---|
186 | "Is this integer odd?" |
---|
187 | (case (typecode n) |
---|
188 | (#.target::tag-fixnum (logbitp 0 (the fixnum n))) |
---|
189 | (#.target::subtag-bignum (%bignum-oddp n)) |
---|
190 | (t (report-bad-arg n 'integer)))) |
---|
191 | |
---|
192 | (defun evenp (n) |
---|
193 | "Is this integer even?" |
---|
194 | (case (typecode n) |
---|
195 | (#.target::tag-fixnum (not (logbitp 0 (the fixnum n)))) |
---|
196 | (#.target::subtag-bignum (not (%bignum-oddp n))) |
---|
197 | (t (report-bad-arg n 'integer)))) |
---|
198 | |
---|
199 | ;; expansion slightly changed |
---|
200 | (defun =-2 (x y) |
---|
201 | (number-case x |
---|
202 | (fixnum (number-case y |
---|
203 | (fixnum (eq x y)) |
---|
204 | (double-float (eq 0 (fixnum-dfloat-compare x y))) |
---|
205 | (short-float (eq 0 (fixnum-sfloat-compare x y))) |
---|
206 | ((bignum ratio) nil) |
---|
207 | (complex (and (zerop (%imagpart y)) (= x (%realpart y)))))) |
---|
208 | (double-float (number-case y ; x |
---|
209 | (double-float (= (the double-float x)(the double-float y))) ;x |
---|
210 | (short-float (with-stack-double-floats ((dy y)) |
---|
211 | (= (the double-float x) (the double-float dy)))) |
---|
212 | (fixnum (eq 0 (fixnum-dfloat-compare y x))) |
---|
213 | (bignum (eq 0 (bignum-dfloat-compare y x))) |
---|
214 | (ratio (= (rational x) y)) |
---|
215 | (complex (and (zerop (%imagpart y)) (= x (%realpart y)))))) |
---|
216 | (short-float (number-case y |
---|
217 | (short-float (= (the short-float x) (the short-float y))) |
---|
218 | (double-float (with-stack-double-floats ((dx x)) |
---|
219 | (= (the double-float dx) (the double-float y)))) |
---|
220 | (fixnum (eq 0 (fixnum-sfloat-compare y x))) |
---|
221 | (bignum (eq 0 (bignum-sfloat-compare y x))) |
---|
222 | (ratio (= (rational x) y)) |
---|
223 | (complex (and (zerop (%imagpart y)) (= x (%realpart y)))))) |
---|
224 | (bignum (number-case y |
---|
225 | (bignum (eq 0 (bignum-compare x y))) |
---|
226 | ((fixnum ratio) nil) |
---|
227 | (double-float (eq 0 (bignum-dfloat-compare x y))) |
---|
228 | (short-float (eq 0 (bignum-sfloat-compare x y))) |
---|
229 | (complex (and (zerop (%imagpart y)) (= x (%realpart y)))))) |
---|
230 | (ratio (number-case y |
---|
231 | (integer nil) |
---|
232 | (ratio |
---|
233 | (and (eql (%numerator x) (%numerator y)) |
---|
234 | (eql (%denominator x) (%denominator y)))) |
---|
235 | (float (= x (rational y))) |
---|
236 | (complex (and (zerop (%imagpart y)) (= x (%realpart y)))))) |
---|
237 | (complex (number-case y |
---|
238 | (complex (and (= (%realpart x) (%realpart y)) |
---|
239 | (= (%imagpart x) (%imagpart y)))) |
---|
240 | ((float rational) |
---|
241 | (and (zerop (%imagpart x)) (= (%realpart x) y))))))) |
---|
242 | |
---|
243 | (defun /=-2 (x y) |
---|
244 | (declare (notinline =-2)) |
---|
245 | (not (= x y))) |
---|
246 | |
---|
247 | |
---|
248 | ; true iff (< x y) is false. |
---|
249 | (defun >=-2 (x y) |
---|
250 | (declare (notinline <-2)) |
---|
251 | (not (< x y))) |
---|
252 | |
---|
253 | |
---|
254 | |
---|
255 | (defun <=-2 (x y) |
---|
256 | (declare (notinline >-2)) |
---|
257 | (not (> x y))) |
---|
258 | |
---|
259 | (defun <-2 (x y) |
---|
260 | (number-case x |
---|
261 | (fixnum (number-case y |
---|
262 | (fixnum (< (the fixnum x) (the fixnum y))) |
---|
263 | (double-float (eq -1 (fixnum-dfloat-compare x y))) |
---|
264 | (short-float (eq -1 (fixnum-sfloat-compare x y))) |
---|
265 | (bignum (bignum-plusp y)) |
---|
266 | (ratio (< x (ceiling (%numerator y)(%denominator y)))))) |
---|
267 | (double-float (number-case y ; x |
---|
268 | (double-float (< (the double-float x)(the double-float y))) ;x |
---|
269 | (short-float (with-stack-double-floats ((dy y)) |
---|
270 | (< (the double-float x) (the double-float dy)))) |
---|
271 | (fixnum (eq 1 (fixnum-dfloat-compare y x))) |
---|
272 | (bignum (eq 1 (bignum-dfloat-compare y x))) |
---|
273 | (ratio (< (rational x) y)))) |
---|
274 | (short-float (number-case y |
---|
275 | (short-float (< (the short-float x) (the short-float y))) |
---|
276 | (double-float (with-stack-double-floats ((dx x)) |
---|
277 | (< (the double-float dx) (the double-float y)))) |
---|
278 | (fixnum (eq 1 (fixnum-sfloat-compare y x))) |
---|
279 | (bignum (eq 1 (bignum-sfloat-compare y x))) |
---|
280 | (ratio (< (rational x) y)))) |
---|
281 | (bignum (number-case y |
---|
282 | (bignum (EQ -1 (bignum-compare x y))) |
---|
283 | (fixnum (not (bignum-plusp x))) |
---|
284 | (ratio (< x (ceiling (%numerator y)(%denominator y)))) |
---|
285 | (double-float (eq -1 (bignum-dfloat-compare x y))) |
---|
286 | (short-float (eq -1 (bignum-sfloat-compare x y))))) |
---|
287 | (ratio (number-case y |
---|
288 | (integer (< (floor (%numerator x)(%denominator x)) y)) |
---|
289 | (ratio |
---|
290 | (< (* (%numerator (the ratio x)) |
---|
291 | (%denominator (the ratio y))) |
---|
292 | (* (%numerator (the ratio y)) |
---|
293 | (%denominator (the ratio x))))) |
---|
294 | (float (< x (rational y))))))) |
---|
295 | |
---|
296 | |
---|
297 | |
---|
298 | (defun >-2 (x y) |
---|
299 | ;(declare (optimize (speed 3)(safety 0))) |
---|
300 | (number-case x |
---|
301 | (fixnum (number-case y |
---|
302 | (fixnum (> (the fixnum x) (the fixnum y))) |
---|
303 | (bignum (not (bignum-plusp y))) |
---|
304 | (double-float (eq 1 (fixnum-dfloat-compare x y))) |
---|
305 | (short-float (eq 1 (fixnum-sfloat-compare x y))) |
---|
306 | ;; or (> (* x denom) num) ? |
---|
307 | (ratio (> x (floor (%numerator y) (%denominator y)))))) |
---|
308 | (double-float (number-case y |
---|
309 | (double-float (> (the double-float x) (the double-float y))) |
---|
310 | (short-float (with-stack-double-floats ((dy y)) |
---|
311 | (> (the double-float x) (the double-float dy)))) |
---|
312 | (fixnum (eq -1 (fixnum-dfloat-compare y x))) |
---|
313 | (bignum (eq -1 (bignum-dfloat-compare y x))) |
---|
314 | (ratio (> (rational x) y)))) |
---|
315 | (short-float (number-case y |
---|
316 | (short-float (> (the short-float x) (the short-float y))) |
---|
317 | (double-float (with-stack-double-floats ((dx x)) |
---|
318 | (> (the double-float dx) (the double-float y)))) |
---|
319 | (fixnum (eq -1 (fixnum-sfloat-compare y x))) |
---|
320 | (bignum (eq -1 (bignum-sfloat-compare y x))) |
---|
321 | (ratio (> (rational x) y)))) |
---|
322 | (bignum (number-case y |
---|
323 | (fixnum (bignum-plusp x)) |
---|
324 | (bignum (eq 1 (bignum-compare x y))) |
---|
325 | ;; or (> (* x demon) num) |
---|
326 | (ratio (> x (floor (%numerator y) (%denominator y)))) |
---|
327 | (double-float (eq 1 (bignum-dfloat-compare x y))) |
---|
328 | (short-float (eq 1 (bignum-sfloat-compare x y))))) |
---|
329 | (ratio (number-case y |
---|
330 | ;; or (> num (* y denom)) |
---|
331 | (integer (> (ceiling (%numerator x) (%denominator x)) y)) |
---|
332 | (ratio |
---|
333 | (> (* (%numerator (the ratio x)) |
---|
334 | (%denominator (the ratio y))) |
---|
335 | (* (%numerator (the ratio y)) |
---|
336 | (%denominator (the ratio x))))) |
---|
337 | (float (> x (rational y))))))) |
---|
338 | |
---|
339 | |
---|
340 | ; t if any bits set after exp (unbiased) |
---|
341 | (defun hi-lo-fraction-p (hi lo exp) |
---|
342 | (declare (fixnum hi lo exp)) |
---|
343 | (if (> exp 24) |
---|
344 | (not (eql 0 (%ilogand lo (%ilsr (- exp 25) #xfffffff)))) |
---|
345 | (or (not (zerop lo))(not (eql 0 (%ilogand hi (%ilsr exp #x1ffffff))))))) |
---|
346 | |
---|
347 | |
---|
348 | |
---|
349 | (defun negate-hi-lo (hi lo) |
---|
350 | (setq hi (%ilogxor hi #x3ffffff)) |
---|
351 | (if (eq 0 lo) |
---|
352 | (setq hi (+ hi 1)) |
---|
353 | (setq lo (+ (%ilogxor lo #xfffffff) 1))) |
---|
354 | (values hi lo)) |
---|
355 | |
---|
356 | |
---|
357 | |
---|
358 | (defun fixnum-dfloat-compare (int dfloat) |
---|
359 | (declare (double-float dfloat) (fixnum int)) |
---|
360 | (if (and (eq int 0)(= dfloat 0.0d0)) |
---|
361 | 0 |
---|
362 | (with-stack-double-floats ((d1 int)) |
---|
363 | (locally (declare (double-float d1)) |
---|
364 | (if (eq int (%truncate-double-float->fixnum d1)) |
---|
365 | (cond ((< d1 dfloat) -1) |
---|
366 | ((= d1 dfloat) 0) |
---|
367 | (t 1)) |
---|
368 | ;; Whatever we do here should have the effect |
---|
369 | ;; of comparing the integer to the result of calling |
---|
370 | ;; RATIONAL on the float. We could probably |
---|
371 | ;; skip the call to RATIONAL in more cases, |
---|
372 | ;; but at least check the obvious ones here |
---|
373 | ;; (e.g. different signs) |
---|
374 | (multiple-value-bind (mantissa exponent sign) |
---|
375 | (integer-decode-double-float dfloat) |
---|
376 | (declare (type (integer -1 1) sign) |
---|
377 | (fixnum exponent)) |
---|
378 | (cond ((zerop int) |
---|
379 | (- sign)) |
---|
380 | ((and (< int 0) (eql sign 1)) -1) |
---|
381 | ((and (> int 0) (eql sign -1)) 1) |
---|
382 | (t |
---|
383 | ;; See RATIONAL. Can probably avoid this if |
---|
384 | ;; magnitudes are clearly dissimilar. |
---|
385 | (if (= sign -1) (setq mantissa (- mantissa))) |
---|
386 | (let* ((rat (if (< exponent 0) |
---|
387 | (/ mantissa (ash 1 (the fixnum (- exponent)))) |
---|
388 | (ash mantissa exponent)))) |
---|
389 | (if (< int rat) |
---|
390 | -1 |
---|
391 | (if (eq int rat) |
---|
392 | 0 |
---|
393 | 1))))))))))) |
---|
394 | |
---|
395 | |
---|
396 | |
---|
397 | (defun fixnum-sfloat-compare (int sfloat) |
---|
398 | (declare (short-float sfloat) (fixnum int)) |
---|
399 | (if (and (eq int 0)(= sfloat 0.0s0)) |
---|
400 | 0 |
---|
401 | (#+32-bit-target target::with-stack-short-floats #+32-bit-target ((s1 int)) |
---|
402 | #-32-bit-target let* #-32-bit-target ((s1 (%int-to-sfloat int))) |
---|
403 | (locally |
---|
404 | (declare (short-float s1)) |
---|
405 | (if (eq (%truncate-short-float->fixnum s1) int) |
---|
406 | (cond ((< s1 sfloat) -1) |
---|
407 | ((= s1 sfloat) 0) |
---|
408 | (t 1)) |
---|
409 | ;; Whatever we do here should have the effect |
---|
410 | ;; of comparing the integer to the result of calling |
---|
411 | ;; RATIONAL on the float. We could probably |
---|
412 | ;; skip the call to RATIONAL in more cases, |
---|
413 | ;; but at least check the obvious ones here |
---|
414 | ;; (e.g. different signs) |
---|
415 | (multiple-value-bind (mantissa exponent sign) |
---|
416 | (integer-decode-short-float sfloat) |
---|
417 | (declare (type (integer -1 1) sign) |
---|
418 | (fixnum exponent)) |
---|
419 | (cond ((zerop int) |
---|
420 | (- sign)) |
---|
421 | ((and (< int 0) (eql sign 1)) -1) |
---|
422 | ((and (> int 0) (eql sign -1)) 1) |
---|
423 | (t |
---|
424 | ;; See RATIONAL. Can probably avoid this if |
---|
425 | ;; magnitudes are clearly dissimilar. |
---|
426 | (if (= sign -1) (setq mantissa (- mantissa))) |
---|
427 | (let* ((rat (if (< exponent 0) |
---|
428 | (/ mantissa (ash 1 (the fixnum (- exponent)))) |
---|
429 | (ash mantissa exponent)))) |
---|
430 | (if (< int rat) |
---|
431 | -1 |
---|
432 | (if (eq int rat) |
---|
433 | 0 |
---|
434 | 1))))))))))) |
---|
435 | |
---|
436 | |
---|
437 | |
---|
438 | ;;; lotta stuff to avoid making a rational from a float |
---|
439 | ;;; returns -1 less, 0 equal, 1 greater |
---|
440 | (defun bignum-dfloat-compare (int float) |
---|
441 | (cond |
---|
442 | ((and (eq int 0)(= float 0.0d0)) 0) |
---|
443 | (t |
---|
444 | (let* ((fminus (%double-float-minusp float)) |
---|
445 | (iminus (minusp int)) |
---|
446 | (gt (if iminus -1 1))) |
---|
447 | (declare (fixnum gt)) |
---|
448 | (if (neq fminus iminus) |
---|
449 | gt ; if different signs, done |
---|
450 | (let ((intlen (integer-length int)) |
---|
451 | (exp (- (the fixnum (%double-float-exp float)) 1022))) |
---|
452 | (declare (fixnum intlen exp)) |
---|
453 | (cond |
---|
454 | ((and (not fminus) (< intlen exp)) -1) |
---|
455 | ((> intlen exp) gt) ; if different exp, done |
---|
456 | ((and fminus (or (< (1+ intlen) exp) |
---|
457 | (and (= (1+ intlen) exp) |
---|
458 | (neq (one-bignum-factor-of-two int) intlen)))) |
---|
459 | ;(print 'zow) |
---|
460 | (the fixnum (- gt))) ; ; integer-length is strange for neg powers of 2 |
---|
461 | (t (multiple-value-bind (hi lo)(fixnum-decode-float float) |
---|
462 | (declare (fixnum hi lo)) |
---|
463 | (when fminus (multiple-value-setq (hi lo)(negate-hi-lo hi lo))) |
---|
464 | (let* ((sz 26) ; exp > 28 always |
---|
465 | (pos (- exp 25)) |
---|
466 | (big-bits (%ldb-fixnum-from-bignum int sz pos))) |
---|
467 | (declare (fixnum pos big-bits sz)) |
---|
468 | ;(print (list big-bits hi sz pos)) |
---|
469 | (cond |
---|
470 | ((< big-bits hi) -1) |
---|
471 | ((> big-bits hi) 1) |
---|
472 | (t (let* ((sz (min (- exp 25) 28)) |
---|
473 | (pos (- exp 25 sz)) ; ? |
---|
474 | (ilo (if (< exp 53) (ash lo (- exp 53)) lo)) |
---|
475 | (big-bits (%ldb-fixnum-from-bignum int sz pos))) |
---|
476 | (declare (fixnum pos sz ilo big-bits)) |
---|
477 | ;(PRINT (list big-bits ilo)) |
---|
478 | (cond |
---|
479 | ((< big-bits ilo) -1) |
---|
480 | ((> big-bits ilo) 1) |
---|
481 | ((eq exp 53) 0) |
---|
482 | ((< exp 53) |
---|
483 | (if (not (hi-lo-fraction-p hi lo exp)) 0 -1)) ; -1 if pos |
---|
484 | (t (if (%i< (one-bignum-factor-of-two int) (- exp 53)) 1 0))))))) |
---|
485 | ))))))))) |
---|
486 | |
---|
487 | |
---|
488 | |
---|
489 | ;;; I don't know if it's worth doing a more "real" version of this. |
---|
490 | (defun bignum-sfloat-compare (int float) |
---|
491 | (with-stack-double-floats ((df float)) |
---|
492 | (bignum-dfloat-compare int df))) |
---|
493 | |
---|
494 | ;;;; Canonicalization utilities: |
---|
495 | |
---|
496 | ;;; CANONICAL-COMPLEX -- Internal |
---|
497 | ;;; |
---|
498 | ;;; If imagpart is 0, return realpart, otherwise make a complex. This is |
---|
499 | ;;; used when we know that realpart and imagpart are the same type, but |
---|
500 | ;;; rational canonicalization might still need to be done. |
---|
501 | ;;; |
---|
502 | |
---|
503 | (defun canonical-complex (realpart imagpart) |
---|
504 | (if (eql imagpart 0) |
---|
505 | realpart |
---|
506 | (%make-complex realpart imagpart))) |
---|
507 | |
---|
508 | |
---|
509 | |
---|
510 | |
---|
511 | (defun +-2 (x y) |
---|
512 | (number-case x |
---|
513 | (fixnum (number-case y |
---|
514 | (fixnum (+ (the fixnum x) (the fixnum y))) |
---|
515 | (double-float (rat-dfloat + x y)) |
---|
516 | (short-float (rat-sfloat + x y)) |
---|
517 | (bignum (add-bignum-and-fixnum y x)) |
---|
518 | (complex (complex (+ x (%realpart y)) |
---|
519 | (%imagpart y))) |
---|
520 | (ratio (let* ((dy (%denominator y)) |
---|
521 | (n (+ (* x dy) (%numerator y)))) |
---|
522 | (%make-ratio n dy))))) |
---|
523 | (double-float (number-case y |
---|
524 | (double-float (+ (the double-float x) (the double-float y))) |
---|
525 | (short-float (with-stack-double-floats ((dy y)) |
---|
526 | (+ (the double-float x) (the double-float dy)))) |
---|
527 | (rational (dfloat-rat + x y)) |
---|
528 | (complex (complex (+ x (%realpart y)) |
---|
529 | (%imagpart y))))) |
---|
530 | (short-float (number-case y |
---|
531 | (short-float (+ (the short-float x) (the short-float y))) |
---|
532 | (double-float (with-stack-double-floats ((dx x)) |
---|
533 | (+ (the double-float dx) (the double-float y)))) |
---|
534 | (rational (sfloat-rat + x y)) |
---|
535 | (complex (complex (+ x (%realpart y)) |
---|
536 | (%imagpart y))))) |
---|
537 | (bignum (number-case y |
---|
538 | (bignum (add-bignums x y)) |
---|
539 | (fixnum (add-bignum-and-fixnum x y)) |
---|
540 | (double-float (rat-dfloat + x y)) |
---|
541 | (short-float (rat-sfloat + x y)) |
---|
542 | (complex (complex (+ x (realpart y)) |
---|
543 | (%imagpart y))) |
---|
544 | (ratio |
---|
545 | (let* ((dy (%denominator y)) |
---|
546 | (n (+ (* x dy) (%numerator y)))) |
---|
547 | (%make-ratio n dy))))) |
---|
548 | (complex (number-case y |
---|
549 | (complex (canonical-complex (+ (%realpart x) (%realpart y)) |
---|
550 | (+ (%imagpart x) (%imagpart y)))) |
---|
551 | ((rational float) (complex (+ (%realpart x) y) (%imagpart x))))) |
---|
552 | (ratio (number-case y |
---|
553 | (ratio |
---|
554 | (let* ((nx (%numerator x)) |
---|
555 | (dx (%denominator x)) |
---|
556 | (ny (%numerator y)) |
---|
557 | (dy (%denominator y)) |
---|
558 | (g1 (gcd dx dy))) |
---|
559 | (if (eql g1 1) |
---|
560 | (%make-ratio (+ (* nx dy) (* dx ny)) (* dx dy)) |
---|
561 | (let* ((t1 (+ (* nx (truncate dy g1)) (* (truncate dx g1) ny))) |
---|
562 | (g2 (gcd t1 g1)) |
---|
563 | (t2 (truncate dx g1))) |
---|
564 | (cond ((eql t1 0) 0) |
---|
565 | ((eql g2 1) (%make-ratio t1 (* t2 dy))) |
---|
566 | (t |
---|
567 | (let* ((nn (truncate t1 g2)) |
---|
568 | (t3 (truncate dy g2)) |
---|
569 | (nd (if (eql t2 1) t3 (* t2 t3)))) |
---|
570 | (if (eql nd 1) nn (%make-ratio nn nd))))))))) |
---|
571 | (integer |
---|
572 | (let* ((dx (%denominator x)) (n (+ (%numerator x) (* y dx)))) |
---|
573 | (%make-ratio n dx))) |
---|
574 | (double-float (rat-dfloat + x y)) |
---|
575 | (short-float (rat-sfloat + x y)) |
---|
576 | (complex (complex (+ x (%realpart y)) |
---|
577 | (%imagpart y))))))) |
---|
578 | |
---|
579 | (defun --2 (x y) |
---|
580 | (number-case x |
---|
581 | (fixnum (number-case y |
---|
582 | (fixnum (- (the fixnum x) (the fixnum y))) |
---|
583 | (double-float (rat-dfloat - x y)) |
---|
584 | (short-float (rat-sfloat - x y)) |
---|
585 | (bignum |
---|
586 | (with-small-bignum-buffers ((bx x)) |
---|
587 | (subtract-bignum bx y))) |
---|
588 | (complex (complex (- x (%realpart y)) |
---|
589 | (- (%imagpart y)))) |
---|
590 | (ratio (let* ((dy (%denominator y)) |
---|
591 | (n (- (* x dy) (%numerator y)))) |
---|
592 | (%make-ratio n dy))))) |
---|
593 | (double-float (number-case y |
---|
594 | (double-float (- (the double-float x) (the double-float y))) |
---|
595 | (short-float (with-stack-double-floats ((dy y)) |
---|
596 | (- (the double-float x) (the double-float dy)))) |
---|
597 | (rational (dfloat-rat - x y)) |
---|
598 | (complex (complex (- x (%realpart y)) |
---|
599 | (- (%imagpart y)))))) |
---|
600 | (short-float (number-case y |
---|
601 | (short-float (- (the short-float x) (the short-float y))) |
---|
602 | (double-float (with-stack-double-floats ((dx x)) |
---|
603 | (- (the double-float dx) (the double-float y)))) |
---|
604 | (rational (sfloat-rat - x y)) |
---|
605 | (complex (complex (- x (%realpart y)) |
---|
606 | (- (%imagpart y)))))) |
---|
607 | (bignum (number-case y |
---|
608 | (bignum (subtract-bignum x y)) |
---|
609 | (fixnum (if (eql y target::target-most-negative-fixnum) |
---|
610 | (with-small-bignum-buffers ((by y)) |
---|
611 | (subtract-bignum x by)) |
---|
612 | (add-bignum-and-fixnum x (- y)))) |
---|
613 | (double-float (rat-dfloat - x y)) |
---|
614 | (short-float (rat-sfloat - x y)) |
---|
615 | (complex (complex (- x (realpart y)) |
---|
616 | (- (%imagpart y)))) |
---|
617 | (ratio |
---|
618 | (let* ((dy (%denominator y)) |
---|
619 | (n (- (* x dy) (%numerator y)))) |
---|
620 | (%make-ratio n dy))))) |
---|
621 | (complex (number-case y |
---|
622 | (complex (canonical-complex (- (%realpart x) (%realpart y)) |
---|
623 | (- (%imagpart x) (%imagpart y)))) |
---|
624 | ((rational float) (complex (- (%realpart x) y) (%imagpart x))))) |
---|
625 | (ratio (number-case y |
---|
626 | (ratio |
---|
627 | (let* ((nx (%numerator x)) |
---|
628 | (dx (%denominator x)) |
---|
629 | (ny (%numerator y)) |
---|
630 | (dy (%denominator y)) |
---|
631 | (g1 (gcd dx dy))) |
---|
632 | (if (eql g1 1) |
---|
633 | (%make-ratio (- (* nx dy) (* dx ny)) (* dx dy)) |
---|
634 | (let* ((t1 (- (* nx (truncate dy g1)) (* (truncate dx g1) ny))) |
---|
635 | (g2 (gcd t1 g1)) |
---|
636 | (t2 (truncate dx g1))) |
---|
637 | (cond ((eql t1 0) 0) |
---|
638 | ((eql g2 1) (%make-ratio t1 (* t2 dy))) |
---|
639 | (t |
---|
640 | (let* ((nn (truncate t1 g2)) |
---|
641 | (t3 (truncate dy g2)) |
---|
642 | (nd (if (eql t2 1) t3 (* t2 t3)))) |
---|
643 | (if (eql nd 1) nn (%make-ratio nn nd))))))))) |
---|
644 | (integer |
---|
645 | (let* ((dx (%denominator x)) (n (- (%numerator x) (* y dx)))) |
---|
646 | (%make-ratio n dx))) |
---|
647 | (double-float (rat-dfloat - x y)) |
---|
648 | (short-float (rat-sfloat - x y)) |
---|
649 | (complex (complex (- x (%realpart y)) |
---|
650 | (- (%imagpart y)))))))) |
---|
651 | |
---|
652 | |
---|
653 | ;;; BUILD-RATIO -- Internal |
---|
654 | ;;; |
---|
655 | ;;; Given a numerator and denominator with the GCD already divided out, make |
---|
656 | ;;; a canonical rational. We make the denominator positive, and check whether |
---|
657 | ;;; it is 1. |
---|
658 | ;;; |
---|
659 | |
---|
660 | (defun build-ratio (num den) |
---|
661 | (if (minusp den) (setq num (- num) den (- den))) |
---|
662 | (case den |
---|
663 | (0 (divide-by-zero-error 'build-ratio num den)) |
---|
664 | (1 num) |
---|
665 | (t (%make-ratio num den)))) |
---|
666 | |
---|
667 | |
---|
668 | |
---|
669 | |
---|
670 | ;;; MAYBE-TRUNCATE -- Internal |
---|
671 | ;;; |
---|
672 | ;;; Truncate X and Y, but bum the case where Y is 1. |
---|
673 | ;;; |
---|
674 | |
---|
675 | |
---|
676 | (defun maybe-truncate (x y) |
---|
677 | (if (eql y 1) |
---|
678 | x |
---|
679 | (truncate x y))) |
---|
680 | |
---|
681 | |
---|
682 | (defun *-2 (x y) |
---|
683 | ;(declare (optimize (speed 3)(safety 0))) |
---|
684 | (flet ((integer*ratio (x y) |
---|
685 | (if (eql x 0) 0 |
---|
686 | (let* ((ny (%numerator y)) |
---|
687 | (dy (%denominator y)) |
---|
688 | (gcd (gcd x dy))) |
---|
689 | (if (eql gcd 1) |
---|
690 | (%make-ratio (* x ny) dy) |
---|
691 | (let ((nn (* (truncate x gcd) ny)) |
---|
692 | (nd (truncate dy gcd))) |
---|
693 | (if (eql nd 1) |
---|
694 | nn |
---|
695 | (%make-ratio nn nd))))))) |
---|
696 | (complex*real (x y) |
---|
697 | (canonical-complex (* (%realpart x) y) (* (%imagpart x) y)))) |
---|
698 | (number-case x |
---|
699 | (double-float (number-case y |
---|
700 | (double-float (* (the double-float x)(the double-float y))) |
---|
701 | (short-float (with-stack-double-floats ((dy y)) |
---|
702 | (* (the double-float x) (the double-float dy)))) |
---|
703 | (rational (dfloat-rat * x y)) |
---|
704 | (complex (complex*real y x)))) |
---|
705 | (short-float (number-case y |
---|
706 | (double-float (with-stack-double-floats ((dx x)) |
---|
707 | (* (the double-float dx) (the double-float y)))) |
---|
708 | (short-float (* (the short-float x) (the short-float y))) |
---|
709 | (rational (sfloat-rat * x y)) |
---|
710 | (complex (complex*real y x)))) |
---|
711 | (bignum (number-case y |
---|
712 | (fixnum |
---|
713 | (if (eql y target::target-most-negative-fixnum) |
---|
714 | (with-small-bignum-buffers ((by y)) |
---|
715 | (multiply-bignums x by)) |
---|
716 | (multiply-bignum-and-fixnum x y))) |
---|
717 | (bignum (multiply-bignums x y)) |
---|
718 | (double-float (dfloat-rat * y x)) |
---|
719 | (short-float (sfloat-rat * y x)) |
---|
720 | (ratio (integer*ratio x y)) |
---|
721 | (complex (complex*real y x)))) |
---|
722 | (fixnum (number-case y |
---|
723 | (bignum (if (eql x target::target-most-negative-fixnum) |
---|
724 | (with-small-bignum-buffers ((bx x)) |
---|
725 | (multiply-bignums y bx)) |
---|
726 | (multiply-bignum-and-fixnum y x))) |
---|
727 | (fixnum (multiply-fixnums (the fixnum x) (the fixnum y))) |
---|
728 | (short-float (sfloat-rat * y x)) |
---|
729 | (double-float (dfloat-rat * y x)) |
---|
730 | (ratio (integer*ratio x y)) |
---|
731 | (complex (complex*real y x)))) |
---|
732 | (complex (number-case y |
---|
733 | (complex (let* ((rx (%realpart x)) |
---|
734 | (ix (%imagpart x)) |
---|
735 | (ry (%realpart y)) |
---|
736 | (iy (%imagpart y))) |
---|
737 | (canonical-complex (- (* rx ry) (* ix iy)) (+ (* rx iy) (* ix ry))))) |
---|
738 | (real (complex*real x y)))) |
---|
739 | (ratio (number-case y |
---|
740 | (ratio (let* ((nx (%numerator x)) |
---|
741 | (dx (%denominator x)) |
---|
742 | (ny (%numerator y)) |
---|
743 | (dy (%denominator y)) |
---|
744 | (g1 (gcd nx dy)) |
---|
745 | (g2 (gcd dx ny))) |
---|
746 | (build-ratio (* (maybe-truncate nx g1) |
---|
747 | (maybe-truncate ny g2)) |
---|
748 | (* (maybe-truncate dx g2) |
---|
749 | (maybe-truncate dy g1))))) |
---|
750 | (integer (integer*ratio y x)) |
---|
751 | (double-float (rat-dfloat * x y)) |
---|
752 | (short-float (rat-sfloat * x y)) |
---|
753 | (complex (complex*real y x))))))) |
---|
754 | |
---|
755 | |
---|
756 | |
---|
757 | (defun integer*integer (x y &optional res) |
---|
758 | (declare (ignore res)) |
---|
759 | (number-case x |
---|
760 | (fixnum (number-case y |
---|
761 | (fixnum (* (the fixnum x) (the fixnum y))) |
---|
762 | (t (multiply-bignum-and-fixnum y x)))) |
---|
763 | (bignum (number-case y |
---|
764 | (fixnum (multiply-bignum-and-fixnum x y)) |
---|
765 | (t (multiply-bignums x y)))))) |
---|
766 | |
---|
767 | |
---|
768 | |
---|
769 | |
---|
770 | |
---|
771 | ;;; INTEGER-/-INTEGER -- Internal |
---|
772 | ;;; |
---|
773 | ;;; Divide two integers, producing a canonical rational. If a fixnum, we |
---|
774 | ;;; see if they divide evenly before trying the GCD. In the bignum case, we |
---|
775 | ;;; don't bother, since bignum division is expensive, and the test is not very |
---|
776 | ;;; likely to suceed. |
---|
777 | ;;; |
---|
778 | (defun integer-/-integer (x y) |
---|
779 | (if (and (typep x 'fixnum) (typep y 'fixnum)) |
---|
780 | (multiple-value-bind (quo rem) (%fixnum-truncate x y) |
---|
781 | (if (eql 0 rem) |
---|
782 | quo |
---|
783 | (let ((gcd (gcd x y))) |
---|
784 | (declare (fixnum gcd)) |
---|
785 | (if (eql gcd 1) |
---|
786 | (build-ratio x y) |
---|
787 | (build-ratio (%fixnum-truncate x gcd) (%fixnum-truncate y gcd)))))) |
---|
788 | (let ((gcd (gcd x y))) |
---|
789 | (if (eql gcd 1) |
---|
790 | (build-ratio x y) |
---|
791 | (build-ratio (truncate x gcd) (truncate y gcd)))))) |
---|
792 | |
---|
793 | |
---|
794 | |
---|
795 | (defun /-2 (x y) |
---|
796 | (macrolet ((real-complex-/ (x y) |
---|
797 | (let ((ry (gensym)) |
---|
798 | (iy (gensym)) |
---|
799 | (r (gensym)) |
---|
800 | (dn (gensym))) |
---|
801 | `(let* ((,ry (%realpart ,y)) |
---|
802 | (,iy (%imagpart ,y))) |
---|
803 | (if (> (abs ,ry) (abs ,iy)) |
---|
804 | (let* ((,r (/ ,iy ,ry)) |
---|
805 | (,dn (* ,ry (+ 1 (* ,r ,r))))) |
---|
806 | (canonical-complex (/ ,x ,dn) |
---|
807 | (/ (- (* ,x ,r)) ,dn))) |
---|
808 | (let* ((,r (/ ,ry ,iy)) |
---|
809 | (,dn (* ,iy (+ 1 (* ,r ,r))))) |
---|
810 | (canonical-complex (/ (* ,x ,r) ,dn) |
---|
811 | (/ (- ,x) ,dn)))))))) |
---|
812 | (number-case x |
---|
813 | (double-float (number-case y |
---|
814 | (double-float (/ (the double-float x) (the double-float y))) |
---|
815 | (short-float (with-stack-double-floats ((dy y)) |
---|
816 | (/ (the double-float x) |
---|
817 | (the double-float dy)))) |
---|
818 | (rational (dfloat-rat / x y)) |
---|
819 | (complex (real-complex-/ x y)))) |
---|
820 | (short-float (number-case y |
---|
821 | (short-float (/ (the short-float x) (the short-float y))) |
---|
822 | (double-float (with-stack-double-floats ((dx x)) |
---|
823 | (/ (the double-float dx) |
---|
824 | (the double-float y)))) |
---|
825 | (rational (sfloat-rat / x y)) |
---|
826 | (complex (real-complex-/ x y)))) |
---|
827 | (integer (number-case y |
---|
828 | (double-float (rat-dfloat / x y)) |
---|
829 | (short-float (rat-sfloat / x y)) |
---|
830 | (integer (integer-/-integer x y)) |
---|
831 | (complex (real-complex-/ x y)) |
---|
832 | (ratio |
---|
833 | (if (eql 0 x) |
---|
834 | 0 |
---|
835 | (let* ((ny (%numerator y)) |
---|
836 | (dy (%denominator y)) |
---|
837 | (gcd (gcd x ny))) |
---|
838 | (build-ratio (* (maybe-truncate x gcd) dy) |
---|
839 | (maybe-truncate ny gcd))))))) |
---|
840 | (complex (number-case y |
---|
841 | (complex (let* ((rx (%realpart x)) |
---|
842 | (ix (%imagpart x)) |
---|
843 | (ry (%realpart y)) |
---|
844 | (iy (%imagpart y))) |
---|
845 | (if (> (abs ry) (abs iy)) |
---|
846 | (let* ((r (/ iy ry)) |
---|
847 | (dn (+ ry (* r iy)))) |
---|
848 | (canonical-complex (/ (+ rx (* ix r)) dn) |
---|
849 | (/ (- ix (* rx r)) dn))) |
---|
850 | (let* ((r (/ ry iy)) |
---|
851 | (dn (+ iy (* r ry)))) |
---|
852 | (canonical-complex (/ (+ (* rx r) ix) dn) |
---|
853 | (/ (- (* ix r) rx) dn)))))) |
---|
854 | ((rational float) |
---|
855 | (canonical-complex (/ (%realpart x) y) (/ (%imagpart x) y))))) |
---|
856 | (ratio (number-case y |
---|
857 | (double-float (rat-dfloat / x y)) |
---|
858 | (short-float (rat-sfloat / x y)) |
---|
859 | (integer |
---|
860 | (when (eql y 0) |
---|
861 | (divide-by-zero-error '/ x y)) |
---|
862 | (let* ((nx (%numerator x)) (gcd (gcd nx y))) |
---|
863 | (build-ratio (maybe-truncate nx gcd) |
---|
864 | (* (maybe-truncate y gcd) (%denominator x))))) |
---|
865 | (complex (real-complex-/ x y)) |
---|
866 | (ratio |
---|
867 | (let* ((nx (%numerator x)) |
---|
868 | (dx (%denominator x)) |
---|
869 | (ny (%numerator y)) |
---|
870 | (dy (%denominator y)) |
---|
871 | (g1 (gcd nx ny)) |
---|
872 | (g2 (gcd dx dy))) |
---|
873 | (build-ratio (* (maybe-truncate nx g1) |
---|
874 | (maybe-truncate dy g2)) |
---|
875 | (* (maybe-truncate dx g2) |
---|
876 | (maybe-truncate ny g1)))))))))) |
---|
877 | |
---|
878 | |
---|
879 | |
---|
880 | (defun divide-by-zero-error (operation &rest operands) |
---|
881 | (error (make-condition 'division-by-zero |
---|
882 | :operation operation |
---|
883 | :operands operands))) |
---|
884 | |
---|
885 | |
---|
886 | (defun 1+ (number) |
---|
887 | "Returns NUMBER + 1." |
---|
888 | (+-2 number 1)) |
---|
889 | |
---|
890 | (defun 1- (number) |
---|
891 | "Returns NUMBER - 1." |
---|
892 | (--2 number 1)) |
---|
893 | |
---|
894 | |
---|
895 | |
---|
896 | |
---|
897 | (defun conjugate (number) |
---|
898 | "Return the complex conjugate of NUMBER. For non-complex numbers, this is |
---|
899 | an identity." |
---|
900 | (number-case number |
---|
901 | (complex (complex (%realpart number) (- (%imagpart number)))) |
---|
902 | (number number))) |
---|
903 | |
---|
904 | (defun numerator (rational) |
---|
905 | "Return the numerator of NUMBER, which must be rational." |
---|
906 | (number-case rational |
---|
907 | (ratio (%numerator rational)) |
---|
908 | (integer rational))) |
---|
909 | |
---|
910 | (defun denominator (rational) |
---|
911 | "Return the denominator of NUMBER, which must be rational." |
---|
912 | (number-case rational |
---|
913 | (ratio (%denominator rational)) |
---|
914 | (integer 1))) |
---|
915 | |
---|
916 | |
---|
917 | |
---|
918 | (defun abs (number) |
---|
919 | "Return the absolute value of the number." |
---|
920 | (number-case number |
---|
921 | (fixnum |
---|
922 | (locally (declare (fixnum number)) |
---|
923 | (if (minusp number) (- number) number))) |
---|
924 | (double-float |
---|
925 | (%double-float-abs number)) |
---|
926 | (short-float |
---|
927 | (%short-float-abs number)) |
---|
928 | (bignum |
---|
929 | (if (bignum-minusp number)(negate-bignum number) number)) |
---|
930 | (ratio |
---|
931 | (if (minusp number) (- number) number)) |
---|
932 | (complex |
---|
933 | (let ((rx (%realpart number)) |
---|
934 | (ix (%imagpart number))) |
---|
935 | (number-case rx |
---|
936 | (rational |
---|
937 | (sqrt (+ (* rx rx) (* ix ix)))) |
---|
938 | (short-float |
---|
939 | (%short-float (%double-float-hypot (%double-float rx) |
---|
940 | (%double-float ix)))) |
---|
941 | (double-float |
---|
942 | (%double-float-hypot rx ix))))))) |
---|
943 | |
---|
944 | |
---|
945 | |
---|
946 | (defun phase (number) |
---|
947 | "Return the angle part of the polar representation of a complex number. |
---|
948 | For complex numbers, this is (atan (imagpart number) (realpart number)). |
---|
949 | For non-complex positive numbers, this is 0. For non-complex negative |
---|
950 | numbers this is PI." |
---|
951 | (number-case number |
---|
952 | (rational |
---|
953 | (if (minusp number) |
---|
954 | (%short-float pi) |
---|
955 | 0.0f0)) |
---|
956 | (double-float |
---|
957 | (if (minusp number) |
---|
958 | (%double-float pi) |
---|
959 | 0.0d0)) |
---|
960 | (complex |
---|
961 | (atan (%imagpart number) (%realpart number))) |
---|
962 | (short-float |
---|
963 | (if (minusp number) |
---|
964 | (%short-float pi) |
---|
965 | 0.0s0)))) |
---|
966 | |
---|
967 | |
---|
968 | |
---|
969 | ; from Lib;numbers.lisp, sort of |
---|
970 | (defun float (number &optional other) |
---|
971 | "Converts any REAL to a float. If OTHER is not provided, it returns a |
---|
972 | SINGLE-FLOAT if NUMBER is not already a FLOAT. If OTHER is provided, the |
---|
973 | result is the same float format as OTHER." |
---|
974 | (if (null other) |
---|
975 | (if (typep number 'float) |
---|
976 | number |
---|
977 | (%short-float number)) |
---|
978 | (if (typep other 'double-float) |
---|
979 | (%double-float number) |
---|
980 | (if (typep other 'short-float) |
---|
981 | (%short-float number) |
---|
982 | (float number (require-type other 'float)))))) |
---|
983 | |
---|
984 | |
---|
985 | |
---|
986 | |
---|
987 | |
---|
988 | ;;; If the numbers do not divide exactly and the result of (/ number divisor) |
---|
989 | ;;; would be negative then decrement the quotient and augment the remainder by |
---|
990 | ;;; the divisor. |
---|
991 | ;;; |
---|
992 | (defun floor (number &optional divisor) |
---|
993 | "Return the greatest integer not greater than number, or number/divisor. |
---|
994 | The second returned value is (mod number divisor)." |
---|
995 | (if (null divisor)(setq divisor 1)) |
---|
996 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
997 | (if (and (not (zerop rem)) |
---|
998 | (if (minusp divisor) |
---|
999 | (plusp number) |
---|
1000 | (minusp number))) |
---|
1001 | (if (called-for-mv-p) |
---|
1002 | (values (1- tru) (+ rem divisor)) |
---|
1003 | (1- tru)) |
---|
1004 | (values tru rem)))) |
---|
1005 | |
---|
1006 | |
---|
1007 | |
---|
1008 | (defun %fixnum-floor (number divisor) |
---|
1009 | (declare (fixnum number divisor)) |
---|
1010 | (if (eq divisor 1) |
---|
1011 | (values number 0) |
---|
1012 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
1013 | (if (eq rem 0) |
---|
1014 | (values tru 0) |
---|
1015 | (locally (declare (fixnum tru rem)) |
---|
1016 | (if (and ;(not (zerop rem)) |
---|
1017 | (if (minusp divisor) |
---|
1018 | (plusp number) |
---|
1019 | (minusp number))) |
---|
1020 | (values (the fixnum (1- tru)) (the fixnum (+ rem divisor))) |
---|
1021 | (values tru rem))))))) |
---|
1022 | |
---|
1023 | |
---|
1024 | |
---|
1025 | ;;; If the numbers do not divide exactly and the result of (/ number divisor) |
---|
1026 | ;;; would be positive then increment the quotient and decrement the remainder by |
---|
1027 | ;;; the divisor. |
---|
1028 | ;;; |
---|
1029 | (defun ceiling (number &optional divisor) |
---|
1030 | "Return the smallest integer not less than number, or number/divisor. |
---|
1031 | The second returned value is the remainder." |
---|
1032 | (if (null divisor)(setq divisor 1)) |
---|
1033 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
1034 | (if (and (not (zerop rem)) |
---|
1035 | (if (minusp divisor) |
---|
1036 | (minusp number) |
---|
1037 | (plusp number))) |
---|
1038 | (if (called-for-mv-p) |
---|
1039 | (values (+ tru 1) (- rem divisor)) |
---|
1040 | (+ tru 1)) |
---|
1041 | (values tru rem)))) |
---|
1042 | |
---|
1043 | |
---|
1044 | |
---|
1045 | (defun %fixnum-ceiling (number divisor) |
---|
1046 | "Returns the smallest integer not less than number, or number/divisor. |
---|
1047 | The second returned value is the remainder." |
---|
1048 | (declare (fixnum number divisor)) |
---|
1049 | (multiple-value-bind (tru rem) (%fixnum-truncate number divisor) |
---|
1050 | (if (eq 0 rem) |
---|
1051 | (values tru 0) |
---|
1052 | (locally (declare (fixnum tru rem)) |
---|
1053 | (if (and ;(not (zerop rem)) |
---|
1054 | (if (minusp divisor) |
---|
1055 | (minusp number) |
---|
1056 | (plusp number))) |
---|
1057 | (values (the fixnum (+ tru 1))(the fixnum (- rem divisor))) |
---|
1058 | (values tru rem)))))) |
---|
1059 | |
---|
1060 | |
---|
1061 | |
---|
1062 | (defun integer-decode-denorm-short-float (mantissa sign) |
---|
1063 | (declare (fixnum mantissa sign)) |
---|
1064 | (do* ((bias 0 (1+ bias)) |
---|
1065 | (sig mantissa (ash sig 1))) |
---|
1066 | ((logbitp 23 sig) |
---|
1067 | (values sig |
---|
1068 | (- (- IEEE-single-float-bias) |
---|
1069 | IEEE-single-float-digits |
---|
1070 | bias) |
---|
1071 | sign)))) |
---|
1072 | |
---|
1073 | |
---|
1074 | (defun integer-decode-short-float (sfloat) |
---|
1075 | (multiple-value-bind (mantissa exp sign)(fixnum-decode-short-float sfloat) |
---|
1076 | (let* ((biased (- exp IEEE-single-float-bias IEEE-single-float-digits))) |
---|
1077 | (setq sign (if (eql 0 sign) 1 -1)) |
---|
1078 | (if (eq exp 255) |
---|
1079 | (error "Can't decode NAN/Inf: ~s" sfloat)) |
---|
1080 | (if (eql 0 exp) |
---|
1081 | (if (eql 0 mantissa) |
---|
1082 | (values 0 biased sign) |
---|
1083 | (integer-decode-denorm-short-float (ash mantissa 1) sign)) |
---|
1084 | (values (logior #x800000 mantissa) biased sign))))) |
---|
1085 | |
---|
1086 | |
---|
1087 | |
---|
1088 | |
---|
1089 | ;;; INTEGER-DECODE-FLOAT -- Public |
---|
1090 | ;;; |
---|
1091 | ;;; Dispatch to the correct type-specific i-d-f function. |
---|
1092 | ;;; |
---|
1093 | (defun integer-decode-float (x) |
---|
1094 | "Returns three values: |
---|
1095 | 1) an integer representation of the significand. |
---|
1096 | 2) the exponent for the power of 2 that the significand must be multiplied |
---|
1097 | by to get the actual value. This differs from the DECODE-FLOAT exponent |
---|
1098 | by FLOAT-DIGITS, since the significand has been scaled to have all its |
---|
1099 | digits before the radix point. |
---|
1100 | 3) -1 or 1 (i.e. the sign of the argument.)" |
---|
1101 | (number-case x |
---|
1102 | (short-float |
---|
1103 | (integer-decode-short-float x)) |
---|
1104 | (double-float |
---|
1105 | (integer-decode-double-float x)))) |
---|
1106 | |
---|
1107 | |
---|
1108 | ;;; %UNARY-TRUNCATE -- Interface |
---|
1109 | ;;; |
---|
1110 | ;;; This function is called when we are doing a truncate without any funky |
---|
1111 | ;;; divisor, i.e. converting a float or ratio to an integer. Note that we do |
---|
1112 | ;;; *not* return the second value of truncate, so it must be computed by the |
---|
1113 | ;;; caller if needed. |
---|
1114 | ;;; |
---|
1115 | ;;; In the float case, we pick off small arguments so that compiler can use |
---|
1116 | ;;; special-case operations. We use an exclusive test, since (due to round-off |
---|
1117 | ;;; error), (float most-positive-fixnum) may be greater than |
---|
1118 | ;;; most-positive-fixnum. |
---|
1119 | ;;; |
---|
1120 | (defun %unary-truncate (number) |
---|
1121 | (number-case number |
---|
1122 | (integer number) |
---|
1123 | (ratio (truncate-no-rem (%numerator number) (%denominator number))) |
---|
1124 | (double-float |
---|
1125 | (if (and (< (the double-float number) |
---|
1126 | (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 0.0d0)) |
---|
1127 | (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 0.0d0) |
---|
1128 | (the double-float number))) |
---|
1129 | (%truncate-double-float->fixnum number) |
---|
1130 | (%truncate-double-float number))) |
---|
1131 | (short-float |
---|
1132 | (if (and (< (the short-float number) |
---|
1133 | (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 0.0s0)) |
---|
1134 | (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 0.0s0) |
---|
1135 | (the short-float number))) |
---|
1136 | (%truncate-short-float->fixnum number) |
---|
1137 | (%truncate-short-float number))))) |
---|
1138 | |
---|
1139 | |
---|
1140 | |
---|
1141 | ; cmucl:compiler:float-tran.lisp |
---|
1142 | (defun xform-truncate (x) |
---|
1143 | (let ((res (%unary-truncate x))) |
---|
1144 | (values res (- x res)))) |
---|
1145 | |
---|
1146 | |
---|
1147 | |
---|
1148 | (defun truncate (number &optional divisor) |
---|
1149 | "Returns number (or number/divisor) as an integer, rounded toward 0. |
---|
1150 | The second returned value is the remainder." |
---|
1151 | (if (null divisor)(setq divisor 1)) |
---|
1152 | (when (not (called-for-mv-p)) |
---|
1153 | (return-from truncate (truncate-no-rem number divisor))) |
---|
1154 | (macrolet |
---|
1155 | ((truncate-rat-dfloat (number divisor) |
---|
1156 | `(with-stack-double-floats ((fnum ,number) |
---|
1157 | (f2)) |
---|
1158 | (let ((res (%unary-truncate (%double-float/-2! fnum ,divisor f2)))) |
---|
1159 | (values res |
---|
1160 | (%double-float--2 fnum (%double-float*-2! (%double-float res f2) ,divisor f2)))))) |
---|
1161 | (truncate-rat-sfloat (number divisor) |
---|
1162 | #+32-bit-target |
---|
1163 | `(target::with-stack-short-floats ((fnum ,number) |
---|
1164 | (f2)) |
---|
1165 | (let ((res (%unary-truncate (%short-float/-2! fnum ,divisor f2)))) |
---|
1166 | (values res |
---|
1167 | (%short-float--2 fnum (%short-float*-2! (%short-float res f2) ,divisor f2))))) |
---|
1168 | #+64-bit-target |
---|
1169 | `(let* ((temp (%short-float ,number)) |
---|
1170 | (res (%unary-truncate (/ (the short-float temp) |
---|
1171 | (the short-float ,divisor))))) |
---|
1172 | (values res |
---|
1173 | (- (the short-float temp) |
---|
1174 | (the short-float (* (the short-float (%short-float res)) |
---|
1175 | (the short-float ,divisor))))))) |
---|
1176 | ) |
---|
1177 | (number-case number |
---|
1178 | (fixnum |
---|
1179 | (number-case divisor |
---|
1180 | (fixnum (if (eq divisor 1) (values number 0) (%fixnum-truncate number divisor))) |
---|
1181 | (bignum (if (eq number target::target-most-negative-fixnum) |
---|
1182 | (with-small-bignum-buffers ((bn number)) |
---|
1183 | (multiple-value-bind (q r) (bignum-truncate bn divisor) |
---|
1184 | (if (eq r bn) |
---|
1185 | (values q number) |
---|
1186 | (values q r)))) |
---|
1187 | (values 0 number))) |
---|
1188 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1189 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1190 | (ratio (let ((q (truncate (* number (%denominator divisor)) ; this was wrong |
---|
1191 | (%numerator divisor)))) |
---|
1192 | (values q (- number (* q divisor))))))) |
---|
1193 | (bignum (number-case divisor |
---|
1194 | (fixnum (if (eq divisor 1) |
---|
1195 | (values number 0) |
---|
1196 | (if (eq divisor target::target-most-negative-fixnum);; << aargh |
---|
1197 | (with-small-bignum-buffers ((bd divisor)) |
---|
1198 | (bignum-truncate number bd)) |
---|
1199 | (bignum-truncate-by-fixnum number divisor)))) |
---|
1200 | (bignum (bignum-truncate number divisor)) |
---|
1201 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1202 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1203 | (ratio (let ((q (truncate (* number (%denominator divisor)) ; so was this |
---|
1204 | (%numerator divisor)))) |
---|
1205 | (values q (- number (* q divisor))))))) |
---|
1206 | (short-float (if (eql divisor 1) |
---|
1207 | (let* ((res (%unary-truncate number))) |
---|
1208 | (values res (- number res))) |
---|
1209 | (number-case divisor |
---|
1210 | (short-float |
---|
1211 | #+32-bit-target |
---|
1212 | (target::with-stack-short-floats ((f2)) |
---|
1213 | (let ((res (%unary-truncate (%short-float/-2! number divisor f2)))) |
---|
1214 | (values res |
---|
1215 | (%short-float--2 |
---|
1216 | number |
---|
1217 | (%short-float*-2! (%short-float res f2) divisor f2))))) |
---|
1218 | #+64-bit-target |
---|
1219 | (let ((res (%unary-truncate |
---|
1220 | (/ (the short-float number) |
---|
1221 | (the short-float divisor))))) |
---|
1222 | (values res |
---|
1223 | (- (the short-float number) |
---|
1224 | (* (the short-float (%short-float res)) |
---|
1225 | (the short-float divisor)))))) |
---|
1226 | ((fixnum bignum ratio) |
---|
1227 | #+32-bit-target |
---|
1228 | (target::with-stack-short-floats ((fdiv divisor) |
---|
1229 | (f2)) |
---|
1230 | (let ((res (%unary-truncate (%short-float/-2! number fdiv f2)))) |
---|
1231 | (values res |
---|
1232 | (%short-float--2 |
---|
1233 | number |
---|
1234 | (%short-float*-2! (%short-float res f2) fdiv f2))))) |
---|
1235 | #+64-bit-target |
---|
1236 | (let* ((fdiv (%short-float divisor)) |
---|
1237 | (res (%unary-truncate |
---|
1238 | (/ (the short-float number) |
---|
1239 | (the short-float fdiv))))) |
---|
1240 | (values res (- number (* res fdiv)))) |
---|
1241 | |
---|
1242 | ) |
---|
1243 | (double-float |
---|
1244 | (with-stack-double-floats ((fnum number) |
---|
1245 | (f2)) |
---|
1246 | (let* ((res (%unary-truncate (%double-float/-2! fnum divisor f2)))) |
---|
1247 | (values res |
---|
1248 | (%double-float--2 |
---|
1249 | fnum |
---|
1250 | (%double-float*-2! (%double-float res f2) divisor f2))))))))) |
---|
1251 | (double-float (if (eql divisor 1) |
---|
1252 | (let ((res (%unary-truncate number))) |
---|
1253 | (values res (- number res))) |
---|
1254 | (number-case divisor |
---|
1255 | ((fixnum bignum ratio short-float) |
---|
1256 | (with-stack-double-floats ((fdiv divisor) |
---|
1257 | (f2)) |
---|
1258 | (let ((res (%unary-truncate (%double-float/-2! number fdiv f2)))) |
---|
1259 | (values res |
---|
1260 | (%double-float--2 |
---|
1261 | number |
---|
1262 | (%double-float*-2! (%double-float res f2) fdiv f2)))))) |
---|
1263 | (double-float |
---|
1264 | (with-stack-double-floats ((f2)) |
---|
1265 | (let ((res (%unary-truncate (%double-float/-2! number divisor f2)))) |
---|
1266 | (values res |
---|
1267 | (%double-float--2 |
---|
1268 | number |
---|
1269 | (%double-float*-2! (%double-float res f2) divisor f2))))))))) |
---|
1270 | (ratio (number-case divisor |
---|
1271 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1272 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1273 | (rational |
---|
1274 | (let ((q (truncate (%numerator number) |
---|
1275 | (* (%denominator number) divisor)))) |
---|
1276 | (values q (- number (* q divisor)))))))))) |
---|
1277 | |
---|
1278 | (defun truncate-no-rem (number divisor) |
---|
1279 | "Returns number (or number/divisor) as an integer, rounded toward 0." |
---|
1280 | (macrolet |
---|
1281 | ((truncate-rat-dfloat (number divisor) |
---|
1282 | `(with-stack-double-floats ((fnum ,number) |
---|
1283 | (f2)) |
---|
1284 | (%unary-truncate (%double-float/-2! fnum ,divisor f2)))) |
---|
1285 | (truncate-rat-sfloat (number divisor) |
---|
1286 | #+32-bit-target |
---|
1287 | `(target::with-stack-short-floats ((fnum ,number) |
---|
1288 | (f2)) |
---|
1289 | (%unary-truncate (%short-float/-2! fnum ,divisor f2))) |
---|
1290 | #+64-bit-target |
---|
1291 | `(let ((fnum (%short-float ,number))) |
---|
1292 | (%unary-truncate (/ (the short-float fnum) |
---|
1293 | (the short-float ,divisor)))))) |
---|
1294 | (number-case number |
---|
1295 | (fixnum |
---|
1296 | (if (eql number target::target-most-negative-fixnum) |
---|
1297 | (if (zerop divisor) |
---|
1298 | (error 'division-by-zero :operation 'truncate :operands (list number divisor)) |
---|
1299 | (with-small-bignum-buffers ((bn number)) |
---|
1300 | (let* ((result (truncate-no-rem bn divisor))) |
---|
1301 | (if (eq result bn) |
---|
1302 | number |
---|
1303 | result)))) |
---|
1304 | (number-case divisor |
---|
1305 | (fixnum (if (eq divisor 1) number (values (%fixnum-truncate number divisor)))) |
---|
1306 | (bignum 0) |
---|
1307 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1308 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1309 | (ratio (let ((q (truncate (* number (%denominator divisor)) |
---|
1310 | (%numerator divisor)))) |
---|
1311 | q))))) |
---|
1312 | (bignum (number-case divisor |
---|
1313 | (fixnum (if (eq divisor 1) number |
---|
1314 | (if (eq divisor target::target-most-negative-fixnum) |
---|
1315 | (with-small-bignum-buffers ((bd divisor)) |
---|
1316 | (bignum-truncate number bd :no-rem)) |
---|
1317 | (bignum-truncate-by-fixnum number divisor)))) |
---|
1318 | (bignum (bignum-truncate number divisor :no-rem)) |
---|
1319 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1320 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1321 | (ratio (let ((q (truncate (* number (%denominator divisor)) |
---|
1322 | (%numerator divisor)))) |
---|
1323 | Q)))) |
---|
1324 | (double-float (if (eql divisor 1) |
---|
1325 | (let ((res (%unary-truncate number))) |
---|
1326 | RES) |
---|
1327 | (number-case divisor |
---|
1328 | ((fixnum bignum ratio) |
---|
1329 | (with-stack-double-floats ((fdiv divisor) |
---|
1330 | (f2)) |
---|
1331 | (let ((res (%unary-truncate (%double-float/-2! number fdiv f2)))) |
---|
1332 | RES))) |
---|
1333 | (short-float |
---|
1334 | (with-stack-double-floats ((ddiv divisor) |
---|
1335 | (f2)) |
---|
1336 | (%unary-truncate (%double-float/-2! number ddiv f2)))) |
---|
1337 | (double-float |
---|
1338 | (with-stack-double-floats ((f2)) |
---|
1339 | (%unary-truncate (%double-float/-2! number divisor f2))))))) |
---|
1340 | (short-float (if (eql divisor 1) |
---|
1341 | (let ((res (%unary-truncate number))) |
---|
1342 | RES) |
---|
1343 | (number-case divisor |
---|
1344 | ((fixnum bignum ratio) |
---|
1345 | #+32-bit-target |
---|
1346 | (target::with-stack-short-floats ((fdiv divisor) |
---|
1347 | (f2)) |
---|
1348 | (let ((res (%unary-truncate (%short-float/-2! number fdiv f2)))) |
---|
1349 | RES)) |
---|
1350 | #+64-bit-target |
---|
1351 | (%unary-truncate (/ (the short-float number) |
---|
1352 | (the short-float (%short-float divisor))))) |
---|
1353 | (short-float |
---|
1354 | #+32-bit-target |
---|
1355 | (target::with-stack-short-floats ((ddiv divisor) |
---|
1356 | (f2)) |
---|
1357 | (%unary-truncate (%short-float/-2! number ddiv f2))) |
---|
1358 | #+64-bit-target |
---|
1359 | (%unary-truncate (/ (the short-float number) |
---|
1360 | (the short-float (%short-float divisor))))) |
---|
1361 | (double-float |
---|
1362 | (with-stack-double-floats ((n2 number) |
---|
1363 | (f2)) |
---|
1364 | (%unary-truncate (%double-float/-2! n2 divisor f2))))))) |
---|
1365 | (ratio (number-case divisor |
---|
1366 | (double-float (truncate-rat-dfloat number divisor)) |
---|
1367 | (short-float (truncate-rat-sfloat number divisor)) |
---|
1368 | (rational |
---|
1369 | (let ((q (truncate (%numerator number) |
---|
1370 | (* (%denominator number) divisor)))) |
---|
1371 | Q))))))) |
---|
1372 | |
---|
1373 | |
---|
1374 | ;;; %UNARY-ROUND -- Interface |
---|
1375 | ;;; |
---|
1376 | ;;; Similar to %UNARY-TRUNCATE, but rounds to the nearest integer. If we |
---|
1377 | ;;; can't use the round primitive, then we do our own round-to-nearest on the |
---|
1378 | ;;; result of i-d-f. [Note that this rounding will really only happen with |
---|
1379 | ;;; double floats, since the whole single-float fraction will fit in a fixnum, |
---|
1380 | ;;; so all single-floats larger than most-positive-fixnum can be precisely |
---|
1381 | ;;; represented by an integer.] |
---|
1382 | ;;; |
---|
1383 | ;;; returns both values today |
---|
1384 | |
---|
1385 | (defun %unary-round (number) |
---|
1386 | (number-case number |
---|
1387 | (integer (values number 0)) |
---|
1388 | (ratio (let ((q (round (%numerator number) (%denominator number)))) |
---|
1389 | (values q (- number q)))) |
---|
1390 | (double-float |
---|
1391 | (if (and (< (the double-float number) |
---|
1392 | (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 1.0d0)) |
---|
1393 | (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 1.0d0) |
---|
1394 | (the double-float number))) |
---|
1395 | (let ((round (%unary-round-to-fixnum number))) |
---|
1396 | (values round (- number round))) |
---|
1397 | (multiple-value-bind (trunc rem) (truncate number) |
---|
1398 | (if (not (%double-float-minusp number)) |
---|
1399 | (if (or (> rem 0.5d0)(and (= rem 0.5d0) (oddp trunc))) |
---|
1400 | (values (+ trunc 1) (- rem 1.0d0)) |
---|
1401 | (values trunc rem)) |
---|
1402 | (if (or (> rem -0.5d0)(and (evenp trunc)(= rem -0.5d0))) |
---|
1403 | (values trunc rem) |
---|
1404 | (values (1- trunc) (+ 1.0d0 rem))))))) |
---|
1405 | (short-float |
---|
1406 | (if (and (< (the short-float number) |
---|
1407 | (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 1.0s0)) |
---|
1408 | (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 1.0s0) |
---|
1409 | (the double-float number))) |
---|
1410 | (let ((round (%unary-round-to-fixnum number))) |
---|
1411 | (values round (- number round))) |
---|
1412 | (multiple-value-bind (trunc rem) (truncate number) |
---|
1413 | (if (not (%short-float-minusp number)) |
---|
1414 | (if (or (> rem 0.5s0)(and (= rem 0.5s0) (oddp trunc))) |
---|
1415 | (values (+ trunc 1) (- rem 1.0s0)) |
---|
1416 | (values trunc rem)) |
---|
1417 | (if (or (> rem -0.5s0)(and (evenp trunc)(= rem -0.5s0))) |
---|
1418 | (values trunc rem) |
---|
1419 | (values (1- trunc) (+ 1.0s0 rem))))))))) |
---|
1420 | |
---|
1421 | (defun %unary-round-to-fixnum (number) |
---|
1422 | (number-case number |
---|
1423 | (double-float |
---|
1424 | (%round-nearest-double-float->fixnum number)) |
---|
1425 | (short-float |
---|
1426 | (%round-nearest-short-float->fixnum number)))) |
---|
1427 | |
---|
1428 | |
---|
1429 | |
---|
1430 | |
---|
1431 | ; cmucl:compiler:float-tran.lisp |
---|
1432 | #| |
---|
1433 | (defun xform-round (x) |
---|
1434 | (let ((res (%unary-round x))) |
---|
1435 | (values res (- x res)))) |
---|
1436 | |# |
---|
1437 | |
---|
1438 | #| |
---|
1439 | (defun round (number &optional divisor) |
---|
1440 | "Rounds number (or number/divisor) to nearest integer. |
---|
1441 | The second returned value is the remainder." |
---|
1442 | (if (null divisor)(setq divisor 1)) |
---|
1443 | (if (eql divisor 1) |
---|
1444 | (xform-round number) |
---|
1445 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
1446 | (let ((thresh (if (integerp divisor) (ash (abs divisor) -1)(/ (abs divisor) 2)))) ; does this need to be a ratio? |
---|
1447 | (cond ((or (> rem thresh) |
---|
1448 | (and (= rem thresh) (oddp tru))) |
---|
1449 | (if (minusp divisor) |
---|
1450 | (values (- tru 1) (+ rem divisor)) |
---|
1451 | (values (+ tru 1) (- rem divisor)))) |
---|
1452 | ((let ((-thresh (- thresh))) |
---|
1453 | (or (< rem -thresh) |
---|
1454 | (and (= rem -thresh) (oddp tru)))) |
---|
1455 | (if (minusp divisor) |
---|
1456 | (values (+ tru 1) (- rem divisor)) |
---|
1457 | (values (- tru 1) (+ rem divisor)))) |
---|
1458 | (t (values tru rem))))))) |
---|
1459 | |# |
---|
1460 | |
---|
1461 | |
---|
1462 | (defun %fixnum-round (number divisor) |
---|
1463 | (declare (fixnum number divisor)) |
---|
1464 | (multiple-value-bind (quo rem)(truncate number divisor) ; should => %fixnum-truncate |
---|
1465 | (if (= 0 rem) |
---|
1466 | (values quo rem) |
---|
1467 | (locally (declare (fixnum quo rem)) |
---|
1468 | (let* ((minusp-num (minusp number)) |
---|
1469 | (minusp-div (minusp divisor)) |
---|
1470 | (2rem (* rem (if (neq minusp-num minusp-div) -2 2)))) |
---|
1471 | ;(declare (fixnum 2rem)) ; no way jose |
---|
1472 | ;(truncate (1- most-positive-fixnum) most-positive-fixnum) |
---|
1473 | ; 2rem has same sign as divisor |
---|
1474 | (cond (minusp-div |
---|
1475 | (if (or (< 2rem divisor) |
---|
1476 | (and (= 2rem divisor)(logbitp 0 quo))) |
---|
1477 | (if minusp-num |
---|
1478 | (values (the fixnum (+ quo 1))(the fixnum (- rem divisor))) |
---|
1479 | (values (the fixnum (- quo 1))(the fixnum (+ rem divisor)))) |
---|
1480 | (values quo rem))) |
---|
1481 | (t (if (or (> 2rem divisor) |
---|
1482 | (and (= 2rem divisor)(logbitp 0 quo))) |
---|
1483 | (if minusp-num |
---|
1484 | (values (the fixnum (- quo 1))(the fixnum (+ rem divisor))) |
---|
1485 | (values (the fixnum (+ quo 1))(the fixnum (- rem divisor)))) |
---|
1486 | (values quo rem))))))))) |
---|
1487 | #| |
---|
1488 | ; + + => + + |
---|
1489 | ; + - => - + |
---|
1490 | ; - + => - - |
---|
1491 | ; - - => + - |
---|
1492 | (defun %fixnum-round (number divisor) |
---|
1493 | (declare (fixnum number divisor)) |
---|
1494 | "Rounds number (or number/divisor) to nearest integer. |
---|
1495 | The second returned value is the remainder." |
---|
1496 | (if (eq divisor 1) |
---|
1497 | (values number 0) |
---|
1498 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
1499 | (if (= 0 rem) |
---|
1500 | (values tru rem) |
---|
1501 | (locally (declare (fixnum tru rem)) |
---|
1502 | (let* ((minusp-num (minusp number)) |
---|
1503 | (minusp-div (minusp divisor)) |
---|
1504 | (half-div (ash (if minusp-div (- divisor) divisor) -1)) |
---|
1505 | (abs-rem (if minusp-num (- rem) rem))) |
---|
1506 | (declare (fixnum half-div abs-rem)) ; true of abs-rem? |
---|
1507 | (if (or (> abs-rem half-div) |
---|
1508 | (and |
---|
1509 | (not (logbitp 0 divisor)) |
---|
1510 | (logbitp 0 tru) ; oddp |
---|
1511 | (= abs-rem half-div))) |
---|
1512 | (if (eq minusp-num minusp-div) |
---|
1513 | (values (the fixnum (+ tru 1))(the fixnum (- rem divisor))) |
---|
1514 | (values (the fixnum (- tru 1))(the fixnum (+ rem divisor)))) |
---|
1515 | (values tru rem)))))))) |
---|
1516 | |# |
---|
1517 | |
---|
1518 | |
---|
1519 | |
---|
1520 | ;; makes 1 piece of garbage instead of average of 2 |
---|
1521 | (defun round (number &optional divisor) |
---|
1522 | "Rounds number (or number/divisor) to nearest integer. |
---|
1523 | The second returned value is the remainder." |
---|
1524 | (if (null divisor)(setq divisor 1)) |
---|
1525 | (if (eql divisor 1) |
---|
1526 | (%unary-round number) |
---|
1527 | (multiple-value-bind (tru rem) (truncate number divisor) |
---|
1528 | (if (= 0 rem) |
---|
1529 | (values tru rem) |
---|
1530 | (let* ((mv-p (called-for-mv-p)) |
---|
1531 | (minusp-num (minusp number)) |
---|
1532 | (minusp-div (minusp divisor)) |
---|
1533 | (2rem (* rem (if (neq minusp-num minusp-div) -2 2)))) |
---|
1534 | ; 2rem has same sign as divisor |
---|
1535 | (cond (minusp-div |
---|
1536 | (if (or (< 2rem divisor) |
---|
1537 | (and (= 2rem divisor)(oddp tru))) |
---|
1538 | (if mv-p |
---|
1539 | (if minusp-num |
---|
1540 | (values (+ tru 1)(- rem divisor)) |
---|
1541 | (values (- tru 1)(+ rem divisor))) |
---|
1542 | (if minusp-num (+ tru 1)(- tru 1))) |
---|
1543 | (values tru rem))) |
---|
1544 | (t (if (or (> 2rem divisor) |
---|
1545 | (and (= 2rem divisor)(oddp tru))) |
---|
1546 | (if mv-p |
---|
1547 | (if minusp-num |
---|
1548 | (values (- tru 1)(+ rem divisor)) |
---|
1549 | (values (+ tru 1)(- rem divisor))) |
---|
1550 | (if minusp-num (- tru 1)(+ tru 1))) |
---|
1551 | (values tru rem))))))))) |
---|
1552 | |
---|
1553 | |
---|
1554 | ;; #-PPC IN L1-NUMBERS.LISP (or implement %%numdiv) |
---|
1555 | ;; Anyone caught implementing %%numdiv will be summarily executed. |
---|
1556 | (defun rem (number divisor) |
---|
1557 | "Returns second result of TRUNCATE." |
---|
1558 | (number-case number |
---|
1559 | (fixnum |
---|
1560 | (number-case divisor |
---|
1561 | (fixnum (nth-value 1 (%fixnum-truncate number divisor))) |
---|
1562 | (bignum |
---|
1563 | (if (and (eql number target::target-most-negative-fixnum) |
---|
1564 | (eql divisor (- target::target-most-negative-fixnum))) |
---|
1565 | 0 |
---|
1566 | number)) |
---|
1567 | (t (nth-value 1 (truncate number divisor))))) |
---|
1568 | (bignum |
---|
1569 | (number-case divisor |
---|
1570 | (fixnum |
---|
1571 | (if (eq divisor target::target-most-negative-fixnum) |
---|
1572 | (nth-value 1 (truncate number divisor)) |
---|
1573 | (bignum-truncate-by-fixnum-no-quo number divisor))) |
---|
1574 | (bignum |
---|
1575 | (bignum-rem number divisor)) |
---|
1576 | (t (nth-value 1 (truncate number divisor))))) |
---|
1577 | (t (nth-value 1 (truncate number divisor))))) |
---|
1578 | |
---|
1579 | ;; #-PPC IN L1-NUMBERS.LISP (or implement %%numdiv) |
---|
1580 | ;; See above. |
---|
1581 | (defun mod (number divisor) |
---|
1582 | "Returns second result of FLOOR." |
---|
1583 | (let ((rem (rem number divisor))) |
---|
1584 | (if (and (not (zerop rem)) |
---|
1585 | (if (minusp divisor) |
---|
1586 | (plusp number) |
---|
1587 | (minusp number))) |
---|
1588 | (+ rem divisor) |
---|
1589 | rem))) |
---|
1590 | |
---|
1591 | (defun cis (theta) |
---|
1592 | "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)." |
---|
1593 | (cond ((complexp theta) |
---|
1594 | (error "Argument to CIS is complex: ~S" theta)) |
---|
1595 | ((or (typep theta 'ratio) |
---|
1596 | (> (abs theta) #.(ash 1 23))) |
---|
1597 | (if (typep theta 'double-float) |
---|
1598 | (%extended-cis theta) |
---|
1599 | (coerce (%extended-cis theta) '(complex single-float)))) |
---|
1600 | (t |
---|
1601 | (complex (cos theta) (sin theta))))) |
---|
1602 | |
---|
1603 | |
---|
1604 | (defun complex (realpart &optional (imagpart 0)) |
---|
1605 | "Return a complex number with the specified real and imaginary components." |
---|
1606 | (number-case realpart |
---|
1607 | (short-float |
---|
1608 | (number-case imagpart |
---|
1609 | (short-float (canonical-complex realpart imagpart)) |
---|
1610 | (double-float (canonical-complex (%double-float realpart) imagpart)) |
---|
1611 | (rational (canonical-complex realpart (%short-float imagpart))))) |
---|
1612 | (double-float |
---|
1613 | (number-case imagpart |
---|
1614 | (double-float (canonical-complex |
---|
1615 | (the double-float realpart) |
---|
1616 | (the double-float imagpart))) |
---|
1617 | (short-float (canonical-complex realpart (%double-float imagpart))) |
---|
1618 | (rational (canonical-complex |
---|
1619 | (the double-float realpart) |
---|
1620 | (the double-float (%double-float imagpart)))))) |
---|
1621 | (rational (number-case imagpart |
---|
1622 | (double-float (canonical-complex |
---|
1623 | (the double-float (%double-float realpart)) |
---|
1624 | (the double-float imagpart))) |
---|
1625 | (short-float (canonical-complex (%short-float realpart) imagpart)) |
---|
1626 | (rational (canonical-complex realpart imagpart)))))) |
---|
1627 | |
---|
1628 | ;; #-PPC IN L1-NUMBERS.LISP |
---|
1629 | (defun realpart (number) |
---|
1630 | "Extract the real part of a number." |
---|
1631 | (number-case number |
---|
1632 | (complex (%realpart number)) |
---|
1633 | (number number))) |
---|
1634 | |
---|
1635 | ;; #-PPC IN L1-NUMBERS.LISP |
---|
1636 | (defun imagpart (number) |
---|
1637 | "Extract the imaginary part of a number." |
---|
1638 | (number-case number |
---|
1639 | (complex (%imagpart number)) |
---|
1640 | (float (* 0 number)) |
---|
1641 | (rational 0))) |
---|
1642 | |
---|
1643 | (defun logand-2 (x y) |
---|
1644 | (number-case x |
---|
1645 | (fixnum (number-case y |
---|
1646 | (fixnum |
---|
1647 | (%ilogand (the fixnum x)(the fixnum y))) |
---|
1648 | (bignum (fix-big-logand x y)))) |
---|
1649 | (bignum (number-case y |
---|
1650 | (fixnum (fix-big-logand y x)) |
---|
1651 | (bignum (bignum-logical-and x y)))))) |
---|
1652 | |
---|
1653 | (defun logior-2 (x y) |
---|
1654 | (number-case x |
---|
1655 | (fixnum (number-case y |
---|
1656 | (fixnum (%ilogior2 x y)) |
---|
1657 | (bignum |
---|
1658 | (if (zerop x) |
---|
1659 | y |
---|
1660 | (with-small-bignum-buffers ((bx x)) |
---|
1661 | (bignum-logical-ior bx y)))))) |
---|
1662 | (bignum (number-case y |
---|
1663 | (fixnum (if (zerop y) |
---|
1664 | x |
---|
1665 | (with-small-bignum-buffers ((by y)) |
---|
1666 | (bignum-logical-ior x by)))) |
---|
1667 | (bignum (bignum-logical-ior x y)))))) |
---|
1668 | |
---|
1669 | (defun logxor-2 (x y) |
---|
1670 | (number-case x |
---|
1671 | (fixnum (number-case y |
---|
1672 | (fixnum (%ilogxor2 x y)) |
---|
1673 | (bignum |
---|
1674 | (with-small-bignum-buffers ((bx x)) |
---|
1675 | (bignum-logical-xor bx y))))) |
---|
1676 | (bignum (number-case y |
---|
1677 | (fixnum (with-small-bignum-buffers ((by y)) |
---|
1678 | (bignum-logical-xor x by))) |
---|
1679 | (bignum (bignum-logical-xor x y)))))) |
---|
1680 | |
---|
1681 | |
---|
1682 | |
---|
1683 | ; see cmucl:compiler:srctran.lisp for transforms |
---|
1684 | |
---|
1685 | (defun lognand (integer1 integer2) |
---|
1686 | "Complement the logical AND of INTEGER1 and INTEGER2." |
---|
1687 | (lognot (logand integer1 integer2))) |
---|
1688 | |
---|
1689 | (defun lognor (integer1 integer2) |
---|
1690 | "Complement the logical AND of INTEGER1 and INTEGER2." |
---|
1691 | (lognot (logior integer1 integer2))) |
---|
1692 | |
---|
1693 | (defun logandc1 (x y) |
---|
1694 | "Return the logical AND of (LOGNOT integer1) and integer2." |
---|
1695 | (number-case x |
---|
1696 | (fixnum (number-case y |
---|
1697 | (fixnum (%ilogand (%ilognot x) y)) |
---|
1698 | (bignum (fix-big-logandc1 x y)))) ; (%ilogand-fix-big (%ilognot x) y)))) |
---|
1699 | (bignum (number-case y |
---|
1700 | (fixnum (fix-big-logandc2 y x)) ; (%ilogandc2-fix-big y x)) |
---|
1701 | (bignum (bignum-logandc2 y x)))))) ;(bignum-logical-and (bignum-logical-not x) y)))))) |
---|
1702 | |
---|
1703 | |
---|
1704 | #| ; its in numbers |
---|
1705 | (defun logandc2 (integer1 integer2) |
---|
1706 | "Returns the logical AND of integer1 and (LOGNOT integer2)." |
---|
1707 | (logand integer1 (lognot integer2))) |
---|
1708 | |# |
---|
1709 | |
---|
1710 | (defun logorc1 (integer1 integer2) |
---|
1711 | "Return the logical OR of (LOGNOT integer1) and integer2." |
---|
1712 | (logior (lognot integer1) integer2)) |
---|
1713 | |
---|
1714 | #| |
---|
1715 | (defun logorc2 (integer1 integer2) |
---|
1716 | "Returns the logical OR of integer1 and (LOGNOT integer2)." |
---|
1717 | (logior integer1 (lognot integer2))) |
---|
1718 | |# |
---|
1719 | |
---|
1720 | (defun logtest (integer1 integer2) |
---|
1721 | "Predicate which returns T if logand of integer1 and integer2 is not zero." |
---|
1722 | ; (not (zerop (logand integer1 integer2))) |
---|
1723 | (number-case integer1 |
---|
1724 | (fixnum (number-case integer2 |
---|
1725 | (fixnum (not (= 0 (%ilogand integer1 integer2)))) |
---|
1726 | (bignum (logtest-fix-big integer1 integer2)))) |
---|
1727 | (bignum (number-case integer2 |
---|
1728 | (fixnum (logtest-fix-big integer2 integer1)) |
---|
1729 | (bignum (bignum-logtest integer1 integer2)))))) |
---|
1730 | |
---|
1731 | |
---|
1732 | |
---|
1733 | (defun lognot (number) |
---|
1734 | "Return the bit-wise logical not of integer." |
---|
1735 | (number-case number |
---|
1736 | (fixnum (%ilognot number)) |
---|
1737 | (bignum (bignum-logical-not number)))) |
---|
1738 | |
---|
1739 | (defun logcount (integer) |
---|
1740 | "Count the number of 1 bits if INTEGER is positive, and the number of 0 bits |
---|
1741 | if INTEGER is negative." |
---|
1742 | (number-case integer |
---|
1743 | (fixnum |
---|
1744 | (%ilogcount (if (minusp (the fixnum integer)) |
---|
1745 | (%ilognot integer) |
---|
1746 | integer))) |
---|
1747 | (bignum |
---|
1748 | (bignum-logcount integer)))) |
---|
1749 | |
---|
1750 | |
---|
1751 | |
---|
1752 | (defun ash (integer count) |
---|
1753 | "Shifts integer left by count places preserving sign. - count shifts right." |
---|
1754 | (etypecase integer |
---|
1755 | (fixnum |
---|
1756 | (etypecase count |
---|
1757 | (fixnum |
---|
1758 | (if (eql integer 0) |
---|
1759 | 0 |
---|
1760 | (if (eql count 0) |
---|
1761 | integer |
---|
1762 | (let ((length (integer-length (the fixnum integer)))) |
---|
1763 | (declare (fixnum length count)) |
---|
1764 | (cond ((and (plusp count) |
---|
1765 | (> (+ length count) |
---|
1766 | (- (1- target::nbits-in-word) target::fixnumshift))) |
---|
1767 | (with-small-bignum-buffers ((bi integer)) |
---|
1768 | (bignum-ashift-left bi count))) |
---|
1769 | ((and (minusp count) (< count (- (1- target::nbits-in-word)))) |
---|
1770 | (if (minusp integer) -1 0)) |
---|
1771 | (t (%iash (the fixnum integer) count))))))) |
---|
1772 | (bignum |
---|
1773 | (if (minusp count) |
---|
1774 | (if (minusp integer) -1 0) |
---|
1775 | (error "Count ~s too large for ASH" count))))) |
---|
1776 | (bignum |
---|
1777 | (etypecase count |
---|
1778 | (fixnum |
---|
1779 | (if (eql count 0) |
---|
1780 | integer |
---|
1781 | (if (plusp count) |
---|
1782 | (bignum-ashift-left integer count) |
---|
1783 | (bignum-ashift-right integer (- count))))) |
---|
1784 | (bignum |
---|
1785 | (if (minusp count) |
---|
1786 | (if (minusp integer) -1 0) |
---|
1787 | (error "Count ~s too large for ASH" count))))))) |
---|
1788 | |
---|
1789 | (defun integer-length (integer) |
---|
1790 | "Return the number of significant bits in the absolute value of integer." |
---|
1791 | (number-case integer |
---|
1792 | (fixnum |
---|
1793 | (%fixnum-intlen (the fixnum integer))) |
---|
1794 | (bignum |
---|
1795 | (bignum-integer-length integer)))) |
---|
1796 | |
---|
1797 | |
---|
1798 | ; not CL, used below |
---|
1799 | (defun byte-mask (size) |
---|
1800 | (1- (ash 1 (the fixnum size)))) |
---|
1801 | |
---|
1802 | (defun byte-position (bytespec) |
---|
1803 | "Return the position part of the byte specifier bytespec." |
---|
1804 | (if (> bytespec 0) |
---|
1805 | (- (integer-length bytespec) (logcount bytespec)) |
---|
1806 | (- bytespec))) |
---|
1807 | |
---|
1808 | |
---|
1809 | ; CMU CL returns T. |
---|
1810 | (defun upgraded-complex-part-type (type) |
---|
1811 | "Return the element type of the most specialized COMPLEX number type that |
---|
1812 | can hold parts of type SPEC." |
---|
1813 | (declare (ignore type)) |
---|
1814 | 'real) |
---|
1815 | |
---|
1816 | ;;; This is the MRG31k3p random number generator described in |
---|
1817 | ;;; P. L'Ecuyer and R. Touzin, "Fast Combined Multiple Recursive |
---|
1818 | ;;; Generators with Multipliers of the form a = +/- 2^d +/- 2^e", |
---|
1819 | ;;; Proceedings of the 2000 Winter Simulation Conference, Dec. 2000, |
---|
1820 | ;;; 683--689. |
---|
1821 | ;;; |
---|
1822 | ;;; A link to the paper is available on L'Ecuyer's web site: |
---|
1823 | ;;; http://www.iro.umontreal.ca/~lecuyer/papers.html. |
---|
1824 | ;;; |
---|
1825 | ;;; This generator has a period of about 2^185. It produces values in |
---|
1826 | ;;; in the half-open interval [0, 2^31 - 1). |
---|
1827 | ;;; |
---|
1828 | ;;; It uses 6 words of state. |
---|
1829 | |
---|
1830 | (defconstant mrg31k3p-m1 #.(- (expt 2 31) 1)) |
---|
1831 | (defconstant mrg31k3p-m2 #.(- (expt 2 31) 21069)) |
---|
1832 | (defconstant mrg31k3p-limit #.(1- (expt 2 31)) |
---|
1833 | "Exclusive upper bound on values returned by %mrg31k3p.") |
---|
1834 | |
---|
1835 | |
---|
1836 | ;;; This is a portable version of the MRG31k3p generator. It's not |
---|
1837 | ;;; too bad in a 64-bit CCL, but the generator pretty much has to be |
---|
1838 | ;;; in LAP for 32-bit ports. |
---|
1839 | #-(or x8632-target ppc32-target x8664-target ppc64-target arm-target) |
---|
1840 | (defun %mrg31k3p (state) |
---|
1841 | (let* ((v (random.mrg31k3p-state state))) |
---|
1842 | (declare (type (simple-array (unsigned-byte 32) (*)) v) |
---|
1843 | (optimize speed)) |
---|
1844 | (let ((y1 (+ (+ (ash (logand (aref v 1) #x1ff) 22) |
---|
1845 | (ash (aref v 1) -9)) |
---|
1846 | (+ (ash (logand (aref v 2) #xffffff) 7) |
---|
1847 | (ash (aref v 2) -24))))) |
---|
1848 | (declare (type (unsigned-byte 32) y1)) |
---|
1849 | (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1)) |
---|
1850 | (incf y1 (aref v 2)) |
---|
1851 | (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1)) |
---|
1852 | (setf (aref v 2) (aref v 1) |
---|
1853 | (aref v 1) (aref v 0) |
---|
1854 | (aref v 0) y1)) |
---|
1855 | (let ((y1 (+ (ash (logand (aref v 3) #xffff) 15) |
---|
1856 | (* 21069 (ash (aref v 3) -16)))) |
---|
1857 | (y2 (+ (ash (logand (aref v 5) #xffff) 15) |
---|
1858 | (* 21069 (ash (aref v 5) -16))))) |
---|
1859 | (declare (type (unsigned-byte 32) y1 y2)) |
---|
1860 | (if (>= y1 mrg31k3p-m2) (decf y1 mrg31k3p-m2)) |
---|
1861 | (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2)) |
---|
1862 | (incf y2 (aref v 5)) |
---|
1863 | (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2)) |
---|
1864 | (incf y2 y1) |
---|
1865 | (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2)) |
---|
1866 | (setf (aref v 5) (aref v 4) |
---|
1867 | (aref v 4) (aref v 3) |
---|
1868 | (aref v 3) y2)) |
---|
1869 | (let* ((x10 (aref v 0)) |
---|
1870 | (x20 (aref v 3))) |
---|
1871 | (if (<= x10 x20) |
---|
1872 | (+ (- x10 x20) mrg31k3p-m1) |
---|
1873 | (- x10 x20))))) |
---|
1874 | |
---|
1875 | (eval-when (:compile-toplevel :execute) |
---|
1876 | (declaim (inline %16-random-bits))) |
---|
1877 | |
---|
1878 | (defun %16-random-bits (state) |
---|
1879 | (logand #xffff (the fixnum (%mrg31k3p state)))) |
---|
1880 | |
---|
1881 | #+64-bit-target |
---|
1882 | (defun %big-fixnum-random (number state) |
---|
1883 | (declare (fixnum number) |
---|
1884 | (ftype (function (random-state) fixnum) %mrg31k3p)) |
---|
1885 | (let ((low (ldb (byte 30 0) (%mrg31k3p state))) |
---|
1886 | (high (ldb (byte 30 0) (%mrg31k3p state)))) |
---|
1887 | (declare (fixnum low high)) |
---|
1888 | (fast-mod (logior low (the fixnum (ash high 30))) |
---|
1889 | number))) |
---|
1890 | |
---|
1891 | ;;; When using a dead simple random number generator, it's reasonable |
---|
1892 | ;;; to take 16 bits of the output and discard the rest. With a more |
---|
1893 | ;;; expensive generator, however, it may be worthwhile to do more bit |
---|
1894 | ;;; fiddling here here so that we can use all of the random bits |
---|
1895 | ;;; produced by %mrg31k2p. |
---|
1896 | #+32-bit-target |
---|
1897 | (defun %bignum-random (number state) |
---|
1898 | (let* ((bits (+ (integer-length number) 8)) |
---|
1899 | (half-words (ash (the fixnum (+ bits 15)) -4)) |
---|
1900 | (long-words (ash (+ half-words 1) -1)) |
---|
1901 | (dividend (%alloc-misc long-words target::subtag-bignum)) |
---|
1902 | (16-bit-dividend dividend) |
---|
1903 | (index 1)) |
---|
1904 | (declare (fixnum long-words index bits) |
---|
1905 | (dynamic-extent dividend) |
---|
1906 | (type (simple-array (unsigned-byte 16) (*)) 16-bit-dividend) ;lie |
---|
1907 | (optimize (speed 3) (safety 0))) |
---|
1908 | (loop |
---|
1909 | ;; This had better inline due to the lie above, or it will error |
---|
1910 | #+big-endian-target |
---|
1911 | (setf (aref 16-bit-dividend index) (%16-random-bits state)) |
---|
1912 | #+little-endian-target |
---|
1913 | (setf (aref 16-bit-dividend (the fixnum (1- index))) |
---|
1914 | (%16-random-bits state)) |
---|
1915 | (decf half-words) |
---|
1916 | (when (<= half-words 0) (return)) |
---|
1917 | #+big-endian-target |
---|
1918 | (setf (aref 16-bit-dividend (the fixnum (1- index))) |
---|
1919 | (%16-random-bits state)) |
---|
1920 | #+little-endian-target |
---|
1921 | (setf (aref 16-bit-dividend index) (%16-random-bits state)) |
---|
1922 | (decf half-words) |
---|
1923 | (when (<= half-words 0) (return)) |
---|
1924 | (incf index 2)) |
---|
1925 | ;; The bignum code expects normalized bignums |
---|
1926 | (let* ((result (mod dividend number))) |
---|
1927 | (if (eq dividend result) |
---|
1928 | (copy-uvector result) |
---|
1929 | result)))) |
---|
1930 | |
---|
1931 | (defun %float-random (number state) |
---|
1932 | (let ((ratio (gvector :ratio (random target::target-most-positive-fixnum state) target::target-most-positive-fixnum))) |
---|
1933 | (declare (dynamic-extent ratio)) |
---|
1934 | (* number ratio))) |
---|
1935 | |
---|
1936 | (defun random (number &optional (state *random-state*)) |
---|
1937 | (if (not (typep state 'random-state)) (report-bad-arg state 'random-state)) |
---|
1938 | (cond |
---|
1939 | ((and (fixnump number) (> (the fixnum number) 0)) |
---|
1940 | #+32-bit-target |
---|
1941 | (fast-mod (%mrg31k3p state) number) |
---|
1942 | #+64-bit-target |
---|
1943 | (if (< number mrg31k3p-limit) |
---|
1944 | (fast-mod (%mrg31k3p state) number) |
---|
1945 | (%big-fixnum-random number state))) |
---|
1946 | ((and (typep number 'double-float) (> (the double-float number) 0.0)) |
---|
1947 | (%float-random number state)) |
---|
1948 | ((and (typep number 'short-float) (> (the short-float number) 0.0s0)) |
---|
1949 | (%float-random number state)) |
---|
1950 | ((and (bignump number) (> number 0)) |
---|
1951 | (%bignum-random number state)) |
---|
1952 | (t (report-bad-arg number '(or (integer (0)) (float (0.0))))))) |
---|
1953 | |
---|
1954 | (eval-when (:compile-toplevel :execute) |
---|
1955 | (defmacro bignum-abs (nexp) |
---|
1956 | (let ((n (gensym))) |
---|
1957 | `(let ((,n ,nexp)) |
---|
1958 | (if (bignum-minusp ,n) (negate-bignum ,n) ,n)))) |
---|
1959 | |
---|
1960 | (defmacro fixnum-abs (nexp) |
---|
1961 | (let ((n (gensym))) |
---|
1962 | `(let ((,n ,nexp)) |
---|
1963 | (if (minusp (the fixnum ,n)) |
---|
1964 | (if (eq ,n target::target-most-negative-fixnum) |
---|
1965 | (- ,n) |
---|
1966 | (the fixnum (- (the fixnum ,n)))) |
---|
1967 | ,n)))) |
---|
1968 | ) |
---|
1969 | |
---|
1970 | |
---|
1971 | ;;; TWO-ARG-GCD -- Internal |
---|
1972 | ;;; |
---|
1973 | ;;; Do the GCD of two integer arguments. With fixnum arguments, we use the |
---|
1974 | ;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly |
---|
1975 | ;;; structurified), otherwise we call BIGNUM-GCD. We pick off the special case |
---|
1976 | ;;; of 0 before the dispatch so that the bignum code doesn't have to worry |
---|
1977 | ;;; about "small bignum" zeros. |
---|
1978 | ;;; |
---|
1979 | (defun gcd-2 (n1 n2) |
---|
1980 | ;(declare (optimize (speed 3)(safety 0))) |
---|
1981 | (cond |
---|
1982 | ((eql n1 0) (%integer-abs n2)) |
---|
1983 | ((eql n2 0) (%integer-abs n1)) |
---|
1984 | (t (number-case n1 |
---|
1985 | (fixnum |
---|
1986 | (number-case n2 |
---|
1987 | (fixnum |
---|
1988 | (if (eql n1 target::target-most-negative-fixnum) |
---|
1989 | (if (eql n2 target::target-most-negative-fixnum) |
---|
1990 | (- target::target-most-negative-fixnum) |
---|
1991 | (bignum-fixnum-gcd (- target::target-most-negative-fixnum) (abs n2))) |
---|
1992 | (if (eql n2 target::target-most-negative-fixnum) |
---|
1993 | (bignum-fixnum-gcd (- target::target-most-negative-fixnum) (abs n1)) |
---|
1994 | (locally |
---|
1995 | (declare (optimize (speed 3) (safety 0)) |
---|
1996 | (fixnum n1 n2)) |
---|
1997 | (if (minusp n1)(setq n1 (the fixnum (- n1)))) |
---|
1998 | (if (minusp n2)(setq n2 (the fixnum (- n2)))) |
---|
1999 | (%fixnum-gcd n1 n2))))) |
---|
2000 | (bignum (if (eql n1 target::target-most-negative-fixnum) |
---|
2001 | (%bignum-bignum-gcd n2 (- target::target-most-negative-fixnum)) |
---|
2002 | (bignum-fixnum-gcd (bignum-abs n2)(fixnum-abs n1)))))) |
---|
2003 | (bignum |
---|
2004 | (number-case n2 |
---|
2005 | (fixnum |
---|
2006 | (if (eql n2 target::target-most-negative-fixnum) |
---|
2007 | (%bignum-bignum-gcd (bignum-abs n1)(fixnum-abs n2)) |
---|
2008 | (bignum-fixnum-gcd (bignum-abs n1)(fixnum-abs n2)))) |
---|
2009 | (bignum (%bignum-bignum-gcd n1 n2)))))))) |
---|
2010 | |
---|
2011 | #| |
---|
2012 | (defun fixnum-gcd (n1 n2) |
---|
2013 | (declare (optimize (speed 3) (safety 0)) |
---|
2014 | (fixnum n1 n2)) |
---|
2015 | (do* ((k 0 (%i+ 1 k)) |
---|
2016 | (n1 n1 (%iasr 1 n1)) |
---|
2017 | (n2 n2 (%iasr 1 n2))) |
---|
2018 | ((oddp (logior n1 n2)) |
---|
2019 | (do ((temp (if (oddp n1) (the fixnum (- n2)) (%iasr 1 n1)) |
---|
2020 | (%iasr 1 temp))) |
---|
2021 | (nil) |
---|
2022 | (declare (fixnum temp)) |
---|
2023 | (when (oddp temp) |
---|
2024 | (if (plusp temp) |
---|
2025 | (setq n1 temp) |
---|
2026 | (setq n2 (- temp))) |
---|
2027 | (setq temp (the fixnum (- n1 n2))) |
---|
2028 | (when (zerop temp) |
---|
2029 | (let ((res (%ilsl k n1))) |
---|
2030 | (return res)))))) |
---|
2031 | (declare (fixnum n1 n2 k)))) |
---|
2032 | |# |
---|
2033 | |
---|
2034 | |
---|
2035 | |
---|
2036 | (defun %quo-1 (n) |
---|
2037 | (/ 1 n)) |
---|
2038 | |
---|
2039 | ;; Compute (sqrt (+ (* x x) (* y y))), but |
---|
2040 | ;; try to be a little more careful about it. |
---|
2041 | ;; Both x and y must be double-floats. |
---|
2042 | (defun %double-float-hypot (x y) |
---|
2043 | (with-stack-double-floats ((a) (b) (c)) |
---|
2044 | (%%double-float-abs! x a) |
---|
2045 | (%%double-float-abs! y b) |
---|
2046 | (when (> a b) |
---|
2047 | (psetq a b b a)) |
---|
2048 | (if (= b 0d0) |
---|
2049 | 0d0 |
---|
2050 | (progn |
---|
2051 | (%double-float/-2! a b c) |
---|
2052 | (* b (fsqrt (+ 1d0 (* c c)))))))) |
---|
2053 | |
---|