0.7.12.4:

	Fix bug in LOG for rational numbers near 1.
	... also do neato things for e.g. (log <bignum> 3), so that
		answers can be computed even though <bignum> is bigger
		than MOST-POSITIVE-DOUBLE-FLOAT.
	... regression tests will be along shortly, I hope

diff --git a/NEWS b/NEWS
index c4182e0..408b56f 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1502,6 +1502,9 @@ changes in sbcl-0.7.12 relative to sbcl-0.7.11:
     ** CONSTANTP now returns true for all self-evaluating objects.
 
 changes in sbcl-0.7.13 relative to sbcl-0.7.12:
+  * fixed a bug in LOG, so that LOG of a rational argument near 1 now
+    gives a closer approximation to the right answer than previously.
+    (thanks to Raymond Toy)
   * fixed some bugs revealed by Paul Dietz' test suite:
     ** ARRAY-IN-BOUNDS-P now allows arbitrary integers as arguments,
        not just nonnegative fixnums;

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 9acba10..b221854 100644 (file)
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -269,18 +269,53 @@
             (* base power)
             (exp (* power (log base)))))))))
 
+;;; FIXME: Maybe rename this so that it's clearer that it only works
+;;; on integers?
+(defun log2 (x)
+  (declare (type integer x))
+  ;; CMUCL comment:
+  ;;
+  ;;   Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n +
+  ;;   log2(f).  So we grab the top few bits of x and scale that
+  ;;   appropriately, take the log of it and add it to n.
+  ;;
+  ;; Motivated by an attempt to get LOG to work better on bignums.
+  (let ((n (integer-length x)))
+    (if (< n sb!vm:double-float-digits)
+       (log (coerce x 'double-float) 2.0d0)
+       (let ((f (ldb (byte sb!vm:double-float-digits
+                           (- n sb!vm:double-float-digits))
+                     x)))
+         (+ n (log (scale-float (coerce f 'double-float)
+                                (- sb!vm:double-float-digits))
+                   2.0d0))))))
+
 (defun log (number &optional (base nil base-p))
   #!+sb-doc
   "Return the logarithm of NUMBER in the base BASE, which defaults to e."
   (if base-p
-      (if (zerop base)
-         base                          ; ANSI spec
-         (/ (log number) (log base)))
+      (cond
+       ((zerop base) base) ; ANSI spec
+       ((and (typep number '(integer (0) *))
+             (typep base '(integer (0) *)))
+        (coerce (/ (log2 number) (log2 base)) 'single-float))
+       (t (/ (log number) (log base))))
       (number-dispatch ((number number))
-       (((foreach fixnum bignum ratio))
+       (((foreach fixnum bignum))
+        (if (minusp number)
+            (complex (log (- number)) (coerce pi 'single-float))
+            (coerce (/ (log2 number) (log (exp 1.0d0) 2.0d0)) 'single-float)))
+       ((ratio)
         (if (minusp number)
             (complex (log (- number)) (coerce pi 'single-float))
-            (coerce (%log (coerce number 'double-float)) 'single-float)))
+            (let ((numerator (numerator number))
+                  (denominator (denominator number)))
+              (if (= (integer-length numerator)
+                     (integer-length denominator))
+                  (coerce (%log1p (coerce (- number 1) 'double-float))
+                          'single-float)
+                  (coerce (- (log numerator) (log denominator))
+                             'single-float)))))
        (((foreach single-float double-float))
         ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
         ;; Since this doesn't seem to be an implementation issue

diff --git a/version.lisp-expr b/version.lisp-expr
index 4c4af7d..ec7e398 100644 (file)
--- a/version.lisp-expr
+++ b/version.lisp-expr
@@ -18,4 +18,4 @@
 ;;; versions, especially for internal versions off the main CVS
 ;;; branch, it gets hairier, e.g. "0.pre7.14.flaky4.13".)
 
-"0.7.12.3"
+"0.7.12.4"