Re: Bug in numeric multiplication - Mailing list pgsql-hackers

I wrote:
> I pushed this, but while looking at it, my attention was drawn to this
> bit down near the end of the loop:

>         /*
>          * The dividend digit we are about to replace might still be nonzero.
>          * Fold it into the next digit position.  We don't need to worry about
>          * overflow here since this should nearly cancel with the subtraction
>          * of the divisor.
>          */
>         div[qi + 1] += div[qi] * NBASE;

> In the first place, I'm not feeling especially convinced by that handwave
> about there being no risk of overflow.  In the second place, this is
> indisputably failing to consider whether maxdiv might need to increase.
> If I add
>         Assert(Abs(div[qi + 1]) <= (maxdiv * (NBASE-1)));
> right after this, the regression tests blow up.  So it looks to me like
> we've got some more work to do.

After thinking some more about what this is doing, it seems to me that
we could avoid changing div[qi + 1] at all here, and instead deal with
leftover dividend digits by shoving them into the floating-point part of
the calculation.  All that we really need to have happen as a consequence
of the above is to inject an additional "div[qi] * NBASE" term into future
calculations of fdividend.  So we can do that out-of-band and avoid
mucking up the maxdiv invariant.

Also, I realized that this bit:

                /*
                 * All the div[] digits except possibly div[qi] are now in the
                 * range 0..NBASE-1.
                 */
                maxdiv = Abs(newdig) / (NBASE - 1);
                maxdiv = Max(maxdiv, 1);

is unduly conservative, because, while indeed div[qi] might be out of
range, there is no need to consider it in maxdiv anymore, because the new
maxdiv value only determines what will happen in future loop iterations
where the current div[qi] will be irrelevant.  So we should just reset
maxdiv to 1 unconditionally here.

So that led me to the attached patch, which passes regression tests fine.
I'm not sure that it's provably okay though.  The loose ends are to show
that div[qi] can't overflow an int during the divisor-subtraction step
and that "outercarry" remains bounded.   Testing suggests that outercarry
can't exceed INT_MAX/NBASE, but I don't see how to prove that.

            regards, tom lane

diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index dcdc5cf..03a646b 100644
*** a/src/backend/utils/adt/numeric.c
--- b/src/backend/utils/adt/numeric.c
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6186,6192 ****
      double        fdividend,
                  fdivisor,
                  fdivisorinverse,
!                 fquotient;
      int            qi;
      int            i;

--- 6186,6193 ----
      double        fdividend,
                  fdivisor,
                  fdivisorinverse,
!                 fquotient,
!                 outercarry;
      int            qi;
      int            i;

*************** div_var_fast(NumericVar *var1, NumericVa
*** 6274,6289 ****
       * To avoid overflow in maxdiv itself, it represents the max absolute
       * value divided by NBASE-1, ie, at the top of the loop it is known that
       * no div[] entry has an absolute value exceeding maxdiv * (NBASE-1).
       */
      maxdiv = 1;

      /*
!      * Outer loop computes next quotient digit, which will go into div[qi]
       */
      for (qi = 0; qi < div_ndigits; qi++)
      {
          /* Approximate the current dividend value */
!         fdividend = (double) div[qi];
          for (i = 1; i < 4; i++)
          {
              fdividend *= NBASE;
--- 6275,6297 ----
       * To avoid overflow in maxdiv itself, it represents the max absolute
       * value divided by NBASE-1, ie, at the top of the loop it is known that
       * no div[] entry has an absolute value exceeding maxdiv * (NBASE-1).
+      *
+      * Note that maxdiv only includes digit positions that are still part of
+      * the dividend, ie, strictly speaking the above holds only for div[i]
+      * where i >= qi, at the top of the loop.
       */
      maxdiv = 1;

      /*
!      * Outer loop computes next quotient digit, which will go into div[qi].
!      * "outercarry" represents high-order dividend digits carried across
!      * iterations.
       */
+     outercarry = 0;
      for (qi = 0; qi < div_ndigits; qi++)
      {
          /* Approximate the current dividend value */
!         fdividend = outercarry * NBASE + (double) div[qi];
          for (i = 1; i < 4; i++)
          {
              fdividend *= NBASE;
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6320,6340 ****
                          carry = 0;
                      div[i] = newdig;
                  }
!                 newdig = div[qi] + carry;
!                 div[qi] = newdig;

                  /*
                   * All the div[] digits except possibly div[qi] are now in the
!                  * range 0..NBASE-1.
                   */
!                 maxdiv = Abs(newdig) / (NBASE - 1);
!                 maxdiv = Max(maxdiv, 1);

                  /*
                   * Recompute the quotient digit since new info may have
                   * propagated into the top four dividend digits
                   */
!                 fdividend = (double) div[qi];
                  for (i = 1; i < 4; i++)
                  {
                      fdividend *= NBASE;
--- 6328,6347 ----
                          carry = 0;
                      div[i] = newdig;
                  }
!                 div[qi] += carry;

                  /*
                   * All the div[] digits except possibly div[qi] are now in the
!                  * range 0..NBASE-1; and we don't care about including div[qi]
!                  * anymore, so we can reset maxdiv to 1.
                   */
!                 maxdiv = 1;

                  /*
                   * Recompute the quotient digit since new info may have
                   * propagated into the top four dividend digits
                   */
!                 fdividend = outercarry * NBASE + (double) div[qi];
                  for (i = 1; i < 4; i++)
                  {
                      fdividend *= NBASE;
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6360,6370 ****

          /*
           * The dividend digit we are about to replace might still be nonzero.
!          * Fold it into the next digit position.  We don't need to worry about
!          * overflow here since this should nearly cancel with the subtraction
!          * of the divisor.
           */
!         div[qi + 1] += div[qi] * NBASE;

          div[qi] = qdigit;
      }
--- 6367,6382 ----

          /*
           * The dividend digit we are about to replace might still be nonzero.
!          * Carry it forward to future division steps in "outercarry".
!          *
!          * Although it might appear that outercarry would soon grow large
!          * enough to lose precision, this is not actually possible because
!          * each divisor-subtraction step will drive the accumulated value
!          * towards zero.  We could probably store it in an int without
!          * overflow problems, but since the value is only used in floating
!          * calculations, we might as well store it as a double.
           */
!         outercarry = outercarry * NBASE + (double) div[qi];

          div[qi] = qdigit;
      }
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6372,6378 ****
      /*
       * Approximate and store the last quotient digit (div[div_ndigits])
       */
!     fdividend = (double) div[qi];
      for (i = 1; i < 4; i++)
          fdividend *= NBASE;
      fquotient = fdividend * fdivisorinverse;
--- 6384,6390 ----
      /*
       * Approximate and store the last quotient digit (div[div_ndigits])
       */
!     fdividend = outercarry * NBASE + (double) div[qi];
      for (i = 1; i < 4; i++)
          fdividend *= NBASE;
      fquotient = fdividend * fdivisorinverse;