[Gcl-devel] Re: [Maxima] float to bfloat

[Camm Maguire]

Greetings!

"Stavros Macrakis" <address@hidden> writes:

> I contributed an improved rationalize function to GCL last year -- I
> don't know what its status is.  Sigh.  One of these days I'll have to
> take the plunge into building and running the latest and greatest
> releases with CVS.
> 

Thanks for resurrecting this thread!  You can see the initial
discussion at 

http://mail.gnu.org/archive/html/gcl-devel/2003-01/msg00016.html

with your contribution at

http://mail.gnu.org/archive/html/gcl-devel/2003-01/msg00029.html

or

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))))))



We were concerned about a comment in gcl_numlib.lsp:

;; although the following is correct code in that it approximates the
;; x to within eps, it does not preserve (eql (float (rationalize x) x) x)
;; since the test for eql is more strict than the float-epsilon

to which you remarked

I agree that the spec for rationalize is ambiguous, and making it behave
better is not high priority.  But the comment in numlib.lsp doesn't make
sense, because (float (rationalize x)) should be *precisely* equal to
x -- I'm afraid I can't look at the code in more detail right now.
Anyway, this is not a high priority.

and I said


4) As you already realize, I'm most happy to put in a better
   rationalize if (eq (float (rationalize x)) x) 


We left things there, I believe.  I've just written a little test for
your float-rationalize which shows that it does not preserve (eql x
(float (rationalize x))) -- perhaps you could look into this if you're
still interested?

(defun test-f (x i) (multiple-value-bind (q r s) (integer-decode-float
x) (format t "~S~%" x) (and (eql x (float (* q (expt 2 r)))) (eql
(float (float-rationalize x)) x )(if (> i 0) (test-f (float (* (+ q 1)
(expt 2 r))) (- i 1)) t))))

(test-f 1.0 10)
1.0
1.0000000000000002
NIL

(float (float-rationalize 1.0000000000000002))

1.0000000000000004

Take care,

>         -s
> 
> 
> _______________________________________________
> Maxima mailing list
> address@hidden
> http://www.math.utexas.edu/mailman/listinfo/maxima
> 
> 
> 

-- 
Camm Maguire                                            address@hidden
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah