Towards finding "the worst" input cases

[Tim Peters]

Jon raises some good points, and I have no argument with what he says.
I'm reduced again to saying (& apologizing (half-heartedly <grin>) for
this) that I just can't make time to spell out what these various cases
do-- and don't --really show.  They're to be used only by consenting
adults who've played with this stuff enough to fill in the (many)
missing details and qualifications for themselves.

Some minor glosses, then, on Jon's points:

>  ... although more careful and considered comments will have to wait
>  until much later.

That's the difference between us:  if I waited until I had the time I'd
never say anything <grin -- but this is the *last* thing my employer
wants me looking at, & I can't afford to take it out of "sleep time"
anymore>.

>  [description of similar kind of search based on looking for powers of
>   10 close to powers of 2, starting from a good value for log2(10)]

It's been so long since I've seen Matula's article (I know the one you
mean, though) that I can't remember how he approached it.  Continued
fractions are the "natural" way, & I both assume that's how you & he did
it and that the appearance of continued fractions in what I did is more
than coincidence ...

>  ...
>  But it seems to me that the "worst case" number you turned up, however
>  interesting in its own right, doesn't quite support the conclusion:

Absolutely right, that particular <grin> number doesn't.  These numbers
are "interesting" in that they're close to, but not exactly at, the
halfway point between adjacent representable numbers in a binary fp
format with 53 bits of precision (note that some of the cases, as noted
obliquely in a previous msg, aren't even relevant to the IEEE double
format per se).  Not all are relevant to nearest/even rounding, because
just as you note nearest/even is certain to get some right "by
predictable accident" (only the ones like xxx0.1000... and xxx1.0111...
are hard for nearest/even); similarly not all are relevant to nearest/up
("add .5 & chop") rounding either; etc.  What to *do* with these takes
some thought I haven't said much about, and truly sorry to have left
everyone to do it themselves.

On the other hand, don't undervalue the examples either.  E.g., suppose
that what I suspect "the worst" (in the sense above) case to be (for,
e.g., 17-digits-or-less input) is in fact the worst case (I'm sure now I
know *how* to find "the worst" case efficiently, not sure I wrote the
program-- which turned out to be pretty hairy --correctly).  Then if you
fake an fp format with enough extra precision to handle that case, you
will necessarily have faked enough extra precision to handle the cases
that are *actually* troublesome for IEEE nearest/even too -- that is,
while some of the "bad cases" I've trotted out are not actually relevant
to nearest/even rounding (although some are), the cases that *are*
actually troublesome to nearest/even rounding cannot be worse.

This shouldn't be sneezed at too quickly:  the worst case for 17-digit-
or-less input I found (and it may in fact-- or may not <grin> --be *the*
worst case> "only" required another 66 bits of precision.  If that's so,
the Clinger- and Gay-like escapes to much larger integer precisions
aren't really necessary (for 17-digit-or-less decimal -> IEEE double
input).

Extending this shaky <grin> chain, the benefit to *that* is that the
worst-case time could be dramatically reduced, and there are in turn
several benefits from that, two of which are subtle:

1) In parallel environments it's a Bad Thing to have operations that
   take time wildly dependent on the particular value of the inputs
   (e.g., that makes static load-balancing impossible from the word
   go).

2) As someone darkly hinted here long ago, the nastier secure sites some
   of us deal with have problems with system algorithms whose runtimes
   can be used in straightforward ways to deduce things about the
   algorithms' inputs ("hey, it took 50x longer this time, so they must
   be printing a number with an exponent > 200 ... that's one whale of a
   bomb" <0.7 grin>).

>  ...
>  Maybe this got discussed in subsequent discussions, so my apologies if
>  I'm just rehashing old comments.

No, surprisingly (to me), it didn't -- so I suspect everyone understood
that.  Or, that they didn't <grin -- ya just never know>.  I'm delighted
you brought it up, 'cuz I wouldn't have taken the time without a push.

>  Ideally, I'd like to program up a programmable-precision input
>  converter and give some experimental whacks at it; given that only a
>  few runs would ever be made by it, I'm sure that Lisp's bignums are
>  fully adequate in performance.

Addressing this to everyone: If you have a language with bignums, you
can whip up an arbitrary-precision fp converter in 15 minutes tops
(although whenever a change in systems or languages (or employers ...)
forces me to do this, I usually spend two more *days* grafting on all
the formatting options everyone wants ...).  The thrust of late has been
so much toward "fast" routines that it might help someone to sketch the
*easy* way to do it (skip this if you already know how):

1) Express the input as a rational number, top/bot; this is very easy
   with bignums.  I'll assume everything is positive, and that top isn't
   0, etc -- the end cases are obvious.

