Subject: Re: Optimize numeric multiplication for one and two base-NBASE digit multiplicands.

From: "Joel Jacobson" <[email-masked]>
Date: Mon, 01 Jul 2024 15:14:58 +0200

On Mon, Jul 1, 2024, at 15:11, Joel Jacobson wrote:
> Not really sure why. Maybe the code I tried can be optimized further:

If anyone want experiment with the int128 version,
here is a patch that adds a separate numeric_mul_patched() function,
so it's easier to benchmark against the unmodified numeric_mul().

/Joel

---

From: "Joel Jacobson" <[email-masked]>
Date: Mon, 01 Jul 2024 21:56:23 +0200

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

---

From: Dean Rasheed <[email-masked]>
Date: Mon, 1 Jul 2024 23:19:36 +0100

On Mon, 1 Jul 2024 at 20:56, Joel Jacobson <[email-masked]> wrote:
>
> Below is a more realistic benchmark
>

I had a play with this, and came up with a slightly different way of
doing it that works for var2 of any size, as long as var1 is just 1 or
2 digits.

Repeating your benchmark where both numbers have up to 2 NBASE-digits,
this new approach was slightly faster:

SELECT SUM(num1*num2) FROM bench_mul_var; -- HEAD
Time: 4762.990 ms (00:04.763)
Time: 4332.166 ms (00:04.332)
Time: 4276.211 ms (00:04.276)
Time: 4247.321 ms (00:04.247)
Time: 4166.738 ms (00:04.167)

SELECT SUM(num1*num2) FROM bench_mul_var; -- v2 patch
Time: 4398.812 ms (00:04.399)
Time: 3672.668 ms (00:03.673)
Time: 3650.227 ms (00:03.650)
Time: 3611.420 ms (00:03.611)
Time: 3534.218 ms (00:03.534)

SELECT SUM(num1*num2) FROM bench_mul_var; -- this patch
Time: 3350.596 ms (00:03.351)
Time: 3336.224 ms (00:03.336)
Time: 3335.599 ms (00:03.336)
Time: 3336.990 ms (00:03.337)
Time: 3351.453 ms (00:03.351)

(This was on an older Intel Core i9-9900K, so I'm not sure why all the
timings are faster. What compiler settings are you using?)

The approach taken in this patch only uses 32-bit integers, so in
theory it could be extended to work for var1ndigits = 3, 4, or even
more, but the code would get increasingly complex, and I ran out of
steam at 2 digits. It might be worth trying though.

