Skip to content

Commit 70cb8dd

Browse files
cmputeJacob Zhongclaude
authored
Various improvements on integer multiplication (including NTT) (#74)
* Add plan for NTT implementation * WIP: implemented ntt mul * WIP: tidy up * WIP: some param tuning * WIP: further tune b_pack * WIP: minor test improvements Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Some minor improvements Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Fix CI * Fix CI again * Change solinas to proth for NTT * Tidy up num-modular usage * Tune the NTT threshold * Tidy up * Remove todos * Fix 32-bit and clippy CI failures in NTT module Use native Word/Lane types throughout pack.rs (instead of u64/u32) so the same source compiles cleanly on both 32-bit and 64-bit targets. Resolve the remaining clippy warnings (unnecessary_cast, let_and_return, too_many_arguments, needless_range_loop, type_complexity) that were failing the -D warnings CI run. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * CI: run clippy on 32-bit Word target Adds a second clippy step to the existing Clippy job that runs with --cfg force_bits="32", catching Word-width-dependent lints like the unnecessary_cast warnings recently fixed in the NTT module. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Fix clippy warnings on 32-bit Word in dashu-float Extract radix once in from_str_native to avoid repeated B as u32 casts, and silence the identity try_into() in num_traits::Num::from_str_radix on 32-bit Word targets. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Remove development-only NTT tests Tests removed are fully subsumed by the schoolbook-comparison and roundtrip tests added later in development. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Fix fmt --------- Co-authored-by: Jacob Zhong <jacob@rimbot.com> Co-authored-by: Claude Opus 4.7 <noreply@anthropic.com>
1 parent a4d1cc0 commit 70cb8dd

23 files changed

Lines changed: 2054 additions & 94 deletions

File tree

.github/workflows/tests.yml

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -162,4 +162,9 @@ jobs:
162162
with:
163163
toolchain: stable
164164
components: clippy
165-
- run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings
165+
- name: Clippy (default / 64-bit Word)
166+
run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings
167+
- name: Clippy (32-bit Word)
168+
env:
169+
RUSTFLAGS: --cfg force_bits="32"
170+
run: cargo clippy --all-features --all-targets --workspace --exclude dashu-python -- -D warnings

Cargo.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@ std = ["dashu-base/std", "dashu-int/std", "dashu-float/std", "dashu-ratio/std"]
3636
# stable features
3737
serde = ["dashu-int/serde", "dashu-float/serde", "dashu-ratio/serde"]
3838
num-order = ["dashu-int/num-order", "dashu-float/num-order", "dashu-ratio/num-order"]
39+
tuning = ["dashu-int/tuning"]
3940
zeroize = ["dashu-int/zeroize", "dashu-float/zeroize", "dashu-ratio/zeroize"]
4041

4142
# unstable features

TODO.md

Lines changed: 0 additions & 69 deletions
This file was deleted.

float/src/parse.rs

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,12 @@ impl<const B: Word> Repr<B> {
2727
pub fn from_str_native(mut src: &str) -> Result<(Self, usize), ParseError> {
2828
assert!(MIN_RADIX as Word <= B && B <= MAX_RADIX as Word);
2929

30+
// B is guaranteed to be in 2..=36 by the assert above; the cast to u32
31+
// is needed because `from_str_radix` takes a u32 radix. On 32-bit Word
32+
// targets the cast is a no-op.
33+
#[allow(clippy::unnecessary_cast)]
34+
let radix: u32 = B as u32;
35+
3036
// parse and remove the sign
3137
let sign = match src.strip_prefix('-') {
3238
Some(s) => {
@@ -100,14 +106,14 @@ impl<const B: Word> Repr<B> {
100106
return Err(ParseError::UnsupportedRadix);
101107
} else {
102108
let digits = int_str.len() - int_str.matches('_').count();
103-
(UBig::from_str_radix(&src[..dot], B as u32)?, digits, B as u32)
109+
(UBig::from_str_radix(&src[..dot], radix)?, digits, radix)
104110
}
105111
} else {
106112
if pmarker {
107113
// prefix is required for using `p` as scale marker
108114
return Err(ParseError::UnsupportedRadix);
109115
}
110-
(UBig::ZERO, 0, B as u32)
116+
(UBig::ZERO, 0, radix)
111117
};
112118

113119
// parse fractional part
@@ -139,7 +145,7 @@ impl<const B: Word> Repr<B> {
139145
return Err(ParseError::UnsupportedRadix);
140146
} else {
141147
ndigits = src.len() - src.matches('_').count();
142-
UBig::from_str_radix(src, B as u32)?
148+
UBig::from_str_radix(src, radix)?
143149
}
144150
};
145151

float/src/third_party/num_traits.rs

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -133,7 +133,7 @@ impl<R: Round, const B: Word> num_traits::Num for FBig<R, B> {
133133
#[inline]
134134
fn from_str_radix(s: &str, radix: u32) -> Result<Self, Self::FromStrRadixErr> {
135135
// the conversion might a fail with 16-bit words.
136-
#[allow(clippy::unnecessary_fallible_conversions)]
136+
#[allow(clippy::unnecessary_fallible_conversions, clippy::useless_conversion)]
137137
let r: Word = radix.try_into().map_err(|_| ParseError::UnsupportedRadix)?;
138138
if r == B {
139139
#[allow(deprecated)] // TODO(v0.5): remove after from_str_native is made private.

integer/CHANGELOG.md

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,15 +3,29 @@
33
## Unreleased
44

55
### Add
6+
- NTT-based multiplication using Proth primes (`K·2^N + 1`), combined via Garner CRT. Supports 64-bit and 32-bit Word targets. Threshold at 4 000 words (~256 kbits).
7+
- Asymmetric NTT chunking: when one operand is much larger than the other, the shorter operand is forward-transformed once and reused across chunks.
68
- `UBig::from_u64` and `IBig::from_i64`, const on 32-bit and 64-bit targets.
79

810
### Improve
911
- Basecase (schoolbook) multiplication now uses an dword mult inner kernel (two multiplier words per sweep over the accumulator, mirroring GMP's `mpn_addmul_2` and `mpn_submul_2`), roughly halving accumulator memory traffic.
1012
- Addition and subtraction carry/borrow propagation now uses `Word` (u64/u32) instead of `bool` throughout the architecture-specific `add_with_carry` and `sub_with_borrow` functions, eliminating `bool`↔Word conversions in the inner loops.
13+
- Lowered the Karatsuba→Toom-3 multiplication threshold from 192 to 96 words, giving Toom-Cook-3 at ~6000 bits instead of ~12000 bits — closes the gap with malachite at ~10000-bit sizes.
14+
- NTT coefficient width increased from 16 to 64 bits (K_eff=3 for 64-bit, K_eff=2 otherwise), roughly halving the transform length at each step.
15+
- NTT multiplication auto-selects `K_eff = 2` primes when headroom allows, skipping the third prime.
16+
- Multiplication thresholds can be overridden at runtime via `DASHU_THRESHOLD_SIMPLE`, `DASHU_THRESHOLD_KARATSUBA`, and `DASHU_THRESHOLD_NTT` environment variables (requires `tuning` feature).
1117

12-
### Improve
13-
- Logarithm for very large values uses power-sequence decomposition, replacing iterative single-step multiplication.
14-
- Improve power-of-two base formatting ([#3](https://github.com/cmpute/dashu/pull/3))
18+
### Change
19+
- NTT multiplication now uses Proth primes (`K·2^N + 1`) instead of Solinas primes, improving modular reduction speed.
20+
- NTT threshold lowered from 40 000 to 4 000 words.
21+
- NTT enabled for 32-bit Word targets.
22+
- Arch-specific NTT prime definitions under `arch/generic_{32,64}_bit/ntt.rs`.
23+
24+
### Fix
25+
- `pack.rs` test used 64-bit literals that overflowed `Word` (`u32`) on 32-bit targets, breaking the test build.
26+
- `pack.rs` now uses native `Word`/`Lane` types throughout instead of `u64`/`u32`, fixing clippy `unnecessary_cast` warnings on 64-bit.
27+
- `test_unpack_carry_propagation` had a hardcoded 64-bit shift assumption; now derived from `Word::BITS` so it works on 32-bit.
28+
- Various clippy warnings (`let_and_return`, `too_many_arguments`, `needless_range_loop`, `type_complexity`) resolved across the NTT module.
1529

1630
## 0.4.2
1731

integer/Cargo.toml

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ all-features = true
1919
[features]
2020
default = ["std", "num-order"]
2121
std = ["dashu-base/std"]
22+
tuning = ["std"]
2223

2324
# unstable dependencies
2425
rand = ["rand_v08"]
@@ -30,7 +31,7 @@ dashu-base = { version = "0.4.1", default-features = false, path = "../base" }
3031
cfg-if = { version = "1.0.0" }
3132
static_assertions = { version = "1.1" }
3233
rustversion = { version = "1.0.0" }
33-
num-modular = { version = "0.6.1" }
34+
num-modular = { version = "0.6.4" }
3435

3536
# stable dependencies
3637
num-order = { optional = true, version = "1.2.0", default-features = false }

integer/benches/primitive.rs

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -150,6 +150,31 @@ fn ubig_ilog_large(criterion: &mut Criterion) {
150150
group.finish();
151151
}
152152

153+
fn ubig_mul_asymmetric(criterion: &mut Criterion) {
154+
let mut rng = StdRng::seed_from_u64(SEED);
155+
let mut group = criterion.benchmark_group("ubig_mul_asymmetric");
156+
group.plot_config(PlotConfiguration::default().summary_scale(AxisScale::Logarithmic));
157+
158+
// b just above the NTT threshold (4 000 words = 256 kbits → use 500 kbits).
159+
let b_bits = 500_000;
160+
let b = random_ubig(b_bits, &mut rng);
161+
162+
// a ranges from 1 kbit (below Karatsuba threshold) to heavily
163+
// asymmetric (10×), exercising all chunked-mul code paths.
164+
for &a_bits in &[
165+
1_000, 10_000, 100_000, 500_000, 1_000_000, 2_000_000, 5_000_000,
166+
] {
167+
let a = random_ubig(a_bits, &mut rng);
168+
group.bench_with_input(
169+
BenchmarkId::from_parameter(format!("{a_bits}/{b_bits}")),
170+
&(a, &b),
171+
|bencher, (ta, tb)| bencher.iter(|| ta * *tb),
172+
);
173+
}
174+
175+
group.finish();
176+
}
177+
153178
criterion_group!(
154179
benches,
155180
ubig_add,
@@ -163,6 +188,7 @@ criterion_group!(
163188
ubig_modulo_pow,
164189
ubig_pow_large_base,
165190
ubig_ilog_large,
191+
ubig_mul_asymmetric,
166192
);
167193

168194
criterion_main!(benches);

integer/src/arch/generic_32_bit/mod.rs

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,4 +4,5 @@ pub(crate) mod add;
44
#[path = "../generic/digits.rs"]
55
pub(crate) mod digits;
66

7+
pub(crate) mod ntt;
78
pub(crate) mod word;
Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
//! NTT primes and constants for 32-bit Word targets.
2+
//!
3+
//! Uses Proth primes of the form `K * 2^N + 1`.
4+
//! All constants computed by `integer/src/mul/ntt/compute_constants.py`.
5+
6+
use num_modular::FixedProth32;
7+
8+
// Proth reducer instances — each with a different (N, K) pair.
9+
pub type Rp0 = FixedProth32<26, 7>;
10+
pub type Rp1 = FixedProth32<27, 15>;
11+
pub type Rp2 = FixedProth32<27, 17>;
12+
13+
pub const P0: Rp0 = FixedProth32::<26, 7>;
14+
pub const P1: Rp1 = FixedProth32::<27, 15>;
15+
pub const P2: Rp2 = FixedProth32::<27, 17>;
16+
17+
pub const K: usize = 3;
18+
pub const MAX_LOG_N: u32 = 26;
19+
pub const B_PACK_MIN: u32 = 8;
20+
pub const B_PACK_CANDIDATES: &[u32] = &[32, 16, 8];
21+
22+
pub type Lane = u32;
23+
24+
/// Primitive `MAX_LOG_N`-th roots of unity for each prime.
25+
pub const OMEGA_MAX: [Lane; K] = [
26+
0x0000088b, // P0
27+
0x3a26eef8, // P1
28+
0x1aa0ab5e, // P2
29+
];
30+
31+
pub const CRT_INV_IJ: [[Lane; K]; K] = [[0, 0x4e42c85b, 0x5fb425ef], [0, 0, 0x44000009], [0, 0, 0]];
32+
33+
/// Prime moduli indexed by PI.
34+
pub const MODULI: [Lane; K] = [Rp0::MODULUS, Rp1::MODULUS, Rp2::MODULUS];
35+
36+
#[cfg(test)]
37+
mod tests {
38+
use super::*;
39+
use num_modular::Reducer;
40+
41+
type ReducerFns = (fn(Lane) -> Lane, fn(Lane) -> Lane, fn(Lane) -> Lane);
42+
43+
#[test]
44+
fn test_primes_proth_form() {
45+
assert_eq!(MODULI[0], 7u32 * (1u32 << 26) + 1);
46+
assert_eq!(MODULI[1], 15u32 * (1u32 << 27) + 1);
47+
assert_eq!(MODULI[2], 17u32 * (1u32 << 27) + 1);
48+
}
49+
50+
#[test]
51+
fn test_primes_v2() {
52+
for &p in &MODULI {
53+
let v2 = (p - 1).trailing_zeros();
54+
assert!(v2 >= MAX_LOG_N, "v2(p-1) = {v2} < MAX_LOG_N");
55+
}
56+
}
57+
58+
#[test]
59+
fn test_omega_order() {
60+
for (pi, &omega_max) in OMEGA_MAX.iter().enumerate() {
61+
let p = MODULI[pi];
62+
let (sqr, to_m, from_m): ReducerFns = match pi {
63+
0 => {
64+
(|w| P0.reduce((w as u64) * (w as u64)), |v| P0.transform(v), |v| P0.residue(v))
65+
}
66+
1 => {
67+
(|w| P1.reduce((w as u64) * (w as u64)), |v| P1.transform(v), |v| P1.residue(v))
68+
}
69+
2 => {
70+
(|w| P2.reduce((w as u64) * (w as u64)), |v| P2.transform(v), |v| P2.residue(v))
71+
}
72+
_ => unreachable!(),
73+
};
74+
75+
let mut w = to_m(omega_max);
76+
for _ in 0..MAX_LOG_N - 1 {
77+
w = sqr(w);
78+
}
79+
assert_eq!(from_m(w), p - 1, "omega^(2^(MAX_LOG_N-1)) != -1 mod p for prime {pi}");
80+
w = sqr(w);
81+
assert_eq!(from_m(w), 1, "omega^(2^MAX_LOG_N) != 1 mod p for prime {pi}");
82+
}
83+
}
84+
}

0 commit comments

Comments
 (0)