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

Last change on this file since 10152 was 10152, checked in by rme, 12 years ago

32-bit %BIGNUM-RANDOM: account for endianness (so that random produces
the same results on ppc32 and x8632).

FIXNUM-SFLOAT-COMPARE, TRUNCATE, TRUNCATE-NO-REM: conditionalize on
32-bit-target rather than ppc32-target.

most-positive-fixnum => target::target-most-positive-fixnum
most-negative-fixnum => target::target-most-negative-fixnum

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