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.

[cuviper]

I want to first acknowledge this with thanks -- it's quite impressive to be near GMP performance!

However to set expectations: this is also a large PR, and I will need some studying to understand what's going on, so it may take me a while to find time to review this.

[byeongkeunahn]

Thanks. I'd like to note a few changes in the latest commit, which brings 10-15% performance gains and reduced memory footprint:

• The updated code uses decimation-in-frequency (DIF) in forward transform and decimation-in-time (DIT) in inverse transform.
• The prime numbers are replaced with bigger primes to enable tighter bit-packing in some cases that occur periodically.
  • P1 = 14_259_017_916_245_606_401, Max NTT length = 2^22 * 3^21 * 5^2 = 1_096_847_532_018_892_800
  • P2 = 17_984_575_660_032_000_001, Max NTT length = 2^19 * 3^17 * 5^6 = 1_057_916_215_296_000_000
  • P3 = 17_995_154_822_184_960_001, Max NTT length = 2^17 * 3^22 * 5^4 = 2_570_736_403_169_280_000
  • P2 and P3 are used for the two-prime NTT, whereas the three-prime NTT uses all three.

[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...

[byeongkeunahn]

I’ve made the suggested improvements:

• Cleaned up the commit history.
  The original branch remains available here for reference.
• Rebased the branch onto the latest main.
• Replaced the unsafe routines with safe equivalents. Benchmarks showed no performance regression.
• Added comments explaining why the NTT routine requires its own Montgomery arithmetic routines.

I also made a small improvement to the planner based on exhaustive Criterion.rs benchmarks covering operand sizes of up to 300 million bits each. This reduced execution time by 1.2% on average.

Some of this work was completed with assistance from ChatGPT Codex.

Thanks again for the detailed review. Please let me know if there are any remaining issues.