The Table-Maker's Dilemma and correctly-rounded transcendental functions

[David Hough]

Jon White relates:

> Sometime in the late 1970's, I had a chance to talk about these ideas
> with Kahan at Berkeley; I asked him if he could show me a "closed
> form" for the relation between the significand and exponent sizes and
> the amount of extra precision required.  He conjectured that there was
> no simple, mathematical formula for this.  Since every implementation
> of the IEEE standard will have only a finite number of representable
> numbers, then in theory you could just find out the answer by exhaustion.
> I think I did so for the PDP10 format (23-bit significand, but not IEEE
> format), but don't recall the sharp result; it was either 6 or 7 or 8
> extra bits needed.

A lot of people, including me, have been confused about the
implications of the table-maker's dilemma, which strictly speaking
applies to transcendental functions but as a practical matter affects
people's thinking about algebraic and even rational functions like base
conversion.

The table-maker's dilemma derives from this fact about the usual
transcendental functions f:  given a machine-representable x, the
amount of computation necessary to determine a correctly-rounded f(x)
can't be bounded in advance independent of x.  Or in White's framework,
the amount of extra precision necessary can't be determined
analytically.  This is because f(x) may in general be arbitrarily close
to a rounding boundary between representable numbers.

So the table-maker's dilemma, in general, when he appears to be
computing an f(x) that's exactly between two representable numbers, is
when to quit computing to higher precision and start trying to prove
that he's discovered a previously unknown rational x for which f(x) is
rational.  I think for the commonest transcendental f's you can prove
that there aren't any such undiscovered cases, but the issue is
probably open for higher transcendental functions.

The correct fact, and the correct dilemma, give rise to the following
commonly-held but incorrect inference:

	It is impossible to compute correctly-rounded elementary
	transcendental functions.

or in another form:

	Any algorithm that computes correctly-rounded elementary
	transcendental functions must be hopelessly expensive.

These inferences are wrong, because they overlook the fact that White
remembered in the context of base conversion: there are only finitely
many representable numbers.  For the common transcendentals, all the
f(x) except the obvious ones are not rational, hence not on any
boundary between representable numbers.  For each representable number,
there will be some high enough precision in which you can conclude that
f(x) is definitively on one side or another of a boundary.  Thus there
is a worst case, and an algorithm that computed in that high enough
precision would suffice.  That takes care of the first form of the
incorrect inference.

The second form of the incorrect inference is trickier, but again
White's experience with base conversion points out the RISCy
resolution: handle the common case cheaply and quickly; fall back to a
more expensive method only when the cheap computation is inconclusive.
For instance, correctly-rounded single-precision elementary
transcendental functions can be computed relatively cheaply in average
cost by using a double-precision algorithm that provides an error bound
that tells you whether the single- precision result is correctly
rounded to single precision.  In the rare cases that fails, you can use
integer arithmetic or probably even table lookup to resolve the issue.

Acceptably fast correctly-rounded double-precision transcendentals are
still a research area unless high-quality quadruple-precision hardware
is available; that just defers the problem to correctly-rounded
quadruple-precision transcendentals.  I think the general approach for
building correctly-rounded working-precision transcendentals out of
only working-precision variables and registers has to depend on
appropriate use of doubleD-precision arithmetic, as I described a while
back, to produce the unrounded working-precision f(x) as a sum of two
working-precision variables u + v and an error bound z, so that 
	|f(x) - (u+v)| <= z.  
Then if (u+v+z) rounded to working precision in the
positive direction, and (u+v-z) rounded to working precision in the
negative direction, both produce the same result, you're done.
(Three-operand addition was mentioned in the same note as doubleD
precision.) If that fails, and you had analytical reasons to expect
failures to be extremely rare, the first thing to try (during
development) is to print out a diagnostic, with the intent of building
up a table of exceptional values that have to be specially handled.  If
the table is known in advance or discovered empirically to be too big,
then you would have to implement another algorithm in higher-precision
integer arithmetic to decide the case.

Customers will disagree about whether this is too high a price to pay
for getting the same numerical results on different IEEE 754 machines
in a networked environment, particularly since that goal also requires
constraints on expression evaluation that eliminate some of the
potential benefit of micro-tasking parallel processors.