source: trunk/ccl/lib/numbers.lisp @ 5682

Last change on this file since 5682 was 5682, checked in by gb, 14 years ago

Explicitly require NUMBER-CASE-MACRO.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.8 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;; Lib;numbers.lisp - Lisp arithmetic code.
18
19(in-package "CCL")
20
21(eval-when (:compile-toplevel :execute)
22 (require :number-macros)
23 (require :number-case-macro)
24 #+(and cross-compiling 64-bit-target)
25 (declaim (ftype function %single-float-atanh %single-float-acosh
26                 %single-float-asinh %single-float-tanh
27                 %single-float-cosh %single-float-sinh)))
28
29
30
31(defconstant double-float-positive-infinity
32  #.(let* ((division-by-zero (get-fpu-mode  :division-by-zero)))
33      (declare (notinline /))
34      (unwind-protect
35           (progn
36             (ccl:set-fpu-mode :division-by-zero nil)
37             (/ 0d0))
38        (ccl:set-fpu-mode :division-by-zero division-by-zero))))
39
40(defconstant double-float-negative-infinity
41  #.(let* ((division-by-zero (get-fpu-mode  :division-by-zero)))
42      (declare (notinline /))
43      (unwind-protect
44           (progn
45             (ccl:set-fpu-mode :division-by-zero nil)
46             (/ -0d0))
47        (ccl:set-fpu-mode :division-by-zero division-by-zero))))
48
49
50
51(defun parse-float (str len off) 
52  ; we cant assume this really is a float but dont call with eg s1 or e1
53  (let ((integer 0)(expt 0)(sign 0)(done 0)(digits 0) point-pos type) 
54    (setq integer
55          (do ((n off (1+ n))
56               (first t nil)
57               (maxn  (+ off len)))
58              ((>= n maxn) integer)
59            (declare (fixnum n maxn))
60            (let ((c (%schar str n)))
61              (cond ((eq c #\.)
62                     (setq point-pos digits))
63                    ((and first (eq c #\+)))
64                    ((and first (eq c #\-))
65                     (setq sign -1))
66                    ((memq c '(#\s #\f #\S #\F))
67                     (setq type 'short-float)
68                     (return integer))
69                    ((memq c '(#\d #\l  #\D  #\L))
70                     (setq type 'double-float)
71                     (return integer))
72                    ((memq c '(#\e #\E))
73                     (return integer))
74                    ((setq c (digit-char-p c))
75                     (setq digits (1+ digits))
76                     (setq integer (+ c (* 10 integer))))                 
77                    (t (return-from parse-float nil)))
78              (setq done (1+ done)))))
79    (when point-pos
80      (setq expt  (%i- point-pos digits)))
81    (when (null type)
82      (setq type *read-default-float-format*))
83    (when (> len done)
84      (let ((eexp nil) (inf nil) (nan nil) (esign 1) c (xsign-n -1))
85        (do ((n (%i+ off done 1) (1+ n))
86             (first t nil))
87            ((>= n (+ off len)))
88          (declare (fixnum n))
89          (setq c (%schar str n))
90          (cond ((and first (or (eq c #\+)(eq c #\-)))
91                 (when (eq c #\-)(setq esign -1))
92                 (setq xsign-n (1+ n)))
93                ((and (= n xsign-n)
94                      (or (eq c #\+)(eq c #\-)))
95                 (if (eq c #\-)
96                     (setq nan t)
97                     (setq inf t)))
98                ((setq c (digit-char-p c))
99                 (setq eexp (+ c (* (or eexp 0) 10))))
100                (t (return-from parse-float nil))))
101        (when (not eexp)(return-from parse-float nil))
102        (cond 
103         (inf 
104          (return-from parse-float
105            (coerce (if (minusp sign)
106                        double-float-negative-infinity
107                        double-float-positive-infinity)
108                    type)))
109         (nan 
110          (return-from parse-float
111            (let* ((invalid (ccl:get-fpu-mode :invalid)))
112            (unwind-protect
113                (progn
114                  (ccl:set-fpu-mode :invalid nil)
115                  (coerce
116                   ;; we could also have used a double-float-nan
117                   ;; variable binding here:
118                   (+ double-float-positive-infinity
119                      double-float-positive-infinity)
120                   type))
121              (ccl:set-fpu-mode :invalid invalid)))))
122         (expt (setq expt (%i+ expt (* esign eexp))))
123         (t (return-from parse-float nil)))))
124    (fide sign integer expt (subtypep type 'short-float))))
125
126
127;; an interesting test case: 1.448997445238699
128;; The correct result is 6525704354437805 x 2^-52
129;; Incorrect is          6525704354437806 x 2^-52
130;; (from Will Clinger, "How to Read Floating Point Numbers Accurately",
131;;  ACM SIGPLAN'90 Conference on Programming Language Design and Implementation")
132;; Doug Curries numbers 214748.3646, 1073741823/5000
133
134
135;; Sane read losers
136;; 15871904747836473438871.0e-8
137;; 3123927307537977993905.0-13
138;; 17209940865514936528.0e-6
139;; "13.60447536e132" => adds some gratuitous drech
140;; "94824331561426550.889e182"
141;; "1166694.64175277e-150" => 1.1666946417527701E-144
142;; "3109973217844.55680988601e-173"
143;; "817332.e-184" => 8.173320000000001E-179
144;; "2695.13e-180" => 2.6951300000000002E-177
145;; "89.85345789e-183" => 8.985345789000001E-182
146;; "0864813880.29e140" => 8.648138802899999E+148
147;; "5221.e-193" => 5.2209999999999995E-190
148;; "7.15628e-175" => 7.156280000000001E-175
149
150(defparameter float-powers-of-5  nil)
151(defparameter integer-powers-of-5 nil)
152
153(defun 5-to-e (e)
154  (declare (fixnum e)(optimize (speed 3)(safety 0)))
155  (if (> e 335)
156    (* (5-to-e 335)(5-to-e (- e 335))) ; for the dude who types 200 digits and e-500
157    (if (< e 12)
158      (svref integer-powers-of-5 e)
159      (multiple-value-bind (q r) (truncate e 12) ; was floor
160        (declare (fixnum q r))       
161        (if (eql r 0)
162          (svref integer-powers-of-5 (%i+ q 11))
163          (* (svref integer-powers-of-5 r)
164             (svref integer-powers-of-5 (%i+ q 11))))))))
165
166(defun float-5-to-e (e)
167  (if (> e 22)  ; shouldnt happen
168    (expt 5.0d0 e)
169    (svref float-powers-of-5 e)))
170
171(defparameter a-short-float nil)
172
173(eval-when (:compile-toplevel :execute)
174  ; number of bits for mantissa before rounding
175  (defconstant *short-float-extended-precision* 28)
176  (defconstant *double-float-extended-precision* 60)
177  ; number of mantissa bits including hidden bit
178  (defconstant *double-float-precision* (1+ IEEE-double-float-mantissa-width))
179  (defconstant *short-float-precision* (1+ IEEE-single-float-mantissa-width))
180  (defconstant *double-float-bias* IEEE-double-float-bias)
181  (defconstant *double-float-max-exponent* (1+ IEEE-double-float-normal-exponent-max))
182  (defconstant *double-float-max-exact-power-of-5* 23)
183  ;(defconstant *short-float-max-exact-integer-length* 24)
184  (defconstant *double-float-max-exact-integer-length* 53)
185)
186
187
188
189
190(eval-when (:compile-toplevel :execute)
191  (defconstant *short-float-max-exact-power-of-5* 10)
192  (defconstant *short-float-bias* IEEE-single-float-bias)
193  (defconstant *short-float-max-exact-integer-length* 24)
194  (defconstant *short-float-max-exponent* (1+ IEEE-single-float-normal-exponent-max))
195)
196
197 
198;; this stuff  could be in a shared file
199
200(defun fide #|float-integer-with-decimal-exponent|# (sign integer power-of-10 &optional short)
201  ;; take care of the zero case
202  (when (zerop integer)
203    (return-from fide ;float-integer-with-decimal-exponent
204       (if short
205         (if (minusp sign) -0.0s0 0.0s0)
206         (if (minusp sign) -0.0d0 0.0d0))))
207  (let ((abs-power (abs power-of-10))
208        (integer-length (integer-length integer)))
209    ;; this doesn't work for the above example, so arithmetic must be done wrong
210    ;; This does work if set FPCR precision to double
211    ;; now see if the conversion can be done simply:
212    ;; if both the integer and the power of 10 can be floated exactly, then
213    ;; correct rounding can be done by the multiply or divide
214    (when (or;short
215           (and (<= integer-length 
216                    ;; was (if short 17 53) why 17? see above
217                    (if short *short-float-max-exact-integer-length* *double-float-max-exact-integer-length*)) 
218                ;; (integer-length (expt 5 23)) => 54
219                ;; was (if short 5 23)
220                (< abs-power  (if short 
221                                *short-float-max-exact-power-of-5*
222                                *double-float-max-exact-power-of-5*)))) ; we mean < 23 not <=
223      ;; if you care about consing, this could be done in assembly language or whatever,
224      ;; since all integers fit in 53 bits
225      (return-from fide ;float-integer-with-decimal-exponent
226        (let* ((signed-integer (prog1 (if (minusp sign) (- integer) integer)))
227               (float (float signed-integer (if short 0.0s0 0.0d0)))
228               (10-to-power (scale-float (float-5-to-e abs-power) abs-power)))
229          ;; coerce to short-float does not whine about undeflow, but does re overflow
230          (when short (setq 10-to-power (coerce 10-to-power 'short-float)))
231          (if (zerop abs-power)
232            float
233            (if (minusp power-of-10)
234              (/ float  10-to-power)
235              (* float  10-to-power))))))
236    (try-harder sign integer power-of-10 short)))
237
238
239(defun try-harder (sign integer power-of-10 short)
240  (flet ((ovf (&optional under)
241           (if under
242             (if (get-fpu-mode :underflow)
243               (error 'floating-point-underflow
244                      :operation 'scale
245                      :operands (list sign integer power-of-10)))
246             (if (get-fpu-mode :overflow)
247               (error 'floating-point-overflow
248                      :operation 'scale
249                      :operands (list sign integer power-of-10))))
250           (return-from try-harder
251             (if under
252               (if short
253                 (if (minusp sign) -0.0s0 0.0s0)                 
254                 (if (minusp sign) 0.0d0 0.0d0))
255               (if short
256                 (if (minusp sign) most-negative-short-float most-positive-short-float)             
257                 (if (minusp sign) most-negative-double-float most-positive-double-float))))))
258  (let* ((integer-length (integer-length integer)) new-int power-of-2)
259    (if (minusp power-of-10)
260      (progn 
261        ;; avoid creating enormous integers with 5-to-e only to error later
262        (when (< power-of-10 -335)
263          (let ((poo (+ (round integer-length 3.2) power-of-10)))
264            ;; overestimate digits in integer
265            (when (< poo -335) (ovf t))
266            ;; this case occurs if 600+ digits
267            (when (> poo 335) (ovf))))
268        (let* ((divisor (5-to-e (- power-of-10)))
269               ;; make sure we will have enough bits in the quotient
270               ;; (and a couple extra for rounding)
271               (shift-factor (+ (- (integer-length divisor) integer-length)
272                                (if short *short-float-extended-precision* *double-float-extended-precision*)))
273               (scaled-integer integer))
274          (if (plusp shift-factor)
275            (setq scaled-integer (ash integer shift-factor))
276            (setq divisor (ash divisor (- shift-factor))))
277          (multiple-value-bind (quotient remainder)(floor scaled-integer divisor)
278            (unless (zerop remainder) ; whats this - tells us there's junk below
279              (setq quotient (logior quotient 1)))
280            (setq new-int quotient)
281            (setq power-of-2  (- power-of-10 shift-factor)))))
282      (progn
283        (when (> power-of-10 335)(ovf))
284        (setq new-int (* integer (5-to-e power-of-10)))
285        (setq power-of-2 power-of-10)))
286    (float-and-scale-and-round sign new-int power-of-2 short))))
287
288
289(defun float-and-scale-and-round (sign integer power-of-2 short &optional result)
290  (let* ((length (integer-length integer))
291         (lowbits 0)
292         (prec (if short *short-float-precision* *double-float-precision*))
293         (ep (if short *short-float-extended-precision* *double-float-extended-precision*)))
294    (when (<= length prec)
295      ;; float can be done exactly, so do it the easy way
296      (return-from float-and-scale-and-round
297        (scale-float (float (if (minusp sign) (- integer) integer) (if short a-short-float))
298                     power-of-2)))   
299    (let* ((exponent (+ length power-of-2))
300           (biased-exponent (+ exponent (if short *short-float-bias* *double-float-bias*)))
301           (sticky-residue nil))
302      (cond
303       ((<= biased-exponent 0)
304        ;; denormalize the number
305        (setf sticky-residue (not (zerop (ldb integer (byte (- 1 biased-exponent) 0)))))
306        (setf integer (ash integer (- biased-exponent 1)))
307        (setf biased-exponent 0)))
308      (let ((lowest (min ep length)))
309        (when (and (> length ep)(not (zerop (ldb (byte (- length ep) 0) integer))))
310          (setq integer (logior integer (ash 1 (- length ep)))))
311        ; somewhere between 1 and (- ep prec) bits
312        (setq lowbits (ash (ldb (byte (- lowest prec) (- length lowest)) integer) (- ep lowest))))
313      (let* ((significand (ldb (byte (1- prec) (- length prec)) integer)))
314        (when (and (not (zerop (ldb (byte 1 (- length (1+ prec))) integer)))   ; round bit
315                   (or sticky-residue (oddp significand)
316                       (not (zerop (ldb (byte (- ep prec 1) 0) lowbits)))))
317          ;; round up
318          (setf significand (ldb (byte (1- prec) 0) (+ significand 1)))
319          (when (zerop significand)
320            (incf biased-exponent)))
321        (cond ((and (zerop biased-exponent)
322                    (zerop significand)
323                    (get-fpu-mode :underflow))
324               (error 'floating-point-underflow
325                      :operation 'scale
326                      :operands (list sign integer power-of-2)))
327              ((>= biased-exponent (if short *short-float-max-exponent* *double-float-max-exponent*))
328               (cond 
329                     (t
330                      (if (get-fpu-mode :overflow)
331                        (error 'floating-point-overflow
332                               :operation 'scale
333                               :operands (list sign integer power-of-2)))
334                      (setf significand 0)                     
335                      (setq biased-exponent (if short *short-float-max-exponent* *double-float-max-exponent*))))))
336        (values
337         (if short 
338           (make-short-float-from-fixnums (ldb (byte 23 0) significand)
339                                          biased-exponent
340                                          sign #-64-bit-target result)
341           (make-float-from-fixnums (ldb (byte 24 28) significand)
342                                    (ldb (byte 28 0) significand)
343                                    biased-exponent
344                                    sign result))
345         lowbits)))))
346
347
348
349
350(defparameter a-short-float 1.0s0)
351
352
353(defmethod print-object ((rs random-state) stream)
354  (let* ((byte (byte 16 (- 16 target::fixnum-shift))))
355    (format stream "#.(~S ~S ~S)"   ;>> #.GAG!!!
356            'ccl::initialize-random-state
357            (ldb byte (%svref rs 1))
358            (ldb byte (%svref rs 2)))))
359
360
361
362(defun float-radix (float)
363  "Return (as an integer) the radix b of its floating-point argument."
364  (require-type float 'float)
365  2)
366
367(defun float-digits (float)
368  (if (typep (require-type float 'float) 'short-float)
369    IEEE-single-float-digits
370    IEEE-double-float-digits))
371
372(defun number-arg (arg)
373  (if (numberp arg) arg (%badarg arg 'number)))
374
375
376
377
378
379;==> Needs a transform...
380(defun logandc2 (integer1 integer2)
381  "Bitwise AND INTEGER1 with (LOGNOT INTEGER2)."
382  (logandc1 integer2 integer1))
383
384(defun logorc2 (integer1 integer2)
385  "Bitwise OR INTEGER1 with (LOGNOT INTEGER2)."
386  (logorc1 integer2 integer1))
387
388
389
390; Figure that the common (2-arg) case is caught by a compiler transform anyway.
391(defun gcd (&lexpr numbers)
392  "Return the greatest common divisor of the arguments, which must be
393  integers. Gcd with no arguments is defined to be 0."
394  (let* ((count (%lexpr-count numbers)))
395    (declare (fixnum count))   
396    (if (zerop count)
397      0
398      (let* ((n0 (%lexpr-ref numbers count 0)))
399        (if (= count 1)
400          (%integer-abs n0)
401          (do* ((i 1 (1+ i)))
402               ((= i count) n0)
403            (declare (fixnum i))
404            (setq n0 (gcd-2 n0 (%lexpr-ref numbers count i)))))))))
405
406(defun lcm-2 (n0 n1)
407  (or (typep n0 'integer) (report-bad-arg n0 'integer))
408  (or (typep n1 'integer) (report-bad-arg n1 'integer))
409  (locally (declare (integer n0 n1))
410    (if (zerop n0)
411      0
412      (if (zerop n1)
413        0
414        (let* ((small (if (< n0 n1) n0 n1))
415               (large (if (eq small n0) n1 n0)))
416          (values (truncate (abs (* n0 n1)) (gcd large small))))))))
417
418(defun lcm (&lexpr numbers)
419  "Return the least common multiple of one or more integers. LCM of no
420  arguments is defined to be 1."
421  (let* ((count (%lexpr-count numbers)))
422    (declare (fixnum count))   
423    (if (zerop count)
424      1
425      (let* ((n0 (%lexpr-ref numbers count 0)))
426        (if (= count 1)
427          (%integer-abs n0)
428          (if (= count 2)
429            (lcm-2 n0 (%lexpr-ref numbers count 1))
430            (do* ((i 1 (1+ i)))
431                 ((= i count) n0)
432              (declare (fixnum i))
433              (setq n0 (lcm-2 n0 (%lexpr-ref numbers count i))))))))))
434
435
436#|
437(defun rationalize (x)
438  (etypecase x
439    (rational x)
440    (real
441     (cond ((minusp x) (- (rationalize (- x))))
442           ((zerop x) 0)
443           (t
444            (let ((eps (etypecase x
445                         (single-float single-float-epsilon)
446                         (double-float double-float-epsilon)))
447                  (y ())
448                  (a ()))
449              (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
450                   (num (setq a (truncate x))
451                        (+ (* (setq a (truncate y)) num) onum))
452                   (den 1 (+ (* a den) oden))
453                   (onum 1 num)
454                   (oden 0 den))
455                  ((and (not (zerop den))
456                        (not (> (abs (/ (- x (/ (float num x)
457                                                (float den x)))
458                                        x))
459                                eps)))
460                   (integer-/-integer num den)))))))))
461|#
462
463(defun rationalize (number)
464  "Converts any REAL to a RATIONAL.  Floats are converted to a simple rational
465  representation exploiting the assumption that floats are only accurate to
466  their precision.  RATIONALIZE (and also RATIONAL) preserve the invariant:
467      (= x (float (rationalize x) x))"
468  (if (floatp number)
469    (labels ((simpler-rational (less-predicate lonum loden hinum hiden
470                                               &aux (trunc (if (eql less-predicate #'<=)
471                                                             #'ceiling
472                                                             #'(lambda (n d) (1+ (floor n d)))))
473                                               (term (funcall trunc lonum loden)))
474               ;(pprint (list lonum loden hinum hiden))
475               (if (funcall less-predicate (* term hiden) hinum)
476                 (values term 1)
477                 (multiple-value-bind 
478                   (num den)
479                   (simpler-rational less-predicate hiden (- hinum (* (1- term) hiden))
480                                     loden (- lonum (* (1- term) loden)))
481                   (values (+ den (* (1- term) num)) num)))))                           
482      (multiple-value-bind (fraction exponent sign) (integer-decode-float number)
483        ;; the first 2 tests may be unnecessary - I think the check
484        ;; for denormalized is compensating for a bug in 3.0 re
485        ;; floating a rational (in order to pass tests in
486        ;; ppc-test-arith).
487        (if (or (and (typep number 'double-float)  ; is it denormalized
488                     (eq exponent #.(nth-value 1 (integer-decode-float least-positive-double-float)))) ; aka -1074))
489                (eq exponent #.(nth-value 1 (integer-decode-float least-positive-short-float))) ; aka -149))
490                (zerop (logand fraction (1- fraction)))) ; or a power of two
491          (rational number)
492          (if (minusp exponent)
493            ;;less than 1
494            (let ((num (ash fraction 2))
495                  (den (ash 1 (- 2 exponent))))
496              (multiple-value-bind 
497                (n d)
498                (simpler-rational (if (evenp fraction) #'<= #'<)
499                                  (- num 2) ;(if (zerop (logand fraction (1- fraction))) 1 2))
500                                  den  (+ num 2) den)
501                (when (minusp sign)
502                  (setq n (- n)))
503                (/ n d)))
504            ;;greater than 1
505            (ash (if (minusp number) (- fraction) fraction) exponent)))))
506    (rational number)))
507#|
508(defun testrat (&optional (n 1000))
509  (dotimes (i n)
510    (let* (( numerator (random (ash 1 63)))
511          (denominator (random (ash 1 63)))
512          (sign  (if (zerop (random 2)) 1 -1))
513          (trial (float (/ (* sign numerator) denominator)))
514          (rat (rationalize trial)))
515      (when (not (= (float rat) trial))
516        (error "Rationalize failed. Input ~s Rational ~s Float ~s" trial rat (float rat))))))
517
518; smallest fails in 3.0 - powers of 2 - works here but we cheat a bit
519(defun testrat2 ()
520  (let ((f least-positive-double-float))
521    (dotimes (i 100)
522      (when (not (= (float (rationalize f)) f))
523        (cerror "a" "rat failed ~s ~s" f i))
524      (setq f (* f 2)))))
525
526; fails a lot in 3.0 - not powers of 2 - works here
527(defun testrat3 ()
528  (let ((f least-positive-double-float))
529    (dotimes (i 1000)
530      (let ((f2 (* (+ i i 1) f)))
531        (when (not (= (float (rationalize f2)) f2))
532          (cerror "a" "rat failed ~s ~s" f2 i)))))
533  (let ((f least-negative-double-float))
534    (dotimes (i 1000)
535      (let ((f2 (* (+ i i 1) f)))
536        (when (not (= (float (rationalize f2)) f2))
537          (cerror "a" "rat failed ~s ~s" f2 i))))))
538
539(defun testrat31 ()
540  (let ((f least-positive-short-float))
541    (dotimes (i 1000)
542      (let ((f2 (* (+ i i 1) f)))
543        (when (not (= (float (rationalize f2)) f2))
544          (cerror "a" "rat failed ~s ~s" f2 i)))))
545  (let ((f least-negative-short-float))
546    (dotimes (i 1000)
547      (let ((f2 (* (+ i i 1) f)))
548        (when (not (= (float (rationalize f2)) f2))
549          (cerror "a" "rat failed ~s ~s" f2 i))))))
550
551; works in 3.0 - and here
552(defun testrat4 ()
553  (let ((f least-positive-normalized-double-float))
554    (dotimes (i 1000)
555      (let ((f2 (* (+ i i 1) f)))
556        (when (not (= (float (rationalize f2)) f2))
557          (cerror "a" "rat failed ~s ~s" f2 i)))))
558  (let ((f least-negative-normalized-double-float))
559    (dotimes (i 100)
560      (let ((f2 (* (+ i i 1) f)))
561        (when (not (= (float (rationalize f2)) f2))
562          (cerror "a" "rat failed ~s ~s" f2 i))))))
563       
564   
565|#
566
567#| now in l1-numbers.lisp
568(defun logeqv (&lexpr numbers)
569  (let* ((count (%lexpr-count numbers)))
570    (declare (fixnum count))
571    (if (zerop count)
572      -1
573      (let* ((n0 (%lisp-word-ref numbers count)))
574        (if (= count 1)
575          (require-type n0 'integer)
576          (do* ((i 1 (1+ i)))
577               ((= i count) n0)
578            (declare (fixnum i))
579            (declare (optimize (speed 3) (safety 0)))
580            (setq n0 (logeqv-2 (%lexpr-ref numbers count i) n0))))))))
581|#
582
583
584(defparameter *boole-ops* 
585  (vector
586   #'(lambda (i1 i2) (declare (ignore i1 i2)) 0)
587   #'(lambda (i1 i2) (declare (ignore i1 i2)) -1)
588   #'(lambda (i1 i2) (declare (ignore i2)) i1)
589   #'(lambda (i1 i2) (declare (ignore i1)) i2)
590   #'(lambda (i1 i2) (declare (ignore i2)) (lognot i1))
591   #'(lambda (i1 i2) (declare (ignore i1)) (lognot i2))
592   #'(lambda (i1 i2) (logand i1 i2))
593   #'(lambda (i1 i2) (logior i1 i2))
594   #'(lambda (i1 i2) (logxor i1 i2))
595   #'(lambda (i1 i2) (logeqv i1 i2))
596   #'(lambda (i1 i2) (lognand i1 i2))
597   #'(lambda (i1 i2) (lognor i1 i2))
598   #'(lambda (i1 i2) (logandc1 i1 i2))
599   #'(lambda (i1 i2) (logandc2 i1 i2))
600   #'(lambda (i1 i2) (logorc1 i1 i2))
601   #'(lambda (i1 i2) (logorc2 i1 i2))))
602 
603
604
605;===> Change these constants to match maclisp!!
606(defun boole (op integer1 integer2)
607  "Bit-wise boolean function on two integers. Function chosen by OP:
608        0       BOOLE-CLR
609        1       BOOLE-SET
610        2       BOOLE-1
611        3       BOOLE-2
612        4       BOOLE-C1
613        5       BOOLE-C2
614        6       BOOLE-AND
615        7       BOOLE-IOR
616        8       BOOLE-XOR
617        9       BOOLE-EQV
618        10      BOOLE-NAND
619        11      BOOLE-NOR
620        12      BOOLE-ANDC1
621        13      BOOLE-ANDC2
622        14      BOOLE-ORC1
623        15      BOOLE-ORC2"
624  (unless (and (typep op 'fixnum)
625               (locally (declare (fixnum op))
626                 (and (>= op 0)
627                      (<= op 15))))
628    (report-bad-arg op '(integer 0 15)))
629  (funcall (%svref *boole-ops* op)
630           (require-type integer1 'integer)
631           (require-type integer2 'integer)))
632
633
634(defun %integer-power (b e)
635  (declare (type unsigned-byte e))
636  (if (zerop e)
637    (+ 1 (* b 0))
638    (if (eql b 2)
639      (ash 1 e)
640      (do* ((next (ash e -1) (ash e -1))
641            (oddexp (oddp e) (oddp e))
642            (total (if oddexp b 1) (if oddexp (* b total) total)))
643           ((zerop next) total)
644        (declare (type unsigned-byte next))
645        (setq b (* b b) e next)))))
646
647(defun signum (x)
648  "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
649  (cond ((complexp x) (if (zerop x) x (/ x (abs x))))
650        ((rationalp x) (if (plusp x) 1 (if (zerop x) 0 -1)))
651        ((zerop x) (float 0.0 x))
652        (t (float-sign x))))
653
654
655
656; Thanks to d34676@tansei.cc.u-tokyo.ac.jp (Akira KURIHARA)
657(defun isqrt (n &aux n-len-quarter n-half n-half-isqrt
658                init-value iterated-value)
659  "Return the root of the nearest integer less than n which is a perfect
660   square.  Argument n must be a non-negative integer"
661  (cond
662   ((eql n 0) 0)
663   ; this fails sometimes - do we care? 70851992595801818865024053174 or #x80000000
664   ; maybe we do - its used by dotimes
665   ;((not (int>0-p n)) (report-bad-arg n '(integer 0))) ;'unsigned-byte)) ; Huh?
666   ((or (not (integerp n))(minusp n))(report-bad-arg n '(integer 0)))
667   ((> n 24)            ; theoretically (> n 7) ,i.e., n-len-quarter > 0
668    (setq n-len-quarter (ash (integer-length n) -2))
669    (setq n-half (ash n (- (ash n-len-quarter 1))))
670    (setq n-half-isqrt (isqrt n-half))
671    (setq init-value (ash (1+ n-half-isqrt) n-len-quarter))
672    (loop
673      (setq iterated-value (ash (+ init-value (floor n init-value)) -1))
674      (if (not (< iterated-value init-value))
675        (return init-value)
676        (setq init-value iterated-value))))
677   ((> n 15) 4)
678   ((> n  8) 3)
679   ((> n  3) 2)
680   (t 1)))
681
682
683(defun sinh (x)
684  "Return the hyperbolic sine of NUMBER."
685  (if (complexp x) 
686    (/ (- (exp x) (exp (- x))) 2)
687    (if (typep x 'double-float)
688      (%double-float-sinh! x (%make-dfloat))
689      #+32-bit-target
690      (ppc32::with-stack-short-floats ((sx x))
691        (%single-float-sinh! sx (%make-sfloat)))
692      #+64-bit-target
693        (%single-float-sinh (%short-float x)))))
694
695
696(defun cosh (x)
697  "Return the hyperbolic cosine of NUMBER."
698  (if (complexp x) 
699    (/ (+ (exp x) (exp (- x))) 2)
700    (if (typep x 'double-float)
701      (%double-float-cosh! x (%make-dfloat))
702      #+32-bit-target
703      (ppc32::with-stack-short-floats ((sx x))
704        (%single-float-cosh! sx (%make-sfloat)))
705      #+64-bit-target
706      (%single-float-cosh (%short-float x)))))
707
708(defun tanh (x)
709  "Return the hyperbolic tangent of NUMBER."
710  (if (complexp x) 
711    (/ (sinh x) (cosh x))
712    (if (typep x 'double-float)
713      (%double-float-tanh! x (%make-dfloat))
714      #+32-bit-target
715      (ppc32::with-stack-short-floats ((sx x))
716        (%single-float-tanh! sx (%make-sfloat)))
717      #+64-bit-target
718      (%single-float-tanh (%short-float x)))))
719
720(defun asinh (x)
721  "Return the hyperbolic arc sine of NUMBER."
722  (if (complexp x) 
723    (log (+ x (sqrt (+ 1 (* x x)))))
724    (if (typep x 'double-float)
725      (%double-float-asinh! x (%make-dfloat))
726      #+32-bit-target
727      (ppc32::with-stack-short-floats ((sx x))
728        (%single-float-asinh! sx (%make-sfloat)))
729      #+64-bit-target
730      (%single-float-asinh (%short-float x)))))
731
732(defun acosh (x)
733  "Return the hyperbolic arc cosine of NUMBER."
734  (if (and (realp x) (<= 1.0 x))
735    (if (typep x 'double-float)
736      (%double-float-acosh! x (%make-dfloat))
737      #+32-bit-target
738      (ppc32::with-stack-short-floats ((sx x))
739        (%single-float-acosh! sx (%make-sfloat)))
740      #+64-bit-target
741      (%single-float-acosh (%short-float x)))
742    (* 2 (log (+ (sqrt (/ (1+ x) 2)) (sqrt (/ (1- x) 2)))))))
743
744(defun atanh (x)
745  "Return the hyperbolic arc tangent of NUMBER."
746  (if (and (realp x) (<= -1.0 (setq x (float x)) 1.0))
747    (if (typep x 'double-float)
748      (%double-float-atanh! x (%make-dfloat))
749      #+32-bit-target
750      (%single-float-atanh! x (%make-sfloat))
751      #+64-bit-target
752      (%single-float-atanh x))
753    (/ (log (/ (+ 1 x) (- 1 x))) 2)))
754
755
756(defun ffloor (number &optional divisor)
757  "Same as FLOOR, but returns first value as a float."
758  (multiple-value-bind (q r) (floor number divisor)
759    (values (float q (if (floatp r) r 0.0)) r)))
760
761(defun fceiling (number &optional divisor)
762  "Same as CEILING, but returns first value as a float."
763  (multiple-value-bind (q r) (ceiling number divisor)
764    (values (float q (if (floatp r) r 0.0)) r)))
765
766(defun ftruncate (number &optional divisor)
767  "Same as TRUNCATE, but returns first value as a float."
768  (multiple-value-bind (q r) (truncate number divisor)
769    (values (float q (if (floatp r) r 0.0)) r)))
770
771(defun fround (number &optional divisor)
772  "Same as ROUND, but returns first value as a float."
773  (multiple-value-bind (q r) (round number divisor)
774    (values (float q (if (floatp r) r 0.0)) r)))
775
776(defun rational (number)
777  "RATIONAL produces a rational number for any real numeric argument. This is
778  more efficient than RATIONALIZE, but it assumes that floating-point is
779  completely accurate, giving a result that isn't as pretty."
780  (if (floatp number)
781    (multiple-value-bind (s e sign)
782        (number-case number
783          (short-float
784           (integer-decode-short-float number))
785          (double-float
786           (integer-decode-double-float number)))
787      (if (eq sign -1) (setq s (- s)))
788      (if (%iminusp e)
789        (/ s (ash 1 (%i- 0 e)))
790        (ash s e)))
791    (if (rationalp number)
792      number
793      (report-bad-arg number 'real))))
794
795; make power tables for floating point reader
796(progn
797  (setq float-powers-of-5 (make-array 23))
798  (let ((array float-powers-of-5))
799    (dotimes (i 23)
800      (setf (svref array i)  (float (expt 5 i) 0.0d0))))
801  (setq integer-powers-of-5 (make-array (+ 12 (floor 324 12))))
802  (let ((array integer-powers-of-5))
803    (dotimes (i 12)
804      (setf (svref array i)  (expt 5 i)))
805    (dotimes (i (floor 324 12))
806      (setf (svref array (+ i 12)) (expt 5 (* 12 (1+ i)))))))
807
808
809(provide 'numbers)
810
Note: See TracBrowser for help on using the repository browser.