Regards,
Dean
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 5510a20..adbfd5c
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -8748,6 +8748,74 @@ mul_var(const NumericVar *var1, const Nu
 	}
 
 	/*
+	 * Simplified fast-path computation, if var1 has just one or two digits.
+	 * This is significantly faster, since it avoids allocating a separate
+	 * digit array, making multiple passes over var2, and having separate
+	 * carry-propagation passes.
+	 */
+	if (var1ndigits <= 2)
+	{
+		NumericDigit *res_buf;
+
+		/* 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 directly, in one pass, propagating the
+		 * carry up as we go.
+		 */
+		switch (var1ndigits)
+		{
+			case 1:
+				carry = 0;
+				for (i = res_ndigits - 3; i >= 0; i--)
+				{
+					newdig = (int) var1digits[0] * var2digits[i] + carry;
+					res_digits[i + 1] = (NumericDigit) (newdig % NBASE);
+					carry = newdig / NBASE;
+				}
+				res_digits[0] = (NumericDigit) carry;
+				break;
+
+			case 2:
+				newdig = (int) var1digits[1] * var2digits[res_ndigits - 4];
+				res_digits[res_ndigits - 2] = (NumericDigit) (newdig % NBASE);
+				carry = newdig / NBASE;
+
+				for (i = res_ndigits - 4; i > 0; i--)
+				{
+					newdig = (int) var1digits[0] * var2digits[i] +
+						(int) var1digits[1] * var2digits[i - 1] + carry;
+					res_digits[i + 1] = (NumericDigit) (newdig % NBASE);
+					carry = newdig / NBASE;
+				}
+
+				newdig = (int) var1digits[0] * var2digits[0] + carry;
+				res_digits[1] = (NumericDigit) (newdig % NBASE);
+				res_digits[0] = (NumericDigit) (newdig / 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 */
+		strip_var(result);
+
+		return;
+	}
+
+	/*
 	 * 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.

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 09:49:49 +0200

> this new approach was slightly faster:
...
> timings are faster. What compiler settings are you using?)

Strange. I just did `./configure` with a --prefix.

Compiler settings on my Intel Core i9-14900K machine:

$ pg_config | grep -E '^(CC|CFLAGS|CPPFLAGS|LDFLAGS)'
CC = gcc
CPPFLAGS = -D_GNU_SOURCE
CFLAGS = -Wall -Wmissing-prototypes -Wpointer-arith -Wdeclaration-after-statement -Werror=vla -Wendif-labels -Wmissing-format-attribute -Wimplicit-fallthrough=3 -Wcast-function-type -Wshadow=compatible-local -Wformat-security -fno-strict-aliasing -fwrapv -fexcess-precision=standard -Wno-format-truncation -Wno-stringop-truncation -O2
CFLAGS_SL = -fPIC
LDFLAGS = -Wl,--as-needed -Wl,-rpath,'/home/joel/pg-dev/lib',--enable-new-dtags
LDFLAGS_EX =
LDFLAGS_SL =

> Attachments:
> * optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.txt

Really nice!

I've benchmarked your patch on my three machines with great results.
I added a setseed() step, to make the benchmarks reproducible,
shouldn't matter much since it should statistically average out, but I thought why not.

SELECT setseed(0.12345);
\timing

/*
 * Apple M3 Max
 */

SELECT SUM(num1*num2) FROM bench_mul_var; -- HEAD
Time: 3622.342 ms (00:03.622)
Time: 3029.786 ms (00:03.030)
Time: 3046.497 ms (00:03.046)
Time: 3035.910 ms (00:03.036)
Time: 3034.073 ms (00:03.034)

SELECT SUM(num1*num2) FROM bench_mul_var; -- optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.txt
Time: 2484.685 ms (00:02.485)
Time: 2478.341 ms (00:02.478)
Time: 2494.397 ms (00:02.494)
Time: 2470.987 ms (00:02.471)
Time: 2490.215 ms (00:02.490)

/*
 * Intel Core i9-14900K
 */

SELECT SUM(num1*num2) FROM bench_mul_var; -- HEAD
Time: 2555.569 ms (00:02.556)
Time: 2523.145 ms (00:02.523)
Time: 2518.671 ms (00:02.519)
Time: 2514.501 ms (00:02.515)
Time: 2516.919 ms (00:02.517)

SELECT SUM(num1*num2) FROM bench_mul_var; -- optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.txt
Time: 2246.441 ms (00:02.246)
Time: 2243.900 ms (00:02.244)
Time: 2245.350 ms (00:02.245)
Time: 2245.080 ms (00:02.245)
Time: 2247.856 ms (00:02.248)

/*
 * AMD Ryzen 9 7950X3D
 */

SELECT SUM(num1*num2) FROM bench_mul_var; -- HEAD
Time: 3037.497 ms (00:03.037)
Time: 3010.037 ms (00:03.010)
Time: 3000.956 ms (00:03.001)
Time: 2989.424 ms (00:02.989)
Time: 2984.911 ms (00:02.985)

SELECT SUM(num1*num2) FROM bench_mul_var; -- optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.txt
Time: 2645.530 ms (00:02.646)
Time: 2640.472 ms (00:02.640)
Time: 2638.613 ms (00:02.639)
Time: 2637.889 ms (00:02.638)
Time: 2638.054 ms (00:02.638)

/Joel

---

From: Dean Rasheed <[email-masked]>
Date: Tue, 2 Jul 2024 09:22:15 +0100

> > * optimize-numeric-mul_var-small-var1-arbitrary-var2.patch.txt

Shortly after posting that, I realised that there was a small bug. This bit:

            case 2:
                newdig = (int) var1digits[1] * var2digits[res_ndigits - 4];

isn't quite right in the case where rscale is less than the full
result. In that case, the least significant result digit has a
contribution from one more var2 digit, so it needs to be:

                newdig = (int) var1digits[1] * var2digits[res_ndigits - 4];
                if (res_ndigits - 3 < var2ndigits)
                    newdig += (int) var1digits[0] * var2digits[res_ndigits - 3];

That doesn't noticeably affect the performance though. Update attached.

> shouldn't matter much since it should statistically average out, but I thought why not.

Nice. The results on the Apple M3 Max look particularly impressive.

I think it'd probably be worth trying to extend this to 3 or maybe 4
var1 digits, since that would cover a lot of "everyday" sized numeric
values that a lot of people might be using. I don't think we should go
beyond that though.

Regards,
Dean
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
new file mode 100644
index 5510a20..f77ad8b
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -8748,6 +8748,76 @@ mul_var(const NumericVar *var1, const Nu
 	}
 
 	/*
+				if (res_ndigits - 3 < var2ndigits)
+					newdig += (int) var1digits[0] * var2digits[res_ndigits - 3];

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 11:05:03 +0200

> That doesn't noticeably affect the performance though. Update attached.

Nice catch. Could we add a test somehow that would test mul_var()
with rscale less than the full result, that would catch bugs like this one
and others?

I created a new benchmark, that specifically tests different var1ndigits.
I've only run it on Apple M3 Max yet. More to come.

\timing
SELECT setseed(0.12345);
CREATE TABLE bench_mul_var_var1ndigits_1 (var1 numeric, var2 numeric);
INSERT INTO bench_mul_var_var1ndigits_1 (var1, var2)
SELECT random(0::numeric,9999::numeric), random(10000000::numeric,1e32::numeric) FROM generate_series(1,1e8);
CREATE TABLE bench_mul_var_var1ndigits_2 (var1 numeric, var2 numeric);
INSERT INTO bench_mul_var_var1ndigits_2 (var1, var2)
SELECT random(10000000::numeric,99999999::numeric), random(100000000000::numeric,1e32::numeric) FROM generate_series(1,1e8);
CREATE TABLE bench_mul_var_var1ndigits_3 (var1 numeric, var2 numeric);
INSERT INTO bench_mul_var_var1ndigits_3 (var1, var2)
SELECT random(100000000000::numeric,999999999999::numeric), random(1000000000000000::numeric,1e32::numeric) FROM generate_series(1,1e8);
CREATE TABLE bench_mul_var_var1ndigits_4 (var1 numeric, var2 numeric);
INSERT INTO bench_mul_var_var1ndigits_4 (var1, var2)
SELECT random(1000000000000000::numeric,9999999999999999::numeric), random(10000000000000000000::numeric,1e32::numeric) FROM generate_series(1,1e8);

/*
 * Apple M3 Max
 */

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_1; -- HEAD
Time: 2986.952 ms (00:02.987)
Time: 2991.765 ms (00:02.992)
Time: 2987.253 ms (00:02.987)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_1; -- v2-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 2874.841 ms (00:02.875)
Time: 2883.070 ms (00:02.883)
Time: 2899.973 ms (00:02.900)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_2; -- HEAD
Time: 3459.556 ms (00:03.460)
Time: 3304.983 ms (00:03.305)
Time: 3299.728 ms (00:03.300)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_2; -- v2-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3053.140 ms (00:03.053)
Time: 3065.227 ms (00:03.065)
Time: 3069.995 ms (00:03.070)

/*
 * Just for completeness, also testing var1ndigits 3 and 4,
 * although no change is expected since not yet implemented.
 */

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- HEAD
Time: 3809.005 ms (00:03.809)
Time: 3438.260 ms (00:03.438)
Time: 3453.920 ms (00:03.454)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- v2-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3437.592 ms (00:03.438)
Time: 3457.586 ms (00:03.458)
Time: 3474.344 ms (00:03.474)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_4; -- HEAD
Time: 4133.193 ms (00:04.133)
Time: 3554.477 ms (00:03.554)
Time: 3560.855 ms (00:03.561)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_4; -- v2-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3508.540 ms (00:03.509)
Time: 3566.721 ms (00:03.567)
Time: 3524.083 ms (00:03.524)

> values that a lot of people might be using. I don't think we should go

I think so too. I'm working on var1ndigits=3 now.

/Joel

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 12:23:36 +0200

>>                     newdig += (int) var1digits[0] * var2digits[res_ndigits - 3];

Just to be able to verify mul_var() is working as expected when
rscale is less than the full result, I've added a numeric_mul_patched()
function that takes a third rscale_adjustment parameter:

I've tried to get a different result with and without the fix,
but I'm failing to hit the bug.

Can you think of an example that should trigger the bug?

Example usage:

SELECT numeric_mul_patched(9999.999,2,0);
var1ndigits=1 var2ndigits=2 rscale=3
 numeric_mul_patched
---------------------
           19999.998
(1 row)

SELECT numeric_mul_patched(9999.999,2,1);
var1ndigits=1 var2ndigits=2 rscale=4
 numeric_mul_patched
---------------------
          19999.9980
(1 row)

SELECT numeric_mul_patched(9999.999,2,-1);
var1ndigits=1 var2ndigits=2 rscale=2
 numeric_mul_patched
---------------------
            20000.00
(1 row)

SELECT numeric_mul_patched(9999.999,2,-2);
var1ndigits=1 var2ndigits=2 rscale=1
 numeric_mul_patched
---------------------
             20000.0
(1 row)

/Joel

---

From: Dean Rasheed <[email-masked]>
Date: Tue, 2 Jul 2024 12:44:18 +0100

> function that takes a third rscale_adjustment parameter:

Yeah, we could expose such a function, and maybe it would be generally
useful, not just for testing, but that's really a separate proposal.
In general, mul_var() with reduced rscale doesn't guarantee correctly
rounded results though, which might make it less generally acceptable.

For this patch though, the aim is just to ensure the results are the
same as before.

> Can you think of an example that should trigger the bug?

9999.0001 * 5000.9999_9999 with rscale = 0 triggers it (returned
50004999 instead of 50005000).

Regards,
Dean

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 18:20:35 +0200

> 50004999 instead of 50005000).

Thanks, helpful.

Attached patch adds the var1ndigits=3 case.

Benchmark:

/*
 * Apple M3 Max
 */

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- v3-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3230.017 ms (00:03.230)
Time: 3244.094 ms (00:03.244)
Time: 3246.370 ms (00:03.246)

/*
 * Intel Core i9-14900K
 */

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- HEAD
Time: 4082.759 ms (00:04.083)
Time: 4077.372 ms (00:04.077)
Time: 4080.830 ms (00:04.081)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- v3-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3989.021 ms (00:03.989)
Time: 3978.263 ms (00:03.978)
Time: 3955.331 ms (00:03.955)

/*
 * AMD Ryzen 9 7950X3D
 */

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- HEAD
Time: 3791.271 ms (00:03.791)
Time: 3760.138 ms (00:03.760)
Time: 3758.773 ms (00:03.759)

SELECT SUM(var1*var2) FROM bench_mul_var_var1ndigits_3; -- v3-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch
Time: 3400.824 ms (00:03.401)
Time: 3373.714 ms (00:03.374)
Time: 3345.505 ms (00:03.346)

Regards,
Joel

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 20:53:23 +0200

On Tue, Jul 2, 2024, at 18:20, Joel Jacobson wrote:
> * v3-optimize-numeric-mul_var-small-var1-arbitrary-var2.patch

Hmm, v3 contains a bug which I haven't been able to solve yet.
Reporting now to avoid time waste reviewing it since it's buggy.

The attached patch is how I tested and found the bug.
It contains a file test-mul-var.sql, which tests mul_var()
with varying rscale, using the SQL-callable numeric_mul_patched(),
which third argument is the rscale_adjustment.

Out of 2481600 random tests, all passed except 4:

SELECT
     result = numeric_mul_patched(var1,var2,rscale_adjustment),
     COUNT(*)
FROM test_numeric_mul_patched
GROUP BY 1
ORDER BY 1;
 ?column? |  count
----------+---------
 f        |       4
 t        | 2481596
(2 rows)

SELECT
    var1,
    var2,
    var1*var2 AS full_resolution,
    rscale_adjustment,
    result AS expected,
    numeric_mul_patched(var1,var2,rscale_adjustment) AS actual
FROM test_numeric_mul_patched
WHERE result <> numeric_mul_patched(var1,var2,rscale_adjustment);
       var1        |      var2      |      full_resolution      | rscale_adjustment | expected | actual
-------------------+----------------+---------------------------+-------------------+----------+--------
     676.797214075 | 0.308068877759 | 208.500158210502929257925 |               -21 |      209 |    208
         555.07029 | 0.381033094735 |     211.50015039415392315 |               -17 |      212 |    211
    0.476637718921 |      66.088276 |     31.500165120061470196 |               -18 |       32 |     31
 0.060569165063082 |      998.85933 |   60.50007563356949425506 |               -20 |       61 |     60
(4 rows)

Trying to wrap my head around what could cause this.

It's rounding down instead of up, and these cases all end with decimal .500XXXX.

Regards,
Joel

From 074aeb223ab496f23c2075eabd35e6e76241d1d8 Mon Sep 17 00:00:00 2001
From: Joel Jakobsson <[email-masked]>
Date: Mon, 1 Jul 2024 07:17:50 +0200
Subject: [PATCH] Add SQL-callable numeric_mul_patched() to bench Simplified
 fast-path computation

---
 src/backend/utils/adt/numeric.c | 438 ++++++++++++++++++++++++++++++++
 src/include/catalog/pg_proc.dat |   3 +
 src/include/utils/numeric.h     |   2 +
 test-mul-var.sql                |  48 ++++
 4 files changed, 491 insertions(+)
 create mode 100644 test-mul-var.sql

diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 5510a203b0..8f5d553f15 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,131 @@ 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);
+	int32		rscale_adjustment = PG_GETARG_INT32(2);
+	Numeric		res;
+
+	res = numeric_mul_patched_opt_error(num1, num2, rscale_adjustment, 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, int32 rscale_adjustment, 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 + rscale_adjustment);
+
+	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 +8992,316 @@ 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;
+	}
+
+	/*
+	 * Simplified fast-path computation, if var1 has just one or two digits.
+	 * This is significantly faster, since it avoids allocating a separate
+	 * digit array, making multiple passes over var2, and having separate
+	 * carry-propagation passes.
+	 */
+	if (var1ndigits <= 3)
+	{
+		NumericDigit *res_buf;
+
+		/* 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 directly, in one pass, propagating the
+		 * carry up as we go.
+		 */
+		switch (var1ndigits)
+		{
+			case 1:
+				carry = 0;
+				for (i = res_ndigits - 3; i >= 0; i--)
+				{
+					newdig = (int) var1digits[0] * var2digits[i] + carry;
+					res_digits[i + 1] = (NumericDigit) (newdig % NBASE);
+					carry = newdig / NBASE;
+				}
+				res_digits[0] = (NumericDigit) carry;
+				break;
+
+			case 2:
+				newdig = (int) var1digits[1] * var2digits[res_ndigits - 4];
+				if (res_ndigits - 3 < var2ndigits)
+					newdig += (int) var1digits[0] * var2digits[res_ndigits - 3];
+				res_digits[res_ndigits - 2] = (NumericDigit) (newdig % NBASE);
+				carry = newdig / NBASE;
+				for (i = res_ndigits - 4; i >= 1; i--)
+				{
+					newdig = (int) var1digits[0] * var2digits[i] +
+						(int) var1digits[1] * var2digits[i - 1] + carry;
+					res_digits[i + 1] = (NumericDigit) (newdig % NBASE);
+					carry = newdig / NBASE;
+				}
+				newdig = (int) var1digits[0] * var2digits[0] + carry;
+				res_digits[1] = (NumericDigit) (newdig % NBASE);
+				res_digits[0] = (NumericDigit) (newdig / NBASE);
+				break;
+
+			case 3:
+				newdig = (int) var1digits[2] * var2digits[res_ndigits - 5];
+				if (res_ndigits - 4 < var2ndigits)
+					newdig += (int) var1digits[1] * var2digits[res_ndigits - 4];
+				if (res_ndigits - 3 < var2ndigits)
+					newdig += (int) var1digits[0] * var2digits[res_ndigits - 3];
+				res_digits[res_ndigits - 2] = (NumericDigit) (newdig % NBASE);
+				carry = newdig / NBASE;
+				for (i = res_ndigits - 4; i >= 2; i--)
+				{
+					newdig = carry;
+					if (i < var2ndigits)
+						newdig += (int) var1digits[0] * var2digits[i];
+					if (i - 1 >= 0 && i - 1 < var2ndigits)
+						newdig += (int) var1digits[1] * var2digits[i - 1];
+					if (i - 2 >= 0 && i - 2 < var2ndigits)
+						newdig += (int) var1digits[2] * var2digits[i - 2];
+					res_digits[i + 1] = (NumericDigit) (newdig % NBASE);
+					carry = newdig / NBASE;
+				}
+				newdig = carry;
+				if (var2ndigits > 1)
+					newdig += (int) var1digits[0] * var2digits[1];
+				if (var2ndigits > 0)
+					newdig += (int) var1digits[1] * var2digits[0];
+				res_digits[2] = (NumericDigit) (newdig % NBASE);
+				carry = newdig / NBASE;
+				newdig = (int) var1digits[0] * var2digits[0] + carry;
+				res_digits[1] = (NumericDigit) (newdig % NBASE);
+				res_digits[0] = (NumericDigit) (newdig / 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 */
+		strip_var(result);
+
+		return;
+	}
+
+	/*
+	 * 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 d4ac578ae6..5b3024cb6d 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 int4', 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..454a56da9a 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,
+									 int32 rscale_adjustment, 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,
diff --git a/test-mul-var.sql b/test-mul-var.sql
new file mode 100644
index 0000000000..ee7e3855bc
--- /dev/null
+++ b/test-mul-var.sql
@@ -0,0 +1,48 @@
+CREATE TABLE test_numeric_mul_patched (
+  var1 numeric,
+  var2 numeric,
+  rscale_adjustment int,
+  result numeric
+);
+
+DO $$
+DECLARE
+var1 numeric;
+var2 numeric;
+BEGIN
+  FOR i IN 1..100 LOOP
+    RAISE NOTICE '%', i;
+    FOR var1ndigits IN 1..4 LOOP
+      FOR var2ndigits IN 1..4 LOOP
+        FOR var1dscale IN 0..(var1ndigits*4) LOOP
+          FOR var2dscale IN 0..(var2ndigits*4) LOOP
+            FOR rscale_adjustment IN 0..(var1dscale+var2dscale) LOOP
+              var1 := round(random(
+                format('1%s',repeat('0',(var1ndigits-1)*4-1))::numeric,
+                format('%s',repeat('9',var1ndigits*4))::numeric
+              ) / 10::numeric^var1dscale, var1dscale);
+              var2 := round(random(
+                format('1%s',repeat('0',(var2ndigits-1)*4-1))::numeric,
+                format('%s',repeat('9',var2ndigits*4))::numeric
+              ) / 10::numeric^var2dscale, var2dscale);
+              INSERT INTO test_numeric_mul_patched
+                (var1, var2, rscale_adjustment)
+              VALUES
+                (var1, var2, -rscale_adjustment);
+            END LOOP;
+          END LOOP;
+        END LOOP;
+      END LOOP;
+    END LOOP;
+  END LOOP;
+END $$;
+
+-- First, set result with a numeric_mul_patched() version where
+-- the Simplified fast-path computation code has been commented out.
+UPDATE test_numeric_mul_patched SET result = numeric_mul_patched(var1, var2, rscale_adjustment);
+
+-- Then, recompile with the Simplified fast-path computation code,
+-- and check if any differences can be found:
+SELECT *, numeric_mul_patched(var1,var2,rscale_adjustment)
+FROM test_numeric_mul_patched
+WHERE result IS DISTINCT FROM numeric_mul_patched(var1,var2,rscale_adjustment);
-- 
2.45.1

---

From: "Joel Jacobson" <[email-masked]>
Date: Tue, 02 Jul 2024 21:55:19 +0200

> It's rounding down instead of up, and these cases all end with decimal .500XXXX.

Interesting, I actually think there is a bug in the normal mul_var() code.
Found a case that rounds down, when it should round up:

Calling mul_var() with:
var1=51.2945442386599
var2=0.828548712212
rscale=0

returns 42, but I think it should return 43,
since 51.2945442386599*0.828548712212=42.5000285724431241296446988

But maybe this is expected and OK, having to do with MUL_GUARD_DIGITS?

/Joel

---

From: Dean Rasheed <[email-masked]>
Date: Wed, 3 Jul 2024 11:41:58 +0100

> But maybe this is expected and OK, having to do with MUL_GUARD_DIGITS?

No, that's not a bug. It's to be expected. The point of using only
MUL_GUARD_DIGITS extra digits is to get the correctly rounded result
*most of the time*, while saving a lot of effort by not computing the
full result.

The callers of mul_var() that ask for reduced rscale results are the
transcendental functions like ln_var() and exp_var(), which are
working to some limited precision intended to compute the final result
reasonably accurately to a particular rscale.

Regards,
Dean