Re: floating point

Below is a concrete example of a useful procedure whose correctness
depends upon the predictability of IEEE floating point arithmetic. It
probably won't work, and may not even terminate, if the compiler uses
extended precision for some operations but double precision for
others.


I wrote this code over a decade ago, to counter pro-flaky-FP sentiment
within the Scheme community. You can translate it into your favorite
language quite easily. You don't even have to remove the tail
recursion, because its depth is bounded by the floating point
precision + maximum binary exponent.


The basic problem is that the IEEE standard was conceived as a
standard for hardware, and says scarcely a word about high level
languages or compilers. IEEE Std 754-1985 goes to great lengths to
ensure that IEEE floating point arithmetic is predictable at the
hardware level, but Kahan himself has urged compiler writers to use
extended precision for intermediate results, without seeming to
appreciate how this leads to unpredictable floating point arithmetic
at the language level, where almost all programmers live. See Kahan's
1997 John von Neumann lecture on The Baleful Effect of Computer
Languages and Benchmarks upon Applied Mathematics, Physics, and
Chemistry, available at
http://http.cs.berkeley.edu/~wkahan/SIAMjvnl.ps.


Compare Kahan's attitude with that of David Goldberg's, who touted the
predictability of IEEE floating point in his article on What Every
Computer Scientist Should Know About Floating Point. Doug Priest's
Appendix D brings Goldberg's article back to earth by explaining the
real-world consequences of letting the compiler choose the precision
for intermediate results. The Goldberg article and Priest's appendix
are available at http://www.validgh.com/.


Will
--------


; Given a procedure f representing a continuous real function,
; and an interval [x0, x1] such that (f x0) and (f x1) are of
; opposite signs, returns an inexact number approximating a
; root of f in that interval.
;
; This Scheme code works in systems that represent inexact
; numbers as IEEE floating point numbers with a fixed precision
; (single, double, extended). The result x is a best
; approximation to a root of the mathematical function obtained
; by linear interpolation from the values of the procedure f on
; the floating point numbers representable within the interval.
;
; This code may not work in systems that allow the precision of
; floating point computations to vary at the whim of the compiler.


(define (find-root f x0 x1)
    (let ((x0 (if (exact? x0) (exact->inexact x0) x0))
                (x1 (if (exact? x1) (exact->inexact x1) x1)))
        (define (loop lo hi)
            (let ((mid (/ (+ lo hi) 2)))
                (cond ((or (= mid lo) (= mid hi))
                              (choose lo hi))
                            (else (let ((fmid (f mid)))
                                            (cond ((zero? fmid) mid)
                                                        ((positive? fmid)
                                                          (loop lo mid))
                                                        ((negative? fmid)
                                                          (loop mid hi))
                                                        (else ???)))))))
        (define (choose x y)
            (let ((fx (f x))
                        (fy (f y)))
                (if (< (abs fx) (abs fy))
                        x
                        y)))
        (let ((fx0 (f x0))
                    (fx1 (f x1)))
            (cond ((zero? fx0) x0)
                        ((zero? fx1) x1)
                        ((and (positive? fx0) (negative? fx1))
                          (loop x1 x0)))
                        ((and (positive? fx1) (negative? fx0))
                          (loop x0 x1))
                        (else (error "(f x0) and (f x1) have the same sign" f x0 x1)))))

Return to the comp.compilers page.
Search the comp.compilers archives again.