From 3f063954c98d21ea8a95388d01db96a1e056c34d Mon Sep 17 00:00:00 2001
From: Raymond Toy <toy.raymond@gmail.com>
Date: Mon, 24 Nov 2014 13:41:00 -0800
Subject: [PATCH] Fix bug in dd-%log2 where (dd-%log2 100w0) returned 6.16
 instead of 6.64.

Tests were already added to tests/irrat.lisp.
---
 src/code/irrat-dd.lisp | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp
index 438339fa8..45b422e76 100644
--- a/src/code/irrat-dd.lisp
+++ b/src/code/irrat-dd.lisp
@@ -1243,7 +1243,17 @@ pi/4    11001001000011111101101010100010001000010110100011000 010001101001100010
 
 ;;; dd-%log2
 ;;; Base 2 logarithm.
-
+;;;
+;;; The argument is separated into its exponent and fractional
+;;; parts.  If the exponent is between -1 and +1, the (natural)
+;;; logarithm of the fraction is approximated by
+;;;
+;;;     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
+;;;
+;;; Otherwise, setting  z = 2(x-1)/x+1),
+;;; 
+;;;     log(x) = z + z**3 R(z)/S(z).
+;;;
 (let ((P (make-array 13 :element-type 'double-double-float
 		     :initial-contents
 		     '(
@@ -1326,7 +1336,7 @@ pi/4    11001001000011111101101010100010001000010110100011000 010001101001100010
 		      ;; 2*(x-1)/(x+1)
 		      (setf z (- x 0.5w0))
 		      (decf z 0.5w0)
-		      (setf y (+ (* 0.5w0 z) 0.5w0))))
+		      (setf y (+ (* 0.5w0 x) 0.5w0))))
 	       (setf x (/ z y))
 	       (setf z (* x x))
 	       (setf y (* x (/ (* z (poly-eval z r))
-- 
GitLab