Various improvements on integer multiplication (including NTT)

.github/workflows/tests.yml

@@ -162,4 +162,9 @@ jobs:
         with:
           toolchain: stable
           components: clippy
-      - run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings
+      - name: Clippy (default / 64-bit Word)
+        run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings
+      - name: Clippy (32-bit Word)
+        env:
+          RUSTFLAGS: --cfg force_bits="32"
+        run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings

Cargo.toml

@@ -36,6 +36,7 @@ std = ["dashu-base/std", "dashu-int/std", "dashu-float/std", "dashu-ratio/std"]
 # stable features
 serde = ["dashu-int/serde", "dashu-float/serde", "dashu-ratio/serde"]
 num-order = ["dashu-int/num-order", "dashu-float/num-order", "dashu-ratio/num-order"]
+tuning = ["dashu-int/tuning"]
 zeroize = ["dashu-int/zeroize", "dashu-float/zeroize", "dashu-ratio/zeroize"]
 
 # unstable features

float/src/parse.rs

@@ -27,6 +27,12 @@ impl<const B: Word> Repr<B> {
     pub fn from_str_native(mut src: &str) -> Result<(Self, usize), ParseError> {
         assert!(MIN_RADIX as Word <= B && B <= MAX_RADIX as Word);
 
+        // B is guaranteed to be in 2..=36 by the assert above; the cast to u32
+        // is needed because `from_str_radix` takes a u32 radix. On 32-bit Word
+        // targets the cast is a no-op.
+        #[allow(clippy::unnecessary_cast)]
+        let radix: u32 = B as u32;
+
         // parse and remove the sign
         let sign = match src.strip_prefix('-') {
             Some(s) => {
@@ -100,14 +106,14 @@ impl<const B: Word> Repr<B> {
                     return Err(ParseError::UnsupportedRadix);
                 } else {
                     let digits = int_str.len() - int_str.matches('_').count();
-                    (UBig::from_str_radix(&src[..dot], B as u32)?, digits, B as u32)
+                    (UBig::from_str_radix(&src[..dot], radix)?, digits, radix)
                 }
             } else {
                 if pmarker {
                     // prefix is required for using `p` as scale marker
                     return Err(ParseError::UnsupportedRadix);
                 }
-                (UBig::ZERO, 0, B as u32)
+                (UBig::ZERO, 0, radix)
             };
 
             // parse fractional part
@@ -139,7 +145,7 @@ impl<const B: Word> Repr<B> {
                 return Err(ParseError::UnsupportedRadix);
             } else {
                 ndigits = src.len() - src.matches('_').count();
-                UBig::from_str_radix(src, B as u32)?
+                UBig::from_str_radix(src, radix)?
             }
         };
 

float/src/third_party/num_traits.rs

@@ -133,7 +133,7 @@ impl<R: Round, const B: Word> num_traits::Num for FBig<R, B> {
     #[inline]
     fn from_str_radix(s: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
         // the conversion might a fail with 16-bit words.
-        #[allow(clippy::unnecessary_fallible_conversions)]
+        #[allow(clippy::unnecessary_fallible_conversions, clippy::useless_conversion)]
         let r: Word = radix.try_into().map_err(|_| ParseError::UnsupportedRadix)?;
         if r == B {
             #[allow(deprecated)] // TODO(v0.5): remove after from_str_native is made private.

integer/CHANGELOG.md

@@ -3,15 +3,29 @@
 ## Unreleased
 
 ### 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.
 - `UBig::from_u64` and `IBig::from_i64`, const on 32-bit and 64-bit targets.
 
 ### Improve
 - Basecase (schoolbook) multiplication now uses an dword mult inner kernel (two multiplier words per sweep over the accumulator, mirroring GMP's `mpn_addmul_2` and `mpn_submul_2`), roughly halving accumulator memory traffic.
 - Addition and subtraction carry/borrow propagation now uses `Word` (u64/u32) instead of `bool` throughout the architecture-specific `add_with_carry` and `sub_with_borrow` functions, eliminating `bool`↔Word conversions in the inner loops.
+- 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).
 
-### Improve
-- Logarithm for very large values uses power-sequence decomposition, replacing iterative single-step multiplication.
-- Improve power-of-two base formatting ([#3](https://github.com/cmpute/dashu/pull/3))
+### Change
+- NTT multiplication now uses Proth primes (`K·2^N + 1`) instead of Solinas primes, improving modular reduction speed.
+- NTT threshold lowered from 40 000 to 4 000 words.
+- NTT enabled for 32-bit Word targets.
+- Arch-specific NTT prime definitions under `arch/generic_{32,64}_bit/ntt.rs`.
+
+### Fix
+- `pack.rs` test used 64-bit literals that overflowed `Word` (`u32`) on 32-bit targets, breaking the test build.
+- `pack.rs` now uses native `Word`/`Lane` types throughout instead of `u64`/`u32`, fixing clippy `unnecessary_cast` warnings on 64-bit.
+- `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.
 
 ## 0.4.2
 

integer/Cargo.toml

@@ -19,6 +19,7 @@ all-features = true
 [features]
 default = ["std", "num-order"]
 std = ["dashu-base/std"]
+tuning = ["std"]
 
 # unstable dependencies
 rand = ["rand_v08"]
@@ -30,7 +31,7 @@ dashu-base = { version = "0.4.1", default-features = false, path = "../base" }
 cfg-if = { version = "1.0.0" }
 static_assertions = { version = "1.1" }
 rustversion = { version = "1.0.0" }
-num-modular = { version = "0.6.1" }
+num-modular = { version = "0.6.4" }
 
 # stable dependencies
 num-order = { optional = true, version = "1.2.0", default-features = false }

integer/benches/primitive.rs

@@ -150,6 +150,31 @@ fn ubig_ilog_large(criterion: &mut Criterion) {
     group.finish();
 }
 
+fn ubig_mul_asymmetric(criterion: &mut Criterion) {
+    let mut rng = StdRng::seed_from_u64(SEED);
+    let mut group = criterion.benchmark_group("ubig_mul_asymmetric");
+    group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
+
+    // b just above the NTT threshold (4 000 words = 256 kbits → use 500 kbits).
+    let b_bits = 500_000;
+    let b = random_ubig(b_bits, &mut rng);
+
+    // a ranges from 1 kbit (below Karatsuba threshold) to heavily
+    // asymmetric (10×), exercising all chunked-mul code paths.
+    for &a_bits in &[
+        1_000, 10_000, 100_000, 500_000, 1_000_000, 2_000_000, 5_000_000,
+    ] {
+        let a = random_ubig(a_bits, &mut rng);
+        group.bench_with_input(
+            BenchmarkId::from_parameter(format!("{a_bits}/{b_bits}")),
+            &(a, &b),
+            |bencher, (ta, tb)| bencher.iter(|| ta * *tb),
+        );
+    }
+
+    group.finish();
+}
+
 criterion_group!(
     benches,
     ubig_add,
@@ -163,6 +188,7 @@ criterion_group!(
     ubig_modulo_pow,
     ubig_pow_large_base,
     ubig_ilog_large,
+    ubig_mul_asymmetric,
 );
 
 criterion_main!(benches);

integer/src/arch/generic_32_bit/mod.rs

@@ -4,4 +4,5 @@ pub(crate) mod add;
 #[path = "../generic/digits.rs"]
 pub(crate) mod digits;
 
+pub(crate) mod ntt;
 pub(crate) mod word;

integer/src/arch/generic_32_bit/ntt.rs

@@ -0,0 +1,84 @@
+//! NTT primes and constants for 32-bit Word targets.
+//!
+//! Uses Proth primes of the form `K * 2^N + 1`.
+//! All constants computed by `integer/src/mul/ntt/compute_constants.py`.
+
+use num_modular::FixedProth32;
+
+// Proth reducer instances — each with a different (N, K) pair.
+pub type Rp0 = FixedProth32<26, 7>;
+pub type Rp1 = FixedProth32<27, 15>;
+pub type Rp2 = FixedProth32<27, 17>;
+
+pub const P0: Rp0 = FixedProth32::<26, 7>;
+pub const P1: Rp1 = FixedProth32::<27, 15>;
+pub const P2: Rp2 = FixedProth32::<27, 17>;
+
+pub const K: usize = 3;
+pub const MAX_LOG_N: u32 = 26;
+pub const B_PACK_MIN: u32 = 8;
+pub const B_PACK_CANDIDATES: &[u32] = &[32, 16, 8];
+
+pub type Lane = u32;
+
+/// Primitive `MAX_LOG_N`-th roots of unity for each prime.
+pub const OMEGA_MAX: [Lane; K] = [
+    0x0000088b, // P0
+    0x3a26eef8, // P1
+    0x1aa0ab5e, // P2
+];
+
+pub const CRT_INV_IJ: [[Lane; K]; K] = [[0, 0x4e42c85b, 0x5fb425ef], [0, 0, 0x44000009], [0, 0, 0]];
+
+/// Prime moduli indexed by PI.
+pub const MODULI: [Lane; K] = [Rp0::MODULUS, Rp1::MODULUS, Rp2::MODULUS];
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use num_modular::Reducer;
+
+    type ReducerFns = (fn(Lane) -> Lane, fn(Lane) -> Lane, fn(Lane) -> Lane);
+
+    #[test]
+    fn test_primes_proth_form() {
+        assert_eq!(MODULI[0], 7u32 * (1u32 << 26) + 1);
+        assert_eq!(MODULI[1], 15u32 * (1u32 << 27) + 1);
+        assert_eq!(MODULI[2], 17u32 * (1u32 << 27) + 1);
+    }
+
+    #[test]
+    fn test_primes_v2() {
+        for &p in &MODULI {
+            let v2 = (p - 1).trailing_zeros();
+            assert!(v2 >= MAX_LOG_N, "v2(p-1) = {v2} < MAX_LOG_N");
+        }
+    }
+
+    #[test]
+    fn test_omega_order() {
+        for (pi, &omega_max) in OMEGA_MAX.iter().enumerate() {
+            let p = MODULI[pi];
+            let (sqr, to_m, from_m): ReducerFns = match pi {
+                0 => {
+                    (|w| P0.reduce((w as u64) * (w as u64)), |v| P0.transform(v), |v| P0.residue(v))
+                }
+                1 => {
+                    (|w| P1.reduce((w as u64) * (w as u64)), |v| P1.transform(v), |v| P1.residue(v))
+                }
+                2 => {
+                    (|w| P2.reduce((w as u64) * (w as u64)), |v| P2.transform(v), |v| P2.residue(v))
+                }
+                _ => unreachable!(),
+            };
+
+            let mut w = to_m(omega_max);
+            for _ in 0..MAX_LOG_N - 1 {
+                w = sqr(w);
+            }
+            assert_eq!(from_m(w), p - 1, "omega^(2^(MAX_LOG_N-1)) != -1 mod p for prime {pi}");
+            w = sqr(w);
+            assert_eq!(from_m(w), 1, "omega^(2^MAX_LOG_N) != 1 mod p for prime {pi}");
+        }
+    }
+}

integer/src/arch/generic_32_bit/word.rs

@@ -9,3 +9,7 @@ pub type DoubleWord = u64;
 
 /// Signed double machine word.
 pub type SignedDoubleWord = i64;
+
+/// Accumulator for the product of three primes (3 × 2^32 ≈ 2^96).
+#[derive(Clone, Copy, Debug, Default)]
+pub struct TripleWord(pub [u32; 3]);

integer/src/arch/generic_64_bit/mod.rs

@@ -4,4 +4,5 @@ pub(crate) mod add;
 #[path = "../generic/digits.rs"]
 pub(crate) mod digits;
 
+pub(crate) mod ntt;
 pub(crate) mod word;