2) Assuming the output base is ob, and you want exactly nob digits in
   the significand (this happens to be Icon code -- straightforward,
   though -- "op:=" is augmented assignment, like C's "op=" -- pow(x,y)
   is x^y, but for "global" efficiency reasons it almost always pays to
   write a wrapper (like the "pow" here) that caches computed results
   for reuse):

   # "Normalize", maintaining the invariant
   #    top/bot  * ob^exp  =  original_top / original_bot
   exp := 0
   # 1:  Assure 1 <= top/bot
   until bot <= top do { top *:= ob; exp -:= 1 }
   # 2:  Assure top/bot < ob
   until top < ob*bot do { bot *:= ob; exp +:= 1 }

   # get the leading nob base-ob digits
   top *:= pow( ob, nob-1 )
   q := top / bot

That's basically it -- at this point q is an integer that takes exactly
nob base-ob digits to express (& I'm certainly not going to explain how
to convert a stinking integer to base ob <grin>), and the truncated
base-ob value is ("||" means string concatentation, and use your
imagination for the rest ...):

first_digit_in_q || "." || remaining_digits_in_q || "* ob^" || exp

You can stick in any rounding rule you like by computing

   remainder := top - q * bot

Then the excess is exactly remainder/bot ulp; so, e.g., if 2*remainder
is bot exactly, you're exactly halfway between two adjacent
representable numbers; etc.  If you *do* round up (add 1 to q), don't
forget to check for overflow:

  if q = pow( ob, nob ) then {
     q := pow( ob, nob-1 )
     exp +:= 1
  }

This has the advantage of being *very* easy to understand (perhaps it
takes a while to "get it" the first time, though), but runs real slow
because of the normalization loops.  How to speed that part up depends a
lot on your language; in Icon it happens to be extremely cheap to get
an *approximation* to the number of decimal digits in a bignum, and I
use that as a poor-man's log10 function to compute a good guess to exp
at the start:

   temp := -integer( (approx_dec_digs(top) - approx_dec_digs(bot)) *
	   	     log(10.,ob) )  # that's the log of 10. to base ob
   exp := -temp
   if temp > 0 then top *:= pow(ob,temp)
   else if temp < 0 then bot *:= pow(ob,-temp)

followed by the code above.  This usually gets top/bot in such a range
that the "until" loops find nothing to do, else almost always gets them
close enough that one iteration polishes it off.

>  [discussion of Lucid Common Lisp's READer algorithm]
>  ...
>     -- fallback to (horrors) bignums when all else fails.
>  I note that when the fallback is taken, frequently the result from the
>  wide-precision integer case is in fact accurate anyway.  To date, I haven't
>  been able to characterize the frequency or reason for this, but I suspect
>  there is something more lurking there than a simple coding bug (on the
>  other hand, I haven't really looked for the reason.)
>
>  So we have come back to the same old question I couldn't get an answer
>  to in 1977 -- namely, how can we prove (apart from exhaustion) that
>  "wide-integer" precision of exactly N bits is adequate for 64-bit
>  float-double conversions?

Well, I think it's clear we can't, because N bits aren't enough
regardless of how big N is ... if you'll put an upper bound on the
number of significant input digits you'll pay attention to, the question
is answerable (in theory), else since the possible decimal inputs are
dense in the reals the user can get arbitrarily close to any endcase the
rounding discipline has trouble with (hence no N is "big enough").

Once we put a bound on the number of decimal digits we'll look at, I
strongly suspect Kahan was correct in conjecturing that there's no
simple answer.  Several people sent me ways of attacking this based on
superficially different frameworks for analysis, but they all turned out
to be clearly equivalent to the problem of minimizing (in the modular
sense)

	(A + B*i) mod C   for some constant integers A,B,C

where i is constrained to lie in some contiguous range.  It doesn't take
a long time playing with this for small values of A,B,C before you're
struck with how darned *unpredictable* it is -- and since this
functional *form* is the basis of many hashing functions, and related to
venerable methods of pseudo-random number generation, it's probably not
surprising that it doesn't display easy-to-predict behavior <0.2 grin>.

In the analytical formulation of this form I wrote up for the "2^-1022"
cases, note that C had a value something like 2^731 (can't look it up
just now).  The smallest value (A+B*i) mod C can possibly take on
(besides 0) is obviously 1, and a rigorous lower bound follows from that
immediately:  if N is large enough to distinguish 1 part in 2^731 (beyond
the 53 bits computed), it's large enough period.  But that's pretty much
equivalent to saying that if you do the whole thing with exact rational
arithmetic you won't get in trouble <grin/sigh> ... I don't know if
there's a good reason for why truly low values of (A+B*i) mod C don't
show up for i in the relevant ranges, but the bounds on i are so small
compared to C that at this point I'd be genuinely surprised if a "truly
low value" *did* show up in range.

Ah well -- just rambling now.  I do hope someone can find the time to
makes some sense out this; it is, I'm convinced, much more a problem of
number theory than of real analysis (that is, analysis on the reals).

vanishing-into-the-night-ly y'rs  - tim

Tim Peters   Kendall Square Research Corp
timaksr.com,  ksr!timaharvard.harvard.edu


ps:  Another thing inappropriately (I think now) brushed off at the start
     of this series was the problem of recognizing an input exactly
     *at* a nearest/even endpoint.  I had a real slick solution for
     that, but unfortunately the Twinkie wrapper was too greasy to hold
     the ink <smile> ...