[Epic] Faster division in integer and decimal kernels #11421

Is your feature request related to a problem or challenge?

Many kernels divide every value in an array by the same divisor. The divisor is known before the loop starts, but not at compile time. In arrow-rs this includes integer / and % by a scalar (integer_op), decimal rescaling in make_downscaler and cast_decimal_to_integer, the rescaling of each i256 product in multiply_fixed_point (divide_and_round), time unit casts (1, 2, 3), and Display for i256, which divides by 10^38 twice per value.

LLVM turns division by a compile-time constant into a multiply and shifts, but it can't do that when the divisor is only known at run time, so each row pays for a full division. For i128 the call to __divti3 or __udivti3 stays even with a literal divisor such as 10^18 (rustc 1.98.1, aarch64).

The usual fix is to compute a reciprocal of the divisor once and turn each division into a high multiply, a shift, and a small correction. This issue uses the method of Granlund and Montgomery (1994), which needs an N-bit multiplier for N-bit dividends.

Downstream projects are already doing this on their own. @Dandandan added a hand-written StrengthReducedU64 (code) to DataFusion's hash repartitioning in apache/datafusion#21900, which made TPC-H SF10 queries up to 1.16x faster. Comet's WideDecimalBinaryExpr rescales every i256 product with div_round_half_up, which computes value / divisor and value % divisor by the same power of ten. That runs Knuth's Algorithm D (div_rem_knuth) twice per row. With these primitives in arrow-rs, arrow's own kernels could use them, and DataFusion and Comet would share one tested implementation.

Describe the solution you'd like

Precomputed divisors in arrow-buffer, next to i256 and its division code, that arrow's kernels and downstream projects can share. The design below follows from the prototype measurements under "Additional context".

Proposed method: Granlund-Montgomery at every width

The proposal is Granlund-Montgomery at every width, with a mask or shift for powers of two. The width of the high multiply changes, from 32x32 to 64 bits for 32-bit types, which vectorizes, up to 256x256 to 512 bits for i256. Within a width, the form depends on the divisor and the target: the N-bit multiplier of figure 4.2 when it exists, the N + 1-bit multiplier of figure 4.1 or 5.1 otherwise, and for u32 on targets with a one-instruction saturating add, the rounded-down multiplier from Labor of Division (Episode III).

Times are prototype measurements in us per 8192-value batch. "Today" is arrow-arith's rem or div kernel for the integer rows, the current arrow-cast code for the Decimal128 rows, and Comet's current code for the i256 rows. The detailed tables are under "Additional context".

Methods by type

This epic changes division by a scalar. Division by an array, where each row has its own divisor, keeps today's methods, because a precomputed divisor would be used only once.

Rust's compiler-builtins provides __divti3 and __udivti3. On x86_64 they build on the div instruction, which divides 128 bits by 64. On other targets they use a sequence of 64-bit hardware divides.

i256::div_rem picks a method by the size of its operands:

Choosing the variant per batch

For i128 and i256, the kernel would pick the narrowest Granlund-Montgomery variant that fits each batch, from a scan of the values' magnitude. One variant per type isn't enough: the 128-bit variant alone made 18-digit Decimal128 casts 16% slower than today (24.8 vs 21.3 us), while the 64-bit variant made them 1.9x faster. In Comet, the 128-bit variant made products that fit in i128 4.2x faster than today. The scan cost 4% to 8% when it found wide values. The choice can't come from the declared precision, because arrow doesn't check stored values against it (with_precision_and_scale).

Rounding and failures

Callers round the quotient in different ways, but each mode is a sign and remainder correction applied to the truncated (q, r). These could be generic helpers on top of wrapping_div_rem instead of separate division routines:

The half-up test is written as |r| >= |d| - |r| because 2 * |r| can overflow when |d| is close to the type's maximum.

