Implement number theoretic transform for large integer multiplication

benches/bigint.rs

@@ -67,33 +67,48 @@ fn fib2(n: usize) -> BigUint {
 }
 
 #[bench]
-fn multiply_0(b: &mut Bencher) {
+fn multiply_8_8(b: &mut Bencher) {
     multiply_bench(b, 1 << 8, 1 << 8);
 }
 
 #[bench]
-fn multiply_1(b: &mut Bencher) {
+fn multiply_8_16(b: &mut Bencher) {
     multiply_bench(b, 1 << 8, 1 << 16);
 }
 
 #[bench]
-fn multiply_2(b: &mut Bencher) {
-    multiply_bench(b, 1 << 16, 1 << 16);
+fn multiply_12_13(b: &mut Bencher) {
+    multiply_bench(b, 1 << 12, 1 << 13);
 }
 
 #[bench]
-fn multiply_3(b: &mut Bencher) {
-    multiply_bench(b, 1 << 16, 1 << 17);
+fn multiply_12_14(b: &mut Bencher) {
+    multiply_bench(b, 1 << 12, 1 << 14);
 }
 
 #[bench]
-fn multiply_4(b: &mut Bencher) {
-    multiply_bench(b, 1 << 12, 1 << 13);
+fn multiply_13_13(b: &mut Bencher) {
+    multiply_bench(b, 1 << 13, 1 << 13);
 }
 
 #[bench]
-fn multiply_5(b: &mut Bencher) {
-    multiply_bench(b, 1 << 12, 1 << 14);
+fn multiply_14_14(b: &mut Bencher) {
+    multiply_bench(b, 1 << 14, 1 << 14);
+}
+
+#[bench]
+fn multiply_15_15(b: &mut Bencher) {
+    multiply_bench(b, 1 << 15, 1 << 15);
+}
+
+#[bench]
+fn multiply_16_16(b: &mut Bencher) {
+    multiply_bench(b, 1 << 16, 1 << 16);
+}
+
+#[bench]
+fn multiply_16_17(b: &mut Bencher) {
+    multiply_bench(b, 1 << 16, 1 << 17);
 }
 
 #[bench]

src/biguint.rs

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

src/biguint/multiplication.rs

@@ -11,6 +11,8 @@ use core::iter::Product;
 use core::ops::{Mul, MulAssign};
 use num_traits::{CheckedMul, FromPrimitive, Zero};
 
+use super::ntt;
+
 #[inline]
 pub(super) fn mac_with_carry(
     a: BigDigit,
@@ -93,8 +95,9 @@ fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {
     // - If y is at least least twice as long as x, split using Half-Karatsuba.
     // - Next we use Karatsuba multiplication (Toom-2), which we have optimized
     //   to avoid unnecessary allocations for intermediate values.
-    // - For the largest inputs we use Toom-3, which better optimizes the
+    // - Next we use Toom-3, which better optimizes the
     //   number of operations, but uses more temporary allocations.
+    // - For the largest inputs we use number-theoretic transform (NTT).
     //
     // The thresholds are somewhat arbitrary, chosen by evaluating the results
     // of `cargo bench --bench bigint multiply`.
@@ -278,7 +281,7 @@ fn mac3(mut acc: &mut [BigDigit], mut b: &[BigDigit], mut c: &[BigDigit]) {
             }
             NoSign => (),
         }
-    } else {
+    } else if x.len() <= 512 {
         // Toom-3 multiplication:
         //
         // Toom-3 is like Karatsuba above, but dividing the inputs into three parts.
@@ -407,6 +410,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
+        // Theorem (CRT).
+        ntt::mac3(acc, b, c);
     }
 }