/[cmucl]/src/code/irrat.lisp
ViewVC logotype

Diff of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.19.2.2 by dtc, Mon Sep 8 00:22:00 1997 UTC revision 1.19.2.3 by pw, Tue Jun 23 11:22:02 1998 UTC
# Line 77  Line 77 
77  (def-math-rtn "pow" 2)  (def-math-rtn "pow" 2)
78  #-x86 (def-math-rtn "sqrt" 1)  #-x86 (def-math-rtn "sqrt" 1)
79  (def-math-rtn "hypot" 2)  (def-math-rtn "hypot" 2)
80  #-hpux  #-(or hpux x86) (def-math-rtn "log1p" 1)
 (def-math-rtn "log1p" 1)  
81    
82  #+x86 ;; These are needed for use by byte-compiled files.  #+x86 ;; These are needed for use by byte-compiled files.
83  (progn  (progn
# Line 149  Line 148 
148      (declare (double-float x)      (declare (double-float x)
149               (values double-float))               (values double-float))
150      (%logb x))      (%logb x))
151      (defun %log1p (x)
152        (declare (double-float x)
153                 (values double-float))
154        (%log1p x))
155    ) ; progn    ) ; progn
156    
157    
# Line 201  Line 204 
204    "Returns BASE raised to the POWER."    "Returns BASE raised to the POWER."
205    (if (zerop power)    (if (zerop power)
206        (1+ (* base power))        (1+ (* base power))
207        (labels ((real-expt (base power rtype)      (labels (;; determine if the double float is an integer.
208                   (let* ((fbase (coerce base 'double-float))               ;;  0 - not an integer
209                          (fpower (coerce power 'double-float))               ;;  1 - an odd int
210                          (res (coerce (%pow fbase fpower) rtype)))               ;;  2 - an even int
211                     (if (and (zerop res) (minusp fbase))               (isint (ihi lo)
212                         (multiple-value-bind (re im)                 (declare (type (unsigned-byte 31) ihi)
213                                              (complex-pow fbase fpower)                          (type (unsigned-byte 32) lo)
214                           (%make-complex (coerce re rtype) (coerce im rtype)))                          (optimize (speed 3) (safety 0)))
215                         res)))                 (let ((isint 0))
216                 (complex-pow (fbase fpower)                   (declare (type fixnum isint))
217                   (let ((pow (%pow (- fbase) fpower))                   (cond ((>= ihi #x43400000)     ; exponent >= 53
218                         (fpower*pi (* fpower pi)))                          (setq isint 2))
219                     (values (* pow (%cos fpower*pi))                         ((>= ihi #x3ff00000)
220                             (* pow (%sin fpower*pi))))))                          (let ((k (- (ash ihi -20) #x3ff)))      ; exponent
221          (declare (inline real-expt))                            (declare (type (mod 53) k))
222          (number-dispatch ((base number) (power number))                            (cond ((> k 20)
223            (((foreach fixnum (or bignum ratio) (complex rational)) integer)                                   (let* ((shift (- 52 k))
224             (intexp base power))                                          (j (logand (ash lo (- shift))))
225            (((foreach single-float double-float) rational)                                          (j2 (ash j shift)))
226             (real-expt base power '(dispatch-type base)))                                     (declare (type (mod 32) shift)
227            (((foreach fixnum (or bignum ratio) single-float)                                              (type (unsigned-byte 32) j j2))
228              (foreach ratio single-float))                                     (when (= j2 lo)
229             (real-expt base power 'single-float))                                       (setq isint (- 2 (logand j 1))))))
230            (((foreach fixnum (or bignum ratio) single-float double-float)                                  ((= lo 0)
231              double-float)                                   (let* ((shift (- 20 k))
232             (real-expt base power 'double-float))                                          (j (ash ihi (- shift)))
233            ((double-float single-float)                                          (j2 (ash j shift)))
234             (real-expt base power 'double-float))                                     (declare (type (mod 32) shift)
235            (((foreach (complex rational) (complex float)) rational)                                              (type (unsigned-byte 31) j j2))
236             (* (expt (abs base) power)                                     (when (= j2 ihi)
237                (cis (* power (phase base)))))                                       (setq isint (- 2 (logand j 1))))))))))
238            (((foreach fixnum (or bignum ratio) single-float double-float)                   isint))
239              complex)               (real-expt (x y rtype)
240             (exp (* power (log base))))                 (let ((x (coerce x 'double-float))
241            (((foreach (complex float) (complex rational)) number)                       (y (coerce y 'double-float)))
242             (exp (* power (log base))))))))                   (declare (double-float x y))
243                     (let* ((x-hi (kernel:double-float-high-bits x))
244                            (x-lo (kernel:double-float-low-bits x))
245                            (x-ihi (logand x-hi #x7fffffff))
246                            (y-hi (kernel:double-float-high-bits y))
247                            (y-lo (kernel:double-float-low-bits y))
248                            (y-ihi (logand y-hi #x7fffffff)))
249                       (declare (type (signed-byte 32) x-hi y-hi)
250                                (type (unsigned-byte 31) x-ihi y-ihi)
251                                (type (unsigned-byte 32) x-lo y-lo))
252                       ;; y==zero: x**0 = 1
253                       (when (zerop (logior y-ihi y-lo))
254                         (return-from real-expt (coerce 1d0 rtype)))
255                       ;; +-NaN return x+y
256                       (when (or (> x-ihi #x7ff00000)
257                                 (and (= x-ihi #x7ff00000) (/= x-lo 0))
258                                 (> y-ihi #x7ff00000)
259                                 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
260                         (return-from real-expt (coerce (+ x y) rtype)))
261                       (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
262                         (declare (type fixnum yisint))
263                         ;; special value of y
264                         (when (and (zerop y-lo) (= y-ihi #x7ff00000))
265                           ;; y is +-inf
266                           (return-from real-expt
267                             (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
268                                    ;; +-1**inf is NaN
269                                    (coerce (- y y) rtype))
270                                   ((>= x-ihi #x3ff00000)
271                                    ;; (|x|>1)**+-inf = inf,0
272                                    (if (>= y-hi 0)
273                                        (coerce y rtype)
274                                        (coerce 0 rtype)))
275                                   (t
276                                    ;; (|x|<1)**-,+inf = inf,0
277                                    (if (< y-hi 0)
278                                        (coerce (- y) rtype)
279                                        (coerce 0 rtype))))))
280    
281                         (let ((abs-x (abs x)))
282                           (declare (double-float abs-x))
283                           ;; special value of x
284                           (when (and (zerop x-lo)
285                                      (or (= x-ihi #x7ff00000) (zerop x-ihi)
286                                          (= x-ihi #x3ff00000)))
287                             ;; x is +-0,+-inf,+-1
288                             (let ((z (if (< y-hi 0)
289                                          (/ 1 abs-x)       ; z = (1/|x|)
290                                          abs-x)))
291                               (declare (double-float z))
292                               (when (< x-hi 0)
293                                 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
294                                        ;; (-1)**non-int
295                                        (let ((y*pi (* y pi)))
296                                          (declare (double-float y*pi))
297                                          (return-from real-expt
298                                            (complex
299                                             (coerce (%cos y*pi) rtype)
300                                             (coerce (%sin y*pi) rtype)))))
301                                       ((= yisint 1)
302                                        ;; (x<0)**odd = -(|x|**odd)
303                                        (setq z (- z)))))
304                               (return-from real-expt (coerce z rtype))))
305    
306                           (if (>= x-hi 0)
307                               ;; x>0
308                               (coerce (kernel::%pow x y) rtype)
309                               ;; x<0
310                               (let ((pow (kernel::%pow abs-x y)))
311                                 (declare (double-float pow))
312                                 (case yisint
313                                   (1 ; Odd
314                                    (coerce (* -1d0 pow) rtype))
315                                   (2 ; Even
316                                    (coerce pow rtype))
317                                   (t ; Non-integer
318                                    (let ((y*pi (* y pi)))
319                                      (declare (double-float y*pi))
320                                      (complex
321                                       (coerce (* pow (%cos y*pi)) rtype)
322                                       (coerce (* pow (%sin y*pi)) rtype)))))))))))))
323          (declare (inline real-expt))
324          (number-dispatch ((base number) (power number))
325            (((foreach fixnum (or bignum ratio) (complex rational)) integer)
326             (intexp base power))
327            (((foreach single-float double-float) rational)
328             (real-expt base power '(dispatch-type base)))
329            (((foreach fixnum (or bignum ratio) single-float)
330              (foreach ratio single-float))
331             (real-expt base power 'single-float))
332            (((foreach fixnum (or bignum ratio) single-float double-float)
333              double-float)
334             (real-expt base power 'double-float))
335            ((double-float single-float)
336             (real-expt base power 'double-float))
337            (((foreach (complex rational) (complex float)) rational)
338             (* (expt (abs base) power)
339                (cis (* power (phase base)))))
340            (((foreach fixnum (or bignum ratio) single-float double-float)
341              complex)
342             (if (and (zerop base) (plusp (realpart power)))
343                 (* base power)
344                 (exp (* power (log base)))))
345            (((foreach (complex float) (complex rational))
346              (foreach complex double-float single-float))
347             (if (and (zerop base) (plusp (realpart power)))
348                 (* base power)
349                 (exp (* power (log base)))))))))
350    
351  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
352    "Return the logarithm of NUMBER in the base BASE, which defaults to e."    "Return the logarithm of NUMBER in the base BASE, which defaults to e."
# Line 269  Line 379 
379           (complex-sqrt number)           (complex-sqrt number)
380           (coerce (%sqrt (coerce number 'double-float)) 'single-float)))           (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
381      (((foreach single-float double-float))      (((foreach single-float double-float))
382       ;; NOTE there is a problem with (at least x86 NPX) of what result       (if (minusp number)
      ;; should be returned for (sqrt -0.0). The x86 hardware FSQRT  
      ;; instruction returns -0d0. The result is that Python will perhaps  
      ;; note the following test in generic sqrt, non-negatively constrained  
      ;; float types will be passed to FSQRT (or libm on other boxes).  
      ;; So, in the interest of consistency of compiled and interpreted  
      ;; codes, the following test is disabled for now. Maybe the float-sign  
      ;; test could be moved to the optimization codes.  
      (if (< (#+nil float-sign #-nil identity number)  
             (coerce 0 '(dispatch-type number)))  
383           (complex-sqrt number)           (complex-sqrt number)
384           (coerce (%sqrt (coerce number 'double-float))           (coerce (%sqrt (coerce number 'double-float))
385                   '(dispatch-type number))))                   '(dispatch-type number))))
# Line 651  For the special cases, the following val Line 752  For the special cases, the following val
752     0              -infinity     0              -infinity
753  "  "
754    (declare (type double-float x))    (declare (type double-float x))
755    (cond ((float-trapping-nan-p x)    (cond ((float-nan-p x)
756           x)           x)
757          ((float-infinity-p x)          ((float-infinity-p x)
758           #.ext:double-float-positive-infinity)           #.ext:double-float-positive-infinity)
# Line 703  and Y are coerced to single-float." Line 804  and Y are coerced to single-float."
804      ;; signal would have been raised.      ;; signal would have been raised.
805      (with-float-traps-masked (:underflow :overflow)      (with-float-traps-masked (:underflow :overflow)
806        (setf rho (+ (square x) (square y)))        (setf rho (+ (square x) (square y)))
807        (cond ((and (or (float-trapping-nan-p rho)        (cond ((and (or (float-nan-p rho)
808                        (float-infinity-p rho))                        (float-infinity-p rho))
809                    (or (float-infinity-p (abs x))                    (or (float-infinity-p (abs x))
810                        (float-infinity-p (abs y))))                        (float-infinity-p (abs y))))
# Line 737  Z may be any number, but the result is a Line 838  Z may be any number, but the result is a
838            (nu 0d0))            (nu 0d0))
839        (declare (double-float x y eta nu))        (declare (double-float x y eta nu))
840    
841        (if (not (float-trapping-nan-p x))        (if (not (float-nan-p x))
842            (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))            (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
843    
844        (cond ((oddp k)        (cond ((oddp k)

Legend:
Removed from v.1.19.2.2  
changed lines
  Added in v.1.19.2.3

  ViewVC Help
Powered by ViewVC 1.1.5