mbutrovich commented on code in PR #11580:
URL: https://github.com/apache/arrow-rs/pull/11580#discussion_r4328291502
##########
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:
I added a bit more documentation.
--
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]
