From 9d22b244e257d2e4cccc321b7d5ed6d90f5ea3a4 Mon Sep 17 00:00:00 2001 From: Dean Rasheed Date: Thu, 18 Jul 2024 18:32:56 +0100 Subject: [PATCH v3 2/2] Optimise numeric multiplication using base-NBASE^2 arithmetic. Currently mul_var() uses the schoolbook multiplication algorithm, which is O(n^2) in the number of NBASE digits. To improve performance for large inputs, convert the inputs to base NBASE^2 before multiplying, which effectively halves the number of digits in each input, theoretically speeding up the computation by a factor of 4. In practice, the actual speedup for large inputs varies between around 3 and 6 times, depending on the system and compiler used. In turn, this significantly reduces the runtime of the numeric_big regression test. For this to work, 64-bit integers are required for the products of base-NBASE^2 digits, so this works best on 64-bit machines, for which it is faster whenever the shorter input has more than 4 or 5 NBASE digits. On 32-bit machines, the additional overheads, especially during carry propagation and the final conversion back to base-NBASE, are significantly higher, and it is only faster when the shorter input has more than around 30 NBASE digits. Therefore, only use this approach above a platform-dependent threshold. For inputs below the threshold, the original base-NBASE algorithm is used, except that it can be simplified because the threshold is low enough that intermediate carry-propagation passes are not required. Above the threshold, the available headroom in 64-bit integers is much larger than for 32-bit integers, so the frequency of carry-propagation passes is greatly reduced. In addition, unsigned integers are used throughout, further increasing the headroom. --- src/backend/utils/adt/numeric.c | 512 +++++++++++++++++++++++++------- 1 file changed, 401 insertions(+), 111 deletions(-) diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c index 9b9b88662a..c463901428 100644 --- a/src/backend/utils/adt/numeric.c +++ b/src/backend/utils/adt/numeric.c @@ -101,6 +101,29 @@ typedef signed char NumericDigit; typedef int16 NumericDigit; #endif +/* + * Above a certain size threshold, it is faster to multiply numbers by + * converting them to base NBASE^2, and use 64-bit integer arithmetic. This + * threshold is determined empirically, and is necessarily higher on 32-bit + * machines, which are less efficient at performing 64-bit arithmetic. + * + * To simplify the computation below this threshold, it intentially kept below + * the point at which intermediate carry-propagation passes may be necessary. + * Therefore, as explained in mul_var(), it must be no larger than + * (PG_UINT32_MAX - PG_UINT32_MAX / NBASE) / (NBASE - 1)^2, which is 42 when + * NBASE is 10000. + */ +#define MUL_64BIT_THRESHOLD_MAX \ + ((PG_UINT32_MAX - PG_UINT32_MAX / NBASE) / (NBASE - 1) / (NBASE - 1)) + +#if SIZEOF_DATUM < 8 +#define MUL_64BIT_THRESHOLD 30 +#else +#define MUL_64BIT_THRESHOLD 4 +#endif + +#define NBASE_SQR (NBASE * NBASE) + /* * The Numeric type as stored on disk. * @@ -8663,6 +8686,85 @@ sub_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result) } +/* + * div_mod_NBASE_SQR() - + * + * Divide a 64-bit integer "num" by NBASE_SQR, returning the quotient and + * remainder. Technically, the remainder could be a 32-bit integer, but the + * caller actually wants a 64-bit integer, so this is more efficient. + */ +static inline uint64 +div_mod_NBASE_SQR(uint64 num, uint64 *rem) +{ + uint64 quot; + + /* ---------- + * On a 32-bit machine, the compiler does 64-bit division using a builtin + * function such as __udivdi3(), which is very slow. Replace that with a + * multiply-and-shift algorithm, based on the way compilers do it on + * 64-bit machines. Assuming that DEC_DIGITS is 4, and NBASE_SQR = 10^8, + * the multiply-and-shift formula is + * + * quot = num / 10^8 = (num * multiplier) >> 90 + * + * where multiplier = ceil(2^90 / 10^8) = 12379400392853802749. + * + * The 2^90 scaling factor here guarantees correct results for all inputs. + * See "Division by Invariant Integers using Multiplication", Torbjorn + * Granlund and Peter L. Montgomery, PLDI '94: Proceedings of the ACM + * SIGPLAN 1994 conference on Programming language design and + * implementation (https://dl.acm.org/doi/pdf/10.1145/178243.178249). + * + * Since num and multiplier are 64-bit unsigned integers, their product is + * a 128-bit unsigned integer, but we only require the high part. This is + * done by decomposing num and multiplier into high and low 32-bit parts, + * and then computing the upper 64 bits of the full product: + * + * num * multiplier = + * (num_hi * multiplier_hi) << 64 + + * (num_hi * multiplier_lo + num_lo * multiplier_hi) << 32 + + * num_lo * multiplier_lo + * + * We don't bother with this optimization for other NBASE values. + * ---------- + */ +#if SIZEOF_DATUM < 8 && DEC_DIGITS == 4 + const uint64 multiplier = UINT64CONST(12379400392853802749); + + /* high and low 32-bit parts of num and multiplier */ +#define UINT64_HI32(x) ((uint32) ((x) >> 32)) +#define UINT64_LO32(x) ((uint32) (x)) + const uint32 multiplier_hi = UINT64_HI32(multiplier); + const uint32 multiplier_lo = UINT64_LO32(multiplier); + uint32 num_hi = UINT64_HI32(num); + uint32 num_lo = UINT64_LO32(num); + uint64 tmp1, + tmp2, + prod_hi; + + /* high 64-bit part of 128-bit product */ + tmp1 = (uint64) num_hi * multiplier_lo; + tmp2 = (uint64) num_lo * multiplier_lo; + tmp1 += UINT64_HI32(tmp2); + prod_hi = (uint64) num_hi * multiplier_hi; + prod_hi += UINT64_HI32(tmp1); + tmp2 = (uint64) num_lo * multiplier_hi; + tmp2 += UINT64_LO32(tmp1); + prod_hi += UINT64_HI32(tmp2); + + /* quotient is in the top 38 bits */ + quot = prod_hi >> 26; +#else + /* just divide normally */ + quot = num / NBASE_SQR; +#endif + /* remainder */ + *rem = num - quot * NBASE_SQR; + + return quot; +} + + /* * mul_var() - * @@ -8677,10 +8779,6 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result, int res_sign; int res_weight; int maxdigits; - int *dig; - int carry; - int maxdig; - int newdig; int var1ndigits; int var2ndigits; NumericDigit *var1digits; @@ -8688,7 +8786,8 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result, NumericDigit *res_digits; int i, i1, - i2; + i2, + i2limit; /* * Arrange for var1 to be the shorter of the two numbers. This improves @@ -8729,141 +8828,332 @@ mul_var(const NumericVar *var1, const NumericVar *var2, NumericVar *result, return; } - /* Determine result sign and (maximum possible) weight */ + /* Determine result sign */ 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.) + * We do the arithmetic in an array "dig[]" of unsigned 32-bit or 64-bit + * integers, depending on the size of var1. * - * 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 var1 has more than MUL_64BIT_THRESHOLD digits, we convert the inputs + * to base NBASE^2 and multiply using 64-bit integer arithmetic, which is + * much faster, since schoolbook multiplication is O(N^2) in the number of + * input digits, and working in base NBASE^2 effectively halves "N". + * + * Below this threshold, we work with the original base-NBASE numbers, and + * use 32-bit integer arithmetic. To simplify the algorithm, we ensure + * that the threshold is low enough so that the number of products being + * added to any element of dig[] is small enough to avoid integer + * overflow. Furthermore, we need to ensure that overflow doesn't occur + * during the final carry-propagation pass. The carry values can be as + * large as PG_UINT32_MAX / NBASE, and so the values in dig[] must not + * exceed PG_UINT32_MAX - PG_UINT32_MAX / NBASE. Since each product of + * digits is at most (NBASE - 1)^2, the number of products must not exceed + * (PG_UINT32_MAX - PG_UINT32_MAX / NBASE) / (NBASE - 1)^2. + */ + StaticAssertStmt(MUL_64BIT_THRESHOLD <= MUL_64BIT_THRESHOLD_MAX, + "MUL_64BIT_THRESHOLD must not exceed MUL_64BIT_THRESHOLD_MAX"); + + if (var1ndigits <= MUL_64BIT_THRESHOLD) + { + uint32 *dig, + *dig_i1_2; + NumericDigit var1digit; + uint32 carry; + uint32 newdig; - if (res_ndigits < 3) - { - /* All input digits will be ignored; so result is zero */ - zero_var(result); - result->dscale = rscale; - return; - } + /* + * Determine the number of result digits to compute and the (maximum + * possible) result weight. 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 at least one extra output + * digit in case rscale-driven rounding produces a carry out of the + * highest exact digit. + */ + res_ndigits = var1ndigits + var2ndigits + 1; + res_weight = var1->weight + var2->weight + 2; - /* - * 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; + maxdigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS + + MUL_GUARD_DIGITS; + res_ndigits = Min(res_ndigits, maxdigits); - /* - * 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 (res_ndigits < 3) + { + /* All input digits will be ignored; so result is zero */ + zero_var(result); + result->dscale = rscale; + return; + } - if (var1digit == 0) - continue; + /* Allocate dig[] to accumulate the digit products */ + dig = (uint32 *) palloc(res_ndigits * sizeof(uint32)); + + /* + * Start by multiplying var2 by the least significant contributing + * digit of var1, storing the results at the end of dig[], and filling + * the leading slots with zeros. + * + * 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. + * + * The loop here is the same as the inner loop below, except that we + * set the results in dig[], rather than adding to them. This 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. + */ + i1 = Min(var1ndigits - 1, res_ndigits - 3); + var1digit = var1digits[i1]; + + i2limit = Min(var2ndigits, res_ndigits - i1 - 2); + dig_i1_2 = &dig[i1 + 2]; + + memset(dig, 0, (i1 + 2) * sizeof(uint32)); + for (i2 = 0; i2 < i2limit; i2++) + dig_i1_2[i2] = var1digit * var2digits[i2]; - /* Time to normalize? */ - maxdig += var1digit; - if (maxdig > (INT_MAX - INT_MAX / NBASE) / (NBASE - 1)) + /* + * Next, multiply var2 by the remaining digits of var1, adding the + * results to dig[] at the appropriate offsets. + */ + for (i1 = i1 - 1; i1 >= 0; i1--) { - /* Yes, do it */ - carry = 0; - for (i = res_ndigits - 1; i >= 0; i--) + var1digit = var1digits[i1]; + if (var1digit != 0) { - newdig = dig[i] + carry; - if (newdig >= NBASE) - { - carry = newdig / NBASE; - newdig -= carry * NBASE; - } - else - carry = 0; - dig[i] = newdig; + i2limit = Min(var2ndigits, res_ndigits - i1 - 2); + dig_i1_2 = &dig[i1 + 2]; + + for (i2 = 0; i2 < i2limit; i2++) + dig_i1_2[i2] += var1digit * var2digits[i2]; } - Assert(carry == 0); - /* Reset maxdig to indicate new worst-case */ - maxdig = 1 + var1digit; } /* - * Add the appropriate multiple of var2 into the accumulator. + * Finally, construct the result digits by propagating carries up, + * normalizing back to base-NBASE. 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); + } + else + { + int var1ndigitpairs; + int var2ndigitpairs; + int res_ndigitpairs; + int pair_offset; + int maxdigitpairs; + uint64 *dig, + *dig_i1_off; + uint32 *var2digitpairs; + uint32 var1digitpair; + uint64 maxdig; + uint64 carry; + uint64 newdig; + + /* + * As above, determine the number of result digits to compute and the + * (maximum possible) result weight, except that here we will be + * working in base NBASE^2 and so we process the digits of each input + * in pairs. + */ + /* digit pairs in each input */ + var1ndigitpairs = (var1ndigits + 1) / 2; + var2ndigitpairs = (var2ndigits + 1) / 2; + + /* digits in exact result */ + res_ndigits = var1ndigits + var2ndigits; + + /* digit pairs in exact result with at least one extra output digit */ + res_ndigitpairs = res_ndigits / 2 + 1; + + /* pair offset to align result to end of dig[] */ + pair_offset = res_ndigitpairs - var1ndigitpairs - var2ndigitpairs + 1; + + /* maximum possible result weight */ + res_weight = var1->weight + var2->weight + 1 + 2 * res_ndigitpairs - + res_ndigits; + + /* truncate computation based on requested rscale */ + maxdigits = res_weight + 1 + (rscale + DEC_DIGITS - 1) / DEC_DIGITS + + MUL_GUARD_DIGITS; + maxdigitpairs = (maxdigits + 1) / 2; + + res_ndigitpairs = Min(res_ndigitpairs, maxdigitpairs); + res_ndigits = 2 * res_ndigitpairs; + + if (res_ndigitpairs <= pair_offset) + { + /* All input digits will be ignored; so result is zero */ + zero_var(result); + result->dscale = rscale; + return; + } + var1ndigitpairs = Min(var1ndigitpairs, res_ndigitpairs - pair_offset); + var2ndigitpairs = Min(var2ndigitpairs, res_ndigitpairs - pair_offset); + + /* + * Since we will be working in base NBASE^2, we make dig[] an array of + * unsigned 64-bit integers, and since PG_UINT64_MAX is much larger + * than NBASE^4, this gives us a lot of headroom to avoid normalizing + * carries immediately. + * + * maxdig tracks the maximum possible value of any dig[] entry; when + * this threatens to exceed PG_UINT64_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 PG_UINT64_MAX / NBASE^2, so really + * we must normalize when digits threaten to exceed PG_UINT64_MAX - + * PG_UINT64_MAX / NBASE^2. * - * 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. + * To avoid overflow in maxdig itself, it actually represents the + * maximum possible value divided by NBASE^2-1, i.e., at the top of + * the loop it is known that no dig[] entry exceeds maxdig * + * (NBASE^2-1). * - * 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. + * The conversion of var1 to base NBASE^2 is done on the fly, as each + * new digit is required. The digits of var2 are converted upfront, + * and stored at the end of dig[]. To avoid loss of precision, the + * input digits are aligned with the start of digit pair array, + * effectively shifting them up (multiplying by NBASE) if the input + * has an odd number of NBASE digits. + */ + dig = (uint64 *) palloc(res_ndigitpairs * sizeof(uint64) + + var2ndigitpairs * sizeof(uint32)); + + /* convert var2 to base NBASE^2, shifting up if length is odd */ + var2digitpairs = (uint32 *) (dig + res_ndigitpairs); + + for (i2 = 0; i2 < var2ndigitpairs - 1; i2++) + var2digitpairs[i2] = var2digits[2 * i2] * NBASE + var2digits[2 * i2 + 1]; + + if (2 * i2 + 1 < var2ndigits) + var2digitpairs[i2] = var2digits[2 * i2] * NBASE + var2digits[2 * i2 + 1]; + else + var2digitpairs[i2] = var2digits[2 * i2] * NBASE; + + /* + * As above, we start by multiplying var2 by the least significant + * contributing digit pair from var1, storing the results at the end + * of dig[], and filling the leading slots with zeros. + * + * Here, however, digit pair i1 of var1 and digit pair i2 of var2 are + * multiplied and added to digit i1+i2+pair_offset of the accumulator + * array. + */ + i1 = var1ndigitpairs - 1; + if (2 * i1 + 1 < var1ndigits) + var1digitpair = var1digits[2 * i1] * NBASE + var1digits[2 * i1 + 1]; + else + var1digitpair = var1digits[2 * i1] * NBASE; + maxdig = var1digitpair; + + i2limit = Min(var2ndigitpairs, res_ndigitpairs - i1 - pair_offset); + dig_i1_off = &dig[i1 + pair_offset]; + + memset(dig, 0, (i1 + pair_offset) * sizeof(uint64)); + for (i2 = 0; i2 < i2limit; i2++) + dig_i1_off[i2] = (uint64) var1digitpair * var2digitpairs[i2]; + + /* + * Next, multiply var2 by the remaining digit pairs of var1, adding + * the results to dig[] at the appropriate offsets. */ + for (i1 = i1 - 1; i1 >= 0; i1--) { - int i2limit = Min(var2ndigits, res_ndigits - i1 - 2); - int *dig_i1_2 = &dig[i1 + 2]; + var1digitpair = var1digits[2 * i1] * NBASE + var1digits[2 * i1 + 1]; + if (var1digitpair == 0) + continue; + + /* Time to normalize? */ + maxdig += var1digitpair; + if (maxdig > (PG_UINT64_MAX - PG_UINT64_MAX / NBASE_SQR) / (NBASE_SQR - 1)) + { + /* Yes, do it (to base NBASE^2) */ + carry = 0; + for (i = res_ndigitpairs - 1; i >= 0; i--) + { + newdig = dig[i] + carry; + if (newdig >= NBASE_SQR) + carry = div_mod_NBASE_SQR(newdig, &newdig); + else + carry = 0; + dig[i] = newdig; + } + Assert(carry == 0); + /* Reset maxdig to indicate new worst-case */ + maxdig = 1 + var1digitpair; + } + + /* Multiply and add */ + i2limit = Min(var2ndigitpairs, res_ndigitpairs - i1 - pair_offset); + dig_i1_off = &dig[i1 + pair_offset]; for (i2 = 0; i2 < i2limit; i2++) - dig_i1_2[i2] += var1digit * var2digits[i2]; + dig_i1_off[i2] += (uint64) var1digitpair * var2digitpairs[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) + /* + * Now we do a final carry propagation pass to normalize back to base + * NBASE^2, and construct the base-NBASE result digits. + */ + alloc_var(result, res_ndigits); + res_digits = result->digits; + carry = 0; + for (i = res_ndigitpairs - 1; i >= 0; i--) { - carry = newdig / NBASE; - newdig -= carry * NBASE; + newdig = dig[i] + carry; + if (newdig >= NBASE_SQR) + carry = div_mod_NBASE_SQR(newdig, &newdig); + else + carry = 0; + res_digits[2 * i + 1] = (NumericDigit) ((uint32) newdig % NBASE); + res_digits[2 * i] = (NumericDigit) ((uint32) newdig / NBASE); } - else - carry = 0; - res_digits[i] = newdig; - } - Assert(carry == 0); + Assert(carry == 0); - pfree(dig); + pfree(dig); + + /* + * Adjust the result weight, if the inputs were shifted up during base + * conversion (if they had an odd number of NBASE digits). + */ + res_weight -= (var1ndigits & 1) + (var2ndigits & 1); + } /* * Finally, round the result to the requested precision. -- 2.35.3