mbutrovich commented on code in PR #11580:
URL: https://github.com/apache/arrow-rs/pull/11580#discussion_r4328191695
##########
arrow-buffer/src/bigint/div.rs:
##########
@@ -72,6 +93,104 @@ fn div_rem_small<const N: usize>(numerator: &[u64; N],
divisor: u64) -> ([u64; N
(numerator, rem_padded)
}
+/// Like [`div_rem_small_word`], but divides each digit by multiplying with a
+/// reciprocal of the divisor, computed once
+#[cfg(any(
+ test,
+ all(target_pointer_width = "64", any(not(target_arch = "x86_64"), miri))
+))]
+fn div_rem_small_reciprocal<const N: usize>(
+ numerator: &[u64; N],
+ divisor: u64,
+) -> ([u64; N], [u64; N]) {
+ // Shifting both operands so the divisor's top bit is set leaves the
quotient unchanged
+ let shift = divisor.leading_zeros();
+ let divisor = divisor << shift;
+ let reciprocal = reciprocal_word(divisor);
+ let numerator = full_shl(numerator, shift);
+
+ // The top digit is below 2^shift, and so below the divisor
+ let mut rem = numerator[N];
+ let mut quotient = [0; N];
+ for i in (0..N).rev() {
+ let (q, r) = div_rem_word_reciprocal(rem, numerator[i], divisor,
reciprocal);
+ quotient[i] = q;
+ rem = r;
+ }
+
+ let mut rem_padded = [0; N];
+ rem_padded[0] = rem >> shift;
+ (quotient, rem_padded)
+}
+
+/// Returns `floor((2^128 - 1) / divisor) - 2^64` for a divisor with its top
bit set
+///
+/// This is algorithm 2 of [Möller and Granlund, "Improved division by
invariant integers" (2011)](https://gmplib.org/~tege/division-paper.pdf),
+/// which avoids a `u128` division
+#[cfg(any(
+ test,
+ all(target_pointer_width = "64", any(not(target_arch = "x86_64"), miri))
+))]
+fn reciprocal_word(divisor: u64) -> u64 {
+ debug_assert_eq!(divisor >> 63, 1);
+
+ // `floor((2^19 - 3 * 2^8) / d9)` for each value of the top 9 bits `d9`
+ const TABLE: [u16; 256] = {
+ let mut table = [0; 256];
+ let mut i = 0;
+ while i < 256 {
+ table[i] = (((1 << 19) - 3 * (1 << 8)) / (i as u32 + 256)) as u16;
+ i += 1;
+ }
+ table
+ };
+
+ let d0 = divisor & 1;
+ let d9 = divisor >> 55;
+ let d40 = (divisor >> 24) + 1;
+ let d63 = (divisor >> 1) + d0;
+
+ // `d9` is in 256..512, so masking gives `d9 - 256` without a bounds check
+ let v0 = u64::from(TABLE[(d9 & 0xff) as usize]);
+ // Two Newton iterations, each about doubling the bits of precision
+ let v1 = (v0 << 11) - ((v0 * v0 * d40) >> 40) - 1;
+ let v2 = (v1 << 13) + ((v1 * ((1 << 60) - v1 * d40)) >> 47);
+ // `2^96 - v2 * d63 + floor(v2 / 2) * d0`, which the paper shows fits in
64 bits
+ let e = ((v2 >> 1) & d0.wrapping_neg()).wrapping_sub(v2.wrapping_mul(d63));
+ let v3 = (v2 << 31).wrapping_add(((u128::from(v2) * u128::from(e)) >> 65)
as u64);
+ // `floor((v3 + 2^64 + 1) * divisor / 2^64)`, without overflowing a `u128`
+ let p = (((u128::from(v3) + 1) * u128::from(divisor)) >> 64) as u64;
+ v3.wrapping_sub(p.wrapping_add(divisor))
+}
Review Comment:
> @mbutrovich, this is a lot to take in. Could you add some explanation,
particularly for the magic numbers? Apart from the factors of 8, I’m not sure
what some of the bit shifts are doing. Clearer variable names would help too.
I think changing the variable names would hurt readability.
`reciprocal_word` exactly matches Algorithm 2 from the paper, and I think
trying to ascribe our own meaning here would make it more difficult to map to
the source. I can try to bring some of the comments from the Algorithm 2
listing, but at some point we have to rely on the reference and not reproduce
the paper in comment form.
<img width="405" height="313" alt="Screenshot 2026-10-20 at 1 37 55 PM"
src="https://github.com/user-attachments/assets/3e0bc6ac-0033-439c-9230-65e6422efd25";
/>
--
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]
