# Diff of /src/code/irrat.lisp

revision 1.36 by toy, Fri Jan 10 17:19:22 2003 UTC revision 1.37 by toy, Wed Jan 29 18:51:48 2003 UTC
# Line 383  Line 383
383           (if (minusp number)           (if (minusp number)
384               (complex (coerce (log (- number)) 'single-float)               (complex (coerce (log (- number)) 'single-float)
385                        (coerce pi 'single-float))                        (coerce pi 'single-float))
386               (coerce (/ (log2 number) #.(log (exp 1) 2d0)) 'single-float)))               (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
387          ((ratio)          ((ratio)
388           (if (minusp number)           (if (minusp number)
389               (complex (coerce (log (- number)) 'single-float)               (complex (coerce (log (- number)) 'single-float)
# Line 394  Line 394
394                     (bot (denominator number)))                     (bot (denominator number)))
395                 ;; If the number of bits in the numerator and                 ;; If the number of bits in the numerator and
396                 ;; denominator are different, just use the fact                 ;; denominator are different, just use the fact
397                 ;; log(x/y) = log(x) - log(y).  However, if they have                 ;; log(x/y) = log(x) - log(y).  But to preserve
398                 ;; the same number of bits, implying the quotient is                 ;; accuracy, we actually do
399                 ;; near one, we use log1p(x) = log(1+x).  Since the                 ;; (log2(x)-log2(y))/log2(e)).
400                 ;; number is rational, we don't lose precision                 ;;
401                 ;; subtracting 1 from it, and converting it to                 ;; However, if the numerator and denominator have the
402                 ;; double-float is accurate.                 ;; same number of bits, implying the quotient is near
403                   ;; one, we use log1p(x) = log(1+x). Since the number is
404                   ;; rational, we don't lose precision subtracting 1 from
405                   ;; it, and converting it to double-float is accurate.
406                 (if (= (integer-length top)                 (if (= (integer-length top)
407                        (integer-length bot))                        (integer-length bot))
408                     (coerce (%log1p (coerce (- number 1) 'double-float))                     (coerce (%log1p (coerce (- number 1) 'double-float))
409                             'single-float)                             'single-float)
410                     (coerce (- (log top) (log bot)) 'single-float)))))                     (coerce (/ (- (log2 top) (log2 bot))
411                                  #.(log (exp 1d0) 2d0))
412                               'single-float)))))
413          (((foreach single-float double-float))          (((foreach single-float double-float))
414           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?           ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
415           ;; 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.36 changed lines Added in v.1.37