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

Last change on this file was 16685, checked in by rme, 4 years ago

Update copyright/license headers in files.

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