# Diff of /src/code/irrat.lisp

revision 1.50 by rtoy, Wed Jul 19 14:58:52 2006 UTC revision 1.51 by rtoy, Thu Jan 18 16:16:13 2007 UTC
# Line 476  Line 476
476               (* base power)               (* base power)
477               (exp (* power (log base)))))))))               (exp (* power (log base)))))))))
478
479  ;; Compute the base 2 log of an integer  ;; Log base 2 of a real number.  The result is a double-precision
480    ;; number (real or complex, as appropriate)
481  (defun log2 (x)  (defun log2 (x)
482    ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n + log2(f).    (labels ((log2-bignum (bignum)
483    ;;               ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
484    ;; So we grab the top few bits of x and scale that appropriately,               ;; + log2(f).
485    ;; take the log of it and add it to n.               ;;
486    (let ((n (integer-length x)))               ;; So we grab the top few bits of x and scale that
487      (if (< n vm:double-float-digits)               ;; appropriately, take the log of it and add it to n.
488          (log (coerce x 'double-float) 2d0)               ;;
489          (let ((exp (min vm:double-float-digits n))               ;; Return n and log2(f) separately.
490                (f (ldb (byte vm:double-float-digits               (if (minusp bignum)
491                              (max 0 (- n vm:double-float-digits)))                   (multiple-value-bind (n frac)
492                        x)))                       (log2-bignum (abs bignum))
493            (+ n (log (scale-float (float f 1d0) (- exp))                     (values n (complex frac #.(/ pi (log 2d0)))))
494                      2d0))))))                   (let ((n (integer-length bignum)))
495                       (if (< n vm:double-float-digits)
496                           (values 0 (log (coerce bignum 'double-float) 2d0))
497                           (let ((exp (min vm:double-float-digits n))
498                                 (f (ldb (byte vm:double-float-digits
499                                               (max 0 (- n vm:double-float-digits)))
500                                         bignum)))
501                             (values n (log (scale-float (float f 1d0) (- exp))
502                                            2d0))))))))
503        (etypecase x
504          (float
505           (/ (log (float x 1d0)) #.(log 2d0)))
506          (ratio
507           (let ((top (numerator x))
508                 (bot (denominator x)))
509             ;; If the number of bits in the numerator and
510             ;; denominator are different, just use the fact
511             ;; log(x/y) = log(x) - log(y).  But to preserve
512             ;; accuracy, we actually do
513             ;; (log2(x)-log2(y))/log2(e)).
514             ;;
515             ;; However, if the numerator and denominator have the
516             ;; same number of bits, implying the quotient is near
517             ;; one, we use log1p(x) = log(1+x). Since the number is
518             ;; rational, we don't lose precision subtracting 1 from
519             ;; it, and converting it to double-float is accurate.
520             (if (= (integer-length top)
521                    (integer-length bot))
522                 (/ (%log1p (coerce (- x 1) 'double-float))
523                    #.(log 2d0))
524                 (multiple-value-bind (top-n top-frac)
525                     (log2-bignum top)
526                   (multiple-value-bind (bot-n bot-frac)
527                       (log2-bignum bot)
528                     (+ (- top-n bot-n)
529                        (- top-frac bot-frac)))))))
530          (integer
531           (multiple-value-bind (n frac)
532               (log2-bignum x)
533             (+ n frac))))))
534
535    ;; Same as above, except we return double-double-float.
536    ;;
537    ;; FIXME:  Can this be merged with the above?  OAOO.
538    #+double-double
539    (defun log2-dd (x)
540      (labels ((log2-bignum (bignum)
541                 ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
542                 ;; + log2(f).
543                 ;;
544                 ;; So we grab the top few bits of x and scale that
545                 ;; appropriately, take the log of it and add it to n.
546                 ;;
547                 ;; Return n and log2(f) separately.
548                 (if (minusp bignum)
549                     (multiple-value-bind (n frac)
550                         (log2-bignum (abs bignum))
551                       (values n (complex frac #.(/ dd-pi (log 2w0)))))
552                     (let ((n (integer-length bignum)))
553                       (if (< n vm:double-double-float-digits)
554                           (values 0 (log (coerce bignum 'double-double-float) 2w0))
555                           (let ((exp (min vm:double-double-float-digits n))
556                                 (f (ldb (byte vm:double-double-float-digits
557                                               (max 0 (- n vm:double-double-float-digits)))
558                                         bignum)))
559                             (values n (log (scale-float (float f 1w0) (- exp))
560                                            2w0))))))))
561        (etypecase x
562          (float
563           (/ (log (float x 1w0)) #.(log 2w0)))
564          (ratio
565           (let ((top (numerator x))
566                 (bot (denominator x)))
567             ;; If the number of bits in the numerator and
568             ;; denominator are different, just use the fact
569             ;; log(x/y) = log(x) - log(y).  But to preserve
570             ;; accuracy, we actually do
571             ;; (log2(x)-log2(y))/log2(e)).
572             ;;
573             ;; However, if the numerator and denominator have the
574             ;; same number of bits, implying the quotient is near
575             ;; one, we use log1p(x) = log(1+x). Since the number is
576             ;; rational, we don't lose precision subtracting 1 from
577             ;; it, and converting it to double-float is accurate.
578             (if (= (integer-length top)
579                    (integer-length bot))
580                 (/ (dd-%log1p (float (- x 1) 1w0))
581                    #.(log 2w0))
582                 (multiple-value-bind (top-n top-frac)
583                     (log2-bignum top)
584                   (multiple-value-bind (bot-n bot-frac)
585                       (log2-bignum bot)
586                     (+ (- top-n bot-n)
587                        (- top-frac bot-frac)))))))
588          (integer
589           (multiple-value-bind (n frac)
590               (log2-bignum x)
591             (+ n frac))))))
592
593  (defun log (number &optional (base nil base-p))  (defun log (number &optional (base nil base-p))
594    "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 498  Line 596
596        (cond ((zerop base)        (cond ((zerop base)
597               ;; ANSI spec               ;; ANSI spec
598               base)               base)
((and (integerp number) (integerp base)
(plusp number) (plusp base))
;; Let's try to do something nice when both the number
;; and the base are positive integers.  Use the rule that
;; log_b(x) = log_2(x)/log_2(b)
(coerce (/ (log2 number) (log2 base)) 'single-float))
599              ((and (realp number) (realp base))              ((and (realp number) (realp base))
600               ;; CLHS 12.1.4.1 says               ;; CLHS 12.1.4.1 says
601               ;;               ;;
# Line 518  Line 610
610               ;; This makes (log 17 10.0) = (log 17.0 10) and so on.               ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
611               (number-dispatch ((number real) (base real))               (number-dispatch ((number real) (base real))
612                 ((double-float                 ((double-float
613                   (foreach double-float single-float fixnum bignum ratio))                   (foreach double-float single-float))
614                  (/ (log number) (log (coerce base 'double-float))))                  (/ (log2 number) (log2 base)))
615                 (((foreach single-float fixnum bignum ratio)                 (((foreach fixnum bignum ratio)
616                     (foreach fixnum bignum ratio single-float))
617                    (let* ((result (/ (log2 number) (log2 base))))
618                      ;; Figure out the right result type
619                      (if (realp result)
620                          (coerce result 'single-float)
621                          (coerce result '(complex single-float)))))
622                   (((foreach fixnum bignum ratio)
623                   double-float)                   double-float)
624                    (/ (log2 number) (log2 base)))
625                   ((single-float
626                     (foreach fixnum bignum ratio))
627                    (let* ((result (/ (log2 number) (log2 base))))
628                      ;; Figure out the right result type
629                      (if (realp result)
630                          (coerce result 'single-float)
631                          (coerce result '(complex single-float)))))
632                   ((double-float
633                     (foreach fixnum bignum ratio))
634                    (/ (log2 number) (log2 base)))
635                   ((single-float double-float)
636                  (/ (log (coerce number 'double-float)) (log base)))                  (/ (log (coerce number 'double-float)) (log base)))
637                 #+double-double                 #+double-double
638                 ((double-double-float                 ((double-double-float
639                   (foreach double-double-float double-float single-float fixnum bignum ratio))                   (foreach fixnum bignum ratio))
640                    (/ (log2-dd number) (log2-dd base)))
641                   #+double-double
642                   ((double-double-float
643                     (foreach double-double-float double-float single-float))
644                  (/ (log number) (log (coerce base 'double-double-float))))                  (/ (log number) (log (coerce base 'double-double-float))))
645                 #+double-double                 #+double-double
646                 (((foreach double-float single-float fixnum bignum ratio)                 (((foreach fixnum bignum ratio)
647                     double-double-float)
648                    (/ (log2-dd number) (log2-dd base)))
649                   #+double-double
650                   (((foreach double-float single-float)
651                   double-double-float)                   double-double-float)
652                  (/ (log (coerce number 'double-double-float)) (log base)))                  (/ (log (coerce number 'double-double-float)) (log base)))
653                 (((foreach single-float fixnum bignum ratio)                 (((foreach single-float)
654                   (foreach single-float fixnum bignum ratio))                   (foreach single-float))
655                  ;; Converting everything to double-float helps the                  ;; Converting everything to double-float helps the
656                  ;; cases like (log 17 10) = (/ (log 17) (log 10)).                  ;; cases like (log 17 10) = (/ (log 17) (log 10)).
657                  ;; This is usually handled above, but if we compute (/                  ;; This is usually handled above, but if we compute (/

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