/[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.21 by dtc, Sat Aug 30 18:21:36 1997 UTC revision 1.22 by dtc, Sat Aug 30 18:29:21 1997 UTC
# Line 201  Line 201 
201    "Returns BASE raised to the POWER."    "Returns BASE raised to the POWER."
202    (if (zerop power)    (if (zerop power)
203        (1+ (* base power))        (1+ (* base power))
204        (labels ((real-expt (base power rtype)      (labels (;; determine if the double float is an integer.
205                   (let* ((fbase (coerce base 'double-float))               ;;  0 - not an integer
206                          (fpower (coerce power 'double-float))               ;;  1 - an odd int
207                          (res (coerce (%pow fbase fpower) rtype)))               ;;  2 - an even int
208                     (if (and (zerop res) (minusp fbase))               (isint (ihi lo)
209                         (multiple-value-bind (re im)                 (declare (type (unsigned-byte 31) ihi)
210                                              (complex-pow fbase fpower)                          (type (unsigned-byte 32) lo)
211                           (%make-complex (coerce re rtype) (coerce im rtype)))                          (optimize (speed 3) (safety 0)))
212                         res)))                 (let ((isint 0))
213                 (complex-pow (fbase fpower)                   (declare (type fixnum isint))
214                   (let ((pow (%pow (- fbase) fpower))                   (cond ((>= ihi #x43400000)     ; exponent >= 53
215                         (fpower*pi (* fpower pi)))                          (setq isint 2))
216                     (values (* pow (%cos fpower*pi))                         ((>= ihi #x3ff00000)
217                             (* pow (%sin fpower*pi))))))                          (let ((k (- (ash ihi -20) #x3ff)))      ; exponent
218          (declare (inline real-expt))                            (declare (type (mod 53) k))
219          (number-dispatch ((base number) (power number))                            (cond ((> k 20)
220            (((foreach fixnum (or bignum ratio) (complex rational)) integer)                                   (let* ((shift (- 52 k))
221             (intexp base power))                                          (j (logand (ash lo (- shift))))
222            (((foreach single-float double-float) rational)                                          (j2 (ash j shift)))
223             (real-expt base power '(dispatch-type base)))                                     (declare (type (mod 32) shift)
224            (((foreach fixnum (or bignum ratio) single-float)                                              (type (unsigned-byte 32) j j2))
225              (foreach ratio single-float))                                     (when (= j2 lo)
226             (real-expt base power 'single-float))                                       (setq isint (- 2 (logand j 1))))))
227            (((foreach fixnum (or bignum ratio) single-float double-float)                                  ((= lo 0)
228              double-float)                                   (let* ((shift (- 20 k))
229             (real-expt base power 'double-float))                                          (j (ash ihi (- shift)))
230            ((double-float single-float)                                          (j2 (ash j shift)))
231             (real-expt base power 'double-float))                                     (declare (type (mod 32) shift)
232            (((foreach (complex rational) (complex float)) rational)                                              (type (unsigned-byte 31) j j2))
233             (* (expt (abs base) power)                                     (when (= j2 ihi)
234                (cis (* power (phase base)))))                                       (setq isint (- 2 (logand j 1))))))))))
235            (((foreach fixnum (or bignum ratio) single-float double-float)                   isint))
236              complex)               (real-expt (x y rtype)
237             (exp (* power (log base))))                 (let ((x (coerce x 'double-float))
238            (((foreach (complex float) (complex rational)) number)                       (y (coerce y 'double-float)))
239             (exp (* power (log base))))))))                   (declare (double-float x y))
240                     (let* ((x-hi (kernel:double-float-high-bits x))
241                            (x-lo (kernel:double-float-low-bits x))
242                            (x-ihi (logand x-hi #x7fffffff))
243                            (y-hi (kernel:double-float-high-bits y))
244                            (y-lo (kernel:double-float-low-bits y))
245                            (y-ihi (logand y-hi #x7fffffff)))
246                       (declare (type (signed-byte 32) x-hi y-hi)
247                                (type (unsigned-byte 31) x-ihi y-ihi)
248                                (type (unsigned-byte 32) x-lo y-lo))
249                       ;; y==zero: x**0 = 1
250                       (when (zerop (logior y-ihi y-lo))
251                         (return-from real-expt (coerce 1d0 rtype)))
252                       ;; +-NaN return x+y
253                       (when (or (> x-ihi #x7ff00000)
254                                 (and (= x-ihi #x7ff00000) (/= x-lo 0))
255                                 (> y-ihi #x7ff00000)
256                                 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
257                         (return-from real-expt (coerce (+ x y) rtype)))
258                       (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
259                         (declare (type fixnum yisint))
260                         ;; special value of y
261                         (when (and (zerop y-lo) (= y-ihi #x7ff00000))
262                           ;; y is +-inf
263                           (return-from real-expt
264                             (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
265                                    ;; +-1**inf is NaN
266                                    (coerce (- y y) rtype))
267                                   ((>= x-ihi #x3ff00000)
268                                    ;; (|x|>1)**+-inf = inf,0
269                                    (if (>= y-hi 0)
270                                        (coerce y rtype)
271                                        (coerce 0 rtype)))
272                                   (t
273                                    ;; (|x|<1)**-,+inf = inf,0
274                                    (if (< y-hi 0)
275                                        (coerce (- y) rtype)
276                                        (coerce 0 rtype))))))
277    
278                         (let ((abs-x (abs x)))
279                           (declare (double-float abs-x))
280                           ;; special value of x
281                           (when (and (zerop x-lo)
282                                      (or (= x-ihi #x7ff00000) (zerop x-ihi)
283                                          (= x-ihi #x3ff00000)))
284                             ;; x is +-0,+-inf,+-1
285                             (let ((z (if (< y-hi 0)
286                                          (/ 1 abs-x)       ; z = (1/|x|)
287                                          abs-x)))
288                               (declare (double-float z))
289                               (when (< x-hi 0)
290                                 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
291                                        ;; (-1)**non-int
292                                        (let ((y*pi (* y pi)))
293                                          (declare (double-float y*pi))
294                                          (return-from real-expt
295                                            (kernel::%make-complex
296                                             (coerce (%cos y*pi) rtype)
297                                             (coerce (%sin y*pi) rtype)))))
298                                       ((= yisint 1)
299                                        ;; (x<0)**odd = -(|x|**odd)
300                                        (setq z (- z)))))
301                               (return-from real-expt (coerce z rtype))))
302    
303                           (if (>= x-hi 0)
304                               ;; x>0
305                               (coerce (kernel::%pow x y) rtype)
306                               ;; x<0
307                               (let ((pow (kernel::%pow abs-x y)))
308                                 (declare (double-float pow))
309                                 (case yisint
310                                   (1 ; Odd
311                                    (coerce (* -1d0 pow) rtype))
312                                   (2 ; Even
313                                    (coerce pow rtype))
314                                   (t ; Non-integer
315                                    (let ((y*pi (* y pi)))
316                                      (declare (double-float y*pi))
317                                      (kernel::%make-complex
318                                       (coerce (* pow (%cos y*pi)) rtype)
319                                       (coerce (* pow (%sin y*pi)) rtype)))))))))))))
320          (declare (inline real-expt))
321          (number-dispatch ((base number) (power number))
322            (((foreach fixnum (or bignum ratio) (complex rational)) integer)
323             (intexp base power))
324            (((foreach single-float double-float) rational)
325             (real-expt base power '(dispatch-type base)))
326            (((foreach fixnum (or bignum ratio) single-float)
327              (foreach ratio single-float))
328             (real-expt base power 'single-float))
329            (((foreach fixnum (or bignum ratio) single-float double-float)
330              double-float)
331             (real-expt base power 'double-float))
332            ((double-float single-float)
333             (real-expt base power 'double-float))
334            (((foreach (complex rational) (complex float)) rational)
335             (* (expt (abs base) power)
336                (cis (* power (phase base)))))
337            (((foreach fixnum (or bignum ratio) single-float double-float)
338              complex)
339             (if (and (zerop base) (plusp (realpart power)))
340                 (* base power)
341                 (exp (* power (log base)))))
342            (((foreach (complex float) (complex rational))
343              (foreach complex double-float single-float))
344             (if (and (zerop base) (plusp (realpart power)))
345                 (* base power)
346                 (exp (* power (log base)))))))))
347    
348  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
349    "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."

Legend:
Removed from v.1.21  
changed lines
  Added in v.1.22

  ViewVC Help
Powered by ViewVC 1.1.5