On Sat, 28 Jan 2023 at 22:14, Joel Jacobson <j...@compiler.org> wrote: > > I found what appears to be a small harmless error in numeric.c, > that seems worthwhile to fix only because it's currently causes confusion. >
Shrug. Looking at git blame, it's been like that for about 20 years, and I wasn't aware of it causing confusion. > HEAD, patched: > sweight = (arg.weight * DEC_DIGITS) / 2 + 1 > You haven't actually said why this formula is more correct than the current one. I believe that it is when arg.weight >= 0, but I'm not convinced it's correct for arg.weight < 0. Given that this is only an approximate computation, which ensures roughly 16 significant digits in the result, but can't guarantee exactly 16 digits, I'm not convinced of the benefits of changing it. > Note, however, that it's still possible to find examples of when sqrt(numeric) > produce results with different precision for different DEC_DIGITS/NBASE > values, > but in such cases, it's intentional, and due to getting additional precision > for free, since the larger the NBASE, the more decimal digits are produced > at the same time per iteration in the calculation. > > Example: > > HEAD, unpatched > DEC_DIGITS sqrt(102::numeric) > 4 10.09950493836208 > 2 10.099504938362078 > 1 10.0995049383620780 > > HEAD, patched: > DEC_DIGITS sqrt(102::numeric) > 4 10.099504938362078 > 2 10.09950493836208 > 1 10.09950493836208 > > According to the comment in numeric_sqrt(), the goal is to give at least > NUMERIC_MIN_SIG_DIGITS (16) significant digits. > > Since 10.09950493836208 has 16 significant digits, we can see above how > DEC_DIGITS==2 causes an additional unnecessary significant digit to be > computed, > and for DEC_DIGITS==1, two additional unnecessary significant digits are > computed. > > The patched version returns 16 significant digits as expected for > DEC_DIGITS==2 > and DEC_DIGITS==1, and for DEC_DIGITS==4 we get an additional digit for free. > You lost me here. In unpatched HEAD, sqrt(102::numeric) produces 10.099504938362078, not 10.09950493836208 (with DEC_DIGITS = 4). And wasn't your previous point that when DEC_DIGITS = 4, the new formula is the same as the old one? > In conclusion, the proposed patch fixes a harmless problem, but is important > to fix, since otherwise, anyone who want to experiment with different > DEC_DIGITS/NBASE combinations by changing the `#if 0` preprocessor values > in the top of numeric.c will get surprising results from sqrt(). > Anyone changing DEC_DIGITS/NBASE will find hundreds of regression test failures due to changes in result precision all over the place (and failures due to the overflow limit changing). I don't see why numeric_sqrt() should be singled out for fixing. > In passing, also add pow10[] values for DEC_DIGITS==2 and DEC_DIGITS==1, > since otherwise it's not possible to compile such DEC_DIGITS values > due to the assert: > > StaticAssertDecl(lengthof(pow10) == DEC_DIGITS, "mismatch with > DEC_DIGITS"); > That might be worth doing, to ensure that the code still compiles for other DEC_DIGITS/NBASE values. I'm not sure how useful that really is anymore though. As the comment at the top says, it's kept mostly for historical reasons. Regards, Dean