[Python-Dev] Round Bug in Python 1.6?

Tim Peters tim_one@email.msn.com
Sat, 8 Apr 2000 21:26:23 -0400


[Vladimir Marangozov]
> I'm not sure that this will be of any interest to you, number crunchers,
> but a research team in computer arithmetics here reported some major
> results lately: they claim that they "solved" the Table Maker's Dilemma
> for most common functions in IEEE-754 double precision arithmetic.
> (and no, don't ask me what this means ;-)

Back in the old days, some people spent decades making tables of various
function values.  A common way was to laboriously compute high-precision
values over a sparse grid, using e.g. series expansions, then extend that to
a fine grid via relatively simple interpolation formulas between the
high-precision results.

You have to compute the sparse grid to *some* "extra" precision in order to
absorb roundoff errors in the interpolated values.  The "dilemma" is
figuring out how *much* extra precision:  too much and it greatly slows the
calculations, too little and the interpolated values are inaccurate.

The "problem cases" for a function f(x) are those x such that the exact
value of f(x) is very close to being exactly halfway between representable
numbers.  In order to round correctly, you have to figure out which
representable number f(x) is closest to.  How much extra precision do you
need to use to resolve this correctly in all cases?

Suppose you're computing f(x) to 2 significant decimal digits, using 4-digit
arithmetic, and for some specific x0 f(x0) turns out to be 41.49 +- 3.
That's not enough to know whether it *should* round to 41 or 42.  So you
need to try again with more precision.  But how much?  You might try 5
digits next, and might get 41.501 +- 3, and you're still stuck.  Try 6 next?
Might be a waste of effort.  Try 20 next?  Might *still* not be enough -- or
could just as well be that 7 would have been enough and you did 10x the work
you needed to do.

Etc.  It turns out that for most functions there's no general way known to
answer the "how much?" question in advance:  brute force is the best method
known.

For various IEEE double precision functions, so far it's turned out that you
need in the ballpark of 40-60 extra accurate bits (beyond the native 53) in
order to round back correctly to 53 in all cases, but there's no *theory*
supporting that.  It *could* require millions of extra bits.

For those wondering "why bother?", the practical answer is this:  if a std
could require correct rounding, functions would be wholly portable across
machines ("correctly rounded" is precisely defined by purely mathematical
means).  That's where IEEE-754 made its huge break with tradition, by
requiring correct rounding for + - * / and sqrt.  The places it left fuzzy
(like string<->float, and all transcendental functions) are the places your
program produces different results when you port it.

Irritating one:  MS VC++ on Intel platforms generates different code for
exp() depending on the optimization level.  They often differ in the last
bit they compute.  This wholly accounts for why Dragon's speech recognition
software sometimes produces subtly (but very visibly!) different results
depending on how it was compiled.  Before I got tossed into this pit, it was
assumed for a year to be either a -O bug or somebody fetching uninitialized
storage.

that's-what-you-get-when-you-refuse-to-define-results-ly y'rs  - tim