From e4730b128c3bfdcfa8e039922e0b2699a32c58e2 Mon Sep 17 00:00:00 2001
From: dtc <dtc>
Date: Mon, 22 May 2000 17:43:09 +0000
Subject: [PATCH] New implementation of the function rationalize, from Raymond
 Toy based on an algorithm from Bruno Haible.

---
 code/float.lisp | 107 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 106 insertions(+), 1 deletion(-)

diff --git a/code/float.lisp b/code/float.lisp
index 7d6e9bba6..3d1ca41ae 100644
--- a/code/float.lisp
+++ b/code/float.lisp
@@ -5,7 +5,7 @@
 ;;; Carnegie Mellon University, and has been placed in the public domain.
 ;;;
 (ext:file-comment
-  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.21 2000/05/06 17:21:16 dtc Exp $")
+  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/float.lisp,v 1.22 2000/05/22 17:43:09 dtc Exp $")
 ;;;
 ;;; **********************************************************************
 ;;;
@@ -1082,6 +1082,7 @@ rounding modes & do ieee round-to-integer.
     ((rational) x)))
 
 
+#+nil
 (defun rationalize (x)
   "Converts any REAL to a RATIONAL.  Floats are converted to a simple rational
   representation exploiting the assumption that floats are only accurate to
@@ -1121,3 +1122,107 @@ rounding modes & do ieee round-to-integer.
 		   (integer-/-integer num den))
 		(declare ((dispatch-type x) xx)))))))
     ((rational) x)))
+
+;;; RATIONALIZE  --  Public
+;;;
+;;; The algorithm here is the method described in CLISP.  Bruno Haible has
+;;; graciously given permission to use this algorithm.  He says, "You can use
+;;; it, if you present the following explanation of the algorithm."
+;;;
+;;; Algorithm (recursively presented):
+;;;   If x is a rational number, return x.
+;;;   If x = 0.0, return 0.
+;;;   If x < 0.0, return (- (rationalize (- x))).
+;;;   If x > 0.0:
+;;;     Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
+;;;     exponent, sign).
+;;;     If m = 0 or e >= 0: return x = m*2^e.
+;;;     Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
+;;;     with smallest possible numerator and denominator.
+;;;     Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
+;;;       But in this case the result will be x itself anyway, regardless of
+;;;       the choice of a. Therefore we can simply ignore this case.
+;;;     Note 2: At first, we need to consider the closed interval [a,b].
+;;;       but since a and b have the denominator 2^(|e|+1) whereas x itself
+;;;       has a denominator <= 2^|e|, we can restrict the seach to the open
+;;;       interval (a,b).
+;;;     So, for given a and b (0 < a < b) we are searching a rational number
+;;;     y with a <= y <= b.
+;;;     Recursive algorithm fraction_between(a,b):
+;;;       c := (ceiling a)
+;;;       if c < b
+;;;         then return c       ; because a <= c < b, c integer
+;;;         else
+;;;           ; a is not integer (otherwise we would have had c = a < b)
+;;;           k := c-1          ; k = floor(a), k < a < b <= k+1
+;;;           return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
+;;;                             ; note 1 <= 1/(b-k) < 1/(a-k)
+;;;
+;;; You can see that we are actually computing a continued fraction expansion.
+;;;
+;;; Algorithm (iterative):
+;;;   If x is rational, return x.
+;;;   Call (integer-decode-float x). It returns a m,e,s (mantissa,
+;;;     exponent, sign).
+;;;   If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
+;;;   Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
+;;;   (positive and already in lowest terms because the denominator is a
+;;;   power of two and the numerator is odd).
+;;;   Start a continued fraction expansion
+;;;     p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
+;;;   Loop
+;;;     c := (ceiling a)
+;;;     if c >= b
+;;;       then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
+;;;            goto Loop
+;;;   finally partial_quotient(c).
+;;;   Here partial_quotient(c) denotes the iteration
+;;;     i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
+;;;   At the end, return s * (p[i]/q[i]).
+;;;   This rational number is already in lowest terms because
+;;;   p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
+;;;
+(defun rationalize (x)
+  "Converts any REAL to a RATIONAL.  Floats are converted to a simple rational
+  representation exploiting the assumption that floats are only accurate to
+  their precision.  RATIONALIZE (and also RATIONAL) preserve the invariant:
+      (= x (float (rationalize x) x))"
+  (number-dispatch ((x real))
+    (((foreach single-float double-float #+long-float long-float))
+     ;; This is a fairly straigtforward implementation of the iterative
+     ;; algorithm above.
+     (multiple-value-bind (frac expo sign)
+	 (integer-decode-float x)
+       (cond ((or (zerop frac) (>= expo 0))
+	      (if (minusp sign)
+		  (- (ash frac expo))
+		  (ash frac expo)))
+	     (t
+	      ;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e),
+	      ;; so build the fraction up immediately, without having to do
+	      ;; a gcd.
+	      (let ((a (build-ratio (- (* 2 frac) 1) (ash 1 (- 1 expo))))
+		    (b (build-ratio (+ (* 2 frac) 1) (ash 1 (- 1 expo))))
+		    (p0 0)
+		    (q0 1)
+		    (p1 1)
+		    (q1 0))
+		(do ((c (ceiling a) (ceiling a)))
+		    ((< c b)
+		     (let ((top (+ (* c p1) p0))
+			   (bot (+ (* c q1) q0)))
+		       (build-ratio (if (minusp sign)
+					(- top)
+					top)
+				    bot)))
+		  (let* ((k (- c 1))
+			 (p2 (+ (* k p1) p0))
+			 (q2 (+ (* k q1) q0)))
+		    (psetf a (/ (- b k))
+			   b (/ (- a k)))
+		    (setf p0 p1
+			  q0 q1
+			  p1 p2
+			  q1 q2))))))))
+    ((rational) x)))
+
-- 
GitLab