Re: floating point

William D Clinger <will@ccs.neu.edu>
5 Oct 1998 20:50:44 -0400

          From comp.compilers

Related articles
Re: inlining + optimization = nuisance bugs luddy@concmp.com (Luddy Harrison) (1998-09-29)
Re: floating point, was inlining + optimization = nuisance bugs chase@world.std.com (David Chase) (1998-10-04)
Re: floating point will@ccs.neu.edu (William D Clinger) (1998-10-05)
Re: floating point comments@cygnus-software.com (Bruce Dawson) (1998-10-07)
Re: floating point will@ccs.neu.edu (William D Clinger) (1998-10-10)
Re: floating point dmcq@fano.demon.co.uk (David McQuillan) (1998-10-13)
Re: floating point darcy@CS.Berkeley.EDU (Joseph D. Darcy) (1998-10-19)
Re: floating point darcy@usul.CS.Berkeley.EDU (1998-10-24)
Re: floating point comments@cygnus-software.com (Bruce Dawson) (1998-11-01)
[8 later articles]
| List of all articles for this month |

From: William D Clinger <will@ccs.neu.edu>
Newsgroups: comp.compilers
Date: 5 Oct 1998 20:50:44 -0400
Organization: Northeastern University
References: 98-09-164 98-10-018
Keywords: arithmetic

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


Post a followup to this message

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