Diff of /src/code/irrat.lisp

revision 1.51 by rtoy, Thu Jan 18 16:16:13 2007 UTC revision 1.52 by rtoy, Thu Jan 18 17:36:22 2007 UTC
# Line 476  Line 476
476               (* base power)               (* base power)
477               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))
478
479  ;; Log base 2 of a real number.  The result is a double-precision  ;; Log base 2 of a real number.  The result is a either a double-float
480  ;; number (real or complex, as appropriate)  ;; or double-double-float number (real or complex, as appropriate),
481  (defun log2 (x)  ;; depending on the type of FLOAT-TYPE.
482    (labels ((log2-bignum (bignum)  (defun log2 (x &optional (float-type 1d0))
483               ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n    (labels ((log-of-2 (f)
484               ;; + log2(f).               ;; log(2), with the precision specified by the type of F
485               ;;               (number-dispatch ((f real))
486               ;; So we grab the top few bits of x and scale that                 ((double-float)
;; appropriately, take the log of it and add it to n.
;;
;; Return n and log2(f) separately.
(if (minusp bignum)
(multiple-value-bind (n frac)
(log2-bignum (abs bignum))
(values n (complex frac #.(/ pi (log 2d0)))))
(let ((n (integer-length bignum)))
(if (< n vm:double-float-digits)
(values 0 (log (coerce bignum 'double-float) 2d0))
(let ((exp (min vm:double-float-digits n))
(f (ldb (byte vm:double-float-digits
(max 0 (- n vm:double-float-digits)))
bignum)))
(values n (log (scale-float (float f 1d0) (- exp))
2d0))))))))
(etypecase x
(float
(/ (log (float x 1d0)) #.(log 2d0)))
(ratio
(let ((top (numerator x))
(bot (denominator x)))
;; If the number of bits in the numerator and
;; denominator are different, just use the fact
;; log(x/y) = log(x) - log(y).  But to preserve
;; accuracy, we actually do
;; (log2(x)-log2(y))/log2(e)).
;;
;; However, if the numerator and denominator have the
;; same number of bits, implying the quotient is near
;; one, we use log1p(x) = log(1+x). Since the number is
;; rational, we don't lose precision subtracting 1 from
;; it, and converting it to double-float is accurate.
(if (= (integer-length top)
(integer-length bot))
(/ (%log1p (coerce (- x 1) 'double-float))
487                  #.(log 2d0))                  #.(log 2d0))
488               (multiple-value-bind (top-n top-frac)                 #+double-double
489                   (log2-bignum top)                 ((double-double-float)
490                 (multiple-value-bind (bot-n bot-frac)                  #.(log 2w0))))
491                     (log2-bignum bot)             (log-2-pi (f)
492                   (+ (- top-n bot-n)               ;; log(pi), with the precision specified by the type of F
493                      (- top-frac bot-frac)))))))               (number-dispatch ((f real))
494        (integer                 ((double-float)
495         (multiple-value-bind (n frac)                  #.(/ pi (log 2d0)))
496             (log2-bignum x)                 #+double-double
497           (+ n frac))))))                 ((double-double-float)
498                    #.(/ dd-pi (log 2w0)))))
499  ;; Same as above, except we return double-double-float.             (log1p (x)
500  ;;               ;; log(1+x), with the precision specified by the type of
501  ;; FIXME:  Can this be merged with the above?  OAOO.               ;; X
502  #+double-double               (number-dispatch ((x real))
503  (defun log2-dd (x)                 (((foreach single-float double-float))
504    (labels ((log2-bignum (bignum)                  (%log1p (float x 1d0)))
505                   #+double-double
506                   ((double-double-float)
507                    (dd-%log1p x))))
508               (log2-bignum (bignum)
509               ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n               ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
510               ;; + log2(f).               ;; + log2(f).
511               ;;               ;;
# Line 548  Line 516
516               (if (minusp bignum)               (if (minusp bignum)
517                   (multiple-value-bind (n frac)                   (multiple-value-bind (n frac)
518                       (log2-bignum (abs bignum))                       (log2-bignum (abs bignum))
519                     (values n (complex frac #.(/ dd-pi (log 2w0)))))                     (values n (complex frac (log-2-pi float-type))))
520                   (let ((n (integer-length bignum)))                   (let ((n (integer-length bignum))
521                     (if (< n vm:double-double-float-digits)                         (float-bits (float-digits float-type)))
522                         (values 0 (log (coerce bignum 'double-double-float) 2w0))                     (if (< n float-bits)
523                         (let ((exp (min vm:double-double-float-digits n))                         (values 0 (log (float bignum float-type)
524                               (f (ldb (byte vm:double-double-float-digits                                        (float 2 float-type)))
525                                             (max 0 (- n vm:double-double-float-digits)))                         (let ((exp (min float-bits n))
526                                 (f (ldb (byte float-bits
527                                               (max 0 (- n float-bits)))
528                                       bignum)))                                       bignum)))
529                           (values n (log (scale-float (float f 1w0) (- exp))                           (values n (log (scale-float (float f float-type) (- exp))
530                                          2w0))))))))                                          (float 2 float-type)))))))))
531      (etypecase x      (etypecase x
532        (float        (float
533         (/ (log (float x 1w0)) #.(log 2w0)))         (/ (log (float x float-type)) (log-of-2 float-type)))
534        (ratio        (ratio
535         (let ((top (numerator x))         (let ((top (numerator x))
536               (bot (denominator x)))               (bot (denominator x)))
# Line 577  Line 547
547           ;; it, and converting it to double-float is accurate.           ;; it, and converting it to double-float is accurate.
548           (if (= (integer-length top)           (if (= (integer-length top)
549                  (integer-length bot))                  (integer-length bot))
550               (/ (dd-%log1p (float (- x 1) 1w0))               (/ (log1p (float (- x 1) float-type))
551                  #.(log 2w0))                  (log-of-2 float-type))
552               (multiple-value-bind (top-n top-frac)               (multiple-value-bind (top-n top-frac)
553                   (log2-bignum top)                   (log2-bignum top)
554                 (multiple-value-bind (bot-n bot-frac)                 (multiple-value-bind (bot-n bot-frac)
# Line 637  Line 607
607                 #+double-double                 #+double-double
608                 ((double-double-float                 ((double-double-float
609                   (foreach fixnum bignum ratio))                   (foreach fixnum bignum ratio))
610                  (/ (log2-dd number) (log2-dd base)))                  (/ (log2 number 1w0) (log2 base 1w0)))
611                 #+double-double                 #+double-double
612                 ((double-double-float                 ((double-double-float
613                   (foreach double-double-float double-float single-float))                   (foreach double-double-float double-float single-float))
# Line 645  Line 615
615                 #+double-double                 #+double-double
616                 (((foreach fixnum bignum ratio)                 (((foreach fixnum bignum ratio)
617                   double-double-float)                   double-double-float)
618                  (/ (log2-dd number) (log2-dd base)))                  (/ (log2 number 1w0) (log2 base 1w0)))
619                 #+double-double                 #+double-double
620                 (((foreach double-float single-float)                 (((foreach double-float single-float)
621                   double-double-float)                   double-double-float)

Legend:
 Removed from v.1.51 changed lines Added in v.1.52