Free42 possible accuracy flaw

[Albert Chan]

(03-22-2022 10:08 PM)Thomas Okken Wrote:  10^2 = 99.99999999999999999999999999999999

That's just the first one of many where bid128_pow() returns an inexact result for an easy case...

Going thru the source, bid128_pow.c, I may have spotted the issue.

Code:

// Compute dominant exponential term
// We really want exp(l + l_lo)

  BIDECIMAL_CALL1(bid128_exp,res,l);

...

  BIDECIMAL_CALL3(bid128_fma,res,res,l_lo,res);

log of result, is done in excess precision, say, accurate to 128+ bits.
However, it use bid128_exp to recover result, *then* adjust.

exp(l + l_lo) = exp(l) * exp(l_lo) ≈ exp(l) * (1 + l_lo)

But, RHS exp(l) already rounded to 34 decimal digits.
Essentially I think bid128_pow does this with the l_lo correction.


10 LN 2 *       → 4.605170185988091368035982909368728
E^X             → 99.99999999999999999999999999999995
4.152E-34       // Note: no need to use fma for correction
X<>Y * LASTX +  → 99.99999999999999999999999999999999

It should have call a custom exp, to handle exp(l + l_lo), *then* round.

round34(exp(4.6051701859880913680359829093687284152)) = 100.