X-Git-Url: http://repo.macrolet.net/gitweb/?a=blobdiff_plain;f=src%2Fcode%2Ffloat.lisp;h=2d8d467a5fedd61d7c3f8a05a87728fa636424f1;hb=f1acc17f83cdfa9454a53bd0ee9bd0e9b9482817;hp=d5ea0b9c0dee95aa7de75cbf224469767418f150;hpb=2912f5f6c2acb2da3b9fcc0f5afd1ca89782a9f8;p=sbcl.git
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)))