mbutrovich opened a new issue, #11421: URL: https://github.com/apache/arrow-rs/issues/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`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-arith/src/numeric.rs#L344-L371)), decimal rescaling in [`make_downscaler`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-cast/src/cast/decimal.rs#L235-L246) and [`cast_decimal_to_integer`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-cast/src/cast/decimal.rs#L786-L801), time unit casts ([1](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-cast/src/cast/mod.rs#L2218), [2](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-cast/src/cast/mod.rs#L2423), [3](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arr ow-cast/src/cast/mod.rs#L2648)), and [`Display for i256`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-buffer/src/bigint/mod.rs#L83-L84), 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 costs a hardware divide for 64-bit and smaller types, a call to `__divti3` or `__udivti3` for `i128`, and a run of Knuth's Algorithm D ([`div_rem_knuth`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-buffer/src/bigint/div.rs#L83)) for `i256`. For `i128` the call 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)](https://gmplib.org/~tege/divcnst-pldi94.pdf), 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](https://github.com/apache/datafusion/blob/db83fccb9ac10f843be04daa5826501d2338caf4/datafusion/physical-plan/src/repartition/mod.rs#L941-L1002)) 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`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/src/math_funcs/wide_decimal_binary_expr.rs#L123-L146), which computes `value / divisor` and `value % divisor` by the same power of ten. That runs Algorithm D 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". Two questions for maintainers are at the end. #### Proposed method: Granlund-Montgomery at every width The proposal is Granlund-Montgomery at every width, with a mask or shift for powers of two. That's one algorithm, one correctness argument and one set of tests, and only the width of the high multiply changes: 32x32 to 64 bits for 32-bit types, which vectorizes, up to 256x256 to 512 bits for `i256`. Times are 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". <img width="1600" height="1012" alt="Image" src="https://github.com/user-attachments/assets/541e3306-3af9-430d-b8e2-cdacfd7435ad"; /> <details><summary>Numbers behind the chart</summary> | Case | Today | Granlund-Montgomery | Speedup | Variant | |---|---|---|---|---| | `UInt32Array % 13` | 3.87 | 0.97 | 4.0x | 32-bit, vectorized | | `Int32Array / 13` | 3.99 | 1.89 | 2.1x | 32-bit, vectorized | | `UInt64Array % 13` | 4.72 | 2.95 | 1.6x | 64-bit | | `UInt64Array % 16` | 4.77 | 0.62 | 7.7x | Mask | | `Decimal128` downscale, 18-digit values | 21.3 | 11.5 | 1.9x | 64-bit, chosen per batch | | `Decimal128` downscale by 10^8, 38-digit values | 32.4 | 24.7 | 1.3x | 128-bit | | `Decimal128` downscale by 10^18, 38-digit values | 73.7 | 24.7 | 3.0x | 128-bit | | Comet `i256` multiply, products fit in `i128` | 266.5 | 63.4 | 4.2x | 128-bit, chosen per batch | | Comet `i256` multiply, 48-digit products | 398.1 | 147.1 | 2.7x | 256-bit | </details> The 32-bit Granlund-Montgomery loops compile to NEON vector code, because their 32x32 to 64-bit multiply has vector forms. The 64-bit and wider methods need a 64x64 to 128-bit high multiply, which has no vector form on NEON, so those loops stay scalar. #### 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`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-array/src/array/primitive_array.rs#L1709-L1715)). #### Rounding and failures Each divisor type would offer `new`, which returns `None` for a zero divisor, `div_rem`, which returns the same results as `wrapping_div` and `wrapping_rem`, and `checked_div_rem`, which returns `None` only for `MIN / -1`. 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 `div_rem` instead of separate division routines: | Rounding | Correction to `(q, r)` | Callers today | |---|---|---| | Toward zero | None | arrow `div` and `rem`, Comet [decimal to integer casts](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/src/conversion_funcs/numeric.rs#L548) | | Floor (`div_floor`) | Subtract 1 from `q` when `r != 0` and the signs of `n` and `d` differ | Comet [`floor`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/src/math_funcs/floor.rs#L84-L86) | | Ceiling (`div_ceil`) | Add 1 to `q` when `r != 0` and the signs match | Comet [`ceil`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/src/math_funcs/ceil.rs#L84-L86) | | Half away from zero (`HALF_UP` in Java's `RoundingMode`) | Move `q` one step away from zero when `\|r\| >= \|d\| - \|r\|` | arrow-cast [`make_downscaler`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-cast/src/cast/decimal.rs#L235-L246), Comet `div_round_half_up` and [`round`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/src/math_funcs/round.rs#L218-L231), DataFusion [`round`](https://github.com/apache/datafusion/blob/db83fccb9ac10f843be04daa5826501d2338caf4/datafusion/functions/src/math/round.rs#L671-L680) | | Euclidean remainder (`rem_euclid`) | Add `\|d\|` to `r` when `r < 0` | Comet shuffle [`pmod`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/shuffle/src/comet_partitioning.rs#L123-L129) | 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 `checked_div_rem` only needs a per-row compare when the divisor is -1. 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`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-array/src/array/primitive_array.rs#L1018) 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 A sealed trait in `arrow-buffer`, next to [`ArrowNativeType`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-buffer/src/native.rs#L53) and `i256`, would give each integer type and `i256` its divisor, and integer and decimal kernels would take it as a bound. An associated type on [`ArrowNativeTypeOp`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-array/src/arithmetic.rs#L41) would also be non-breaking, because `ArrowNativeType` is sealed, but `ArrowNativeTypeOp` covers floats too, which would need a divisor that does nothing. The kernel builds a divisor for each batch and matches on its variant once, outside the loop. #### Tasks Steps that can start now: - [ ] Make [`i256::div_rem`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow-buffer/src/bigint/mod.rs#L488) public. It is private today, so callers that need both results run Algorithm D twice. - [ ] Add the benchmarks the kernel PRs need: integer division by a scalar ([`arithmetic_kernels`](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow/benches/arithmetic_kernels.rs#L29-L71) covers only `Float32`), 18-digit and 38-digit decimal downscale (the [existing case](https://github.com/apache/arrow-rs/blob/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120/arrow/benches/cast_kernels.rs#L453) uses 10-digit values), and `i256` rescaling. - [ ] Measure the prototype on x86-64, where the cost of the hardware 64-bit divide differs across generations. Find out why the fallible decimal cast gains less with Granlund-Montgomery (1.16x) than with `strength_reduce` (1.23x) when the infallible casts show the opposite. Steps after the API is agreed: - [ ] Add the divisor types for 32-bit and 64-bit integers, with tests against `wrapping_div` and `wrapping_rem` that cover `MIN`, `MAX`, `-1`, powers of two and random inputs. - [ ] Add the divisor types for `i128`, `u128` and `i256`, with the same tests plus `10^k - 1`, `10^k`, `10^k + 1` and values on each side of 2^64. - [ ] Add the rounding helpers, tested against `num::integer::div_floor`, `num::integer::div_ceil`, `rem_euclid`, and a half-up reference built on `/` and `%`. - [ ] Add the sealed trait that gives each integer type and `i256` its divisor. - [ ] Use it in each kernel, one PR per kernel with before and after numbers, in this order: integer `/` and `%` by a scalar, decimal downscale, 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`](https://github.com/apache/datafusion/blob/db83fccb9ac10f843be04daa5826501d2338caf4/datafusion/physical-plan/src/repartition/mod.rs#L961-L980) 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 `i32` divisor in [`pmod`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/shuffle/src/comet_partitioning.rs#L123-L129). 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.1.0, so it would adopt these later than DataFusion. - Two Comet changes need no arrow-rs release. Computing the remainder as `value - quotient * divisor` instead of a second `%` made its multiply kernel 1.1x to 1.5x faster. In `pmod`, `|r| < n`, so `r + n` is already in range and the second `%` can go, which made it 1.3x faster. ### Describe alternatives you've considered [Lemire, Kaser and Kurz (2019)](https://arxiv.org/abs/1902.01961) 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. Depending on [`strength_reduce`](https://github.com/ejmahler/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 numbers are 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. Integer remainder and division by a scalar, from the `batches` group of the `divproto` benchmark. Today is arrow-arith's `rem` or `div` kernel with a scalar divisor, and each method runs through `PrimitiveArray::unary`. A divisor of 200 gave results within 5% of 13. <img width="1600" height="996" alt="Image" src="https://github.com/user-attachments/assets/0677a038-1814-4aa8-ad32-2de80d45718a"; /> <details><summary>Numbers behind the chart, in us per batch</summary> | Operation | Today | Granlund-Montgomery | Mask | |---|---|---|---| | `UInt32Array % 13` | 3.87 | **0.97 (4.0x)** | not applicable | | `UInt32Array % 16` | 3.85 | 0.98 (3.9x) | **0.28 (13.7x)** | | `Int32Array / 13` | 3.99 | **1.89 (2.1x)** | not applicable | | `UInt64Array % 13` | 4.72 | **2.95 (1.6x)** | not applicable | | `UInt64Array % 16` | 4.77 | 2.97 (1.6x) | **0.62 (7.7x)** | </details> Comet's `pmod` on 8192 hash values, from the same 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. <img width="1600" height="708" alt="Image" src="https://github.com/user-attachments/assets/bb76edf8-9f85-4513-989d-633271e37299"; /> <details><summary>Numbers behind the chart, in us per batch</summary> | Partitions | Today | One division | Biased Granlund-Montgomery | Mask | |---|---|---|---|---| | 13 | 5.84 | 4.59 (1.3x) | **1.45 (4.0x)** | not applicable | | 16 | 5.98 | 4.52 (1.3x) | 1.45 (4.1x) | **0.28 (21.4x)** | </details> The arrow-cast `decimal_downscale` benchmark, which calls `cast` on an 8192-value array, with each method in `make_downscaler`. Today is main. <img width="1600" height="1236" alt="Image" src="https://github.com/user-attachments/assets/dddacee5-1d11-470d-83f7-d3678fe2d98b"; /> <details><summary>Numbers behind the chart, in us per batch</summary> | Cast | Today | `strength_reduce` | Granlund-Montgomery | Granlund-Montgomery with the per-batch 64-bit path | |---|---|---|---|---| | `(18, 6)` to `(18, 2)`, 18-digit values | 21.3 | 19.5 (1.1x) | 24.8 (0.9x) | **11.5 (1.9x)** | | `(38, 10)` to `(38, 2)`, 15-digit values | 21.4 | 19.5 (1.1x) | 24.7 (0.9x) | **11.4 (1.9x)** | | `(38, 10)` to `(38, 2)`, 38-digit values | 32.4 | 28.8 (1.1x) | **24.7 (1.3x)** | 26.8 (1.2x) | | `(38, 20)` to `(38, 2)`, 38-digit values | 73.7 | 28.7 (2.6x) | **24.7 (3.0x)** | 26.8 (2.8x) | | `(38, 10)` to `(30, 2)`, 30-digit values, fallible | 42.4 | **34.4 (1.2x)** | 36.4 (1.2x) | 38.2 (1.1x) | </details> 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`, which runs Algorithm D twice per row. <img width="1600" height="884" alt="Image" src="https://github.com/user-attachments/assets/3e26139c-8ea7-4b78-9c99-08dcf5315c15"; /> <details><summary>Numbers behind the chart, in us per batch</summary> | Case | Today | One Algorithm D call, remainder by multiply | Granlund-Montgomery, 256-bit | Granlund-Montgomery, 128-bit when products fit | |---|---|---|---|---| | The existing [`multiply`](https://github.com/apache/datafusion-comet/blob/8d0bb01a5a99924024e07c1722e6eddd9260e167/native/spark-expr/benches/wide_decimal.rs#L156) case, 30-digit products | 266.5 | 233.8 (1.1x) | 148.4 (1.8x) | **63.4 (4.2x)** | | 48-digit products, 34-digit quotients | 398.1 | 301.9 (1.3x) | **147.1 (2.7x)** | 153.2 (2.6x) | | 42-digit products, 8-digit quotients | 351.5 | 234.5 (1.5x) | **137.1 (2.6x)** | 143.6 (2.4x) | </details> The prototype is on a fork ([diff](https://github.com/apache/arrow-rs/compare/6b34163d51aa59f7cac2edcb74e8bd1ddc6b9120...mbutrovich:arrow-rs:ee0358cd48c9bdeb7ca5a4c52d479f949c17f736)), and its [README](https://github.com/mbutrovich/arrow-rs/blob/ee0358cd48c9bdeb7ca5a4c52d479f949c17f736/divproto/README.md) explains how to reproduce each table. [`divproto/`](https://github.com/mbutrovich/arrow-rs/tree/ee0358cd48c9bdeb7ca5a4c52d479f949c17f736/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. Questions for maintainers: - Does a public `arrow_buffer::divisor` module, with a `Divisor<T>` type and the sealed trait, work as the API, or would you prefer a different shape or name? - Can the arrow-rs benchmark runners, or anyone with x86-64 hardware, run `divproto` and the arrow-cast benchmark? All the numbers here are from aarch64, and the 64-bit results depend on the hardware divide. -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
