On 18 November 2015 at 22:19, Tom Lane <t...@sss.pgh.pa.us> wrote:
> A bit of progress on this: I've found a test case, namely
>
> select sqrt(999999...
>
Wow.
> Still, this proves that we are onto a real problem.
>
Agreed. I had a go at turning that example into something using
log(base, num) so that the result would be visible in a pure SQL test,
but I didn't have any luck.
>> Another issue here is that with outercarry added into the qdigit
>> computation, it's no longer clear what the bounds of qdigit itself are,
>
> I concluded that that particular issue is a red herring: qdigit should
> always be a fairly accurate estimate of the next output digit, so it
> cannot fall very far outside the range 0..NBASE-1. Testing confirms that
> the range seen in the regression tests is exactly -1 to 10000, which is
> what you'd expect since there can be some roundoff error.
>
Right, I came to the same conclusion.
> Also, after further thought I've been able to construct an argument why
> outercarry stays bounded. See what you think of the comments in the
> attached patch, which incorporates your ideas about postponing the div[qi]
> calculation.
I think the logic is sound, but I worry that that comment is going to
be very difficult to follow.
I had a go at trying to find a simpler approach and came up with the
attached, which computes the value intended for div[qi+1] near the top
of the loop (using doubles) so that it can detect if it will overflow,
and if so it enters the normalisation block. It still needs some work
to prove that the normalised value for fnextdigit won't overflow, but
it feels like a promising, easier-to-follow approach to me. What do
you think?
Regards,
Dean
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index dcdc5cf..02146f7
*** 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,
! fnextdigit;
int qi;
int i;
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6274,6279 ****
--- 6275,6284 ----
* 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, excluding the first dividend digit, ie, strictly speaking
+ * the above holds only for div[i] where i > qi, at the top of the loop.
*/
maxdiv = 1;
*************** div_var_fast(NumericVar *var1, NumericVa
*** 6295,6371 ****
qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
(((int) fquotient) - 1); /* truncate towards -infinity */
! if (qdigit != 0)
{
! /* Do we need to normalize now? */
! maxdiv += Abs(qdigit);
! if (maxdiv > (INT_MAX - INT_MAX / NBASE - 1) / (NBASE - 1))
{
! /* Yes, do it */
! carry = 0;
! for (i = div_ndigits; i > qi; i--)
{
! newdig = div[i] + carry;
! if (newdig < 0)
! {
! carry = -((-newdig - 1) / NBASE) - 1;
! newdig -= carry * NBASE;
! }
! else if (newdig >= NBASE)
! {
! carry = newdig / NBASE;
! newdig -= carry * NBASE;
! }
! else
! 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;
! if (qi + i <= div_ndigits)
! fdividend += (double) div[qi + i];
}
! /* Compute the (approximate) quotient digit */
! fquotient = fdividend * fdivisorinverse;
! qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
! (((int) fquotient) - 1); /* truncate towards -infinity */
! maxdiv += Abs(qdigit);
}
! /* Subtract off the appropriate multiple of the divisor */
! if (qdigit != 0)
! {
! int istop = Min(var2ndigits, div_ndigits - qi + 1);
! for (i = 0; i < istop; i++)
! div[qi + i] -= qdigit * var2digits[i];
}
}
/*
! * 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;
}
--- 6300,6393 ----
qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
(((int) fquotient) - 1); /* truncate towards -infinity */
! /*
! * Do we need to normalize now?
! *
! * We will subtract a multiple of the divisor from the dividend array
! * causing the absolute value of each entry to increase by as much as
! * qdigit * (NBASE-1), so we need to ensure that that won't overflow.
! *
! * In addition, the remainder from the subtraction applied to div[qi]
! * is folded into the next digit div[qi+1], so we need to be careful
! * to ensure that it doesn't threaten to overflow either. Here
! * fnextdigit is the new value that we intend to assign to div[qi+1].
! */
! maxdiv += Abs(qdigit);
!
! fnextdigit = (double) div[qi + 1] +
! ((double) div[qi] - qdigit * var2digits[0]) * NBASE;
! if (var2ndigits > 1)
! fnextdigit -= qdigit * var2digits[1];
!
! if (maxdiv > (INT_MAX - INT_MAX / NBASE - 1) / (NBASE - 1) ||
! fabs(fnextdigit) > INT_MAX - INT_MAX / NBASE - 1)
{
! /* Yes, do it */
! carry = 0;
! for (i = div_ndigits; i > qi; i--)
{
! newdig = div[i] + carry;
! if (newdig < 0)
{
! carry = -((-newdig - 1) / NBASE) - 1;
! newdig -= carry * NBASE;
}
! else if (newdig >= NBASE)
{
! carry = newdig / NBASE;
! newdig -= carry * NBASE;
}
! else
! 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 = (double) div[qi];
! for (i = 1; i < 4; i++)
! {
! fdividend *= NBASE;
! if (qi + i <= div_ndigits)
! fdividend += (double) div[qi + i];
}
+ /* Compute the (approximate) quotient digit */
+ fquotient = fdividend * fdivisorinverse;
+ qdigit = (fquotient >= 0.0) ? ((int) fquotient) :
+ (((int) fquotient) - 1); /* truncate towards -infinity */
+ maxdiv += Abs(qdigit);
+
+ /* Also recompute the new value for div[qi + 1] */
+ fnextdigit = (double) div[qi + 1] +
+ ((double) div[qi] - qdigit * var2digits[0]) * NBASE;
+ if (var2ndigits > 1)
+ fnextdigit -= qdigit * var2digits[1];
}
/*
! * Subtract off the appropriate multiple of the divisor.
! *
! * We don't need to compute div[qi] or div[qi+1], since they have
! * already been computed and folded into fnextdigit.
*/
! if (qdigit != 0)
! {
! int istop = Min(var2ndigits, div_ndigits - qi + 1);
+ for (i = 2; i < istop; i++)
+ div[qi + i] -= qdigit * var2digits[i];
+ }
+ div[qi + 1] = (int) fnextdigit;
div[qi] = qdigit;
}
--
Sent via pgsql-hackers mailing list (pgsql-hackers@postgresql.org)
To make changes to your subscription:
http://www.postgresql.org/mailpref/pgsql-hackers
