Implement number theoretic transform for large integer multiplication by byeongkeunahn · Pull Request #282 · rust-num/num-bigint

byeongkeunahn · 2023-08-28T14:02:31Z

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 · 2023-09-18T22:28:39Z

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 · 2023-09-19T04:40:43Z

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 · 2024-12-13T08:33:02Z

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

hkalbasi · 2024-12-13T09:09:01Z

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 · 2026-07-02T17:54:48Z

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 · 2026-07-10T00:59:25Z

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 · 2026-07-10T22:13:44Z

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 · 2026-08-02T14:34:56Z

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.

cuviper · 2026-08-05T16:43:44Z

+#![allow(clippy::too_many_arguments)]
+#![allow(clippy::similar_names)]
+
+use crate::biguint::Vec;


It's odd to use an absolute path when it's not an intended re-export.

Suggested change

use crate::biguint::Vec;

use super::Vec;

cuviper · 2026-08-05T16:45:21Z

+    }
+    // Modular inverse: a^-1 mod modulus


Style nit here and throughout: please add a blank line between functions/consts/etc.

Even within code blocks, if there's a comment explaining the following code, then also add a blank line before. (except if it's the first part at a new indentation level)

cuviper · 2026-08-05T17:30:16Z

+            for m5 in 0..=Arith::<P>::factors(5) {
+                for m3 in 0..=Arith::<P>::factors(3) {


Might it be worth forcing const-evaluation on these? i.e. using a const { ... } block, or perhaps adding Arith::FACTORS_5 etc.

cuviper · 2026-08-05T17:31:32Z

+                        break;
+                    }
+                    let (mut len, mut m2) = (len as usize, 0);
+                    while len < min_len && m2 < Arith::<P>::factors(2) {


Also FWIW, factors(2) is just (P-1).trailing_zeros() -- but it doesn't really matter either way if we're doing it at compile time.

cuviper · 2026-08-05T17:36:29Z

+                        let (mut tmp, mut cost) = (len, 0);
+                        let mut g_new = 1;
+
+                        // Length-dependent weights for cost estimation.


These weights and cost formulas look like a whole lot of "magic" constants. I suppose they come from your PR statement:

The NTT length is selected by exhaustively evaluating cost estimates for all allowed lengths within a factor of two.

It would be good to have comments directly in the code to explain this, ideally with some instructions how one might reevaluate these costs. If there's any third-party source as well, please link to that.

cuviper · 2026-08-05T18:04:58Z

+#![allow(clippy::many_single_char_names)]
+#![allow(clippy::needless_range_loop)]
+#![allow(clippy::too_many_arguments)]
+#![allow(clippy::similar_names)]


I'm not interested in maintaining a bunch of pedantic clippy suppressions. AFAICS, the only lint here that we hit by default is too_many_arguments, and that's only on ntt5_kernel and ntt6_kernel, so let's just allow specifically on those.

For too-many-args, it might make sense to switch those to arrays, e.g.

const fn ntt5_kernel<const P: u64, const INV: bool, const TWIDDLE: bool>( [w1, w2, w3, w4]: [u64; 4], [a, mut b, mut c, mut d, mut e]: [u64; 5], ) -> [u64; 5] {

... but I haven't checked whether that affects benchmarks.

cuviper · 2026-08-05T18:09:07Z

+    assert!(!x.is_empty() && x.len() == y.len());
+    let (_n, g, m, last_radix) = (plan.n, plan.g, plan.m, plan.last_radix as u64);
+
+    /* multiply by a constant in advance */


Another style nit: please use // comments throughout, rather than /* */

cuviper · 2026-08-05T18:10:29Z

+// Propagates carry from the beginning to the end of acc,
+//   and returns the resulting carry if it is nonzero.


Suggested change

// Propagates carry from the beginning to the end of acc,

// and returns the resulting carry if it is nonzero.

// Propagates carry from the beginning to the end of acc,

// and returns the resulting carry if it is nonzero.

cuviper · 2026-08-05T18:11:33Z

+    // Computes c as u128 * mreduce(v) as u128,
+    //   using d: u64 = mmulmod(P-1, c).


Suggested change

// Computes c as u128 * mreduce(v) as u128,

// using d: u64 = mmulmod(P-1, c).

// Computes c as u128 * mreduce(v) as u128, using d: u64 = mmulmod(P-1, c).

//

cuviper · 2026-08-05T18:13:04Z

+    // Process remaining carries. The addition carry_acc + bitbuf should not overflow
+    //   since bitbuf is underfilled and carry_acc is always 0 or 1.


Suggested change

// Process remaining carries. The addition carry_acc + bitbuf should not overflow

// since bitbuf is underfilled and carry_acc is always 0 or 1.

// Process remaining carries. The addition carry_acc + bitbuf should not overflow

// since bitbuf is underfilled and carry_acc is always 0 or 1.

byeongkeunahn changed the title ~~Implement number theroetic transform for large integer multiplication~~ Implement number theoretic transform for large integer multiplication Aug 28, 2023

byeongkeunahn force-pushed the master branch from eae5b13 to 66faa24 Compare September 19, 2023 04:24

byeongkeunahn force-pushed the master branch 2 times, most recently from 2639746 to 0e41192 Compare September 19, 2023 05:09

hkalbasi mentioned this pull request Dec 13, 2024

Convert to base 10 is significantly slower than similar libraries #315

Open

byeongkeunahn closed this Jul 27, 2026

byeongkeunahn force-pushed the master branch from 2291466 to 9ec740f Compare July 27, 2026 09:35

byeongkeunahn reopened this Aug 2, 2026

Add multiplication with number-theoretic transform

3302555

byeongkeunahn force-pushed the master branch from 15eb95c to 3302555 Compare August 2, 2026 14:38

cuviper requested changes Aug 5, 2026

View reviewed changes

		for m5 in 0..=Arith::<P>::factors(5) {
		for m3 in 0..=Arith::<P>::factors(3) {

		// Propagates carry from the beginning to the end of acc,
		// and returns the resulting carry if it is nonzero.

		// Computes c as u128 * mreduce(v) as u128,
		// using d: u64 = mmulmod(P-1, c).

		// Process remaining carries. The addition carry_acc + bitbuf should not overflow
		// since bitbuf is underfilled and carry_acc is always 0 or 1.

Uh oh!

Conversation

byeongkeunahn commented Aug 28, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

cuviper commented Sep 18, 2023

Uh oh!

byeongkeunahn commented Sep 19, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hkalbasi commented Dec 13, 2024

Uh oh!

hkalbasi commented Dec 13, 2024

Uh oh!

cuviper commented Jul 2, 2026

Uh oh!

byeongkeunahn commented Jul 10, 2026

Uh oh!

cuviper commented Jul 10, 2026

Uh oh!

byeongkeunahn commented Aug 2, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

byeongkeunahn commented Aug 28, 2023 •

edited

Loading

byeongkeunahn commented Sep 19, 2023 •

edited

Loading