source: trunk/source/lib/numbers.lisp

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

Update copyright/license headers in files.

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