Implement number theoretic transform for large integer multiplication

[byeongkeunahn]

This commit implements number theoretic transform (NTT) for large integer multiplication (issue #169).

• To simplify implementation the Schönhage–Strassen algorithm was not used. Instead, three distinct 64-bit primes were carefully chosen to enable NTT up to ~10^18 64bit integers, which allows multiplication up to ~5 x 10^17 64bit integers. Depending on the input length either two or three primes are used, with the latter only used when the inputs consist of at least 2^40 64bit integers. The convolution results modulo primes are merged using the Chinese Remainder Theorem.
• To reduce padding and the number of cycles for NTT, multiple radices are used (radix-2, radix-3, radix-4, radix-5, radix-6). Radix-8 is desirable but actually slower, presumably due to register spill. Although manual SIMD coding may alleviate this issue, it was not used (i) for maximum portability and (ii) since 64bit SIMD multiply is not widely available even on x86/x64 platforms (AVX512). The prime moduli are carefully chosen to support these radices.
  • The NTT length is selected by exhaustively evaluating cost estimates for all allowed lengths within a factor of two.
• Single-word (u64) Montgomery reduction is used for fast modular multiplication.
• 32bit digits are supported by repacking the digits into u64, running the u64 algorithm, and converting back to u32. This results in 32bit builds being about 3-5x slower compared to 64bit builds, which, however, still is an improvement upon the existing algorithms.
• Unbalanced multiplication is enabled when the cost estimates are favorable.
• Based on experimentation, the following thresholds are chosen:
  • For u64 digits, switch to NTT if the shorter integer has at least 512 digits.
  • For u32 digits, switch to NTT if the shorter integer has at least 2,048 digits.
• The three primes are as follows:
  • P1 = 10_237_243_632_176_332_801, Max NTT length = 2^24 * 3^20 * 5^2 = 1_462_463_376_025_190_400
  • P2 = 13_649_658_176_235_110_401, Max NTT length = 2^26 * 3^19 * 5^2 = 1_949_951_168_033_587_200
  • P3 = 14_259_017_916_245_606_401, Max NTT length = 2^22 * 3^21 * 5^2 = 1_096_847_532_018_892_800
  • P1 and P2 are used for the two-prime NTT, whereas the three-prime NTT uses all three.
• The following MIT-licensed projects are used as reference:
  • http://wwwa.pikara.ne.jp/okojisan/otfft-en/index.html
  • https://github.com/Bubbler-4/cg263748

On Ryzen 7 2700X, 64bit, it takes about 15ms for 2.7Mbits x 2.7Mbits and 170ms for 27Mbits x 27Mbits multiplication. This seems comparable to GMP 6.2.1.

[hkalbasi]

Does the chart shows that current algorithm is faster than GMP? That's impressive.

[hkalbasi]

I ran benchmark fib_hex 100m from https://github.com/tczajka/bigint-benchmark-rs on this PR and it made num-bigint twice faster than malachite, slightly faster than gmp and 12x faster than itself.

[cuviper]

Obviously the performance here is impressive, as I said before. However, as I finally started to try reviewing this, I have a few high level objections.

• The long commit history is not something I would want to merge. If there are discrete logical steps in implementing this feature, we can preserve some of that, but not a bunch of "update", "improve", and merges.
• The Montgomery implementation should share src/biguint/monty.rs, extending that if needed.
• So much unsafe code is uncharacteristic for this crate, and I would much rather stick to safe code as much as possible, even if that means sacrificing some performance. Case in point, miri finds UB here:

test test_mul_overflow_2 ... error: Undefined Behavior: attempting a write access using <5072044> at alloc1162269[0x0], but that tag does not exist in the borrow stack for this location
   --> src/biguint/ntt.rs:230:13
    |
230 |             *out.add(i) = Arith::<P>::mreduce(v);
    |             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ this error occurs as part of an access at alloc1162269[0x0..0x8]
    |
    = help: this indicates a potential bug in the program: it performed an invalid operation, but the Stacked Borrows rules it violated are still experimental
    = help: see https://github.com/rust-lang/unsafe-code-guidelines/blob/master/wip/stacked-borrows.md for further information
help: <5072044> was created by a SharedReadWrite retag at offsets [0x30..0x4830]
   --> src/biguint/ntt.rs:518:27
    |
518 |         conv_base::<P>(g, (&mut x[i..]).as_mut_ptr(), (&mut y[i..]).as_mut_ptr(), tf_current);
    |                           ^^^^^^^^^^^^^^^^^^^^^^^^^^
    = note: this is on thread `test_mul_overfl`
    = note: stack backtrace:
            0: num_bigint::biguint::ntt::conv_base::<17984575660032000001>
                at src/biguint/ntt.rs:230:13: 230:49
            1: num_bigint::biguint::ntt::conv::<17984575660032000001>
                at src/biguint/ntt.rs:518:9: 518:94
            2: num_bigint::biguint::ntt::mac3_two_primes
                at src/biguint/ntt.rs:604:5: 604:98
            3: num_bigint::biguint::ntt::mac3_u64
                at src/biguint/ntt.rs:760:9: 760:41
            4: num_bigint::biguint::ntt::mac3
                at src/biguint/ntt.rs:772:5: 772:24
            5: num_bigint::biguint::multiplication::mac3
                at src/biguint/multiplication.rs:358:9: 358:29
            6: num_bigint::biguint::multiplication::mul3
                at src/biguint/multiplication.rs:366:5: 366:31
            7: num_bigint::biguint::multiplication::<impl std::ops::Mul for &num_bigint::BigUint>::mul
                at src/biguint/multiplication.rs:428:31: 428:41
            8: test_mul_overflow_2
                at tests/biguint.rs:1075:20: 1075:27
            9: test_mul_overflow_2::{closure#0}
                at tests/biguint.rs:1068:25: 1068:25

(that particular one does pass under tree-borrows though)

[byeongkeunahn]

Thanks for your review and comments! I’ll try to address the issues, including removing the unsafe code and cleaning up the commit history, even if that comes at the cost of a small performance regression.

Integrating the Montgomery reduction code with src/monty.rs is a bit more complicated, though:

• The NTT code currently uses a 64-bit modulus regardless of the target architecture, whereas src/monty.rs uses BigDigit. Introducing generics could potentially address this mismatch.
• The Montgomery operations in this PR were hand-optimized by inspecting the LLVM output on Godbolt and benchmarking promising variants. They are tailored specifically to the operations used by the NTT routine, including fused operations and the careful omission of overflow checks where they are mathematically unnecessary. As a result, they serve a somewhat different purpose from the code in src/monty.rs, and I’m not yet sure how they could be cleanly integrated there.

[cuviper]

We can treat the monty-consolidation as "nice to have". However, any specific tailoring needs comments, lest that work be undone by a later contributor or maintainer. I'm skeptical that those optimizations wouldn't be useful to the other monty use as well though...