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

Last change on this file since 11142 was 11142, checked in by gb, 11 years ago

Don't try to define two-arg + and - via macros; use ADD-BIGNUM-AND-FIXNUM
as appropriate.

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