source: release/1.9/source/lib/numbers.lisp @ 15685

Last change on this file since 15685 was 15685, checked in by gb, 6 years ago

Propagate r15683 to trunk.

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