Fix a missing part of the Burnikel-Ziegler algorithm

Cargo.toml

@@ -7,7 +7,7 @@ categories = [ "algorithms", "data-structures", "science" ]
 license = "MIT OR Apache-2.0"
 name = "num-bigint"
 repository = "https://github.com/rust-num/num-bigint"
-version = "0.5.0"
+version = "0.5.1"
 readme = "README.md"
 exclude = ["/ci/*", "/.github/*"]
 edition = "2021"

RELEASES.md

@@ -1,3 +1,11 @@
+# Release 0.5.1 (2026-07-04)
+
+- [Fixed a regression in the B-Z division algorithm.][348]
+
+[348]: https://github.com/rust-num/num-bigint/pull/338
+
+**Contributors**: @cuviper
+
 # Release 0.5.0 (2026-07-02)
 
 - [Upgrade to `rand` v0.10 and/or v0.9, and split `rand_core`][338].

src/biguint/division.rs

@@ -295,19 +295,34 @@ fn div_rem_burnikel_ziegler(u: &BigUint, d: &BigUint) -> (BigUint, BigUint) {
         b2: BigUint,
         level: usize,
     ) -> (BigUint, BigUint) {
-        let (mut q, c) = div_two_digit_by_one(a1, a2, b1.clone(), level);
-        let (mut r, rsub) = (concat_biguint(&c, a3, level), &q * &b2);
+        // Algorithm 2 step 3 distinguishes A₁<B₁ or A₁≥B₁, preserving
+        // div_two_digit_by_one's asserted invariant that `ah < b`.
+        let (mut q, r1, d);
+        if a1 < b1 {
+            (q, r1) = div_two_digit_by_one(a1, a2, b1.clone(), level);
+            d = &q * &b2;
+        } else {
+            // set Q = βⁿ-1 and set R₁ = [A₁;A₂] - [B₁;0] + [0;B₁]
+            //                        (= [A₁;A₂] - QB₁)
+            q = BigUint {
+                data: BigDigits::from_vec(vec![BigDigit::MAX; level]),
+            };
+            r1 = concat_biguint(&(a1 - &b1), a2, level) + &b1;
+            // also Q*B₂ = (βⁿ-1)B₂ = βⁿB₂ - B₂
+            d = (&b2 << (level * usize::from(BITS))) - &b2;
+        };
         // The algorithm guarantees at most two corrections are needed.
-        if r < rsub {
+        let mut r = concat_biguint(&r1, a3, level);
+        if r < d {
             let b = concat_biguint(&b1, b2, level);
             q -= 1u32;
             r += &b;
-            if r < rsub {
+            if r < d {
                 q -= 1u32;
                 r += &b;
             }
         }
-        (q, r - rsub)
+        (q, r - d)
     }
 
     let mut level = 1 << ilog2(u.data.len());

tests/biguint.rs

@@ -946,6 +946,30 @@ fn test_div_rem_burnikel_ziegler_u1_equals_d() {
     assert_eq!(r, remainder, "remainder should be 42");
 }
 
+#[test]
+fn test_div_rem_burnikel_ziegler_a1_equals_b1() {
+    // This triggered a bug for the missing a1 >= b1 guard in div_three_halves_by_two
+    // (Algorithm 2, step 3b of the Burnikel-Ziegler paper).
+    //
+    // d has 65 u64-digits (130 u32-digits). In div_two_digit_by_one, the normalization
+    // shift pads d to fill `level` digits. After div_two_digit_by_one_normalized splits
+    // both ah and the shifted d at the half-level, the second div_three_halves_by_two call
+    // receives a remainder whose upper half equals b1 (the upper half of the shifted d).
+    // It then called div_two_digit_by_one(a1, a2, b1, level) with a1 == b1, violating the
+    // a1 < b1 precondition.
+    //
+    // The quotient ((1 << 4096) - 1) is chosen so the remainder is d - 1, which shares
+    // the same upper half as d after normalization.
+    let d = (BigUint::one() << 4159usize) + BigUint::one();
+    let expected_q = (BigUint::one() << 4096usize) - BigUint::one();
+    let expected_r = &d - BigUint::one();
+    let u = &d * &expected_q + &expected_r;
+
+    let (q, r) = u.div_rem(&d);
+    assert_eq!(q, expected_q, "quotient mismatch");
+    assert_eq!(r, expected_r, "remainder mismatch");
+}
+
 #[test]
 fn test_div_ceil() {
     fn check(a: &BigUint, b: &BigUint, d: &BigUint, m: &BigUint) {

tests/fuzzed.rs

@@ -183,3 +183,44 @@ fn fuzzed_mul_2() {
     assert_eq!(c * a * b, result);
     assert_eq!(c * b * a, result);
 }
+
+#[test]
+fn fuzzed_sqrt_1() {
+    // This input failed an assertion in the Burnikel-Ziegler division algorithm, but only on
+    // 32-bit. See also `test_div_rem_burnikel_ziegler_a1_equals_b1` for a more targeted test
+    // which failed on both 32-bit and 64-bit targets.
+    let hex = "\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff\
+        ffffffffffffffffffffffffffffffff974810cf4c437ba04120a1154166bbef02433931c39106e4\
+        6f2837ef218ee18e59ae8d861c29c5b276b221dfc125c5cff35dcadc6caef8f005fc56a94c280c66\
+        8118da37878e14248f36ae7350949811f342eba025d9ae5164b650727dd1c772504faadc3620bdb7\
+        cfee873de995aae72d18e100a8577eac452fc94350cfdafbe88f23c1aba8d830e517c0b4a76ff158\
+        ece8d7f2ebba856777075fe6a9bd6b1026933a35f287614f2903433e6159b4bba785b6b9b43c78d6\
+        1a93549d0be51d466602ea3d63251ee5b99efcc8e2acf47f9ef0cee82caf2d001295a3d27fec6eba\
+        3995492f907ab2241b1127740ab045796990bf62b6c8b1099dcf291df3184c5b9b4262a4752ece67\
+        d8e470674b27731012d9599c129bfa8dcd77990c747915f919cbb81319ede35b942d1cad3d90e33f\
+        9d4d4bb9a856d74e2bc7dd08a199c12c0f7f7b1832725ec180d23347ea657c027af747e2d103e585\
+        1f8c9fdfe738138d30fa810c59b53d1f8d7fb409e97a46c58049de0ebe7d5ae92fe1fb2e8d762d03\
+        3027158643541a06891dde1258624d714389f68b532406c10692d83b97262078ff66cb4907bd7273\
+        6a614764f22801d93cef8ba49fc2482cea302ac9dc0822fe341350dce3feaecd80651415971d4a1d\
+        cf8c0b3cdb189928321496a0ace3ab8d566b9afb26a812fa1fea858929006cdbd988d7b3ee0fb951\
+        87c8e4e120f63a4fab417604b9b4e86a2158cbc52c4a6f76c6bfee18b4969cba2c5e8408ecc80cb3\
+        b371469d375c372c0000000000000000000000000000000000000000000000000000000000000000\
+        00000000000000000000000000000000000000000000000000000000000000000000000000000000\
+        000000000000000000000000000000";
+
+    let x = BigUint::from_str_radix(hex, 16).unwrap();
+    let r = x.sqrt();
+    assert!(r.pow(2) <= x && (r + 1u32).pow(2) > x);
+}