Precision Tail-off?

[Oscar Benjamin]

On Fri, 17 Feb 2023 at 10:29, Stephen Tucker <stephen_tucker at sil.org> wrote:
>
> Thanks, one and all, for your reponses.
>
> This is a hugely controversial claim, I know, but I would consider this
> behaviour to be a serious deficiency in the IEEE standard.
[snip]
>
> Perhaps this observation should be brought to the attention of the IEEE. I
> would like to know their response to it.

Their response would be that they are well aware of what you are
saying and knew about this already since before writing any standards.
The basic limitation of the IEEE standard in this respect is that it
describes individual operations rather than composite operations. Your
calculation involves composing operations, specifically:

    result = x ** (n / d)

The problem is that there is more than one operation so we have to
evaluate this in two steps:

    e = n / d
    result = x ** e

Now the problem is that although n / d is correctly rounded e has a
small error because the exact value of n / d cannot be represented. In
the second operation taking this slightly off value of e as the
intended input means that the correctly rounded result for x ** e is
not the closest float to the true value of the *compound* operation.
The exponentiation operator in particular is very sensitive to changes
in the exponent when the base is large so the tiny error in e leads to
a more noticeable relative error in x ** e.

The only way to prevent this in full generality is to to have a system
in which no intermediate inexact operations are computed eagerly which
means representing expressions symbolically in some way. That is what
the SymPy code I showed does:

In [6]: from sympy import cbrt

In [7]: e = cbrt(123456789000000000000000000000000000000000)

In [8]: print(e)
100000000000*123456789**(1/3)

In [9]: e.evalf(50)
Out[9]: 49793385921817.447440261250171604380899353243631762

Because the *entire* expression is represented here *exactly* as e it
is then possible to evaluate different parts of the expression
repeatedly with different levels of precision and it is necessary to
do that for full accuracy in this case. Here evalf will use more than
50 digits of precision internally so that at the end you have a result
specified to 50 digits but where the error for the entire expression
is smaller than the final digit. If you give it a more complicated
expression then it will use even more digits internally for deeper
parts of the expression tree because that is what is needed to get a
correctly rounded result for the expression as a whole.

This kind of symbolic evaluation is completely outside the scope of
what the IEEE floating point standards are for. Any system based on
fixed precision and eager evaluation will show the same problem that
you have identified. It is very useful though to have a system with
fixed precision and eager evaluation despite these limitations. The
context for which the IEEE standards are mainly intended (e.g. FPU
instructions) is one in which fixed precision and eager evaluation are
the only option.

--
Oscar