RE: [Gcl-devel] RE: GCL bugs (decode-float, rationalize, expt)

[Stavros Macrakis]

Camm--

Here is code for rationalize which returns the simplest rational.  I
have done some testing, but of course haven't been able to test against
the new integer-decode-float (see comment in code).

      -s

(defun float-rationalize (x)
;;; Given x, returns n/d such that float(n/d)=x and d is minimal.
;;; Assuming integer-decode-float works as stated below, automatically
;;; takes into account the precision of the input float.
  (cond ((= x 0) 0)
        ((< x 0) (- (float-rationalize (- x))))
        (t (multiple-value-bind
               (mantissa exponent sign)
               (integer-decode-float x)
;;; We assume that integer-decode-float reflects the internal
;;; representation exactly.  In particular, the number of digits of
;;; precision in mantissa is exactly the number of digits of
;;; precision in the floating-point number, and not reduced to lowest
;;; terms e.g. (integer-decode-float 1.0) => 2^52 -52 1, and NOT 1 0 1.
;;; Similarly for denormalized numbers: (integer-decode-float
;;; least-positive-double-float) should be 1 -1074 1, and not 2^52
;;; -1126 1.  The CL spec is not clear about all this.
             (simplest-in-range
;;; A float represents a half-open interval, but using the closed
;;; interval is OK because the midpoint is always simpler than
;;; the endpoints.  Similarly, it never matters that the interval for
;;; exact powers of 2 (except least-positive-float) is assymmetric.
              (/ (- (* mantissa 2) 1)
                 (expt 2 (- 1 exponent)))
              (/ (+ (* mantissa 2) 1)
                 (expt 2 (- 1 exponent))) )))))

;;; From Appendix C4 of the IEEE Scheme standard
;;; Alan Bawden's implementation of Hardy and Wright's algorithm.
;;; Converted to CL by Stavros Macrakis
(defun simplest-in-range (bottom top)
  ;; 0 < bottom < top
  (multiple-value-bind (bottom-integer bottom-fraction)
      (truncate bottom)
    (if (= bottom-fraction 0)
        bottom
      (multiple-value-bind (top-integer top-fraction)
          (truncate top)
        (if (= bottom-integer top-integer)
            (+ bottom-integer
               (/ 1 (simplest-in-range
                     (/ 1 top-fraction)
                     (/ 1 bottom-fraction))))
          (+ 1 bottom-integer))))))