Free42 possible accuracy flaw

[Thomas Okken]

(03-24-2022 04:53 PM)Albert Chan Wrote:  Can I assume new implementation is never worse than internal bid128_pow() ?

Yes: it only uses multiplications when the result of that calculation is exact, and bid128_pow() in all other cases.

[J-F Garnier]

(03-24-2022 10:20 AM)J-F Garnier Wrote:  I can't find any similar problem with Free42 Binary (still 3.0.10):

(03-24-2022 11:02 AM)Thomas Okken Wrote:  The binary version simply uses the standard C library's pow() function, without any separate special case handling. So the behavior you're seeing there is whatever the standard pow() provides. (Which may not be the same everywhere!)

FWIW, Free42 1.4.71 running on my ancient Pocket PC (w/ Windows Mobile 5.0, do you remember it?), both in decimal (25-digit incarnation) and binary versions, passes the test without problem, with the last place correctly rounded. Great!

I'm still using my Pocket PC for various reasons with the benefit of having both binary and decimal versions of Free42, plus Emu42.

J-F

[Thomas Okken]

Yes, Hugh Steers' library was pretty good! I switched to the Intel library in order to get IEEE-754-2008 features like denormal numbers and unbiased rounding, but if I had known about the problems with that library in advance, maybe I wouldn't have switched.

Regarding decimal and binary side-by-side: my Android and iOS projects are set up so I can build both, and they are built with different app IDs, so they can be installed side by side. The only thing I haven't dealt with is sharing state files between them the way the desktop versions do. Since there is essentially no demand, I haven't bothered to spend time on that, I have built the binary versions a couple of times when people specifically asked me, but that's it.

Regarding Windows Mobile: funny you should mention that. The Windows Mobile version was the most popular download on my web site for a long time, but I've actually never owned such a device, and never even seen Free42 run on one. I did all my work on it with the simulator in eVC++!

[Werner]

Congatulations, Thomas!
3^729 is now equal to (3^9)^81, a feat never seen before (when using floating-point numbers).
Kahan himself mentioned this particular example somewhere in an article of his, and why the results would not agree.
Cheers, Werner

[J-F Garnier]

(03-25-2022 07:03 AM)Werner Wrote:  Congatulations, Thomas!
3^729 is now equal to (3^9)^81, a feat never seen before (when using floating-point numbers).

Matter of fact, Free42 1.4x (so using Hugh Steers' library) also passed the test!
The congratulations are still valid, of course.

J-F

[Thomas Okken]

Credit for 3^729 being equal to (3^9)^81 should go to Intel then (and Hugh), all I did was get out of the way. My logic will now handle only up to 3^71, that being the highest power of 3 where the exact result fits in 34 digits. Now I am curious what makes that an interesting test case, though!

I guess I should look at Y^X with complex Y and integer X next...

[Werner]

LN(3) and LN(3^9) are both between 1 and 10, so the rounding error is the same. Yet in 3^729, it is multiplied by 729, and in (3^9)^81, only by 81 before the exp().
So the larger the base, the more accurate ^ is.
The 42S gives different results for these.
Werner

[Thomas Okken]

Ah, OK. And on top of that, the relative rounding error in numbers just above 1 is almost 10 times larger than that in numbers just below 10.
But this is just an anecdote, no? I bet examples could be found where it goes the other way.

[Thomas Okken]

Incidentally, Hugh Steers' library uses repeated squaring for integer powers too, even for large exponents, but it also uses double precision internally. And by 'double,' I mean 15 base-10000 digits, so effectively 32 decimal guard digits.

[J-F Garnier]

(03-25-2022 08:53 AM)Werner Wrote:  LN(3) and LN(3^9) are both between 1 and 10, so the rounding error is the same. Yet in 3^729, it is multiplied by 729, and in (3^9)^81, only by 81 before the exp().
So the larger the base, the more accurate ^ is.
The 42S gives different results for these.

On the HP Saturn machines, the 3 internal guard digits are just a bit too short for power operations on large numbers.
I gave an example above, but here is one of the worst cases:
1E44^10.5 --> 9.99999999991E461 (9 ULP error)
Note that LN(1E44) has a mantissa close to 1.

J-F

[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.

[Albert Chan]

(03-25-2022 08:35 AM)Thomas Okken Wrote:  Now I am curious what makes that [3^729] an interesting test case, though!

For 34 decimal digits, 3^729 is close to half-way.

round38(3^729) = 6.6281860542418717610517286421447974859E+347

However, Intel Decimal get this right really by luck. (see previous post)
The flawed step, correction *after* rounded exponential, does not hurt the result.

round38(log3 * 729) = 800.88835843905196502713377771652123869

800.8883584390519650271337777165212
E^X                        → 6.628186054241871761051728642144541e347
3.869E-32
X<>Y * LASTX +      → 6.628186054241871761051728642144797e347

[Werner]

(03-25-2022 07:03 AM)Werner Wrote:  Congatulations, Thomas!
3^729 is now equal to (3^9)^81, a feat never seen before (when using floating-point numbers).
Kahan himself mentioned this particular example somewhere in an article of his, and why the results would not agree.
Cheers, Werner

Found the reference (In 'Mathemathics written in Sand', pg 27), and as usual, I misremembered ;-)
The actual calculations are

729^33.5 vs. 3^201

Where on 10-digit machines the latter has an error of more than one ulp, that could only be avoided using more than three extra digits in the intermediary calculations.

Equivalent numbers for Free42 would be

531441^1072.5 vs 9^6435

The former is off by 2 ulps. (I must confess I'm a bit puzzled by the fact that in this case it is the one with the larger base that is off - which goes against my original reasoning. I'm sure Albert will come to the rescue ;-))

Cheers, Werner

[Albert Chan]

(03-30-2022 08:04 AM)Werner Wrote:  531441^1072.5 vs 9^6435

The former is off by 2 ulps.

Normally, we would not expect 2 ulp error or more.
(except when exponent dropped, say from 10 to 9.999...)

Measure in fractions of ulp, from exact result, we expected EXP part 1/2 ulp or less,
correction part of 1/2 ulp or less. Combined error almost always less than 1 full ulp.

For above example, log of result gives:

14139. 14015 51585 71728 25680 61991 92905 81881 ...

Here, log result had huge integer part, which kills the precision.
If we had 128 good bits, (about 38 digits), say 81881 → 81900, we recovered 33 digits:

81881 → 3.5526 15792 57691 44387 92654 15646 88743 56604E+6140
81900 → 3.5526 15792 57691 44387 92654 15646 88750 31601E+6140

Even with correct recovery, EXP *then* ROUND, we were off 1 ULP.
The flawed recover part added another possible ± 1 ULP