Implement number theoretic transform for large integer multiplication

benches/bigint.rs

@@ -87,6 +87,21 @@ fn multiply_3(b: &mut Bencher) {
     multiply_bench(b, 1 << 16, 1 << 17);
 }
 
+#[bench]
+fn multiply_4(b: &mut Bencher) {
+    multiply_bench(b, 100_000, 1_003_741);
+}
+
+#[bench]
+fn multiply_5(b: &mut Bencher) {
+    multiply_bench(b, 2_718_328, 2_738_633);
+}
+
+#[bench]
+fn multiply_6(b: &mut Bencher) {
+    multiply_bench(b, 27_183_279, 27_386_321);
+}
+
 #[bench]
 fn divide_0(b: &mut Bencher) {
     divide_bench(b, 1 << 8, 1 << 6);

benches/factorial.rs

@@ -16,6 +16,36 @@ fn factorial_mul_biguint(b: &mut Bencher) {
     });
 }
 
+fn factorial_product(l: usize, r: usize) -> BigUint {
+    if l >= r {
+        BigUint::from(l)
+    } else {
+        let m = (l+r)/2;
+        factorial_product(l, m) * factorial_product(m+1, r)
+    }
+}
+
+#[bench]
+fn factorial_mul_biguint_dnc_10k(b: &mut Bencher) {
+    b.iter(|| {
+        factorial_product(1, 10_000);
+    });
+}
+
+#[bench]
+fn factorial_mul_biguint_dnc_100k(b: &mut Bencher) {
+    b.iter(|| {
+        factorial_product(1, 100_000);
+    });
+}
+
+#[bench]
+fn factorial_mul_biguint_dnc_300k(b: &mut Bencher) {
+    b.iter(|| {
+        factorial_product(1, 300_000);
+    });
+}
+
 #[bench]
 fn factorial_mul_u32(b: &mut Bencher) {
     b.iter(|| (1u32..1000).fold(BigUint::one(), Mul::mul));

src/biguint.rs

@@ -22,6 +22,7 @@ mod bits;
 mod convert;
 mod iter;
 mod monty;
+mod ntt;
 mod power;
 mod shift;
 

src/biguint/multiplication.rs

@@ -13,6 +13,8 @@ use core::iter::Product;
 use core::ops::{Mul, MulAssign};
 use num_traits::{CheckedMul, FromPrimitive, One, Zero};
 
+use super::ntt;
+
 #[inline]
 pub(super) fn mac_with_carry(
     a: BigDigit,
@@ -97,7 +99,7 @@ fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {
     //   number of operations, but uses more temporary allocations.
     //
     // The thresholds are somewhat arbitrary, chosen by evaluating the results
-    // of `cargo bench --bench bigint multiply`.
+    // of `cargo bench --bench bigint multiply --features rand`.
 
     if x.len() <= 32 {
         // Long multiplication:
@@ -217,7 +219,7 @@ fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {
             }
             NoSign => (),
         }
-    } else {
+    } else if x.len() <= if cfg!(u64_digit) { 512 } else { 2048 } {
         // Toom-3 multiplication:
         //
         // Toom-3 is like Karatsuba above, but dividing the inputs into three parts.
@@ -346,6 +348,14 @@ fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {
                 NoSign => {}
             }
         }
+    } else {
+        // Number-theoretic transform (NTT) multiplication:
+        //
+        // NTT multiplies two integers by computing the convolution of the arrays
+        // modulo a prime. Since the result may exceed the prime, we use two or three
+        // distinct primes and combine the results using the Chinese Remainder
+        // Theroem (CRT).
+        ntt::mac3(acc, b, c);
     }
 }