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

I wrote:
> Note that div[qi+1], and indeed all the remaining dividend digits, are
> large negative values.  This is what you'd expect if the carry propagation
> step hasn't run for awhile, which is a precondition for div[qi] being
> large enough to cause an issue.  When we compute 218943 * 10000, we will
> indeed get an overflow, and the result will wrap around to some large
> negative value (2^32 less than it should be).  Then we will add that to
> div[qi+1], and we'll get *another* overflow, wrapping what nominally
> should have been a negative sum around to a positive value (2^32 more than
> it should be).  So the two overflows cancel and we get exactly the correct
> new value of div[qi+1].

> I do not know whether it's possible to devise a test case where you don't
> get offsetting overflows.  It may be that there's no operational bug here.
> Still, the code is surely not behaving per face value.

After further thought I've pretty well convinced myself that there is
indeed no observable bug, at least as long as you assume that overflow
within the multiplication will behave as stated above.  The proof goes
like this:

We already know that the divisor-subtraction step cannot cause an overflow
(modulo the question of div[qi] possibly exceeding the maxdiv limit,
which I'll get back to).  So what is at stake is just the possibility of
overflow in the final update that transfers leftover digits from div[qi]
to div[qi+1].

To analyze that, consider just the first two dividend and divisor terms in
the qdigit calculation, that is during each loop iteration approximate
qdigit asqdigit ~= trunc((div[qi]*NBASE + div[qi+1]) /                (var2digits[0]*NBASE + var2digits[1]))
This holds whether or not we do a carry propagation step.  Now write
what the divisor-subtraction step updates div[qi] to asdiv[qi]' = div[qi] - qdigit * var2digits[0]
and the updated div[qi+1] asdiv[qi+1]' = div[qi+1] - qdigit * var2digits[1]
So, if we temporarily disregard the possibility of overflow within the
final calculationdiv[qi+1]'' = div[qi+1]' + div[qi]'*NBASE
we can see that div[qi+1]'' is going to be= div[qi+1]' + div[qi]'*NBASE= div[qi+1] - qdigit * var2digits[1] +  (div[qi]
-qdigit * var2digits[0])*NBASE= div[qi]*NBASE + div[qi+1] -  qdigit * (var2digits[0]*NBASE + var2digits[1])
 
Comparing that to the approximate value of qdigit, we can see that
what we've got here is a modulo calculation, and the value of div[qi+1]''
is going to cancel to zero except for a remainder modulo
var2digits[0]*NBASE + var2digits[1], which of course must fall between 0
and NBASE*NBASE+NBASE (since the truncation is towards minus infinity).

Now of course this is only an approximation.  The main source of error is
that we've omitted the lower-order dividend terms.  Including div[qi+2]
could change the result by at most about INT_MAX/NBASE (which is the most
it could add to the numerator in the qdigit expression above, which will
propagate more or less linearly to div[qi+1]'').  Similarly, adding
div[qi+3] could add at most INT_MAX/NBASE/NBASE.  So we aren't anywhere
near the overflow threshold.  Including the lower-order divisor terms
could change the value of qdigit by a multiplicative factor of no more
than about 1/NBASE.  Lastly, roundoff error in the floating-point part
of the calculation could cause qdigit to be off by one count either
way.  None of these effects are going to let the final div[qi+1] value
get to more than two or three times NBASE squared, which is still
an order of magnitude less than INT_MAX.

Therefore, the final div[qi+1] value cannot overflow an int, and even
though it might be larger than what maxdiv would suggest, it's not going
to be large enough to cause overflow in the next loop iteration either.
It obviously won't cause overflow during carry propagation, and as for
the next divisor-subtraction step, we can do an analysis similar to the
above but approximating qdigit with just these terms:qdigit ~= trunc((div[qi] + div[qi+1]/NBASE) / var2digits[0])
Plugging that intodiv[qi]' = div[qi] - qdigit * var2digits[0]
shows that the updated div[qi] in any loop iteration is going to be just
about -div[qi+1]/NBASE, plus a truncation term that's between 0 and
var2digits[0], plus lower-order terms that aren't going to get you very
much past INT_MAX/NBASE.  So div[qi]' is never an overflow hazard in any
loop iteration.

In short, it's impossible for the end result of the div[qi+1] update
calculation to overflow, or even to get large enough to create a problem
in the next loop iteration.  However, the intermediate result
div[qi]*NBASE can overflow as per the example I showed before.

We could just ignore this on the grounds that it doesn't matter given sane
behavior of integer arithmetic.  Or, if we wanted to be more paranoid, we
could do the multiplication-and-addition in int64 arithmetic; but that
would probably slow things down a bit.

I think all this deserves some code comments, but I'm not sure whether
we need to bother with casting to int64.  Thoughts?
        regards, tom lane