# Diff of /src/code/irrat.lisp

revision 1.33 by toy, Thu Sep 5 16:13:46 2002 UTC revision 1.34 by toy, Tue Dec 31 14:28:08 2002 UTC
# Line 348  Line 348
348               (* base power)               (* base power)
349               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))
350
351    ;; Compute the base 2 log of an integer
352    (defun log2 (x)
353      ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n + log2(f).
354      ;;
355      ;; So we grab the top few bits of x and scale that appropriately,
356      ;; take the log of it and add it to n.
357      (let ((n (integer-length x)))
358        (if (< n vm:double-float-digits)
359            (log (coerce x 'double-float) 2d0)
360            (let ((exp (min vm:double-float-digits n))
361                  (f (ldb (byte vm:double-float-digits
362                                (max 0 (- n vm:double-float-digits)))
363                          x)))
364              (+ n (log (scale-float (float f 1d0) (- exp))
365                        2d0))))))
366
367  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
368    "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."
369    (if base-p    (if base-p
370        (if (zerop base)        (cond ((zerop base)
371            base                          ; ANSI spec               ;; ANSI spec
372            (/ (log number) (log base)))               base)
373                ((and (integerp number) (integerp base)
374                      (plusp number) (plusp base))
375                 ;; Let's try to do something nice when both the number
376                 ;; and the base are positive integers.  Use the rule that
377                 ;; log_b(x) = log_2(x)/log_2(b)
378                 (coerce (/ (log2 number) (log2 base)) 'single-float))
379                (t
380                 (/ (log number) (log base))))
381        (number-dispatch ((number number))        (number-dispatch ((number number))
382          (((foreach fixnum bignum ratio))          (((foreach fixnum bignum))
383           (if (minusp number)           (if (minusp number)
384               (complex (log (- number)) (coerce pi 'single-float))               (complex (coerce (log (- number)) 'single-float)
385               (coerce (%log (coerce number 'double-float)) 'single-float)))                        (coerce pi 'single-float))
386                 (coerce (/ (log2 number) #.(log (exp 1) 2d0)) 'single-float)))
387            ((ratio)
388             (if (minusp number)
389                 (complex (coerce (log (- number)) 'single-float)
390                          (coerce pi 'single-float))
391                 (coerce (- (log (numerator number))
392                            (log (denominator number)))
393                         'single-float)))
394          (((foreach single-float double-float))          (((foreach single-float double-float))
395           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
396           ;; Since this doesn't seem to be an implementation issue           ;; Since this doesn't seem to be an implementation issue

Legend:
 Removed from v.1.33 changed lines Added in v.1.34