[PATCH] Fix old thinko in formula to compute sweight in numeric_sqrt(). - Mailing list pgsql-hackers

From Joel Jacobson
Subject [PATCH] Fix old thinko in formula to compute sweight in numeric_sqrt().
Date
Msg-id 06712c29-98e9-43b3-98da-f234d81c6e49@app.fastmail.com
Whole thread Raw
Responses Re: [PATCH] Fix old thinko in formula to compute sweight in numeric_sqrt().
List pgsql-hackers
Hi,

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.

It hasn't caused any problems, since the incorrect formula happens to
always produce the same result for DEC_DIGITS==4.

However, for other DEC_DIGITS values, it causes an undesired variation in the
precision of the results returned by sqrt().

To understand the problem, let's look at the equivalent formula for sweight,
when replacing DEC_DIGITS with the values 1, 2, 4:

HEAD, unpatched:
    sweight = (arg.weight + 1) * DEC_DIGITS / 2 - 1
Rewritten:
    (arg.weight + 1) * 1 / 2 - 1 <=> arg.weight / 2 - 1 / 2
    (arg.weight + 1) * 2 / 2 - 1 <=> arg.weight
    (arg.weight + 1) * 4 / 2 - 1 <=> 2 * arg.weight + 1

HEAD, patched:
    sweight = (arg.weight * DEC_DIGITS) / 2 + 1
Rewritten:
    (arg.weight * 1) / 2 + 1 <=> arg.weight / 2 + 1
    (arg.weight * 2) / 2 + 1 <=> arg.weight + 1
    (arg.weight * 4) / 2 + 1 <=> 2 * arg.weight + 1

As we can see, the equivalent formula for the patched version is arg.weight
times half the DEC_DIGITS, plus one.

The first part of the formula is the same but note how the patched version
gives a constant addition of `+ 1` regardless of the DEC_DIGITS value,
whereas the unpatched version gives strange subtractions/additions
such as `- 1 / 2` and `+ 3`.

Demonstration of the undesired result digit precision variation effect:

HEAD, unpatched:
DEC_DIGITS  sqrt(2::numeric)
4           1.414213562373095
2           1.4142135623730950
1           1.41421356237309505

HEAD, patched:
DEC_DIGITS  sqrt(2::numeric)
4           1.414213562373095
2           1.414213562373095
1           1.414213562373095

The patched version consistently returns 16 significant digits for sqrt(2::numeric)
when DEC_DIGITS is 1, 2 and 4, whereas the unpatched version surprisingly
gives 18 sig. digits for DEC_DIGITS==1 and 17 sig. digits for DEC_DIGITS==2.

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.

To see why we should get an additional digit for the DEC_DIGITS==4 case,
let's enable NUMERIC_DEBUG and look at the result:

SELECT sqrt(102::numeric);
make_result(): NUMERIC w=0 d=0 POS 0102
make_result(): NUMERIC w=0 d=15 POS 0010 0995 0493 8362 0780
        sqrt
--------------------
10.099504938362078
(1 row)

Since 10.099504938362 has only 14 sig. digits, we need one more NBASE digit
in the result, thus 0780 is computed, and we get an extra decimal digit for
free.

Compare this to DEC_DIGITS==2, which for the patched version correctly
returns 10.09950493836208, since the last produced NBASE digit `08`
is sufficient, i.e. with it, the result has 16 sig. decimal digits,
which is enough, since NUMERIC_MIN_SIG_DIGITS==16.

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().

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");

/Joel

pgsql-hackers by date:

Previous
From: "Hayato Kuroda (Fujitsu)"
Date:
Subject: RE: Assertion failure in SnapBuildInitialSnapshot()
Next
From: Noah Misch
Date:
Subject: Re: run pgindent on a regular basis / scripted manner