Free42 possible accuracy flaw
|
03-24-2022, 05:01 PM
Post: #21
|
|||
|
|||
RE: Free42 possible accuracy flaw | |||
03-24-2022, 07:18 PM
(This post was last modified: 03-24-2022 07:22 PM by J-F Garnier.)
Post: #22
|
|||
|
|||
RE: Free42 possible accuracy flaw
(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 |
|||
03-24-2022, 07:45 PM
Post: #23
|
|||
|
|||
RE: Free42 possible accuracy flaw
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++! |
|||
03-25-2022, 07:03 AM
Post: #24
|
|||
|
|||
RE: Free42 possible accuracy flaw
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
03-25-2022, 08:03 AM
Post: #25
|
|||
|
|||
RE: Free42 possible accuracy flaw | |||
03-25-2022, 08:35 AM
Post: #26
|
|||
|
|||
RE: Free42 possible accuracy flaw
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... |
|||
03-25-2022, 08:53 AM
Post: #27
|
|||
|
|||
RE: Free42 possible accuracy flaw
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
03-25-2022, 09:06 AM
Post: #28
|
|||
|
|||
RE: Free42 possible accuracy flaw
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. |
|||
03-25-2022, 09:18 AM
Post: #29
|
|||
|
|||
RE: Free42 possible accuracy flaw
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.
|
|||
03-25-2022, 11:12 AM
(This post was last modified: 03-25-2022 11:14 AM by J-F Garnier.)
Post: #30
|
|||
|
|||
RE: Free42 possible accuracy flaw
(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(). 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 |
|||
03-25-2022, 02:58 PM
Post: #31
|
|||
|
|||
RE: Free42 possible accuracy flaw
(03-22-2022 10:08 PM)Thomas Okken Wrote: 10^2 = 99.99999999999999999999999999999999 Going thru the source, bid128_pow.c, I may have spotted the issue. Code: // Compute dominant exponential term 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. |
|||
03-27-2022, 10:58 PM
Post: #32
|
|||
|
|||
RE: Free42 possible accuracy flaw
(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 |
|||
03-30-2022, 08:04 AM
Post: #33
|
|||
|
|||
RE: Free42 possible accuracy flaw
(03-25-2022 07:03 AM)Werner Wrote: Congatulations, Thomas! 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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
03-30-2022, 11:51 PM
Post: #34
|
|||
|
|||
RE: Free42 possible accuracy flaw
(03-30-2022 08:04 AM)Werner Wrote: 531441^1072.5 vs 9^6435 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 |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: