source: trunk/source/level-0/l0-numbers.lisp @ 15601

Last change on this file since 15601 was 15560, checked in by gb, 6 years ago

In *-2: careful if multiplying MOST-NEGATIVE-FIXNUM by a BIGNUM.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 78.2 KB
Line 
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                     (bignum-truncate bn divisor))
1184                   (values 0 number)))
1185         (double-float (truncate-rat-dfloat number divisor))
1186         (short-float (truncate-rat-sfloat number divisor))
1187         (ratio (let ((q (truncate (* number (%denominator divisor)) ; this was wrong
1188                                   (%numerator divisor))))
1189                  (values q (- number (* q divisor)))))))
1190      (bignum (number-case divisor
1191                (fixnum (if (eq divisor 1)
1192                          (values number 0)
1193                          (if (eq divisor target::target-most-negative-fixnum);; << aargh
1194                            (with-small-bignum-buffers ((bd divisor))
1195                              (bignum-truncate number bd))
1196                            (bignum-truncate-by-fixnum number divisor))))
1197                (bignum (bignum-truncate number divisor))
1198                (double-float  (truncate-rat-dfloat number divisor))
1199                (short-float (truncate-rat-sfloat number divisor))
1200                (ratio (let ((q (truncate (* number (%denominator divisor)) ; so was this
1201                                          (%numerator divisor))))
1202                         (values q (- number (* q divisor)))))))
1203      (short-float (if (eql divisor 1)
1204                     (let* ((res (%unary-truncate number)))
1205                       (values res (- number res)))
1206                     (number-case divisor
1207                       (short-float
1208                        #+32-bit-target
1209                        (target::with-stack-short-floats ((f2))
1210                          (let ((res (%unary-truncate (%short-float/-2! number divisor f2))))
1211                            (values res 
1212                                    (%short-float--2
1213                                     number 
1214                                     (%short-float*-2! (%short-float res f2) divisor f2)))))
1215                        #+64-bit-target
1216                        (let ((res (%unary-truncate
1217                                    (/ (the short-float number)
1218                                       (the short-float divisor)))))
1219                          (values res
1220                                  (- (the short-float number)
1221                                     (* (the short-float (%short-float res))
1222                                        (the short-float divisor))))))
1223                       ((fixnum bignum ratio)
1224                        #+32-bit-target
1225                        (target::with-stack-short-floats ((fdiv divisor)
1226                                                          (f2))
1227                          (let ((res (%unary-truncate (%short-float/-2! number fdiv f2))))
1228                            (values res 
1229                                    (%short-float--2 
1230                                     number 
1231                                     (%short-float*-2! (%short-float res f2) fdiv f2)))))
1232                        #+64-bit-target
1233                        (let* ((fdiv (%short-float divisor))
1234                               (res (%unary-truncate
1235                                     (/ (the short-float number)
1236                                        (the short-float fdiv)))))
1237                          (values res (- number (* res fdiv))))
1238                                     
1239                        )
1240                       (double-float
1241                        (with-stack-double-floats ((fnum number)
1242                                                   (f2))
1243                          (let* ((res (%unary-truncate (%double-float/-2! fnum divisor f2))))
1244                            (values res
1245                                    (%double-float--2
1246                                     fnum
1247                                     (%double-float*-2! (%double-float res f2) divisor f2)))))))))
1248      (double-float (if (eql divisor 1)
1249                      (let ((res (%unary-truncate number)))
1250                        (values res (- number res)))
1251                      (number-case divisor
1252                        ((fixnum bignum ratio short-float)
1253                         (with-stack-double-floats ((fdiv divisor)
1254                                                    (f2))
1255                           (let ((res (%unary-truncate (%double-float/-2! number fdiv f2))))
1256                             (values res 
1257                                     (%double-float--2 
1258                                      number 
1259                                      (%double-float*-2! (%double-float res f2) fdiv f2))))))                       
1260                        (double-float
1261                         (with-stack-double-floats ((f2))
1262                           (let ((res (%unary-truncate (%double-float/-2! number divisor f2))))
1263                             (values res 
1264                                     (%double-float--2
1265                                      number 
1266                                      (%double-float*-2! (%double-float res f2) divisor f2)))))))))
1267      (ratio (number-case divisor
1268               (double-float (truncate-rat-dfloat number divisor))
1269               (short-float (truncate-rat-sfloat number divisor))
1270               (rational
1271                (let ((q (truncate (%numerator number)
1272                                   (* (%denominator number) divisor))))
1273                  (values q (- number (* q divisor))))))))))
1274
1275(defun truncate-no-rem (number  divisor)
1276  "Returns number (or number/divisor) as an integer, rounded toward 0."
1277  (macrolet 
1278    ((truncate-rat-dfloat (number divisor)
1279       `(with-stack-double-floats ((fnum ,number)
1280                                      (f2))
1281         (%unary-truncate (%double-float/-2! fnum ,divisor f2))))
1282     (truncate-rat-sfloat (number divisor)
1283       #+32-bit-target
1284       `(target::with-stack-short-floats ((fnum ,number)
1285                                      (f2))
1286         (%unary-truncate (%short-float/-2! fnum ,divisor f2)))
1287       #+64-bit-target
1288       `(let ((fnum (%short-float ,number)))
1289         (%unary-truncate (/ (the short-float fnum)
1290                           (the short-float ,divisor))))))
1291    (number-case number
1292    (fixnum
1293     (if (eql number target::target-most-negative-fixnum)
1294       (if (zerop divisor)
1295         (error 'division-by-zero :operation 'truncate :operands (list number divisor))
1296         (with-small-bignum-buffers ((bn number))
1297           (let* ((result (truncate-no-rem bn divisor)))
1298             (if (eq result bn)
1299               number
1300               result))))
1301       (number-case divisor
1302         (fixnum (if (eq divisor 1) number (values (%fixnum-truncate number divisor))))
1303         (bignum 0)
1304         (double-float (truncate-rat-dfloat number divisor))
1305         (short-float (truncate-rat-sfloat number divisor))
1306         (ratio (let ((q (truncate (* number (%denominator divisor))
1307                                   (%numerator divisor))))
1308                  q)))))
1309     (bignum (number-case divisor
1310               (fixnum (if (eq divisor 1) number
1311                         (if (eq divisor target::target-most-negative-fixnum)
1312                           (with-small-bignum-buffers ((bd divisor))
1313                             (bignum-truncate number bd :no-rem))
1314                           (bignum-truncate-by-fixnum number divisor))))
1315               (bignum (bignum-truncate number divisor :no-rem))
1316               (double-float  (truncate-rat-dfloat number divisor))
1317               (short-float (truncate-rat-sfloat number divisor))
1318               (ratio (let ((q (truncate (* number (%denominator divisor))
1319                                         (%numerator divisor))))
1320                        Q))))
1321     (double-float (if (eql divisor 1)
1322                     (let ((res (%unary-truncate number)))
1323                       RES)
1324                     (number-case divisor
1325                       ((fixnum bignum ratio)
1326                        (with-stack-double-floats ((fdiv divisor)
1327                                                   (f2))
1328                          (let ((res (%unary-truncate (%double-float/-2! number fdiv f2))))
1329                            RES)))
1330                       (short-float
1331                        (with-stack-double-floats ((ddiv divisor)
1332                                                   (f2))
1333                          (%unary-truncate (%double-float/-2! number ddiv f2))))
1334                       (double-float
1335                        (with-stack-double-floats ((f2))
1336                          (%unary-truncate (%double-float/-2! number divisor f2)))))))
1337     (short-float (if (eql divisor 1)
1338                    (let ((res (%unary-truncate number)))
1339                      RES)
1340                    (number-case divisor
1341                      ((fixnum bignum ratio)
1342                       #+32-bit-target
1343                       (target::with-stack-short-floats ((fdiv divisor)
1344                                                 (f2))
1345                         (let ((res (%unary-truncate (%short-float/-2! number fdiv f2))))
1346                           RES))
1347                       #+64-bit-target
1348                       (%unary-truncate (/ (the short-float number)
1349                                           (the short-float (%short-float divisor)))))
1350                      (short-float
1351                       #+32-bit-target
1352                       (target::with-stack-short-floats ((ddiv divisor)
1353                                                      (f2))
1354                         (%unary-truncate (%short-float/-2! number ddiv f2)))
1355                       #+64-bit-target
1356                       (%unary-truncate (/ (the short-float number)
1357                                           (the short-float (%short-float divisor)))))
1358                      (double-float
1359                       (with-stack-double-floats ((n2 number)
1360                                                      (f2))
1361                         (%unary-truncate (%double-float/-2! n2 divisor f2)))))))
1362    (ratio (number-case divisor
1363                  (double-float (truncate-rat-dfloat number divisor))
1364                  (short-float (truncate-rat-sfloat number divisor))
1365                  (rational
1366                   (let ((q (truncate (%numerator number)
1367                                      (* (%denominator number) divisor))))
1368                     Q)))))))
1369
1370
1371;;; %UNARY-ROUND  --  Interface
1372;;;
1373;;;    Similar to %UNARY-TRUNCATE, but rounds to the nearest integer.  If we
1374;;; can't use the round primitive, then we do our own round-to-nearest on the
1375;;; result of i-d-f.  [Note that this rounding will really only happen with
1376;;; double floats, since the whole single-float fraction will fit in a fixnum,
1377;;; so all single-floats larger than most-positive-fixnum can be precisely
1378;;; represented by an integer.]
1379;;;
1380;;; returns both values today
1381
1382(defun %unary-round (number)
1383  (number-case number
1384    (integer (values number 0))
1385    (ratio (let ((q (round (%numerator number) (%denominator number))))             
1386             (values q (- number q))))
1387    (double-float
1388     (if (and (< (the double-float number) 
1389                 (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 1.0d0))
1390              (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 1.0d0)
1391                 (the double-float number)))
1392       (let ((round (%unary-round-to-fixnum number)))
1393         (values round (- number round)))
1394       (multiple-value-bind (trunc rem) (truncate number)         
1395         (if (not (%double-float-minusp number))
1396           (if (or (> rem 0.5d0)(and (= rem 0.5d0) (oddp trunc)))
1397             (values (+ trunc 1) (- rem 1.0d0))
1398             (values trunc rem))
1399           (if (or (> rem -0.5d0)(and (evenp trunc)(= rem -0.5d0)))
1400             (values trunc rem)
1401             (values (1- trunc) (+ 1.0d0 rem)))))))
1402    (short-float
1403     (if (and (< (the short-float number) 
1404                 (float (1- (ash 1 (- (1- target::nbits-in-word) target::fixnumshift))) 1.0s0))
1405              (< (float (ash -1 (- (1- target::nbits-in-word) target::fixnumshift)) 1.0s0)
1406                 (the double-float number)))
1407       (let ((round (%unary-round-to-fixnum number)))
1408         (values round (- number round)))
1409       (multiple-value-bind (trunc rem) (truncate number)         
1410         (if (not (%short-float-minusp number))
1411           (if (or (> rem 0.5s0)(and (= rem 0.5s0) (oddp trunc)))
1412             (values (+ trunc 1) (- rem 1.0s0))
1413             (values trunc rem))
1414           (if (or (> rem -0.5s0)(and (evenp trunc)(= rem -0.5s0)))
1415             (values trunc rem)
1416             (values (1- trunc) (+ 1.0s0 rem)))))))))
1417
1418(defun %unary-round-to-fixnum (number)
1419  (number-case number
1420    (double-float
1421     (%round-nearest-double-float->fixnum number))
1422    (short-float
1423     (%round-nearest-short-float->fixnum number))))
1424
1425                         
1426                               
1427         
1428; cmucl:compiler:float-tran.lisp
1429#|
1430(defun xform-round (x)
1431  (let ((res (%unary-round x)))
1432    (values res (- x res))))
1433|#
1434
1435#|
1436(defun round (number &optional divisor)
1437  "Rounds number (or number/divisor) to nearest integer.
1438  The second returned value is the remainder."
1439  (if (null divisor)(setq divisor 1))
1440  (if (eql divisor 1)
1441    (xform-round number)
1442    (multiple-value-bind (tru rem) (truncate number divisor)
1443      (let ((thresh (if (integerp divisor) (ash (abs divisor) -1)(/ (abs divisor) 2)))) ; does this need to be a ratio?
1444        (cond ((or (> rem thresh)
1445                   (and (= rem thresh) (oddp tru)))
1446               (if (minusp divisor)
1447                 (values (- tru 1) (+ rem divisor))
1448                 (values (+ tru 1) (- rem divisor))))
1449              ((let ((-thresh (- thresh)))
1450                 (or (< rem -thresh)
1451                     (and (= rem -thresh) (oddp tru))))
1452               (if (minusp divisor)
1453                 (values (+ tru 1) (- rem divisor))
1454                 (values (- tru 1) (+ rem divisor))))
1455              (t (values tru rem)))))))
1456|#
1457
1458
1459(defun %fixnum-round (number divisor)
1460  (declare (fixnum number divisor))
1461  (multiple-value-bind (quo rem)(truncate number divisor) ; should => %fixnum-truncate
1462    (if (= 0 rem)
1463      (values quo rem)
1464      (locally (declare (fixnum quo rem))
1465        (let* ((minusp-num (minusp number))
1466               (minusp-div (minusp divisor))
1467               (2rem (* rem (if (neq minusp-num minusp-div) -2 2))))
1468          ;(declare (fixnum 2rem)) ; no way jose 
1469          ;(truncate (1- most-positive-fixnum) most-positive-fixnum)
1470          ; 2rem has same sign as divisor
1471          (cond (minusp-div             
1472                 (if (or (< 2rem divisor)
1473                         (and (= 2rem divisor)(logbitp 0 quo)))
1474                   (if minusp-num
1475                     (values (the fixnum (+ quo 1))(the fixnum (- rem divisor)))
1476                     (values (the fixnum (- quo 1))(the fixnum (+ rem divisor))))
1477                   (values quo rem)))
1478                (t (if (or (> 2rem divisor)
1479                           (and (= 2rem divisor)(logbitp 0 quo)))
1480                     (if minusp-num
1481                       (values (the fixnum (- quo 1))(the fixnum (+ rem divisor)))
1482                       (values (the fixnum (+ quo 1))(the fixnum (- rem divisor))))
1483                     (values quo rem)))))))))
1484#|
1485; + + => + +
1486; + - => - +
1487; - + => - -
1488; - - => + -
1489(defun %fixnum-round (number divisor)
1490  (declare (fixnum number divisor))
1491  "Rounds number (or number/divisor) to nearest integer.
1492  The second returned value is the remainder."
1493  (if (eq divisor 1)
1494    (values number 0)
1495    (multiple-value-bind (tru rem) (truncate number divisor)
1496      (if (= 0 rem)
1497        (values tru rem)
1498        (locally (declare (fixnum tru rem))
1499          (let* ((minusp-num (minusp number))
1500                 (minusp-div (minusp divisor))
1501                 (half-div (ash (if minusp-div (- divisor) divisor) -1))
1502                 (abs-rem (if minusp-num (- rem) rem)))           
1503            (declare (fixnum half-div abs-rem)) ; true of abs-rem?
1504            (if (or (> abs-rem half-div)
1505                    (and
1506                     (not (logbitp 0 divisor))
1507                     (logbitp 0 tru) ; oddp
1508                     (= abs-rem half-div)))
1509              (if (eq minusp-num minusp-div)
1510                (values (the fixnum (+ tru 1))(the fixnum (- rem divisor)))
1511                (values (the fixnum (- tru 1))(the fixnum (+ rem divisor))))
1512              (values tru rem))))))))
1513|#
1514
1515
1516
1517;; makes 1 piece of garbage instead of average of 2
1518(defun round (number &optional divisor)
1519  "Rounds number (or number/divisor) to nearest integer.
1520  The second returned value is the remainder."
1521  (if (null divisor)(setq divisor 1))
1522  (if (eql divisor 1)
1523    (%unary-round number)
1524    (multiple-value-bind (tru rem) (truncate number divisor)
1525      (if (= 0 rem)
1526        (values tru rem)
1527        (let* ((mv-p (called-for-mv-p))
1528               (minusp-num (minusp number))
1529               (minusp-div (minusp divisor))
1530               (2rem (* rem (if (neq minusp-num minusp-div) -2 2))))
1531          ; 2rem has same sign as divisor
1532          (cond (minusp-div             
1533                 (if (or (< 2rem divisor)
1534                         (and (= 2rem divisor)(oddp tru)))
1535                   (if mv-p
1536                     (if minusp-num
1537                       (values (+ tru 1)(- rem divisor))
1538                       (values (- tru 1)(+ rem divisor)))
1539                     (if minusp-num (+ tru 1)(- tru 1)))
1540                   (values tru rem)))
1541                (t (if (or (> 2rem divisor)
1542                           (and (= 2rem divisor)(oddp tru)))
1543                     (if mv-p
1544                       (if minusp-num
1545                         (values (- tru 1)(+ rem divisor))
1546                         (values (+ tru 1)(- rem divisor)))
1547                       (if minusp-num (- tru 1)(+ tru 1)))
1548                     (values tru rem)))))))))
1549
1550
1551;; #-PPC IN L1-NUMBERS.LISP (or implement %%numdiv)
1552;; Anyone caught implementing %%numdiv will be summarily executed.
1553(defun rem (number divisor)
1554  "Returns second result of TRUNCATE."
1555  (number-case number
1556    (fixnum
1557     (number-case divisor
1558       (fixnum (nth-value 1 (%fixnum-truncate number divisor)))
1559       (bignum number)
1560       (t (nth-value 1 (truncate number divisor)))))
1561    (bignum
1562     (number-case divisor
1563       (fixnum
1564        (if (eq divisor target::target-most-negative-fixnum)
1565          (nth-value 1 (truncate number divisor))
1566          (bignum-truncate-by-fixnum-no-quo number divisor)))
1567       (bignum
1568        (bignum-rem number divisor))
1569       (t (nth-value 1 (truncate number divisor)))))
1570    (t (nth-value 1 (truncate number divisor)))))
1571
1572;; #-PPC IN L1-NUMBERS.LISP (or implement %%numdiv)
1573;; See above.
1574(defun mod (number divisor)
1575  "Returns second result of FLOOR."
1576  (let ((rem (rem number divisor)))
1577    (if (and (not (zerop rem))
1578             (if (minusp divisor)
1579                 (plusp number)
1580                 (minusp number)))
1581        (+ rem divisor)
1582        rem)))
1583
1584(defun cis (theta)
1585  "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
1586  (if (complexp theta)
1587    (error "Argument to CIS is complex: ~S" theta)
1588    (complex (cos theta) (sin theta))))
1589
1590
1591(defun complex (realpart &optional (imagpart 0))
1592  "Return a complex number with the specified real and imaginary components."
1593  (number-case realpart
1594    (short-float
1595      (number-case imagpart
1596         (short-float (canonical-complex realpart imagpart))
1597         (double-float (canonical-complex (%double-float realpart) imagpart))
1598         (rational (canonical-complex realpart (%short-float imagpart)))))
1599    (double-float 
1600     (number-case imagpart
1601       (double-float (canonical-complex
1602                      (the double-float realpart)
1603                      (the double-float imagpart)))
1604       (short-float (canonical-complex realpart (%double-float imagpart)))
1605       (rational (canonical-complex
1606                              (the double-float realpart)
1607                              (the double-float (%double-float imagpart))))))
1608    (rational (number-case imagpart
1609                (double-float (canonical-complex
1610                               (the double-float (%double-float realpart))
1611                               (the double-float imagpart)))
1612                (short-float (canonical-complex (%short-float realpart) imagpart))
1613                (rational (canonical-complex realpart imagpart)))))) 
1614
1615;; #-PPC IN L1-NUMBERS.LISP
1616(defun realpart (number)
1617  "Extract the real part of a number."
1618  (number-case number
1619    (complex (%realpart number))
1620    (number number)))
1621
1622;; #-PPC IN L1-NUMBERS.LISP
1623(defun imagpart (number)
1624  "Extract the imaginary part of a number."
1625  (number-case number
1626    (complex (%imagpart number))
1627    (float (* 0 number))
1628    (rational 0)))
1629
1630(defun logand-2 (x y) 
1631  (number-case x
1632    (fixnum (number-case y
1633              (fixnum
1634               (%ilogand (the fixnum x)(the fixnum y)))
1635              (bignum (fix-big-logand x y))))
1636    (bignum (number-case y
1637              (fixnum (fix-big-logand y x))
1638              (bignum (bignum-logical-and x y))))))
1639
1640(defun logior-2 (x y)
1641  (number-case x
1642    (fixnum (number-case y
1643              (fixnum (%ilogior2 x y))
1644              (bignum
1645               (if (zerop x)
1646                 y
1647                 (with-small-bignum-buffers ((bx x))
1648                   (bignum-logical-ior bx y))))))
1649    (bignum (number-case y
1650              (fixnum (if (zerop y)
1651                        x
1652                        (with-small-bignum-buffers ((by y))
1653                          (bignum-logical-ior x by))))
1654              (bignum (bignum-logical-ior x y))))))
1655
1656(defun logxor-2 (x y)
1657  (number-case x
1658    (fixnum (number-case y
1659              (fixnum (%ilogxor2 x y))
1660              (bignum
1661               (with-small-bignum-buffers ((bx x))
1662                 (bignum-logical-xor bx y)))))
1663    (bignum (number-case y
1664              (fixnum (with-small-bignum-buffers ((by y))
1665                        (bignum-logical-xor x by)))
1666              (bignum (bignum-logical-xor x y))))))
1667
1668               
1669
1670; see cmucl:compiler:srctran.lisp for transforms
1671
1672(defun lognand (integer1 integer2)
1673  "Complement the logical AND of INTEGER1 and INTEGER2."
1674  (lognot (logand integer1 integer2)))
1675
1676(defun lognor (integer1 integer2)
1677  "Complement the logical AND of INTEGER1 and INTEGER2."
1678  (lognot (logior integer1 integer2)))
1679
1680(defun logandc1 (x y)
1681  "Return the logical AND of (LOGNOT integer1) and integer2."
1682  (number-case x
1683    (fixnum (number-case y               
1684              (fixnum (%ilogand (%ilognot x) y))
1685              (bignum  (fix-big-logandc1 x y))))    ; (%ilogand-fix-big (%ilognot x) y))))
1686    (bignum (number-case y
1687              (fixnum  (fix-big-logandc2 y x))      ; (%ilogandc2-fix-big y x))
1688              (bignum (bignum-logandc2 y x))))))    ;(bignum-logical-and (bignum-logical-not x)  y))))))
1689
1690
1691#| ; its in numbers
1692(defun logandc2 (integer1 integer2)
1693  "Returns the logical AND of integer1 and (LOGNOT integer2)."
1694  (logand integer1 (lognot integer2)))
1695|#
1696
1697(defun logorc1 (integer1 integer2)
1698  "Return the logical OR of (LOGNOT integer1) and integer2."
1699  (logior (lognot integer1) integer2))
1700
1701#|
1702(defun logorc2 (integer1 integer2)
1703  "Returns the logical OR of integer1 and (LOGNOT integer2)."
1704  (logior integer1 (lognot integer2)))
1705|#
1706
1707(defun logtest (integer1 integer2)
1708  "Predicate which returns T if logand of integer1 and integer2 is not zero."
1709 ; (not (zerop (logand integer1 integer2)))
1710  (number-case integer1
1711    (fixnum (number-case integer2
1712              (fixnum (not (= 0 (%ilogand integer1 integer2))))
1713              (bignum (logtest-fix-big integer1 integer2))))
1714    (bignum (number-case integer2
1715              (fixnum (logtest-fix-big integer2 integer1))
1716              (bignum (bignum-logtest integer1 integer2)))))) 
1717
1718
1719
1720(defun lognot (number)
1721  "Return the bit-wise logical not of integer."
1722  (number-case number
1723    (fixnum (%ilognot number))
1724    (bignum (bignum-logical-not number))))
1725
1726(defun logcount (integer)
1727  "Count the number of 1 bits if INTEGER is positive, and the number of 0 bits
1728  if INTEGER is negative."
1729  (number-case integer
1730    (fixnum
1731     (%ilogcount (if (minusp (the fixnum integer))
1732                   (%ilognot integer)
1733                   integer)))
1734    (bignum
1735     (bignum-logcount integer))))
1736
1737
1738
1739(defun ash (integer count)
1740  "Shifts integer left by count places preserving sign. - count shifts right."
1741  (etypecase integer
1742    (fixnum
1743     (etypecase count
1744       (fixnum
1745        (if (eql integer 0)
1746          0
1747          (if (eql count 0)
1748            integer
1749            (let ((length (integer-length (the fixnum integer))))
1750              (declare (fixnum length count))
1751              (cond ((and (plusp count)
1752                          (> (+ length count)
1753                             (- (1- target::nbits-in-word) target::fixnumshift)))
1754                     (with-small-bignum-buffers ((bi integer))
1755                       (bignum-ashift-left bi count)))
1756                    ((and (minusp count) (< count (- (1- target::nbits-in-word))))
1757                     (if (minusp integer) -1 0))
1758                    (t (%iash (the fixnum integer) count)))))))
1759       (bignum
1760        (if (minusp count)
1761          (if (minusp integer) -1 0)         
1762          (error "Count ~s too large for ASH" count)))))
1763    (bignum
1764     (etypecase count
1765       (fixnum
1766        (if (eql count 0) 
1767          integer
1768          (if (plusp count)
1769            (bignum-ashift-left integer count)
1770            (bignum-ashift-right integer (- count)))))
1771       (bignum
1772        (if (minusp count)
1773          (if (minusp integer) -1 0)
1774          (error "Count ~s too large for ASH" count)))))))
1775
1776(defun integer-length (integer)
1777  "Return the number of significant bits in the absolute value of integer."
1778  (number-case integer
1779    (fixnum
1780     (%fixnum-intlen (the fixnum integer)))
1781    (bignum
1782     (bignum-integer-length integer))))
1783
1784
1785; not CL, used below
1786(defun byte-mask (size)
1787  (1- (ash 1 (the fixnum size))))
1788
1789(defun byte-position (bytespec)
1790  "Return the position part of the byte specifier bytespec."
1791  (if (> bytespec 0)
1792    (- (integer-length bytespec) (logcount bytespec))
1793    (- bytespec)))
1794
1795
1796; CMU CL returns T.
1797(defun upgraded-complex-part-type (type)
1798  "Return the element type of the most specialized COMPLEX number type that
1799   can hold parts of type SPEC."
1800  (declare (ignore type))
1801  'real)
1802
1803;;; This is the MRG31k3p random number generator described in
1804;;; P. L'Ecuyer and R. Touzin, "Fast Combined Multiple Recursive
1805;;; Generators with Multipliers of the form a = +/- 2^d +/- 2^e",
1806;;; Proceedings of the 2000 Winter Simulation Conference, Dec. 2000,
1807;;; 683--689.
1808;;;
1809;;; A link to the paper is available on L'Ecuyer's web site:
1810;;; http://www.iro.umontreal.ca/~lecuyer/papers.html.
1811;;;
1812;;; This generator has a period of about 2^185.  It produces values in
1813;;; in the half-open interval [0, 2^31 - 1).
1814;;;
1815;;; It uses 6 words of state.
1816
1817(defconstant mrg31k3p-m1 #.(- (expt 2 31) 1))
1818(defconstant mrg31k3p-m2 #.(- (expt 2 31) 21069))
1819(defconstant mrg31k3p-limit #.(1- (expt 2 31))
1820             "Exclusive upper bound on values returned by %mrg31k3p.")
1821
1822
1823;;; This is a portable version of the MRG31k3p generator.  It's not
1824;;; too bad in a 64-bit CCL, but the generator pretty much has to be
1825;;; in LAP for 32-bit ports.
1826#-(or x8632-target ppc32-target x8664-target ppc64-target arm-target)
1827(defun %mrg31k3p (state)
1828  (let* ((v (random.mrg31k3p-state state)))
1829    (declare (type (simple-array (unsigned-byte 32) (*)) v)
1830             (optimize speed))
1831    (let ((y1 (+ (+ (ash (logand (aref v 1) #x1ff) 22)
1832                    (ash (aref v 1) -9))
1833                 (+ (ash (logand (aref v 2) #xffffff) 7)
1834                    (ash (aref v 2) -24)))))
1835      (declare (type (unsigned-byte 32) y1))
1836      (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1))
1837      (incf y1 (aref v 2))
1838      (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1))
1839      (setf (aref v 2) (aref v 1)
1840            (aref v 1) (aref v 0)
1841            (aref v 0) y1))
1842    (let ((y1 (+ (ash (logand (aref v 3) #xffff) 15)
1843                 (* 21069 (ash (aref v 3) -16))))
1844          (y2 (+ (ash (logand (aref v 5) #xffff) 15)
1845                 (* 21069 (ash (aref v 5) -16)))))
1846      (declare (type (unsigned-byte 32) y1 y2))
1847      (if (>= y1 mrg31k3p-m2) (decf y1 mrg31k3p-m2))
1848      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
1849      (incf y2 (aref v 5))
1850      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
1851      (incf y2 y1)
1852      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
1853      (setf (aref v 5) (aref v 4)
1854            (aref v 4) (aref v 3)
1855            (aref v 3) y2))
1856    (let* ((x10 (aref v 0))
1857           (x20 (aref v 3)))
1858      (if (<= x10 x20)
1859        (+ (- x10 x20) mrg31k3p-m1)
1860        (- x10 x20)))))
1861
1862(eval-when (:compile-toplevel :execute)
1863  (declaim (inline %16-random-bits)))
1864
1865(defun %16-random-bits (state)
1866  (logand #xffff (the fixnum (%mrg31k3p state))))
1867
1868#+64-bit-target
1869(defun %big-fixnum-random (number state)
1870  (declare (fixnum number)
1871           (ftype (function (random-state) fixnum) %mrg31k3p))
1872  (let ((low (ldb (byte 30 0) (%mrg31k3p state)))
1873        (high (ldb (byte 30 0) (%mrg31k3p state))))
1874    (declare (fixnum low high))
1875    (fast-mod (logior low (the fixnum (ash high 30)))
1876              number)))
1877
1878;;; When using a dead simple random number generator, it's reasonable
1879;;; to take 16 bits of the output and discard the rest.  With a more
1880;;; expensive generator, however, it may be worthwhile to do more bit
1881;;; fiddling here here so that we can use all of the random bits
1882;;; produced by %mrg31k2p.
1883#+32-bit-target
1884(defun %bignum-random (number state)
1885  (let* ((bits (+ (integer-length number) 8))
1886         (half-words (ash (the fixnum (+ bits 15)) -4))
1887         (long-words (ash (+ half-words 1) -1))
1888         (dividend (%alloc-misc long-words target::subtag-bignum))
1889         (16-bit-dividend dividend)
1890         (index 1))
1891    (declare (fixnum long-words index bits)
1892             (dynamic-extent dividend)
1893             (type (simple-array (unsigned-byte 16) (*)) 16-bit-dividend) ;lie
1894             (optimize (speed 3) (safety 0)))
1895    (loop
1896       ;; This had better inline due to the lie above, or it will error
1897       #+big-endian-target
1898       (setf (aref 16-bit-dividend index) (%16-random-bits state))
1899       #+little-endian-target
1900       (setf (aref 16-bit-dividend (the fixnum (1- index)))
1901             (%16-random-bits state))
1902       (decf half-words)
1903       (when (<= half-words 0) (return))
1904       #+big-endian-target
1905       (setf (aref 16-bit-dividend (the fixnum (1- index)))
1906             (%16-random-bits state))
1907       #+little-endian-target
1908       (setf (aref 16-bit-dividend index) (%16-random-bits state))
1909       (decf half-words)
1910       (when (<= half-words 0) (return))
1911       (incf index 2))
1912    ;; The bignum code expects normalized bignums
1913    (let* ((result (mod dividend number)))
1914      (if (eq dividend result)
1915        (copy-uvector result)
1916        result))))
1917
1918(defun %float-random (number state)
1919  (let ((ratio (gvector :ratio (random target::target-most-positive-fixnum state) target::target-most-positive-fixnum)))
1920    (declare (dynamic-extent ratio))
1921    (* number ratio)))
1922
1923(defun random (number &optional (state *random-state*))
1924  (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
1925  (cond
1926    ((and (fixnump number) (> (the fixnum number) 0))
1927     #+32-bit-target
1928     (fast-mod (%mrg31k3p state) number)
1929     #+64-bit-target
1930     (if (< number mrg31k3p-limit)
1931       (fast-mod (%mrg31k3p state) number)
1932       (%big-fixnum-random number state)))
1933    ((and (typep number 'double-float) (> (the double-float number) 0.0))
1934     (%float-random number state))
1935    ((and (typep number 'short-float) (> (the short-float number) 0.0s0))
1936     (%float-random number state))
1937    ((and (bignump number) (> number 0))
1938     (%bignum-random number state))
1939    (t (report-bad-arg number '(or (integer (0)) (float (0.0)))))))
1940
1941(eval-when (:compile-toplevel :execute)
1942  (defmacro bignum-abs (nexp)
1943    (let ((n (gensym)))
1944      `(let ((,n ,nexp))
1945         (if  (bignum-minusp ,n) (negate-bignum ,n) ,n))))
1946 
1947  (defmacro fixnum-abs (nexp)
1948    (let ((n (gensym)))
1949      `(let ((,n ,nexp))
1950         (if (minusp (the fixnum ,n))
1951           (if (eq ,n target::target-most-negative-fixnum)
1952             (- ,n)
1953             (the fixnum (- (the fixnum ,n))))
1954           ,n))))
1955  )
1956 
1957
1958;;; TWO-ARG-GCD  --  Internal
1959;;;
1960;;;    Do the GCD of two integer arguments.  With fixnum arguments, we use the
1961;;; binary GCD algorithm from Knuth's seminumerical algorithms (slightly
1962;;; structurified), otherwise we call BIGNUM-GCD.  We pick off the special case
1963;;; of 0 before the dispatch so that the bignum code doesn't have to worry
1964;;; about "small bignum" zeros.
1965;;;
1966(defun gcd-2 (n1 n2)
1967  ;(declare (optimize (speed 3)(safety 0)))
1968  (cond 
1969   ((eql n1 0) (%integer-abs n2))
1970   ((eql n2 0) (%integer-abs n1))
1971   (t (number-case n1
1972        (fixnum 
1973         (number-case n2
1974          (fixnum
1975           (if (eql n1 target::target-most-negative-fixnum)
1976             (if (eql n2 target::target-most-negative-fixnum)
1977               (- target::target-most-negative-fixnum)
1978               (bignum-fixnum-gcd (- target::target-most-negative-fixnum) (abs n2)))
1979             (if (eql n2 target::target-most-negative-fixnum)
1980               (bignum-fixnum-gcd (- target::target-most-negative-fixnum) (abs n1))
1981               (locally
1982                   (declare (optimize (speed 3) (safety 0))
1983                            (fixnum n1 n2))
1984                 (if (minusp n1)(setq n1 (the fixnum (- n1))))
1985                 (if (minusp n2)(setq n2 (the fixnum (- n2))))
1986               (%fixnum-gcd n1 n2)))))
1987           (bignum (if (eql n1 target::target-most-negative-fixnum)
1988                     (%bignum-bignum-gcd n2 (- target::target-most-negative-fixnum))
1989                     (bignum-fixnum-gcd (bignum-abs n2)(fixnum-abs n1))))))
1990        (bignum
1991         (number-case n2
1992           (fixnum
1993            (if (eql n2 target::target-most-negative-fixnum)
1994              (%bignum-bignum-gcd (bignum-abs n1)(fixnum-abs n2))
1995              (bignum-fixnum-gcd (bignum-abs n1)(fixnum-abs n2))))
1996           (bignum (%bignum-bignum-gcd n1 n2))))))))
1997
1998#|
1999(defun fixnum-gcd (n1 n2)
2000  (declare (optimize (speed 3) (safety 0))
2001           (fixnum n1 n2))                   
2002  (do* ((k 0 (%i+ 1 k))
2003        (n1 n1 (%iasr 1 n1))
2004        (n2 n2 (%iasr 1 n2)))
2005       ((oddp (logior n1 n2))
2006        (do ((temp (if (oddp n1) (the fixnum (- n2)) (%iasr 1 n1))
2007                   (%iasr 1 temp)))
2008            (nil)
2009          (declare (fixnum temp))
2010          (when (oddp temp)
2011            (if (plusp temp)
2012              (setq n1 temp)
2013              (setq n2 (- temp)))
2014            (setq temp (the fixnum (- n1 n2)))
2015            (when (zerop temp)
2016              (let ((res (%ilsl k n1)))
2017                (return res))))))
2018    (declare (fixnum n1 n2 k))))
2019|#
2020
2021
2022
2023(defun %quo-1 (n)
2024  (/ 1 n))
2025
2026;; Compute (sqrt (+ (* x x) (* y y))), but
2027;; try to be a little more careful about it.
2028;; Both x and y must be double-floats.
2029(defun %double-float-hypot (x y)
2030  (with-stack-double-floats ((a) (b) (c))
2031    (%%double-float-abs! x a)
2032    (%%double-float-abs! y b)
2033    (when (> a b)
2034      (psetq a b b a))
2035    (if (= b 0d0)
2036      0d0
2037      (progn
2038        (%double-float/-2! a b c)
2039        (* b (fsqrt (+ 1d0 (* c c))))))))
2040                                       
Note: See TracBrowser for help on using the repository browser.