source: branches/working-0711/ccl/lib/numbers.lisp @ 12941

Last change on this file since 12941 was 12941, checked in by gz, 10 years ago

Merge some irrelevant (other platforms, unused) trunk changes

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