Division can fail only on a zero divisor or MIN / -1, and the divisor alone says whether either is possible, so only a divisor of -1 needs a per-row compare. A kernel can test for -1 once and use unary for every other divisor, which matters because try_unary visits valid slots one at a time when there are nulls (measurements). Turning failures into nulls (Spark's legacy mode, casts with safe: true) stays in the kernel, the way unary_opt maps None to null today.

Today [NULL, NULL] / 0, [] / 0 and [NULL, NULL] % 0 return Ok, because try_unary skips null slots. [1, NULL] / 0 returns DivideByZero. Kernels should build the divisor only once there's a valid value to divide, so they keep that behavior.

Where it plugs in

#11467 sets the API that the i128, u128 and i256 types would follow: each divisor type is a #[non_exhaustive] enum whose variants are branch-free types, and the variant set differs by type. DivisorU32 has PowerOfTwo, Up, Down and AddIndicator, DivisorU64 has PowerOfTwo, Up and AddIndicator, and the signed types have PowerOfTwo and AddIndicator. Kernels match on the variants with a macro that lists them per type, so the variant sets can differ without a shared trait. A comparison of the macro with a visitor over a sealed trait found them within 4% of each other. The kernel builds a divisor for each batch, matches on its variant once, outside the loop, and moves the variant into its closure. Each divisor type can have a checked_div_rem, implemented as Some(self.wrapping_div_rem(n)) where the divisor can't overflow.

Tasks

Some of the steps speed up division that doesn't need a precomputed divisor, or remove costs that would hide the divisor's gains in the kernel benchmarks.

Steps that can start now:

• Add i256::wrapping_div_rem and i256::checked_div_rem (feat(arrow-buffer): add i256::wrapping_div_rem and i256::checked_div_rem #11425).
• Move the closure into the iterator in PrimitiveArray::unary (perf(arrow-array): move the closure into the iterator in PrimitiveArray::unary #11525).
• Don't build an error for each decimal div value that overflows when scaled (perf(arrow-arith): don't build an error for each decimal div value that overflows when scaled #11479).
• Let try_unary inline the decimal / and % by a scalar closures (try_unary calls the decimal_op closure out of line for every value #11529).
• Add the benchmarks the kernel PRs need.
  • Decimal downscale, including i256 rescaling, and decimal to integer casts (bench(arrow-cast): add decimal downscale and decimal to integer cast benchmarks #11426).
  • Integer and decimal / and % by a scalar, and multiply_fixed_point (bench(arrow-arith): add integer and decimal division by scalar benchmarks #11427).
• Divide once per row where code needs both the quotient and the remainder.
  • multiply_fixed_point (Divide once per row in multiply_fixed_point #11453, perf(arrow-arith): divide once per row in multiply_fixed_point #11526).
  • make_downscaler and scaled_div, through div_rem_wrapping and div_rem_checked on ArrowNativeTypeOp (Add combined division and remainder to ArrowNativeTypeOp #11454, feat(arrow-array): add div_rem_wrapping and div_rem_checked to ArrowNativeTypeOp #11528).
• Speed up general i256 division (Speed up i256 division #11451).
  • Inline bits, add_assign and sub_assign in i256 division (perf(arrow-buffer): inline bits, add_assign and sub_assign in i256 division #11457).
  • Divide magnitudes below 2^128 with one u128 division (perf(arrow-buffer): divide i256 magnitudes below 2^128 with one u128 division #11530).
  • On 64-bit targets other than x86_64, divide by a single word with a precomputed reciprocal (perf(arrow-buffer): divide i256 by a single word with a reciprocal on 64-bit targets other than x86_64 #11580).

Steps that depend on the divisor API, which its first PR will settle:

• Add the divisor types for 32-bit and 64-bit integers (perf(arrow-buffer): add precomputed divisors for 32-bit and 64-bit integers #11467).
• Try two signed forms from @mightsleep's review of perf(arrow-buffer): add precomputed divisors for 32-bit and 64-bit integers #11467: an unsigned high multiply on |n| with the rounded-up multiplier, and for i32 a quotient from an f64 reciprocal rounded up by one ulp. Both aim at baseline x86-64, where SSE2 has no signed 32-bit widening multiply and i32 gains least from the divisor. Each needs a test of every i32 divisor against wrapping_div at both ends of the last quotient steps first, including n = i32::MIN, whose magnitude doesn't fit in 31 bits.
• Add the divisor types for i128, u128 and i256, with the tests from perf(arrow-buffer): add precomputed divisors for 32-bit and 64-bit integers #11467 plus 10^k - 1, 10^k, 10^k + 1 and values on each side of 2^64. Which variants the wider types get, such as Up for u128, is open.
• Add the rounding helpers, tested against num::integer::div_floor, num::integer::div_ceil, rem_euclid, and a half-up reference built on / and %. For a positive divisor, div_euclid can be one unsigned high multiply on the one's complement of n, from the same review.
• Use it in each kernel, with a macro arm per divisor type, one PR per kernel with before and after numbers, in this order: integer / and % by a scalar, decimal downscale, multiply_fixed_point, decimal / and % by a scalar, decimal to integer, time unit casts, and Display for i256.

Follow-ups in other repos

• DataFusion could replace StrengthReducedU64 once it uses an arrow release that has the u64 divisor. A power-of-two case in the divisor type covers what partition_indices does today. The swap needs a repartition benchmark to show there's no regression.
• Comet could use the i256 divisor in WideDecimalBinaryExpr, with the i128 variant chosen per batch, and the u32 divisor in pmod. Adding 2^31 to the hash lets one unsigned remainder give the result, and with a vectorized 32-bit Granlund-Montgomery divisor, pmod was 4.0x faster than today. Comet depends on arrow 59.2.0 and DataFusion 55.2.0, so it would adopt these later than DataFusion.
• Two Comet changes need no arrow-rs release (perf: divide once per row in pmod and div_round_half_up datafusion-comet#6781).

Describe alternatives you've considered

Lemire, Kaser and Kurz compute the remainder from a 2N-bit reciprocal, as DataFusion's StrengthReducedU64 does. In the prototype it beat Granlund-Montgomery only on the 64-bit path of the Decimal128 rounding (9.7 vs 13.5 us per batch, where Granlund-Montgomery is still 1.6x to 1.8x faster than today), tied on UInt64Array % 13, and lost at 32 bits, where Granlund-Montgomery vectorizes. That one gain isn't worth a second algorithm's code and tests. #11467 does use the paper's observation in section 3.1 that the remainder doesn't depend on the sign of the divisor.

Depending on strength_reduce would avoid carrying the code. It supports unsigned types only and has no 256-bit type. In the arrow-cast measurements below it was slower than the fastest hand-written variant in every case except the fallible cast. Its last release (0.2.4) was in November 2022.

Precomputed tables of the reciprocals of 10^k would remove construction cost for decimal rescaling. In the prototype, building a divisor costs 1.4 to 2.5 ns for u64, 22 ns for i128 and 451 ns for i256. The i256 constructor divides one bit at a time and could use Algorithm D instead. That's at most 0.45 us per 8192-value batch. For scalars and very small batches, plain division is cheaper, so kernels would need a row-count threshold or the tables.

This epic is about how a division is computed, not about decimal result types. Precision-loss decimal arithmetic stays out of scope, as @tustvold decided in #4664.

Additional context

These are prototype measurements against main at 6b34163, kept for the design choices above. Each task's PR reports its own before and after numbers from the benchmark bot. They're from aarch64 only: an Apple M5 Max, release builds and criterion, in us per 8192-value batch, with the fastest method in each table row in bold. The divproto and arrow-cast benchmarks used rustc 1.99.0, and the Comet benchmark used 1.98.1. The arrow-cast and Comet numbers are medians of three runs.

#11467 replaces the prototype's 32-bit and 64-bit divisors, and its benchmark bot runs compare wrapping_div_rem with / and % on x86_64, x86_64 at x86-64-v3 and arm64.

Comet's pmod on 8192 hash values, from the divproto benchmark. Today is Comet's code, which divides twice for negative hashes. "One division" replaces the second % with r + n. The biased method adds 2^31 to the hash so that one unsigned remainder gives the result.

The arrow-cast decimal_downscale benchmark, which calls cast on an 8192-value array, with each method in make_downscaler. Today is main.

The per-batch 64-bit path in this chart used the hardware divide. A separate run with Granlund-Montgomery's 64-bit variant there measured within 2% of it (11.2 vs 11.4 us for 18-digit values), and the prototype now uses that variant.

Comet's wide_decimal benchmark, which evaluates the WideDecimalBinaryExpr multiply on an 8192-row batch, with each method in the kernel. Today is div_round_half_up.

The prototype is on a fork (diff), and its README explains how to reproduce each table. divproto/ is a standalone crate with the divisors, their tests and the microbenchmarks. The arrow-cast change in the same diff is the last column of the downscale table, and divproto/variants/ has patches for the other columns and for the Comet kernel.