From 347cb8c2f9f2a3f9b92cad1b326b2ae0b629dea5 Mon Sep 17 00:00:00 2001 From: Christophe Rhodes Date: Wed, 14 Apr 2004 20:07:46 +0000 Subject: [PATCH] 0.8.9.41: More ANSI-fixes, this time for RATIONALIZE ... new implementation due to Bruno Haible; ... include explanation and helpful references; ... fixes RATIONALIZE.1 and RATIONALIZE.3 --- CREDITS | 7 ++- NEWS | 2 + src/code/float.lisp | 131 ++++++++++++++++++++++++++++++++++++++------------- version.lisp-expr | 2 +- 4 files changed, 108 insertions(+), 34 deletions(-) diff --git a/CREDITS b/CREDITS index 6154a06..649a232 100644 --- a/CREDITS +++ b/CREDITS @@ -553,7 +553,7 @@ Brian Downing: Miles Egan: He creates binary packages of SBCL releases for Red Hat and other - (which?) platforms + (which?) platforms. Nathan Froyd: He has fixed various bugs, and also done a lot of internal @@ -564,6 +564,11 @@ Nathan Froyd: can delete a thousand lines of implement-ITERATE macrology from the codebase.) +Bruno Haible: + He devised an accurate continued-fraction-based implementation of + RATIONALIZE, replacing a less-accurate version inherited from + primordial CMUCL. + Matthias Hoelzl: He reported and fixed COMPILE's misbehavior on macros. diff --git a/NEWS b/NEWS index 93ad34f..1706448 100644 --- a/NEWS +++ b/NEWS @@ -2385,6 +2385,8 @@ changes in sbcl-0.8.10 relative to sbcl-0.8.9: greater than 32 handle EOF correctly. * fixed some bugs revealed by Paul Dietz' test suite: ** READ-SEQUENCE now works on ECHO-STREAMs. + ** RATIONALIZE works more according to its specification. (thanks + to Bruno Haible) planned incompatible changes in 0.8.x: * (not done yet, but planned:) When the profiling interface settles diff --git a/src/code/float.lisp b/src/code/float.lisp index d5ea0b9..2d8d467 100644 --- a/src/code/float.lisp +++ b/src/code/float.lisp @@ -797,41 +797,108 @@ uninterruptibly frob the rounding modes & do ieee round-to-integer. (integer-/-integer (ash int ex) (ash 1 digits))))))) ((rational) x))) +;;; This algorithm for RATIONALIZE, due to Bruno Haible, is included +;;; with permission. +;;; +;;; 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. +;;; +;;; See also +;;; Hardy, Wright: An introduction to number theory +;;; and/or +;;; +;;; + (defun rationalize (x) - #!+sb-doc - "Converts any REAL to a RATIONAL. Floats are converted to a simple rational + "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: + 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)) - ;; Thanks to Kim Fateman, who stole this function rationalize-float from - ;; macsyma's rational. Macsyma'a rationalize was written by the legendary - ;; Gosper (rwg). Guy Steele said about Gosper, "He has been called the - ;; only living 17th century mathematician and is also the best pdp-10 - ;; hacker I know." So, if you can understand or debug this code you win - ;; big. - (cond ((minusp x) (- (rationalize (- x)))) - ((zerop x) 0) - (t - (let ((eps (etypecase x - (single-float single-float-epsilon) - (double-float double-float-epsilon) - #!+long-float - (long-float long-float-epsilon))) - (y ()) - (a ())) - (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x))))) - (num (setq a (truncate x)) - (+ (* (setq a (truncate y)) num) onum)) - (den 1 (+ (* a den) oden)) - (onum 1 num) - (oden 0 den)) - ((and (not (zerop den)) - (not (> (abs (/ (- x (/ (float num x) - (float den x))) - x)) - eps))) - (integer-/-integer num den)) - (declare ((dispatch-type x) xx))))))) + ;; 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))) diff --git a/version.lisp-expr b/version.lisp-expr index ec316cf..5cdc010 100644 --- a/version.lisp-expr +++ b/version.lisp-expr @@ -17,4 +17,4 @@ ;;; checkins which aren't released. (And occasionally for internal ;;; versions, especially for internal versions off the main CVS ;;; branch, it gets hairier, e.g. "0.pre7.14.flaky4.13".) -"0.8.9.40" +"0.8.9.41" -- 1.7.10.4