On Fri, 5 Jul 2024 at 12:56, Joel Jacobson <j...@compiler.org> wrote:
>
> Interesting you got so bad bench results for v6-mul_var_int64.patch
> for var1ndigits=4, that patch is actually the winner on AMD Ryzen 9 7950X3D.
Interesting.
> On Intel Core i9-14900K the winner is
> v6-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.
That must be random noise, since
v6-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch doesn't
invoke mul_var_small() for 4-digit inputs.
> On Apple M3 Max, HEAD is the winner.
Importantly, mul_var_int64() is around 1.25x slower there, and it was
even worse on my machine.
Attached is a v7 mul_var_small() patch adding 4-digit support. For me,
this gives a nice speedup:
SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_4;
Time: 5617.150 ms (00:05.617) -- HEAD
Time: 8203.081 ms (00:08.203) -- v6-mul_var_int64.patch
Time: 4750.212 ms (00:04.750) -- v7-mul_var_small.patch
The other advantage, of course, is that it doesn't require 128-bit
integer support.
Regards,
Dean
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 5510a20..b9497e1
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -551,6 +551,8 @@ static void sub_var(const NumericVar *va
static void mul_var(const NumericVar *var1, const NumericVar *var2,
NumericVar *result,
int rscale);
+static void mul_var_small(const NumericVar *var1, const NumericVar *var2,
+ NumericVar *result, int rscale);
static void div_var(const NumericVar *var1, const NumericVar *var2,
NumericVar *result,
int rscale, bool round);
@@ -8707,7 +8709,7 @@ mul_var(const NumericVar *var1, const Nu
var1digits = var1->digits;
var2digits = var2->digits;
- if (var1ndigits == 0 || var2ndigits == 0)
+ if (var1ndigits == 0)
{
/* one or both inputs is zero; so is result */
zero_var(result);
@@ -8715,6 +8717,17 @@ mul_var(const NumericVar *var1, const Nu
return;
}
+ /*
+ * If var1 has 4 digits or fewer, and we are computing the exact result,
+ * with no rounding, delegate to mul_var_small() which uses a faster short
+ * multiplication algorithm.
+ */
+ if (var1ndigits <= 4 && rscale == var1->dscale + var2->dscale)
+ {
+ mul_var_small(var1, var2, result, rscale);
+ return;
+ }
+
/* Determine result sign and (maximum possible) weight */
if (var1->sign == var2->sign)
res_sign = NUMERIC_POS;
@@ -8858,6 +8871,221 @@ mul_var(const NumericVar *var1, const Nu
result->sign = res_sign;
/* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
+}
+
+
+/*
+ * mul_var_small() -
+ *
+ * This has the same API as mul_var, but it assumes that var1 has no more
+ * than 4 digits and var2 has at least as many digits as var1. For variables
+ * satisfying these conditions, the product can be computed more quickly than
+ * the general algorithm used in mul_var.
+ */
+static void
+mul_var_small(const NumericVar *var1, const NumericVar *var2,
+ NumericVar *result, int rscale)
+{
+ int var1ndigits = var1->ndigits;
+ int var2ndigits = var2->ndigits;
+ NumericDigit *var1digits = var1->digits;
+ NumericDigit *var2digits = var2->digits;
+ int res_sign;
+ int res_weight;
+ int res_ndigits;
+ NumericDigit *res_buf;
+ NumericDigit *res_digits;
+ uint32 carry;
+ uint32 term;
+
+ /* Check preconditions */
+ Assert(var1ndigits <= 4);
+ Assert(var2ndigits >= var1ndigits);
+
+ /* Determine result sign and (maximum possible) weight */
+ if (var1->sign == var2->sign)
+ res_sign = NUMERIC_POS;
+ else
+ res_sign = NUMERIC_NEG;
+ res_weight = var1->weight + var2->weight + 2;
+
+ /* Determine the number of result digits to compute - cf. mul_var() */
+ res_ndigits = var1ndigits + var2ndigits + 1;
+
+ if (res_ndigits < 3)
+ {
+ /* All input digits will be ignored; so result is zero */
+ zero_var(result);
+ result->dscale = rscale;
+ return;
+ }
+
+ /* Allocate result digit array */
+ res_buf = digitbuf_alloc(res_ndigits);
+ res_buf[0] = 0; /* spare digit for later rounding */
+ res_digits = res_buf + 1;
+
+ /*
+ * Compute the result digits in reverse, in one pass, propagating the
+ * carry up as we go.
+ *
+ * This computes res_digits[res_ndigits - 2], ... res_digits[0] by summing
+ * the products var1digits[i1] * var2digits[i2] for which i1 + i2 + 1 is
+ * the result index.
+ */
+ switch (var1ndigits)
+ {
+ case 1:
+ /* ---------
+ * 1-digit case:
+ * var1ndigits = 1
+ * var2ndigits >= 1
+ * res_ndigits = var2ndigits + 2
+ * ----------
+ */
+ carry = 0;
+ for (int i = res_ndigits - 3; i >= 0; i--)
+ {
+ term = (uint32) var1digits[0] * var2digits[i] + carry;
+ res_digits[i + 1] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+ }
+ res_digits[0] = (NumericDigit) carry;
+ break;
+
+ case 2:
+ /* ---------
+ * 2-digit case:
+ * var1ndigits = 2
+ * var2ndigits >= 2
+ * res_ndigits = var2ndigits + 3
+ * ----------
+ */
+ /* last result digit and carry */
+ term = (uint32) var1digits[1] * var2digits[res_ndigits - 4];
+ res_digits[res_ndigits - 2] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ /* remaining digits, except for the first two */
+ for (int i = res_ndigits - 4; i >= 1; i--)
+ {
+ term = (uint32) var1digits[0] * var2digits[i] +
+ (uint32) var1digits[1] * var2digits[i - 1] + carry;
+ res_digits[i + 1] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+ }
+
+ /* first two digits */
+ term = (uint32) var1digits[0] * var2digits[0] + carry;
+ res_digits[1] = (NumericDigit) (term % NBASE);
+ res_digits[0] = (NumericDigit) (term / NBASE);
+ break;
+
+ case 3:
+ /* ---------
+ * 3-digit case:
+ * var1ndigits = 3
+ * var2ndigits >= 3
+ * res_ndigits = var2ndigits + 4
+ * ----------
+ */
+ /* last two result digits */
+ term = (uint32) var1digits[2] * var2digits[res_ndigits - 5];
+ res_digits[res_ndigits - 2] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[1] * var2digits[res_ndigits - 5] +
+ (uint32) var1digits[2] * var2digits[res_ndigits - 6] + carry;
+ res_digits[res_ndigits - 3] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ /* remaining digits, except for the first three */
+ for (int i = res_ndigits - 5; i >= 2; i--)
+ {
+ term = (uint32) var1digits[0] * var2digits[i] +
+ (uint32) var1digits[1] * var2digits[i - 1] +
+ (uint32) var1digits[2] * var2digits[i - 2] + carry;
+ res_digits[i + 1] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+ }
+
+ /* first three digits */
+ term = (uint32) var1digits[0] * var2digits[1] +
+ (uint32) var1digits[1] * var2digits[0] + carry;
+ res_digits[2] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[0] * var2digits[0] + carry;
+ res_digits[1] = (NumericDigit) (term % NBASE);
+ res_digits[0] = (NumericDigit) (term / NBASE);
+ break;
+
+ case 4:
+ /* ---------
+ * 4-digit case:
+ * var1ndigits = 4
+ * var2ndigits >= 4
+ * res_ndigits = var2ndigits + 5
+ * ----------
+ */
+ /* last three result digits */
+ term = (uint32) var1digits[3] * var2digits[res_ndigits - 6];
+ res_digits[res_ndigits - 2] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[2] * var2digits[res_ndigits - 6] +
+ (uint32) var1digits[3] * var2digits[res_ndigits - 7] + carry;
+ res_digits[res_ndigits - 3] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[1] * var2digits[res_ndigits - 6] +
+ (uint32) var1digits[2] * var2digits[res_ndigits - 7] +
+ (uint32) var1digits[3] * var2digits[res_ndigits - 8] + carry;
+ res_digits[res_ndigits - 4] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ /* remaining digits, except for the first four */
+ for (int i = res_ndigits - 6; i >= 3; i--)
+ {
+ term = (uint32) var1digits[0] * var2digits[i] +
+ (uint32) var1digits[1] * var2digits[i - 1] +
+ (uint32) var1digits[2] * var2digits[i - 2] +
+ (uint32) var1digits[3] * var2digits[i - 3] + carry;
+ res_digits[i + 1] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+ }
+
+ /* first four digits */
+ term = (uint32) var1digits[0] * var2digits[2] +
+ (uint32) var1digits[1] * var2digits[1] +
+ (uint32) var1digits[2] * var2digits[0] + carry;
+ res_digits[3] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[0] * var2digits[1] +
+ (uint32) var1digits[1] * var2digits[0] + carry;
+ res_digits[2] = (NumericDigit) (term % NBASE);
+ carry = term / NBASE;
+
+ term = (uint32) var1digits[0] * var2digits[0] + carry;
+ res_digits[1] = (NumericDigit) (term % NBASE);
+ res_digits[0] = (NumericDigit) (term / NBASE);
+ break;
+ }
+
+ /* Store the product in result (minus extra rounding digit) */
+ digitbuf_free(result->buf);
+ result->ndigits = res_ndigits - 1;
+ result->buf = res_buf;
+ result->digits = res_digits;
+ result->weight = res_weight - 1;
+ result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
round_var(result, rscale);
/* Strip leading and trailing zeroes */
