@@ -6169,6 +6169,12 @@ div_var(NumericVar *var1, NumericVar *var2, NumericVar *result,
61696169 *exactly. We compute DIV_GUARD_DIGITS extra digits, but there is
61706170 *no certainty that that's enough. We use this only in the transcendental
61716171 *function calculation routines, where everything is approximate anyway.
6172+ *
6173+ *Although we provide a "round" argument for consistency with div_var,
6174+ *it is unwise to use this function with round=false. In truncation mode
6175+ *it is possible to get a result with no significant digits, for example
6176+ *with rscale=0 we might compute 0.99999... and truncate that to 0 when
6177+ *the correct answer is 1.
61726178 */
61736179static void
61746180div_var_fast (NumericVar * var1 ,NumericVar * var2 ,NumericVar * result ,
@@ -6251,8 +6257,9 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
62516257
62526258/*
62536259 * We estimate each quotient digit using floating-point arithmetic, taking
6254- * the first four digits of the (current) dividend and divisor. This must
6255- * be float to avoid overflow.
6260+ * the first four digits of the (current) dividend and divisor. This must
6261+ * be float to avoid overflow. The quotient digits will generally be off
6262+ * by no more than one from the exact answer.
62566263 */
62576264fdivisor = (double )var2digits [0 ];
62586265for (i = 1 ;i < 4 ;i ++ )
@@ -6274,6 +6281,11 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
62746281 * To avoid overflow in maxdiv itself, it represents the max absolute
62756282 * value divided by NBASE-1, ie, at the top of the loop it is known that
62766283 * no div[] entry has an absolute value exceeding maxdiv * (NBASE-1).
6284+ *
6285+ * Actually, though, that holds good only for div[] entries after div[qi];
6286+ * the adjustment done at the bottom of the loop may cause div[qi + 1] to
6287+ * exceed the maxdiv limit, so that div[qi] in the next iteration is
6288+ * beyond the limit. This does not cause problems, as explained below.
62776289 */
62786290maxdiv = 1 ;
62796291
@@ -6325,10 +6337,10 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
63256337
63266338/*
63276339 * All the div[] digits except possibly div[qi] are now in the
6328- * range 0..NBASE-1.
6340+ * range 0..NBASE-1. We do not need to consider div[qi] in
6341+ * the maxdiv value anymore, so we can reset maxdiv to 1.
63296342 */
6330- maxdiv = Abs (newdig ) / (NBASE - 1 );
6331- maxdiv = Max (maxdiv ,1 );
6343+ maxdiv = 1 ;
63326344
63336345/*
63346346 * Recompute the quotient digit since new info may have
@@ -6348,7 +6360,16 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
63486360maxdiv += Abs (qdigit );
63496361}
63506362
6351- /* Subtract off the appropriate multiple of the divisor */
6363+ /*
6364+ * Subtract off the appropriate multiple of the divisor.
6365+ *
6366+ * The digits beyond div[qi] cannot overflow, because we know they
6367+ * will fall within the maxdiv limit. As for div[qi] itself, note
6368+ * that qdigit is approximately trunc(div[qi] / vardigits[0]),
6369+ * which would make the new value simply div[qi] mod vardigits[0].
6370+ * The lower-order terms in qdigit can change this result by not
6371+ * more than about twice INT_MAX/NBASE, so overflow is impossible.
6372+ */
63526373if (qdigit != 0 )
63536374{
63546375int istop = Min (var2ndigits ,div_ndigits - qi + 1 );
@@ -6360,9 +6381,25 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
63606381
63616382/*
63626383 * The dividend digit we are about to replace might still be nonzero.
6363- * Fold it into the next digit position. We don't need to worry about
6364- * overflow here since this should nearly cancel with the subtraction
6365- * of the divisor.
6384+ * Fold it into the next digit position.
6385+ *
6386+ * There is no risk of overflow here, although proving that requires
6387+ * some care. Much as with the argument for div[qi] not overflowing,
6388+ * if we consider the first two terms in the numerator and denominator
6389+ * of qdigit, we can see that the final value of div[qi + 1] will be
6390+ * approximately a remainder mod (vardigits[0]*NBASE + vardigits[1]).
6391+ * Accounting for the lower-order terms is a bit complicated but ends
6392+ * up adding not much more than INT_MAX/NBASE to the possible range.
6393+ * Thus, div[qi + 1] cannot overflow here, and in its role as div[qi]
6394+ * in the next loop iteration, it can't be large enough to cause
6395+ * overflow in the carry propagation step (if any), either.
6396+ *
6397+ * But having said that: div[qi] can be more than INT_MAX/NBASE, as
6398+ * noted above, which means that the product div[qi] * NBASE *can*
6399+ * overflow. When that happens, adding it to div[qi + 1] will always
6400+ * cause a cancelling overflow so that the end result is correct. We
6401+ * could avoid the intermediate overflow by doing the multiplication
6402+ * and addition in int64 arithmetic, but so far there appears no need.
63666403 */
63676404div [qi + 1 ]+= div [qi ]* NBASE ;
63686405
@@ -6381,9 +6418,12 @@ div_var_fast(NumericVar *var1, NumericVar *var2, NumericVar *result,
63816418div [qi ]= qdigit ;
63826419
63836420/*
6384- * Now we do a final carry propagation pass to normalize the result, which
6385- * we combine with storing the result digits into the output. Note that
6386- * this is still done at full precision w/guard digits.
6421+ * Because the quotient digits might be off by one, some of them might be
6422+ * -1 or NBASE at this point. The represented value is correct in a
6423+ * mathematical sense, but it doesn't look right. We do a final carry
6424+ * propagation pass to normalize the digits, which we combine with storing
6425+ * the result digits into the output. Note that this is still done at
6426+ * full precision w/guard digits.
63876427 */
63886428alloc_var (result ,div_ndigits + 1 );
63896429res_digits = result -> digits ;