5 Oct 1998 20:50:44 -0400

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

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.