From 2e182540e895f16b26224ba11b2f38de93ab2672 Mon Sep 17 00:00:00 2001
From: Raymond Toy <toy.raymond@gmail.com>
Date: Sat, 20 Sep 2014 11:55:27 -0700
Subject: [PATCH] Make (log (scale-float 1d0 k) 2) = k

 * src/code/irrat.lisp
   * Improve accuracy of log2 for the case when x=2^k. There might be
     some loss in accuracy for other values of x, however.
 * tests/float.lisp
   * Add test for log2(2^k) = k.
---
 src/code/irrat.lisp | 12 +++++++++++-
 tests/float.lisp    |  5 +++++
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/code/irrat.lisp b/src/code/irrat.lisp
index 5f24bf9a4..2591907ae 100644
--- a/src/code/irrat.lisp
+++ b/src/code/irrat.lisp
@@ -598,7 +598,17 @@
 					(float 2 float-type)))))))))
     (etypecase x
       (float
-       (/ (log (float x float-type)) (log-of-2 float-type)))
+       (multiple-value-bind (f e s)
+	   (decode-float x)
+	 ;; If x = 2^e*f then log2(x) = e + log2(f) = e +
+	 ;; log(f)/log(2). Accuracy could be improved some since this
+	 ;; has 2 rounding operations instead of just 1 for
+	 ;; log(x)/log(2).
+	 (let ((log2 (+ e (/ (log (float f float-type))
+			     (log-of-2 float-type)))))
+	   (if (minusp s)
+	       (complex log2 (log-2-pi float-type))
+	       log2))))
       (ratio
        (let ((top (numerator x))
 	     (bot (denominator x)))
diff --git a/tests/float.lisp b/tests/float.lisp
index 6fe1f50f7..ff07270a2 100644
--- a/tests/float.lisp
+++ b/tests/float.lisp
@@ -10,3 +10,8 @@
 					 (declare (type (double-float (0d0)) x))
 					 (decode-float x)))
 			1d0)))
+
+(define-test log2
+  (loop for k from -1074 to 1023 do
+    (let ((x (scale-float 1d0 k)))
+      (assert-equal k (log x 2)))))
-- 
GitLab