Add Newton (reciprocal based) division for large quotients

integer/CHANGELOG.md

@@ -5,6 +5,7 @@
 ### Add
 - NTT-based multiplication using Proth primes (`K·2^N + 1`), combined via Garner CRT.  Supports 64-bit and 32-bit Word targets.  Threshold at 4 000 words (~256 kbits).
 - Asymmetric NTT chunking: when one operand is much larger than the other, the shorter operand is forward-transformed once and reused across chunks.
+- Newton (reciprocal based) division: computes an approximate reciprocal of the divisor with Brent–Zimmermann's precision-doubling Newton iteration, then develops each quotient block by multiplication with a bounded correction loop. Dispatched (alongside schoolbook and Burnikel–Ziegler) when both the divisor and the quotient exceed ~6 000 words, where it is up to ~2× faster than Burnikel–Ziegler (see the new `ubig_div_asymmetric` benchmark). The threshold can be overridden at runtime via `DASHU_THRESHOLD_NEWTON_DIV` (requires `tuning` feature).
 - `UBig::from_u64` and `IBig::from_i64`, const on 32-bit and 64-bit targets.
 - Specialized Karatsuba squaring: uses 3 recursive squarings instead of 3 multiplications, with simplified diff handling.
 - Specialized Toom-Cook-3 squaring: evaluates a single polynomial instead of two, 5 recursive squarings instead of multiplications.

integer/benches/primitive.rs

@@ -53,6 +53,38 @@ add_binop_benchmark!(ubig_div, div, 6);
 add_binop_benchmark!(ubig_gcd, gcd, 6);
 add_binop_benchmark!(ubig_gcd_ext, gcd_ext, 5);
 
+/// Division with operands of very different sizes (asymmetric): a large dividend
+/// divided by a much smaller divisor, producing a large quotient.
+///
+/// This is the regime where the Newton (reciprocal based) division pays off:
+/// the cost of computing the divisor's reciprocal is amortized across many
+/// quotient blocks, each produced by plain multiplications. The symmetric
+/// `ubig_div` benchmark above produces only a tiny quotient, so it stays on the
+/// schoolbook path regardless of operand size.
+fn ubig_div_asymmetric(criterion: &mut Criterion) {
+    let mut rng = StdRng::seed_from_u64(SEED);
+    let mut group = criterion.benchmark_group("ubig_div_asymmetric");
+    group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    // divisor bit sizes large enough to reach the Newton regime, crossed with
+    // the dividend/divisor ratio.
+    for &divisor_bits in &[100_000usize, 300_000, 500_000, 1_000_000] {
+        for &ratio in &[4usize, 64] {
+            let divisor = random_ubig(divisor_bits, &mut rng);
+            let dividend = random_ubig(divisor_bits * ratio, &mut rng);
+            // dividend is far larger than the divisor, so the quotient has
+            // ~divisor_bits * (ratio-1) bits.
+            group.bench_with_input(
+                BenchmarkId::from_parameter(format!("{divisor_bits}/{ratio}")),
+                &(dividend, divisor),
+                |bencher, (n, d)| bencher.iter(|| n.div(d)),
+            );
+        }
+    }
+
+    group.finish();
+}
+
 fn ubig_pow(criterion: &mut Criterion) {
     let mut group = criterion.benchmark_group("ubig_pow");
     group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
@@ -202,6 +234,7 @@ criterion_group!(
     ubig_sub,
     ubig_mul,
     ubig_div,
+    ubig_div_asymmetric,
     ubig_gcd,
     ubig_gcd_ext,
     ubig_pow,

integer/src/div/mod.rs

@@ -14,6 +14,7 @@ use crate::{
 use alloc::alloc::Layout;
 
 mod divide_conquer;
+mod newton;
 mod simple;
 pub(crate) use simple::div_rem_highest_word;
 
@@ -237,6 +238,8 @@ pub fn memory_requirement_exact(lhs_len: usize, rhs_len: usize) -> Layout {
     assert!(lhs_len >= rhs_len && rhs_len >= 2);
     if rhs_len <= THRESHOLD_SIMPLE || lhs_len - rhs_len <= THRESHOLD_SIMPLE {
         memory::zero_layout()
+    } else if rhs_len > newton::threshold() && lhs_len - rhs_len > newton::threshold() {
+        newton::memory_requirement_exact(lhs_len, rhs_len)
     } else {
         divide_conquer::memory_requirement_exact(lhs_len, rhs_len)
     }
@@ -262,6 +265,8 @@ pub(crate) fn div_rem_in_place(
 
     if rhs.len() <= THRESHOLD_SIMPLE || lhs.len() - rhs.len() <= THRESHOLD_SIMPLE {
         simple::div_rem_in_place(lhs, rhs, fast_div_rhs_top)
+    } else if rhs.len() > newton::threshold() && lhs.len() - rhs.len() > newton::threshold() {
+        newton::div_rem_in_place(lhs, rhs, fast_div_rhs_top, memory)
     } else {
         divide_conquer::div_rem_in_place(lhs, rhs, fast_div_rhs_top, memory)
     }