(Free42) accuracy of LN1+X

[Albert Chan]

(06-14-2020 01:09 AM)Paul Dale Wrote:  This function is actually fairly easy to implement once its problems are realised.

Probably not.

decimal128 unable to round-trip thru binary128 occurs at all ranges, not just tiny x
We just don't noticed it when x is bigger ...

For example, try this:

1e-12 XEQ "LN1+X"     → 9.999999999995000000000003333333334e-13

We expected result should be rounded(x - x²/2 + x³/3), with last digit = 3

x = 1e-12 converted to binary128, we got x' ≈ (10^33+0.08)/10^45, over-estimated 0.45 ULP

log1p(x') ≈ 9.99999999999500000000000333333333413e-13

This example is lucky that roundings compensated the x'-x error.
You can potientially get into double-rounding errors, all on the same side.

However, even if conversion is exact, we still cannot get correctly rounded result.

>>> from gmpy2 import *
>>> get_context().precision = 113 # binary128
>>> y = log1p(mpfr('1e-12'))
>>> format(y, '.33e')
'9.999999999995000000000003333333334e-13'
>>> format(next_below(y), '.33e')
'9.999999999995000000000003333333332e-13'