public inbox for [email protected]
help / color / mirror / Atom feedFrom: Joel Jacobson <[email protected]>
To: Dagfinn Ilmari Mannsåker <[email protected]>
Cc: pgsql-hackers <[email protected]>
Subject: Re: Optimize numeric multiplication for one and two base-NBASE digit multiplicands.
Date: Mon, 01 Jul 2024 21:56:23 +0200
Message-ID: <[email protected]> (raw)
In-Reply-To: <[email protected]>
References: <[email protected]>
<[email protected]>
<[email protected]>
<[email protected]>
On Mon, Jul 1, 2024, at 15:14, Joel Jacobson wrote:
> * 0001-Optimize-mul_var-for-var2ndigits-4.patch
Found a typo, fixed in new version.
The int128 version is still slower though,
I wonder if there is something that can be done to speed it up further.
Below is a more realistic benchmark than just microbenchmarking mul_var(),
not testing the int128 version, but the code for up to 2 NBASE-digits:
```
CREATE TABLE bench_mul_var (num1 numeric, num2 numeric);
INSERT INTO bench_mul_var (num1, num2)
SELECT random(0::numeric,1e8::numeric), random(0::numeric,1e8::numeric) FROM generate_series(1,1e8);
SELECT SUM(num1*num2) FROM bench_mul_var;
Time: 8331.953 ms (00:08.332)
Time: 7415.241 ms (00:07.415)
Time: 7298.296 ms (00:07.298)
Time: 7314.754 ms (00:07.315)
Time: 7289.560 ms (00:07.290)
SELECT SUM(numeric_mul_patched(num1,num2)) FROM bench_mul_var;
Time: 6403.426 ms (00:06.403)
Time: 6401.797 ms (00:06.402)
Time: 6366.136 ms (00:06.366)
Time: 6376.049 ms (00:06.376)
Time: 6317.282 ms (00:06.317)
``
Benchmarked on a Intel Core i9-14900K machine.
/Joel
Attachments:
[application/octet-stream] v2-0001-Optimize-mul_var-for-var2ndigits-4.patch (14.7K, ../[email protected]/2-v2-0001-Optimize-mul_var-for-var2ndigits-4.patch)
download | inline diff:
From b075e18c2ae94f2ada20572c05ff35730d2c5968 Mon Sep 17 00:00:00 2001
From: Joel Jakobsson <[email protected]>
Date: Mon, 1 Jul 2024 07:17:50 +0200
Subject: [PATCH] Optimize mul_var() for var2ndigits <= 4
---
src/backend/utils/adt/numeric.c | 457 ++++++++++++++++++++++++++++++++
src/include/catalog/pg_proc.dat | 3 +
src/include/utils/numeric.h | 2 +
3 files changed, 462 insertions(+)
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 5510a203b0..b74e03eefa 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -551,6 +551,9 @@ static void sub_var(const NumericVar *var1, const NumericVar *var2,
static void mul_var(const NumericVar *var1, const NumericVar *var2,
NumericVar *result,
int rscale);
+static void mul_var_patched(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);
@@ -3115,6 +3118,130 @@ numeric_mul_opt_error(Numeric num1, Numeric num2, bool *have_error)
}
+/*
+ * numeric_mul_patched() -
+ *
+ * This function multiplies two numeric values using the patched algorithm,
+ * designed for efficient handling of large numbers. It's introduced to allow
+ * direct benchmark comparisons with the standard numeric_mul() function.
+ */
+Datum
+numeric_mul_patched(PG_FUNCTION_ARGS)
+{
+ Numeric num1 = PG_GETARG_NUMERIC(0);
+ Numeric num2 = PG_GETARG_NUMERIC(1);
+ Numeric res;
+
+ res = numeric_mul_patched_opt_error(num1, num2, NULL);
+
+ PG_RETURN_NUMERIC(res);
+}
+
+
+/*
+ * numeric_mul_patched_opt_error() -
+ *
+ * Internal version of numeric_mul_patched().
+ * If "*have_error" flag is provided, on error it's set to true, NULL returned.
+ * This is helpful when caller need to handle errors by itself.
+ */
+Numeric
+numeric_mul_patched_opt_error(Numeric num1, Numeric num2, bool *have_error)
+{
+ NumericVar arg1;
+ NumericVar arg2;
+ NumericVar result;
+ Numeric res;
+
+ /*
+ * Handle NaN and infinities
+ */
+ if (NUMERIC_IS_SPECIAL(num1) || NUMERIC_IS_SPECIAL(num2))
+ {
+ if (NUMERIC_IS_NAN(num1) || NUMERIC_IS_NAN(num2))
+ return make_result(&const_nan);
+ if (NUMERIC_IS_PINF(num1))
+ {
+ switch (numeric_sign_internal(num2))
+ {
+ case 0:
+ return make_result(&const_nan); /* Inf * 0 */
+ case 1:
+ return make_result(&const_pinf);
+ case -1:
+ return make_result(&const_ninf);
+ }
+ Assert(false);
+ }
+ if (NUMERIC_IS_NINF(num1))
+ {
+ switch (numeric_sign_internal(num2))
+ {
+ case 0:
+ return make_result(&const_nan); /* -Inf * 0 */
+ case 1:
+ return make_result(&const_ninf);
+ case -1:
+ return make_result(&const_pinf);
+ }
+ Assert(false);
+ }
+ /* by here, num1 must be finite, so num2 is not */
+ if (NUMERIC_IS_PINF(num2))
+ {
+ switch (numeric_sign_internal(num1))
+ {
+ case 0:
+ return make_result(&const_nan); /* 0 * Inf */
+ case 1:
+ return make_result(&const_pinf);
+ case -1:
+ return make_result(&const_ninf);
+ }
+ Assert(false);
+ }
+ Assert(NUMERIC_IS_NINF(num2));
+ switch (numeric_sign_internal(num1))
+ {
+ case 0:
+ return make_result(&const_nan); /* 0 * -Inf */
+ case 1:
+ return make_result(&const_ninf);
+ case -1:
+ return make_result(&const_pinf);
+ }
+ Assert(false);
+ }
+
+ /*
+ * Unpack the values, let mul_var() compute the result and return it.
+ * Unlike add_var() and sub_var(), mul_var() will round its result. In the
+ * case of numeric_mul(), which is invoked for the * operator on numerics,
+ * we request exact representation for the product (rscale = sum(dscale of
+ * arg1, dscale of arg2)). If the exact result has more digits after the
+ * decimal point than can be stored in a numeric, we round it. Rounding
+ * after computing the exact result ensures that the final result is
+ * correctly rounded (rounding in mul_var() using a truncated product
+ * would not guarantee this).
+ */
+ init_var_from_num(num1, &arg1);
+ init_var_from_num(num2, &arg2);
+
+ init_var(&result);
+
+ mul_var_patched(&arg1, &arg2, &result, arg1.dscale + arg2.dscale);
+
+ if (result.dscale > NUMERIC_DSCALE_MAX)
+ round_var(&result, NUMERIC_DSCALE_MAX);
+
+ res = make_result_opt_error(&result, have_error);
+
+ free_var(&result);
+
+ return res;
+}
+
+
/*
* numeric_div() -
*
@@ -8864,6 +8991,336 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result,
strip_var(result);
}
+/*
+ * mul_var_patched() -
+ *
+ * Implements patched multiplication for large numbers, introduced
+ * alongside the unchanged original mul_var(). This function is part of
+ * an optimization effort, allowing direct benchmark comparisons with
+ * mul_var(). It selects full or half patched based on input size.
+ * This is a temporary measure before considering its replacement of
+ * mul_var() based on benchmark outcomes.
+ */
+static void
+mul_var_patched(const NumericVar *var1, const NumericVar *var2,
+ NumericVar *result, int rscale)
+{
+ int res_ndigits;
+ int res_sign;
+ int res_weight;
+ int maxdigits;
+ int *dig;
+ int carry;
+ int maxdig;
+ int newdig;
+ int var1ndigits;
+ int var2ndigits;
+ NumericDigit *var1digits;
+ NumericDigit *var2digits;
+ NumericDigit *res_digits;
+ int i,
+ i1,
+ i2;
+
+ /*
+ * Arrange for var1 to be the shorter of the two numbers. This improves
+ * performance because the inner multiplication loop is much simpler than
+ * the outer loop, so it's better to have a smaller number of iterations
+ * of the outer loop. This also reduces the number of times that the
+ * accumulator array needs to be normalized.
+ */
+ if (var1->ndigits > var2->ndigits)
+ {
+ const NumericVar *tmp = var1;
+
+ var1 = var2;
+ var2 = tmp;
+ }
+
+ /* copy these values into local vars for speed in inner loop */
+ var1ndigits = var1->ndigits;
+ var2ndigits = var2->ndigits;
+ var1digits = var1->digits;
+ var2digits = var2->digits;
+
+ if (var1ndigits == 0 || var2ndigits == 0)
+ {
+ /* one or both inputs is zero; so is result */
+ zero_var(result);
+ result->dscale = rscale;
+ return;
+ }
+
+ /* 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. If the exact result
+ * would have more than rscale fractional digits, truncate the computation
+ * with MUL_GUARD_DIGITS guard digits, i.e., ignore input digits that
+ * would only contribute to the right of that. (This will give the exact
+ * rounded-to-rscale answer unless carries out of the ignored positions
+ * would have propagated through more than MUL_GUARD_DIGITS digits.)
+ *
+ * Note: an exact computation could not produce more than var1ndigits +
+ * var2ndigits digits, but we allocate one extra output digit in case
+ * rscale-driven rounding produces a carry out of the highest exact digit.
+ */
+ res_ndigits = var1ndigits + var2ndigits + 1;
+ maxdigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS +
+ MUL_GUARD_DIGITS;
+ res_ndigits = Min(res_ndigits, maxdigits);
+
+ if (res_ndigits < 3)
+ {
+ /* All input digits will be ignored; so result is zero */
+ zero_var(result);
+ result->dscale = rscale;
+ return;
+ }
+
+ /*
+ * If var1 and var2 are just one or two NBASE digits, their product will
+ * fit in an int64 can be computed directly, which is significantly faster.
+ */
+ if (var2ndigits <= 2)
+ {
+ int64 product = 0;
+ switch (var1ndigits)
+ {
+ case 1:
+ product = var1digits[0];
+ break;
+ case 2:
+ product = var1digits[0] * NBASE + var1digits[1];
+ break;
+ }
+
+ switch (var2ndigits)
+ {
+ case 1:
+ product *= var2digits[0];
+ break;
+ case 2:
+ product *= var2digits[0] * NBASE + var2digits[1];
+ break;
+ }
+
+ alloc_var(result, res_ndigits);
+ res_digits = result->digits;
+ for (i = res_ndigits - 1; i >= 0; i--)
+ {
+ res_digits[i] = product % NBASE;
+ product /= NBASE;
+ }
+ Assert(product == 0);
+
+ /*
+ * Finally, round the result to the requested precision.
+ */
+ result->weight = res_weight;
+ result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
+
+ return;
+ }
+#ifdef HAVE_INT128
+ /*
+ * If var1 and var2 are up to four digits, their product will fit in
+ * an int128 can be computed directly, which is significantly faster.
+ */
+ if (var2ndigits <= 4)
+ {
+ int128 product = 0;
+
+ switch (var1ndigits)
+ {
+ case 1:
+ product = var1digits[0];
+ break;
+ case 2:
+ product = var1digits[0] * NBASE + var1digits[1];
+ break;
+ case 3:
+ product = ((int128) var1digits[0] * NBASE + var1digits[1])
+ * NBASE + var1digits[2];
+ break;
+ case 4:
+ product = (((int128) var1digits[0] * NBASE + var1digits[1])
+ * NBASE + var1digits[2]) * NBASE + var1digits[3];
+ break;
+ }
+
+ switch (var2ndigits)
+ {
+ case 1:
+ product *= var2digits[0];
+ break;
+ case 2:
+ product *= var2digits[0] * NBASE + var2digits[1];
+ break;
+ case 3:
+ product *= ((int128) var2digits[0] * NBASE + var2digits[1])
+ * NBASE + var2digits[2];
+ break;
+ case 4:
+ product *= (((int128) var2digits[0] * NBASE + var2digits[1])
+ * NBASE + var2digits[2]) * NBASE + var2digits[3];
+ break;
+ }
+
+ alloc_var(result, res_ndigits);
+ res_digits = result->digits;
+ for (i = res_ndigits - 1; i >= 0; i--)
+ {
+ res_digits[i] = product % NBASE;
+ product /= NBASE;
+ }
+ Assert(product == 0);
+
+ /*
+ * Finally, round the result to the requested precision.
+ */
+ result->weight = res_weight;
+ result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
+
+ return;
+ }
+#endif
+
+ /*
+ * We do the arithmetic in an array "dig[]" of signed int's. Since
+ * INT_MAX is noticeably larger than NBASE*NBASE, this gives us headroom
+ * to avoid normalizing carries immediately.
+ *
+ * maxdig tracks the maximum possible value of any dig[] entry; when this
+ * threatens to exceed INT_MAX, we take the time to propagate carries.
+ * Furthermore, we need to ensure that overflow doesn't occur during the
+ * carry propagation passes either. The carry values could be as much as
+ * INT_MAX/NBASE, so really we must normalize when digits threaten to
+ * exceed INT_MAX - INT_MAX/NBASE.
+ *
+ * To avoid overflow in maxdig itself, it actually represents the max
+ * possible value divided by NBASE-1, ie, at the top of the loop it is
+ * known that no dig[] entry exceeds maxdig * (NBASE-1).
+ */
+ dig = (int *) palloc0(res_ndigits * sizeof(int));
+ maxdig = 0;
+
+ /*
+ * The least significant digits of var1 should be ignored if they don't
+ * contribute directly to the first res_ndigits digits of the result that
+ * we are computing.
+ *
+ * Digit i1 of var1 and digit i2 of var2 are multiplied and added to digit
+ * i1+i2+2 of the accumulator array, so we need only consider digits of
+ * var1 for which i1 <= res_ndigits - 3.
+ */
+ for (i1 = Min(var1ndigits - 1, res_ndigits - 3); i1 >= 0; i1--)
+ {
+ NumericDigit var1digit = var1digits[i1];
+
+ if (var1digit == 0)
+ continue;
+
+ /* Time to normalize? */
+ maxdig += var1digit;
+ if (maxdig > (INT_MAX - INT_MAX / NBASE) / (NBASE - 1))
+ {
+ /* Yes, do it */
+ carry = 0;
+ for (i = res_ndigits - 1; i >= 0; i--)
+ {
+ newdig = dig[i] + carry;
+ if (newdig >= NBASE)
+ {
+ carry = newdig / NBASE;
+ newdig -= carry * NBASE;
+ }
+ else
+ carry = 0;
+ dig[i] = newdig;
+ }
+ Assert(carry == 0);
+ /* Reset maxdig to indicate new worst-case */
+ maxdig = 1 + var1digit;
+ }
+
+ /*
+ * Add the appropriate multiple of var2 into the accumulator.
+ *
+ * As above, digits of var2 can be ignored if they don't contribute,
+ * so we only include digits for which i1+i2+2 < res_ndigits.
+ *
+ * This inner loop is the performance bottleneck for multiplication,
+ * so we want to keep it simple enough so that it can be
+ * auto-vectorized. Accordingly, process the digits left-to-right
+ * even though schoolbook multiplication would suggest right-to-left.
+ * Since we aren't propagating carries in this loop, the order does
+ * not matter.
+ */
+ {
+ int i2limit = Min(var2ndigits, res_ndigits - i1 - 2);
+ int *dig_i1_2 = &dig[i1 + 2];
+
+ for (i2 = 0; i2 < i2limit; i2++)
+ dig_i1_2[i2] += var1digit * var2digits[i2];
+ }
+ }
+
+ /*
+ * Now we do a final carry propagation pass to normalize the result, which
+ * we combine with storing the result digits into the output. Note that
+ * this is still done at full precision w/guard digits.
+ */
+ alloc_var(result, res_ndigits);
+ res_digits = result->digits;
+ carry = 0;
+ for (i = res_ndigits - 1; i >= 0; i--)
+ {
+ newdig = dig[i] + carry;
+ if (newdig >= NBASE)
+ {
+ carry = newdig / NBASE;
+ newdig -= carry * NBASE;
+ }
+ else
+ carry = 0;
+ res_digits[i] = newdig;
+ }
+ Assert(carry == 0);
+
+ pfree(dig);
+
+ /*
+ * Finally, round the result to the requested precision.
+ */
+ result->weight = res_weight;
+ result->sign = res_sign;
+
+ /* Round to target rscale (and set result->dscale) */
+ round_var(result, rscale);
+
+ /* Strip leading and trailing zeroes */
+ strip_var(result);
+
+}
+
/*
* div_var() -
diff --git a/src/include/catalog/pg_proc.dat b/src/include/catalog/pg_proc.dat
index 6a5476d3c4..15dde0f3c3 100644
--- a/src/include/catalog/pg_proc.dat
+++ b/src/include/catalog/pg_proc.dat
@@ -4465,6 +4465,9 @@
{ oid => '1726',
proname => 'numeric_mul', prorettype => 'numeric',
proargtypes => 'numeric numeric', prosrc => 'numeric_mul' },
+{ oid => '6347',
+ proname => 'numeric_mul_patched', prorettype => 'numeric',
+ proargtypes => 'numeric numeric', prosrc => 'numeric_mul_patched' },
{ oid => '1727',
proname => 'numeric_div', prorettype => 'numeric',
proargtypes => 'numeric numeric', prosrc => 'numeric_div' },
diff --git a/src/include/utils/numeric.h b/src/include/utils/numeric.h
index 43c75c436f..9036c9db50 100644
--- a/src/include/utils/numeric.h
+++ b/src/include/utils/numeric.h
@@ -97,6 +97,8 @@ extern Numeric numeric_sub_opt_error(Numeric num1, Numeric num2,
bool *have_error);
extern Numeric numeric_mul_opt_error(Numeric num1, Numeric num2,
bool *have_error);
+extern Numeric numeric_mul_patched_opt_error(Numeric num1, Numeric num2,
+ bool *have_error);
extern Numeric numeric_div_opt_error(Numeric num1, Numeric num2,
bool *have_error);
extern Numeric numeric_mod_opt_error(Numeric num1, Numeric num2,
--
2.45.1
view thread (13+ messages) latest in thread
reply
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Reply to all the recipients using the --to and --cc options:
reply via email
To: [email protected]
Cc: [email protected], [email protected]
Subject: Re: Optimize numeric multiplication for one and two base-NBASE digit multiplicands.
In-Reply-To: <[email protected]>
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
This inbox is served by agora; see mirroring instructions
for how to clone and mirror all data and code used for this inbox