Minor improve on division and implement Montgomery reducer

integer/CHANGELOG.md

@@ -6,6 +6,7 @@
 - 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.
 - `UBig::from_u64` and `IBig::from_i64`, const on 32-bit and 64-bit targets.
+- `monty` module: Montgomery modular arithmetic for odd moduli. New [`MontgomeryRepr`](integer::monty::MontgomeryRepr) (precomputed Montgomery constants) and [`Montgomery`](integer::monty::Montgomery) (values in Montgomery form) with multiplication, squaring, addition, subtraction, negation, doubling, exponentiation, and inversion. Single/double-word moduli delegate to `num-modular`; multi-word moduli use word-by-word REDC with a double-word "addmul_2" kernel (the operand product reuses the crate's fast multiply), beating the Barrett division path of `modular::Reduced` for multiplication/squaring/exponentiation at roughly 256–4096 bits. For inverse-heavy computation, `modular::Reduced` remains faster (a Montgomery inverse must convert out of and back into Montgomery form).
 - 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.
 - Specialized NTT squaring: single forward transform instead of two, pointwise square instead of multiply.
@@ -18,7 +19,13 @@
 - Lowered the Karatsuba→Toom-3 multiplication threshold from 192 to 96 words, giving Toom-Cook-3 at ~6000 bits instead of ~12000 bits — closes the gap with malachite at ~10000-bit sizes.
 - NTT coefficient width increased from 16 to 64 bits (K_eff=3 for 64-bit, K_eff=2 otherwise), roughly halving the transform length at each step.
 - NTT multiplication auto-selects `K_eff = 2` primes when headroom allows, skipping the third prime.
+- Multiplication thresholds can be overridden at runtime via `DASHU_THRESHOLD_SIMPLE`, `DASHU_THRESHOLD_KARATSUBA`, and `DASHU_THRESHOLD_NTT` environment variables (requires `tuning` feature).
+- Division threshold (schoolbook ↔ divide-and-conquer crossover) can be overridden at runtime via the `DASHU_THRESHOLD_SIMPLE_DIV` environment variable (requires `tuning` feature). Values below 3 are clamped to 3 to uphold the divide-and-conquer algorithm's `n_lo >= 2` invariant.
+- Montgomery multiplication: the `sqr()` method and `pow_nontrivial` entry point now avoid a redundant clone+overwrite of the multi-word value, saving one `Box<[Word]>` allocation and an `s`-word copy per squaring.
+- Montgomery multiplication: `mul_in_place_large` uses pointer-identity instead of full element comparison to detect self-multiplication (`a *= a`), avoiding an `O(s)` scan on the common distinct-operands path.
+- Montgomery multiplication: `mul_normalized_large` and `sqr_normalized_large` share a common `finish_monty_product` helper for the REDC+canonicalize pipeline tail.
 - Multiplication thresholds can be overridden at runtime via `DASHU_THRESHOLD_SIMPLE_MUL`, `DASHU_THRESHOLD_KARATSUBA_MUL`, and `DASHU_THRESHOLD_NTT_MUL` environment variables (requires `tuning` feature).
+- Montgomery multiplication: `&Montgomery * &Montgomery` for Large operands now builds the result directly from scratch memory instead of cloning one operand then overwriting every word via `mul_in_place_large`, saving an `s`-word copy per multiply. (Add, sub, and neg on references still clone since their in-place operations seed the output buffer from the operand value.)
 
 ### Change
 - Multiplication threshold env vars renamed with `_MUL` suffix: `DASHU_THRESHOLD_SIMPLE_MUL`, `DASHU_THRESHOLD_KARATSUBA_MUL`, `DASHU_THRESHOLD_NTT_MUL` (was without suffix).
@@ -35,6 +42,11 @@
 - `test_unpack_carry_propagation` had a hardcoded 64-bit shift assumption; now derived from `Word::BITS` so it works on 32-bit.
 - Various clippy warnings (`let_and_return`, `too_many_arguments`, `needless_range_loop`, `type_complexity`) resolved across the NTT module.
 
+### Refactor
+- Extracted `simple::MIN_LEN = 3` in the division module (analogous to `mul::karatsuba::MIN_LEN`); the tuning override now clamps against this named constant instead of a magic literal.
+- Moved the word-level multiplication kernels (`add_mul_word_same_len_in_place`, `add_mul_word_in_place`, `sub_mul_word_same_len_in_place`) from `mul/mod.rs` into `mul/simple.rs` alongside the other schoolbook kernels, and widened `add_mul_dword_same_len_in_place` from `pub(crate)` to `pub`. The old `mul::*` paths still work via re-exports.
+- Added `forward_modular_binop_to_assign!`, `impl_modular_commutative_op_for_ref!`, `impl_modular_binop_ref_ref_by_clone!`, and `forward_modular_binop_to_ref_ref!` macros in `helper_macros.rs` to generate the boilerplate `Op<T>`/`Op<&T>`/`OpAssign<T>` impls for any modular type parameterized by a lifetime (taking the target type as a parameter). Both `Montgomery` and `Reduced` Add, Sub, Mul and Div now share these macros instead of repeating the four-by-value/by-ref/assign variants per operator.
+
 ## 0.4.2
 
 - Add `UBig::ones`.

integer/Cargo.toml

@@ -75,3 +75,8 @@ harness = false
 name = "io"
 required-features = ["rand"]
 harness = false
+
+[[bench]]
+name = "modular"
+required-features = ["rand"]
+harness = false

integer/benches/modular.rs

@@ -0,0 +1,152 @@
+//! Benchmarks comparing the two modular-arithmetic backends:
+//! Barrett reduction (`ConstDivisor` / `Reduced`) vs Montgomery reduction
+//! (`MontgomeryRepr` / `Montgomery`).
+//!
+//! Run (full):  `cargo bench -p dashu-int --bench modular --features rand`
+//! Run (quick): `cargo bench -p dashu-int --bench modular --features rand -- --sample-size 10 --warm-up-time 1 --measurement-time 1`
+
+use criterion::{
+    criterion_group, criterion_main, AxisScale, BenchmarkId, Criterion, PlotConfiguration,
+};
+use dashu_int::{fast_div::ConstDivisor, monty::MontgomeryRepr, UBig};
+use rand_v08::prelude::*;
+
+const SEED: u64 = 1;
+
+/// Bit sizes of the moduli benchmarked. On a 64-bit target these span 4, 8, 16, 32, 64,
+/// 128 and 256 words — covering the schoolbook, Karatsuba, Toom-3 and NTT regimes.
+const BITS: &[usize] = &[256, 512, 1024, 2048, 4096, 8192, 16384];
+/// A shorter list for the (expensive) pow benchmark.
+const POW_BITS: &[usize] = &[256, 1024, 4096];
+
+fn random_ubig<R: Rng + ?Sized>(bits: usize, rng: &mut R) -> UBig {
+    rng.gen_range(UBig::ONE << (bits - 1)..UBig::ONE << bits)
+}
+
+/// An odd modulus in the given bit range (Montgomery requires an odd modulus).
+fn random_odd_ubig<R: Rng + ?Sized>(bits: usize, rng: &mut R) -> UBig {
+    random_ubig(bits, rng) | UBig::ONE
+}
+
+/// Benchmark a binary modular operation (`*`, `+`, `-`) for both backends.
+macro_rules! binop_bench {
+    ($group:ident, $op:tt) => {
+        fn $group(c: &mut Criterion) {
+            let mut rng = StdRng::seed_from_u64(SEED);
+            let mut group = c.benchmark_group(stringify!($group));
+            group.plot_config(
+                PlotConfiguration::default().summary_scale(AxisScale::Logarithmic),
+            );
+
+            for &bits in BITS {
+                let m = random_odd_ubig(bits, &mut rng);
+                let barrett = ConstDivisor::new(m.clone());
+                let monty = MontgomeryRepr::new(m.clone());
+                let a = random_ubig(bits, &mut rng);
+                let b = random_ubig(bits, &mut rng);
+                let (ba, bb) = (barrett.reduce(a.clone()), barrett.reduce(b.clone()));
+                let (ma, mb) = (monty.reduce(a.clone()), monty.reduce(b.clone()));
+
+                group.bench_with_input(
+                    BenchmarkId::new("barrett", bits),
+                    &(ba, bb),
+                    |bencher, (a, b)| bencher.iter(|| a $op b),
+                );
+                group.bench_with_input(
+                    BenchmarkId::new("monty", bits),
+                    &(ma, mb),
+                    |bencher, (a, b)| bencher.iter(|| a $op b),
+                );
+            }
+
+            group.finish();
+        }
+    };
+}
+
+binop_bench!(modular_mul, *);
+binop_bench!(modular_add, +);
+binop_bench!(modular_sub, -);
+
+fn modular_sqr(c: &mut Criterion) {
+    let mut rng = StdRng::seed_from_u64(SEED);
+    let mut group = c.benchmark_group("modular_sqr");
+    group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    for &bits in BITS {
+        let m = random_odd_ubig(bits, &mut rng);
+        let barrett = ConstDivisor::new(m.clone());
+        let monty = MontgomeryRepr::new(m.clone());
+        let a = random_ubig(bits, &mut rng);
+        let ba = barrett.reduce(a.clone());
+        let ma = monty.reduce(a.clone());
+
+        group.bench_with_input(BenchmarkId::new("barrett", bits), &ba, |b, a| b.iter(|| a.sqr()));
+        group.bench_with_input(BenchmarkId::new("monty", bits), &ma, |b, a| b.iter(|| a.sqr()));
+    }
+
+    group.finish();
+}
+
+fn modular_pow(c: &mut Criterion) {
+    let mut rng = StdRng::seed_from_u64(SEED);
+    let mut group = c.benchmark_group("modular_pow");
+    group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    for &bits in POW_BITS {
+        if bits >= 4096 {
+            group.sample_size(10);
+        }
+        let m = random_odd_ubig(bits, &mut rng);
+        let barrett = ConstDivisor::new(m.clone());
+        let monty = MontgomeryRepr::new(m.clone());
+        let a = random_ubig(bits, &mut rng);
+        let e = random_ubig(bits, &mut rng);
+        let ba = barrett.reduce(a.clone());
+        let ma = monty.reduce(a.clone());
+
+        group.bench_with_input(BenchmarkId::new("barrett", bits), &(ba, &e), |b, (a, e)| {
+            b.iter(|| a.pow(e))
+        });
+        group.bench_with_input(BenchmarkId::new("monty", bits), &(ma, &e), |b, (a, e)| {
+            b.iter(|| a.pow(e))
+        });
+    }
+
+    group.finish();
+}
+
+fn modular_inv(c: &mut Criterion) {
+    let mut rng = StdRng::seed_from_u64(SEED);
+    let mut group = c.benchmark_group("modular_inv");
+    group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    for &bits in BITS {
+        let m = random_odd_ubig(bits, &mut rng);
+        let barrett = ConstDivisor::new(m.clone());
+        let monty = MontgomeryRepr::new(m.clone());
+        // 2 is always coprime to an odd modulus, so the inverse always exists.
+        let ba = barrett.reduce(2u8);
+        let ma = monty.reduce(2u8);
+
+        group.bench_with_input(BenchmarkId::new("barrett", bits), &ba, |b, a| {
+            b.iter(|| a.clone().inv().unwrap())
+        });
+        group.bench_with_input(BenchmarkId::new("monty", bits), &ma, |b, a| {
+            b.iter(|| a.clone().inv().unwrap())
+        });
+    }
+
+    group.finish();
+}
+
+criterion_group!(
+    benches,
+    modular_mul,
+    modular_sqr,
+    modular_add,
+    modular_sub,
+    modular_pow,
+    modular_inv,
+);
+criterion_main!(benches);

integer/src/div/divide_conquer.rs

@@ -38,12 +38,14 @@ pub(crate) fn div_rem_in_place(
     fast_div_rhs_top: FastDivideNormalized2,
     memory: &mut Memory,
 ) -> bool {
-    assert!(lhs.len() > rhs.len() + div::THRESHOLD_SIMPLE && rhs.len() > div::THRESHOLD_SIMPLE);
+    assert!(
+        lhs.len() > rhs.len() + div::threshold::simple() && rhs.len() > div::threshold::simple()
+    );
 
     let mut overflow = false;
     let n = rhs.len();
     let mut m = lhs.len();
-    assert!(n > div::THRESHOLD_SIMPLE && m >= n);
+    assert!(n > div::threshold::simple() && m >= n);
     while m >= 2 * n {
         let o = div_rem_in_place_same_len(&mut lhs[m - 2 * n..m], rhs, fast_div_rhs_top, memory);
         if o {
@@ -71,9 +73,9 @@ fn div_rem_in_place_same_len(
     memory: &mut Memory,
 ) -> bool {
     let n = rhs.len();
-    assert!(n > div::THRESHOLD_SIMPLE && lhs.len() == 2 * n);
+    assert!(n > div::threshold::simple() && lhs.len() == 2 * n);
     // To guarantee n_lo >= 2.
-    const_assert!(div::THRESHOLD_SIMPLE >= 3);
+    const_assert!(div::THRESHOLD_SIMPLE_DEFAULT >= 3);
     let n_lo = n / 2;
 
     // Divide lhs[n_lo..] by rhs, putting quotient in lhs[n+n_lo..] and remainder in lhs[n_lo..n+n_lo].
@@ -109,7 +111,7 @@ fn div_rem_in_place_small_quotient(
     assert!(n >= 2 && lhs.len() >= n);
     let m = lhs.len() - n;
     assert!(m < n);
-    if m <= div::THRESHOLD_SIMPLE {
+    if m <= div::threshold::simple() {
         return div::simple::div_rem_in_place(lhs, rhs, fast_div_rhs_top);
     }
     // Use top m words of the divisor to get a quotient approximation. It may be too large by at most 2.

integer/src/div/mod.rs

@@ -12,13 +12,35 @@ use crate::{
     shift,
 };
 use alloc::alloc::Layout;
+use static_assertions::const_assert;
 
 mod divide_conquer;
 mod simple;
 pub(crate) use simple::div_rem_highest_word;
 
 /// If divisor or quotient is at most this length, use the simple division algorithm.
-const THRESHOLD_SIMPLE: usize = 32;
+const THRESHOLD_SIMPLE_DEFAULT: usize = 32;
+const_assert!(THRESHOLD_SIMPLE_DEFAULT >= simple::MIN_LEN);
+
+/// Environment-variable override for the division threshold.
+///
+/// When the `tuning` feature is active the user may set
+/// `DASHU_THRESHOLD_SIMPLE_DIV` to override the compile-time default. Values
+/// below [`simple::MIN_LEN`] are clamped to that constant.
+mod threshold {
+    #[inline]
+    pub fn simple() -> usize {
+        #[cfg(feature = "tuning")]
+        {
+            if let Ok(s) = std::env::var("DASHU_THRESHOLD_SIMPLE_DIV") {
+                if let Ok(v) = s.parse::<usize>() {
+                    return v.max(super::simple::MIN_LEN);
+                }
+            }
+        }
+        super::THRESHOLD_SIMPLE_DEFAULT
+    }
+}
 
 /// Normalize a divisor represented as words.
 ///
@@ -235,7 +257,7 @@ pub(crate) const fn fast_rem_by_normalized_dword(
 /// Memory requirement for division.
 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 {
+    if rhs_len <= threshold::simple() || lhs_len - rhs_len <= threshold::simple() {
         memory::zero_layout()
     } else {
         divide_conquer::memory_requirement_exact(lhs_len, rhs_len)
@@ -260,7 +282,7 @@ pub(crate) fn div_rem_in_place(
 ) -> bool {
     debug_assert!(lhs.len() >= rhs.len() && rhs.len() >= 2);
 
-    if rhs.len() <= THRESHOLD_SIMPLE || lhs.len() - rhs.len() <= THRESHOLD_SIMPLE {
+    if rhs.len() <= threshold::simple() || lhs.len() - rhs.len() <= threshold::simple() {
         simple::div_rem_in_place(lhs, rhs, fast_div_rhs_top)
     } else {
         divide_conquer::div_rem_in_place(lhs, rhs, fast_div_rhs_top, memory)

integer/src/div/simple.rs

@@ -9,6 +9,14 @@ use crate::{
     primitive::{double_word, highest_dword, split_dword},
 };
 
+/// Minimum supported operand length for the simple division algorithm.
+///
+/// The divide-and-conquer split (`n_lo = n / 2`) needs at least 2 low words, so the
+/// threshold that selects between simple and divide-and-conquer division must be
+/// `>= 3` (otherwise the recursion produces empty sub-problems). Values below 3 are
+/// silently clamped to this constant.
+pub const MIN_LEN: usize = 3;
+
 /// Division in place using the simple algorithm.
 ///
 /// Divide lhs by rhs, replacing the top words of lhs by the quotient and the