From fa9df23102703764d173d18dae97cce654f3aab7 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 22:35:46 +0400 Subject: [PATCH 01/13] refactor(math): add ceil_log_2, drop malachite CeilingLogBase2 in sync locator Branchless ilog2-based helper for the sync block-locator capacity hint; removes the consensus locator's only malachite-base import. --- consensus/src/processes/sync/mod.rs | 3 +- math/src/lib.rs | 84 +++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 2 deletions(-) diff --git a/consensus/src/processes/sync/mod.rs b/consensus/src/processes/sync/mod.rs index 0d6c73a98b..84de4a4844 100644 --- a/consensus/src/processes/sync/mod.rs +++ b/consensus/src/processes/sync/mod.rs @@ -4,7 +4,6 @@ use itertools::Itertools; use kaspa_consensus_core::errors::sync::{SyncManagerError, SyncManagerResult}; use kaspa_database::prelude::StoreResultExt; use kaspa_hashes::Hash; -use kaspa_math::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; use parking_lot::RwLock; use crate::model::{ @@ -141,7 +140,7 @@ impl< return Err(SyncManagerError::LowHashHigherThanHighHash(low, high)); } - let mut locator = Vec::with_capacity((high_index - low_index).ceiling_log_base_2() as usize); + let mut locator = Vec::with_capacity(kaspa_math::ceil_log_2(high_index - low_index) as usize); let mut step = 1; let mut current_index = high_index; while current_index > low_index { diff --git a/math/src/lib.rs b/math/src/lib.rs index e8c319eea0..a1a450c479 100644 --- a/math/src/lib.rs +++ b/math/src/lib.rs @@ -11,6 +11,16 @@ construct_uint!(Uint256, 4); construct_uint!(Uint320, 5); construct_uint!(Uint3072, 48); +/// Returns the ceiling of the base-2 logarithm of `x`, i.e. the smallest `k` such that `2^k >= x`. +/// +/// # Panics +/// Panics if `x` is 0 (the base-2 logarithm of 0 is undefined). +#[inline] +pub const fn ceil_log_2(x: u64) -> u64 { + // power of two -> floor; not a power of two -> floor + 1 + x.ilog2() as u64 + (!x.is_power_of_two()) as u64 +} + #[derive(thiserror::Error, Debug)] pub enum Error { #[error("{0:?}")] @@ -187,3 +197,77 @@ mod tests { assert_eq!(r / newr, expected); } } + +#[cfg(test)] +mod ceil_log_2_tests { + use crate::ceil_log_2; + + /// Independent reference: the smallest `k` such that `2^k >= x`. Computed in `u128` so it + /// stays correct for `x` near `u64::MAX` (where the answer is 64). + fn oracle(x: u64) -> u64 { + assert!(x != 0); + let mut k = 0u64; + while (1u128 << k) < x as u128 { + k += 1; + } + k + } + + /// A spread of inputs exercising the dense low range, every power-of-2 boundary, and the top. + fn sample_inputs() -> impl Iterator { + let dense = 1..=8192u64; + let boundaries = (0..64u32).flat_map(|k| { + let p = 1u64 << k; + // p-1 (clamped away from 0), p, p+1 + [p.saturating_sub(1).max(1), p, p + 1] + }); + dense.chain(boundaries).chain([u64::MAX]) + } + + #[test] + fn known_values() { + for (x, expected) in [ + (1u64, 0u64), + (2, 1), + (3, 2), + (4, 2), + (5, 3), + (7, 3), + (8, 3), + (9, 4), + (1023, 10), + (1024, 10), + (1025, 11), + (1u64 << 63, 63), + ((1u64 << 63) + 1, 64), + (u64::MAX, 64), + ] { + assert_eq!(ceil_log_2(x), expected, "ceil_log_2({x})"); + } + } + + #[test] + fn correctness_matches_oracle() { + for x in sample_inputs() { + assert_eq!(ceil_log_2(x), oracle(x), "x={x}"); + } + } + + #[test] + fn compatibility_matches_malachite() { + // Direct equivalence with the malachite `CeilingLogBase2` this replaced. Tied to the + // malachite dependency; remove alongside it once `mod_inverse` is migrated off malachite. + // Permanent correctness coverage is provided by `correctness_matches_oracle`. + use crate::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; + for x in sample_inputs() { + assert_eq!(ceil_log_2(x), x.ceiling_log_base_2(), "x={x}"); + } + } + + #[test] + #[should_panic] + fn panics_on_zero() { + // Matches malachite's "Cannot take the base-2 logarithm of 0." panic. + let _ = ceil_log_2(0); + } +} From e73a448a24855c2238cabd07fe001d969503c7f4 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 22:36:06 +0400 Subject: [PATCH 02/13] chore: bump workspace MSRV to 1.95.0 core::hint::cold_path, used by the modular-inverse routine added next, is stable as of 1.95. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 940798afb3..d000e1c5dd 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -70,7 +70,7 @@ members = [ ] [workspace.package] -rust-version = "1.91.0" +rust-version = "1.95.0" version = "2.0.0" authors = ["Kaspa developers"] license = "ISC" From eeff2ee9c03262c4e7eed01dcae772de98f141fb Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 22:37:32 +0400 Subject: [PATCH 03/13] feat(math): in-repo Lehmer 3072-bit modular inverse Stack-allocated variable-time modular inverse for the MuHash element via Lehmer's extended GCD (Knuth TAOCP 4.5.2 Algorithm L): a two-word half-GCD guess peels ~62 bits per 2x2 reduction matrix, with an exact Euclidean step as the fallback. 100% safe Rust with no new dependencies; cross-checked bit-for-bit against the existing malachite mod_inverse on 3072-bit vectors. The optional lehmer-instrument feature exposes per-inverse op counters. --- math/Cargo.toml | 6 + math/src/lehmer.rs | 1116 ++++++++++++++++++++++++++++++++++++++++++++ math/src/lib.rs | 1 + 3 files changed, 1123 insertions(+) create mode 100644 math/src/lehmer.rs diff --git a/math/Cargo.toml b/math/Cargo.toml index 3c40dbe756..c19ec7ea9e 100644 --- a/math/Cargo.toml +++ b/math/Cargo.toml @@ -28,6 +28,12 @@ workflow-wasm.workspace = true rand_chacha.workspace = true criterion.workspace = true +[features] +# Bench/diagnostics only: accumulate per-inverse operation counts in `lehmer` +# (matrices, divides, apply-limb iterations) into atomics. Off by default; has +# zero effect on the production path. +lehmer-instrument = [] + [[bench]] name = "bench" harness = false diff --git a/math/src/lehmer.rs b/math/src/lehmer.rs new file mode 100644 index 0000000000..13ee7e5c20 --- /dev/null +++ b/math/src/lehmer.rs @@ -0,0 +1,1116 @@ +//! Variable-time modular inversion via Lehmer's extended-GCD algorithm, +//! specialized for the 3072-bit MuHash element and fully stack-allocated. +//! +//! This is a *Euclidean* extended GCD (true quotients), not the binary +//! Bernstein-Yang "safegcd" divstep. The distinction matters: safegcd divides by +//! 2 each step, so its cofactors carry a `2^-62k` factor that must be paid down +//! with a full-width modular reduction every round, whereas a Euclidean GCD keeps +//! a single cofactor pair of clean integers that grow gradually (~`n/2` limbs on +//! average) and need only one sign fixup at the end. +//! +//! Outer loop ([`invert`]): each iteration extracts a 2x2 reduction matrix from +//! the aligned top *two* words of the two remainders (the half-GCD guess +//! [`hgcd2`]) and applies it to the full remainders ([`apply_matrix_xy`]) and to +//! the cofactor pair ([`apply_matrix_t`]), peeling ~62 bits per matrix (~50 +//! matrices for a 3072-bit inverse). When the guess cannot confirm even one +//! quotient it falls back to a single exact Euclidean division step. +//! +//! The guess ([`hgcd2`]) works on the top two `u64` words: a quotient of 1 is +//! resolved by one two-word subtraction and a top-word compare, so a hardware +//! divide is needed only for `q >= 2`. A trailing half-limb refinement keeps +//! every matrix entry below `2^63`, so the matrix-apply accumulators never +//! overflow and carry no per-step guard. +//! +//! Sign convention: track only `(t0, t1)`, the cofactor of `value`, kept +//! non-negative and growing; flip a `swapped` flag on each Euclidean swap and +//! recover the answer's sign from it at exit. +//! +//! Algorithm: Knuth, TAOCP Vol. 2, sec. 4.5.2, Algorithm L. Adapted from the +//! MIT-licensed 256-bit implementation by Dean Little (Copyright (c) 2026 Dean +//! Little, , `src/lehmer.rs`), +//! generalized to 48 limbs with a two-word half-GCD guess. Uses native +//! `u128`/`i128` accumulators plus this crate's `Uint3072` for the rare +//! multi-limb-divisor fallback; no external dependency. + +use crate::Uint3072; +use core::cmp::Ordering; + +/// Diagnostics-only per-inverse operation counters (enabled with the +/// `lehmer-instrument` feature). Used by `examples/lehmer_stats.rs` to measure +/// the work distribution. No effect on the production path. +#[cfg(feature = "lehmer-instrument")] +pub mod instrument { + use core::sync::atomic::AtomicU64; + pub static MATRICES: AtomicU64 = AtomicU64::new(0); + pub static FALLBACK_STEPS: AtomicU64 = AtomicU64::new(0); + pub static DIVIDES: AtomicU64 = AtomicU64::new(0); + pub static APPLY_XY_LIMBS: AtomicU64 = AtomicU64::new(0); + pub static APPLY_T_LIMBS: AtomicU64 = AtomicU64::new(0); + pub static SWAPS: AtomicU64 = AtomicU64::new(0); + pub static INVERSES: AtomicU64 = AtomicU64::new(0); + pub static GCD_DIV_NORM: AtomicU64 = AtomicU64::new(0); +} + +macro_rules! count { + ($name:ident += $v:expr) => {{ + #[cfg(feature = "lehmer-instrument")] + instrument::$name.fetch_add($v, core::sync::atomic::Ordering::Relaxed); + }}; +} + +/// Limb count of the operands (3072 bits in base 2^64). +const N: usize = 48; +/// Limb count for the cofactor buffers. The cofactor of `value` is bounded in +/// magnitude by the modulus (< 2^3072, i.e. <= `N` limbs); the extra headroom +/// absorbs carry-out during a matrix application and the multi-limb fallback. +const T: usize = N + 8; +// Matrix entries from `hgcd2` are below `2^63` by construction (the half-limb +// refinement discards the least significant half limb), so the `i128`/`u128` +// accumulators in the matrix applications cannot overflow +// (`a*x + b*y + carry < 2^128` with `a, b < 2^63` and `x, y < 2^64`). + +/// Compute the multiplicative inverse of `value` modulo `modulus` (an odd +/// modulus, e.g. the MuHash prime), via Lehmer's extended GCD. `value`, +/// `modulus` and `out` are little-endian base-2^64 limbs (`N = 48` words), and +/// `value` must already be reduced (`value < modulus`). +/// +/// Returns `true` and writes the inverse into `out` when +/// `gcd(value, modulus) == 1`; returns `false` (leaving `out` unspecified) +/// otherwise. A zero `value` yields `false`. +pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { + count!(INVERSES += 1); + if is_zero(&value) { + core::hint::cold_path(); + return false; + } + + // Euclidean state on (x, y) with x >= y, plus the cofactor of `value`. + // Invariant: x ≡ -t0·value (mod m) and y ≡ t1·value (mod m) when + // `swapped == false`, and the roles flip with each swap. `x`, `y` are owned + // working buffers (they are reduced in place). + let mut x: [u64; N] = modulus; + let mut y: [u64; N] = value; + let mut x_len = top_len(&x); + let mut y_len = top_len(&y); + let mut t0 = [0u64; T]; + let mut t1 = [0u64; T]; + t1[0] = 1; + let mut t0_len = 1usize; // 0 has conventional length 1 + let mut t1_len = 1usize; + let mut swapped = false; + + // Multi-limb phase: run until y collapses to a single limb, then finish + // with a u64 extended-Euclidean tail. + while y_len > 1 { + let (x_hi, y_hi) = highest_two_words_normalized(&x, &y, x_len); + let guess = hgcd2((x_hi >> 64) as u64, x_hi as u64, (y_hi >> 64) as u64, y_hi as u64); + + if let Some((a, b, c, d)) = guess { + count!(MATRICES += 1); + (x_len, y_len) = apply_matrix_xy(&mut x, &mut y, x_len, a, b, c, d); + apply_matrix_t(&mut t0, &mut t1, &mut t0_len, &mut t1_len, a, b, c, d); + + // A partial guess can leave x < y; restore x >= y, toggling the sign. + if x_len < y_len || (x_len == y_len && cmp_prefix(&x, &y, x_len).is_lt()) { + core::mem::swap(&mut x, &mut y); + core::mem::swap(&mut x_len, &mut y_len); + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + } else { + // The guess almost always confirms a quotient, so this exact-division + // path is cold; keep it out of the hot loop body. + core::hint::cold_path(); + count!(FALLBACK_STEPS += 1); + // Guess could not confirm a quotient: one exact Euclidean step. + let (q, r) = div_rem_step(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); + x = y; + x_len = y_len; + y = r; + y_len = top_len(&y); + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + } + + // Single-limb tail. `y` is now <= 1 limb; one u64 extended Euclidean step + // replaces the remaining ~tens of `t_add_qmul` iterations. + if y_len == 1 && y[0] != 0 { + if x[1..N].iter().any(|&w| w != 0) { + let (q, r) = div_rem_step(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); + let old_y = y[0]; + x = [0u64; N]; + x[0] = old_y; + y = r; + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + + let (g, cx, cy) = gcd_ext_u64(x[0], y[0]); + + // Fold the sign of the chosen cofactor into `swapped`, then combine + // magnitudes: new_t0 = |cx|·t0 + |cy|·t1. + if cx < 0 || (cx == 0 && cy > 0) { + swapped = !swapped; + } + let mut new_t0 = [0u64; T]; + let mut new_t0_len = 1usize; + add_mul_word(&mut new_t0, &mut new_t0_len, cx.unsigned_abs(), &t0, t0_len); + add_mul_word(&mut new_t0, &mut new_t0_len, cy.unsigned_abs(), &t1, t1_len); + t0 = new_t0; + + x = [0u64; N]; + x[0] = g; + x_len = if g == 0 { 0 } else { 1 }; + } + + // For an odd prime modulus and `value` in [1, m), the gcd is 1. + if !(x_len == 1 && x[0] == 1) { + return false; + } + + // x = 1 = ±t0·value (mod m). With `swapped`, the inverse is +t0; otherwise + // it is -t0 ≡ m - t0. Reduce t0 into [0, m) first. + let mut inv = reduce_mod(&t0, &modulus); + if !swapped && !is_zero(&inv) { + inv = sub_n(&modulus, &inv); + } + *out = inv; + true +} + +/// Two-word subtract `(h1·2^64 + l1) - (h2·2^64 + l2)`, returning `(hi, lo)`. +/// Callers guarantee the minuend is the larger, so no underflow off the top. +#[inline(always)] +fn sub2(h1: u64, l1: u64, h2: u64, l2: u64) -> (u64, u64) { + let (lo, borrow) = l1.overflowing_sub(l2); + (h1.wrapping_sub(h2).wrapping_sub(borrow as u64), lo) +} + +/// Fold the next half limb up into a fresh full word, `(hi << 32) | (lo >> 32)`, +/// for the handoff from the main hgcd2 loop (working on top *words*) to the +/// half-limb refinement loop. `hi < 2^32` here, so the result fits a `u64`. +#[inline(always)] +fn fold_half(hi: u64, lo: u64) -> u64 { + (hi << 32).wrapping_add(lo >> 32) +} + +/// The half-GCD guess: from the aligned top *two* words of the remainders +/// (`(ah, al) >=~ (bh, bl)`), run Euclidean steps on the 128-bit approximation +/// and return the 2x2 reduction matrix in this module's `(a, b, c, d)` +/// convention, `(x', y') = (a·x - b·y, d·y - c·x)`. Returns `None` when no +/// quotient can be confirmed (the caller then takes one exact division step). +/// +/// A quotient of 1 is resolved by one two-word subtraction and a top-word +/// compare, so [`gcd_div`] (the hardware divide) is reached only for `q >= 2`. +/// The trailing half-limb refinement (the `HALF_LIMIT_2` loop) keeps every +/// matrix entry below `2^63`, which is why the `apply_matrix_*` accumulators +/// need no per-step overflow guard. All matrix arithmetic is `wrapping_*` +/// (entries are `< 2^63`, so it never actually wraps) to avoid the workspace +/// `overflow-checks = true` panic pads. +#[inline] +fn hgcd2(mut ah: u64, mut al: u64, mut bh: u64, mut bl: u64) -> Option<(u64, u64, u64, u64)> { + // Need at least two leading bits of headroom in both top words. The guess + // almost always confirms a quotient, so every `return None` here is cold. + if ah < 2 || bh < 2 { + core::hint::cold_path(); + return None; + } + + // Reduction matrix M = (m00 m01; m10 m11), accumulated by the same column + // updates as the full-width apply; returned remapped to (a,b,c,d). + let (mut m00, mut m01, mut m10, mut m11): (u64, u64, u64, u64); + if ah > bh || (ah == bh && al > bl) { + (ah, al) = sub2(ah, al, bh, bl); // a -= b + if ah < 2 { + return None; + } + (m01, m10) = (1, 0); + } else { + (bh, bl) = sub2(bh, bl, ah, al); // b -= a + if bh < 2 { + return None; + } + (m01, m10) = (0, 1); + } + (m00, m11) = (1, 1); + + const HALF: u32 = 32; + const HALF_LIMIT_1: u64 = 1 << HALF; + + let mut subtract_a = ah < bh; + let mut subtract_a1 = false; + let mut done = false; + loop { + if subtract_a { + subtract_a = false; + } else { + if ah == bh { + done = true; + break; + } + if ah < HALF_LIMIT_1 { + // Top words collapsed to a half limb: fold in the next half from + // the low words and finish in the refinement loop below. + ah = fold_half(ah, al); + bh = fold_half(bh, bl); + break; + } + // Subtract a -= q·b (affects the second column of M). + (ah, al) = sub2(ah, al, bh, bl); + if ah < 2 { + done = true; + break; + } + if ah <= bh { + m01 = m01.wrapping_add(m00); // q = 1: no divide + m11 = m11.wrapping_add(m10); + } else { + let n = ((ah as u128) << 64) | al as u128; + let d = ((bh as u128) << 64) | bl as u128; + let (mut q, r) = gcd_div(n, d); + ah = (r >> 64) as u64; + al = r as u64; + if ah < 2 { + m01 = m01.wrapping_add(q.wrapping_mul(m00)); // q correct, a too small + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + done = true; + break; + } + q = q.wrapping_add(1); // one subtraction already taken above + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + } + } + if ah == bh { + done = true; + break; + } + if bh < HALF_LIMIT_1 { + ah = fold_half(ah, al); + bh = fold_half(bh, bl); + subtract_a1 = true; + break; + } + // Subtract b -= q·a (affects the first column of M). + (bh, bl) = sub2(bh, bl, ah, al); + if bh < 2 { + done = true; + break; + } + if bh <= ah { + m00 = m00.wrapping_add(m01); // q = 1: no divide + m10 = m10.wrapping_add(m11); + } else { + let n = ((bh as u128) << 64) | bl as u128; + let d = ((ah as u128) << 64) | al as u128; + let (mut q, r) = gcd_div(n, d); + bh = (r >> 64) as u64; + bl = r as u64; + if bh < 2 { + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + done = true; + break; + } + q = q.wrapping_add(1); + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + } + } + + // Half-limb refinement: peel a bit more (single-word divides on the folded + // top words) until |a - b| fits in one limb + 1 bit. This is what keeps the + // matrix entries below 2^63, discarding the least significant half limb. + if !done { + const HALF_LIMIT_2: u64 = 1 << (HALF + 1); + loop { + if subtract_a1 { + subtract_a1 = false; + } else { + ah = ah.wrapping_sub(bh); + if ah < HALF_LIMIT_2 { + break; + } + if ah <= bh { + m01 = m01.wrapping_add(m00); + m11 = m11.wrapping_add(m10); + } else { + let q = ah / bh; + ah %= bh; + if ah < HALF_LIMIT_2 { + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + break; + } + let q = q.wrapping_add(1); + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + } + } + bh = bh.wrapping_sub(ah); + if bh < HALF_LIMIT_2 { + break; + } + if ah >= bh { + m00 = m00.wrapping_add(m01); + m10 = m10.wrapping_add(m11); + } else { + let q = bh / ah; + bh %= ah; + if bh < HALF_LIMIT_2 { + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + break; + } + let q = q.wrapping_add(1); + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + } + } + } + + // Remap M to this module's apply convention: + // x' = m11·x - m01·y, y' = m00·y - m10·x. + Some((m11, m01, m10, m00)) +} + +/// `(q, r) = (n / d, n % d)` for 128-bit `n, d` with `d >= 2^64`, computed from +/// hardware `u64` divisions (via [`div_2by1`]) instead of the `u128 / u128` +/// software routine. The quotient is `< 2^64` since `n < 2^128 <= d * 2^64`. +#[inline(always)] +fn gcd_div(n: u128, d: u128) -> (u64, u128) { + count!(DIVIDES += 1); + let (n1, n0) = ((n >> 64) as u64, n as u64); + let (mut d1, mut d0) = ((d >> 64) as u64, d as u64); + debug_assert!(d1 != 0); + let (q, r) = (n1 / d1, n1 % d1); + if q > d1 { + // Non-normalized divisor: when `d1`'s top word is small, the single-word + // estimate `q = n1/d1` can exceed `d1`. Rare but genuinely reachable (a + // b-step can drop the divisor below 2^32 before an a-step divides by it), + // so this is required for correctness, not merely defensive; marked cold + // because it is rare. Normalize `d1` to have its top bit set, then do one + // 128/64 division. + core::hint::cold_path(); + count!(GCD_DIV_NORM += 1); + let c = d1.leading_zeros(); + let wc = 64 - c; + let n2 = n1 >> wc; + let n1n = (n1 << c) | (n0 >> wc); + let n0n = n0 << c; + d1 = (d1 << c) | (d0 >> wc); + d0 <<= c; + let (mut q, rem1) = div_2by1(n2, n1n, d1); + let prod = (q as u128) * (d0 as u128); + let (mut t1, mut t0) = ((prod >> 64) as u64, prod as u64); + if t1 > rem1 || (t1 == rem1 && t0 > n0n) { + q -= 1; + let (nt0, br) = t0.overflowing_sub(d0); + t0 = nt0; + t1 = t1.wrapping_sub(d1).wrapping_sub(br as u64); + } + let (rr0, br) = n0n.overflowing_sub(t0); + let rr1 = rem1.wrapping_sub(t1).wrapping_sub(br as u64); + let rr = ((rr1 as u128) << 64) | rr0 as u128; + (q, rr >> c) // undo normalization + } else { + let mut q = q; + let prod = (q as u128).wrapping_mul(d0 as u128); + let (mut t1, mut t0) = ((prod >> 64) as u64, prod as u64); + if t1 > r || (t1 == r && t0 > n0) { + q -= 1; + let (nt0, br) = t0.overflowing_sub(d0); + t0 = nt0; + t1 = t1.wrapping_sub(d1).wrapping_sub(br as u64); + } + let (rr0, br) = n0.overflowing_sub(t0); + let rr1 = r.wrapping_sub(t1).wrapping_sub(br as u64); + (q, ((rr1 as u128) << 64) | rr0 as u128) + } +} + +/// Multiply-accumulate: `acc + m·w + carry`, returning `(low 64 bits, high 64 +/// bits)`. The high word is the carry into the next limb. A single widening +/// `64x64->128` multiply plus an add-with-carry; the carry is one word. +#[inline(always)] +fn mac(acc: u64, m: u64, w: u64, carry: u64) -> (u64, u64) { + let t = (m as u128).wrapping_mul(w as u128).wrapping_add(acc as u128).wrapping_add(carry as u128); + (t as u64, (t >> 64) as u64) +} + +/// Multiply-subtract-borrow: `acc - m·w - borrow`, returning `(low 64 bits, +/// borrow)`. The borrow is the magnitude subtracted off the next limb. It fits +/// in one word because callers pass `m < 2^63` (the [`hgcd2`] entry bound), so +/// `m·w + borrow < 2^127 + 2^64` and its high word stays below `2^63`. +#[inline(always)] +fn msb(acc: u64, m: u64, w: u64, borrow: u64) -> (u64, u64) { + let sub = (m as u128).wrapping_mul(w as u128).wrapping_add(borrow as u128); + let (lo, b) = acc.overflowing_sub(sub as u64); + (lo, ((sub >> 64) as u64).wrapping_add(b as u64)) +} + +/// `(x, y) <- (a·x - b·y, d·y - c·x)`, over the `len` live limbs only (the +/// current larger length), returning the new `(x_len, y_len)`. Both results are +/// non-negative and fit in `len` limbs for matrices from [`hgcd2`], so +/// the limbs at `[len..N]` (already zero) stay zero. Tracking lengths here lets +/// passes shrink with the remainders instead of always touching all `N` limbs. +/// +/// Each output uses a multiply-accumulate / multiply-subtract structure: a `mac` +/// carry chain for the additive term and an independent `msb` borrow chain for +/// the subtractive one, both single-word. Keeping the two chains independent +/// shortens the loop-carried dependency versus a fused signed-`i128` accumulator, +/// which would serialize a two-register carry plus a per-limb sign-extension +/// every step. +#[inline] +fn apply_matrix_xy(x: &mut [u64; N], y: &mut [u64; N], len: usize, a: u64, b: u64, c: u64, d: u64) -> (usize, usize) { + // A single up-front bound: once the compiler knows `len <= N`, every `[i]` + // with `i < len` below is provably in-bounds, so the per-limb bounds checks + // collapse into this one branch and no `unsafe` is needed. `len <= N` always + // holds here, so this never actually panics. + assert!(len <= N, "apply_matrix_xy: len exceeds N"); + count!(APPLY_XY_LIMBS += len as u64); + // x' = a·x - b·y: `carry_ax` accumulates a·x, `borrow_x` peels off b·y. + let mut carry_ax: u64 = 0; + let mut borrow_x: u64 = 0; + // y' = d·y - c·x. + let mut carry_dy: u64 = 0; + let mut borrow_y: u64 = 0; + + for i in 0..len { + let xi = x[i]; + let yi = y[i]; + + let (px, ncx) = mac(0, a, xi, carry_ax); + let (nx, nbx) = msb(px, b, yi, borrow_x); + carry_ax = ncx; + borrow_x = nbx; + + let (py, ncy) = mac(0, d, yi, carry_dy); + let (ny, nby) = msb(py, c, xi, borrow_y); + carry_dy = ncy; + borrow_y = nby; + + x[i] = nx; + y[i] = ny; + } + // Non-negative result fitting in `len` limbs: the positive carry-out exactly + // cancels the negative borrow-out (both name the same high limb, which is 0). + debug_assert_eq!(carry_ax, borrow_x, "apply_matrix_xy: x went negative or overflowed len limbs"); + debug_assert_eq!(carry_dy, borrow_y, "apply_matrix_xy: y went negative or overflowed len limbs"); + + let nlx = x[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + let nly = y[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + (nlx, nly) +} + +/// `(t0, t1) <- (a·t0 + b·t1, c·t0 + d·t1)`, growing by at most one limb. +/// All-positive (the cofactor matrix uses `+` signs), so `u128` accumulators +/// suffice. Lengths are tracked inline; the cell at `max_len` is always written +/// (carry or zero) to keep `t[len..] == 0` without a separate clearing pass. +#[allow(clippy::too_many_arguments)] +#[inline] +fn apply_matrix_t(t0: &mut [u64; T], t1: &mut [u64; T], t0_len: &mut usize, t1_len: &mut usize, a: u64, b: u64, c: u64, d: u64) { + let max_len = (*t0_len).max(*t1_len); + count!(APPLY_T_LIMBS += max_len as u64); + // Single up-front bound; see `apply_matrix_xy`. With `max_len < T` known, the + // body `[i]` (i < max_len), the carry write `[max_len]`, and the length scan + // (`max_len + 1 <= T`) are all provably in-bounds, collapsing the per-limb + // checks into this one branch. The cofactor magnitude is bounded by the + // modulus with `T = N + 8` limbs of headroom, so `max_len < T` always holds. + assert!(max_len < T, "cofactor exceeded T limbs"); + let (a, b, c, d) = (a as u128, b as u128, c as u128, d as u128); + + let mut c0: u128 = 0; + let mut c1: u128 = 0; + + for i in 0..max_len { + let ti0 = t0[i] as u128; + let ti1 = t1[i] as u128; + + let p_at0 = a.wrapping_mul(ti0); + let p_bt1 = b.wrapping_mul(ti1); + let p_ct0 = c.wrapping_mul(ti0); + let p_dt1 = d.wrapping_mul(ti1); + + let n0 = c0.wrapping_add(p_at0).wrapping_add(p_bt1); + let n1 = c1.wrapping_add(p_ct0).wrapping_add(p_dt1); + + t0[i] = n0 as u64; + t1[i] = n1 as u64; + c0 = n0 >> 64; + c1 = n1 >> 64; + } + // Carry-out (provably < 2^64); also clears any stale cell at `max_len`. + t0[max_len] = c0 as u64; + t1[max_len] = c1 as u64; + + let nl0 = t0[..max_len + 1].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + let nl1 = t1[..max_len + 1].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + *t0_len = nl0.max(1); + *t1_len = nl1.max(1); +} + +/// `t0 += q · t1`, where `q` is up to `N` limbs (the exact Euclidean quotient). +fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T], t1_len: usize) { + for (i, &qi) in q.iter().enumerate() { + if qi == 0 { + continue; + } + let qi = qi as u128; + let mut carry: u64 = 0; + for (j, &t1j) in t1.iter().enumerate().take(t1_len) { + let k = i + j; + if k >= T { + debug_assert_eq!(carry, 0, "t_add_qmul: truncated nonzero limb"); + break; + } + let prod = qi.wrapping_mul(t1j as u128).wrapping_add(t0[k] as u128).wrapping_add(carry as u128); + t0[k] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = i + t1_len; + while carry > 0 && k < T { + let (s, c) = t0[k].overflowing_add(carry); + t0[k] = s; + carry = c as u64; + k += 1; + } + debug_assert!(carry == 0, "t_add_qmul: carry lost off the top"); + } + *t0_len = top_len_t(t0); +} + +/// `out += w · t` (multi-limb unsigned), used by the single-limb tail. +fn add_mul_word(out: &mut [u64; T], out_len: &mut usize, w: u64, t: &[u64; T], t_len: usize) { + if w != 0 { + let w = w as u128; + let mut carry: u64 = 0; + for j in 0..t_len { + // `j < t_len <= T` always, but the bound also lets the compiler clamp + // the trip count and skip the per-limb bounds check. + if j >= T { + break; + } + let prod = w.wrapping_mul(t[j] as u128).wrapping_add(out[j] as u128).wrapping_add(carry as u128); + out[j] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = t_len; + while carry > 0 && k < T { + let (s, c) = out[k].overflowing_add(carry); + out[k] = s; + carry = c as u64; + k += 1; + } + debug_assert!(carry == 0, "add_mul_word: carry lost off the top"); + } + *out_len = top_len_t(out); +} + +/// `(q, r) = (x / y, x % y)`. Fast path for a single-limb divisor (the common +/// case once Lehmer has shrunk `y`); otherwise [`Uint3072::div_rem`]. +fn div_rem_step(x: &[u64; N], y: &[u64; N]) -> ([u64; N], [u64; N]) { + if y[1..N].iter().all(|&w| w == 0) { + let (q, r0) = div_n_by_1(x, y[0]); + let mut r_arr = [0u64; N]; + r_arr[0] = r0; + return (q, r_arr); + } + let (q, r) = Uint3072(*x).div_rem(Uint3072(*y)); + (q.0, r.0) +} + +/// `(q, r) = (x / d, x % d)` for a nonzero single-limb divisor `d` (the common +/// case once Lehmer has collapsed `y`). The divisor is constant across all `N` +/// limbs, so a reciprocal is computed once and each limb costs two multiplies +/// instead of a hardware divide. Precomputed-reciprocal division +/// (Moller-Granlund, 2011). +fn div_n_by_1(x: &[u64; N], d: u64) -> ([u64; N], u64) { + debug_assert!(d != 0); + let mut q = [0u64; N]; + let s = d.leading_zeros(); + if s == 0 { + // `d` already has its top bit set: no normalization shift needed. + let v = invert_limb(d); + let mut r: u64 = 0; + for i in (0..N).rev() { + (q[i], r) = div_2by1_preinv(r, x[i], d, v); + } + (q, r) + } else { + // Normalize: divide `x << s` by `d << s` (same quotient); the bits shifted + // off the top of `x` seed the running remainder, and `x % d = r >> s`. + let d_norm = d << s; + let v = invert_limb(d_norm); + let mut r = x[N - 1] >> (64 - s); + for i in (0..N).rev() { + let lo = if i == 0 { 0 } else { x[i - 1] }; + let cur = (x[i] << s) | (lo >> (64 - s)); + (q[i], r) = div_2by1_preinv(r, cur, d_norm, v); + } + (q, r >> s) + } +} + +/// Reciprocal of a normalized 64-bit divisor `d` (top bit set): +/// `v = floor((2^128 - 1) / d) - 2^64`, the "3/2" reciprocal consumed by +/// [`div_2by1_preinv`] (Moller-Granlund, "Improved division by invariant +/// integers", 2011). The lone wide divide here runs once per [`div_n_by_1`] +/// call (amortized over all `N` limbs), so it is off the per-limb hot path. +#[inline] +const fn invert_limb(d: u64) -> u64 { + debug_assert!(d >> 63 == 1, "invert_limb requires a normalized divisor"); + // v = floor((2^128 - 1) / d) - 2^64 = floor(((!d)·2^64 + !0) / d): subtracting + // 2^64·d from the numerator drops the quotient by exactly 2^64, and that + // numerator's high word `!d` is `< d` for a normalized `d`, so the quotient + // fits a u64 and the software 128/64 [`div_2by1`] computes it -- avoiding the + // slow `u128 / u128` (`__udivti3`). + div_2by1(!d, u64::MAX, d).0 +} + +/// `(nh·2^64 + nl) / d -> (q, r)` via the precomputed reciprocal `di` from +/// [`invert_limb`]. Requires `d` normalized (top bit set) and `nh < d` (so the +/// quotient fits a `u64`); the returned `r` satisfies `r < d`, sustaining that +/// precondition limb to limb. One widening multiply, a branchless first +/// correction, and a rarely-taken second. All `wrapping_*` (the workspace builds +/// with `overflow-checks = true`). +#[inline(always)] +fn div_2by1_preinv(nh: u64, nl: u64, d: u64, di: u64) -> (u64, u64) { + debug_assert!(d >> 63 == 1 && nh < d); + // (qh:ql) = nh·di + (nh + 1)·2^64 + nl + let prod = (nh as u128).wrapping_mul(di as u128); + let (mut qh, ql) = ((prod >> 64) as u64, prod as u64); + let (ql, carry) = ql.overflowing_add(nl); + qh = qh.wrapping_add(nh.wrapping_add(1)).wrapping_add(carry as u64); + + let mut r = nl.wrapping_sub(qh.wrapping_mul(d)); + // First correction: if the estimate `qh` is one too high, `r` wraps above + // `ql`; `mask` is all-ones then, folding the `-1`/`+d` fixup in branch-free. + let mask = 0u64.wrapping_sub((r > ql) as u64); + qh = qh.wrapping_add(mask); + r = r.wrapping_add(mask & d); + // Second correction (rare): a still-too-small quotient leaves `r >= d`. + if r >= d { + r = r.wrapping_sub(d); + qh = qh.wrapping_add(1); + } + (qh, r) +} + +/// `(hi·2^64 + lo) / d -> (q, r)`. Knuth Algorithm D specialized to 128/64. +/// Precondition: `d != 0` and `hi < d` (so the quotient fits in a `u64`). +/// +/// Not on the per-limb hot path (that goes through [`div_n_by_1`]'s reciprocal +/// step, [`div_2by1_preinv`]): this generic divide is reached only once per +/// [`div_n_by_1`] (to seed the reciprocal, via [`invert_limb`]) and from +/// [`gcd_div`]'s rare normalization branch, so a plain software long division is +/// fine. +#[inline] +const fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) { + debug_assert!(d != 0); + debug_assert!(hi < d, "div_2by1 quotient would overflow u64"); + + let s = d.leading_zeros(); + let d = if s == 0 { d } else { d << s }; + let un32 = if s == 0 { hi } else { (hi << s) | (lo >> (64 - s)) }; + let un10 = if s == 0 { lo } else { lo << s }; + + let vn1 = d >> 32; + let vn0 = d & 0xFFFF_FFFF; + let un1 = un10 >> 32; + let un0 = un10 & 0xFFFF_FFFF; + + let mut q1 = un32 / vn1; + let mut rhat = un32 - q1 * vn1; + while q1 >= (1u64 << 32) || q1 * vn0 > (rhat << 32) | un1 { + q1 -= 1; + rhat += vn1; + if rhat >= (1u64 << 32) { + break; + } + } + + let un21 = (un32 << 32).wrapping_add(un1).wrapping_sub(q1.wrapping_mul(d)); + let mut q0 = un21 / vn1; + let mut rhat = un21 - q0 * vn1; + while q0 >= (1u64 << 32) || q0 * vn0 > (rhat << 32) | un0 { + q0 -= 1; + rhat += vn1; + if rhat >= (1u64 << 32) { + break; + } + } + + let r = (un21 << 32).wrapping_add(un0).wrapping_sub(q0.wrapping_mul(d)) >> s; + ((q1 << 32) | q0, r) +} + +/// Extended Euclidean GCD of two `u64`s: `(gcd, cx, cy)` with `cx·x + cy·y = g`. +/// Called once per inversion, so `i128` intermediates are fine. +fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { + let (mut a, mut b, mut c, mut d) = (1i128, 0i128, 0i128, 1i128); + while y != 0 { + let q = (x / y) as i128; + let r = x.wrapping_sub((q as u64).wrapping_mul(y)); + let (nc, nd) = (a.wrapping_sub(q.wrapping_mul(c)), b.wrapping_sub(q.wrapping_mul(d))); + a = c; + b = d; + c = nc; + d = nd; + x = y; + y = r; + } + (x, a as i64, b as i64) +} + +/// Reduce a cofactor (`<= N` significant limbs) into `[0, m)`. The cofactor is +/// bounded by `m` at loop exit, so a few subtractions suffice; the `div_rem` +/// branch is a defensive fallback. +fn reduce_mod(value: &[u64; T], m: &[u64; N]) -> [u64; N] { + debug_assert!(value[N..T].iter().all(|&w| w == 0), "cofactor exceeded N limbs"); + let mut out = [0u64; N]; + out.copy_from_slice(&value[..N]); + for _ in 0..8 { + if cmp_n(&out, m).is_lt() { + return out; + } + out = sub_n(&out, m); + } + (Uint3072(out) % Uint3072(*m)).0 +} + +/// Top 128-bit prefix of `x` and of `y`, both shifted left by the same amount so +/// the prefix of the larger operand `x` has its top bit set (the normalization +/// [`hgcd2`] expects). `y` is read at the same limb positions as `x`, so its +/// prefix is naturally the smaller. Requires `x` the larger operand and +/// `x_len >= 2`. +#[inline] +fn highest_two_words_normalized(x: &[u64; N], y: &[u64; N], x_len: usize) -> (u128, u128) { + debug_assert!(x_len >= 2); + let i = x_len - 1; + let lz = x[i].leading_zeros(); + let lo_idx_ok = i >= 2; + // Combine the limbs at positions (i, i-1, i-2) into the top 128 bits after a + // left shift by `lz`. `arr[i]`'s `lz` leading zeros guarantee no overflow. + let combine = |hi: u64, mid: u64, lo: u64| -> u128 { + let base = ((hi as u128) << 64) | (mid as u128); + if lz == 0 { base } else { (base << lz) | ((lo >> (64 - lz)) as u128) } + }; + let x_lo = if lo_idx_ok { x[i - 2] } else { 0 }; + let y_lo = if lo_idx_ok { y[i - 2] } else { 0 }; + (combine(x[i], x[i - 1], x_lo), combine(y[i], y[i - 1], y_lo)) +} + +#[inline] +fn is_zero(x: &[u64; N]) -> bool { + x.iter().all(|&w| w == 0) +} + +#[inline] +fn top_len(x: &[u64; N]) -> usize { + for i in (0..N).rev() { + if x[i] != 0 { + return i + 1; + } + } + 0 +} + +/// Significant length of a cofactor buffer; 0 maps to 1 by convention. +#[inline] +fn top_len_t(x: &[u64; T]) -> usize { + for i in (0..T).rev() { + if x[i] != 0 { + return i + 1; + } + } + 1 +} + +#[inline] +fn sub_n(a: &[u64; N], b: &[u64; N]) -> [u64; N] { + let mut out = [0u64; N]; + let mut borrow = 0u64; + for i in 0..N { + let (d1, b1) = a[i].overflowing_sub(b[i]); + let (d2, b2) = d1.overflowing_sub(borrow); + out[i] = d2; + borrow = (b1 | b2) as u64; + } + out +} + +#[inline] +fn cmp_n(a: &[u64; N], b: &[u64; N]) -> Ordering { + for i in (0..N).rev() { + match a[i].cmp(&b[i]) { + Ordering::Equal => continue, + other => return other, + } + } + Ordering::Equal +} + +#[inline] +fn cmp_prefix(a: &[u64; N], b: &[u64; N], n: usize) -> Ordering { + for i in (0..n).rev() { + match a[i].cmp(&b[i]) { + Ordering::Equal => continue, + other => return other, + } + } + Ordering::Equal +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Uint3072; + use rand_chacha::{ + ChaCha8Rng, + rand_core::{RngCore, SeedableRng}, + }; + + // MuHash prime: 2^3072 - 1103717. + const MUHASH_PRIME: Uint3072 = { + let mut max = Uint3072::MAX; + max.0[0] -= 1103717 - 1; + max + }; + + fn lehmer_inv(value: [u64; N], modulus: [u64; N]) -> Option<[u64; N]> { + let mut out = [0u64; N]; + invert(value, modulus, &mut out).then_some(out) + } + + fn malachite_inv(value: [u64; N], modulus: [u64; N]) -> Option<[u64; N]> { + Uint3072(value).mod_inverse(Uint3072(modulus)).map(|inv| inv.0) + } + + fn addmod(a: Uint3072, b: Uint3072, m: Uint3072) -> Uint3072 { + let (res, overflow) = a.overflowing_add(b); + if overflow || res >= m { res.overflowing_sub(m).0 } else { res } + } + + /// `(a * b) mod m` via binary double-and-add (overflow-safe oracle). + fn mulmod(a: Uint3072, b: Uint3072, m: Uint3072) -> Uint3072 { + let mut result = Uint3072::ZERO; + let mut base = a % m; + let mut exp = b; + while !exp.is_zero() { + if exp.0[0] & 1 == 1 { + result = addmod(result, base, m); + } + base = addmod(base, base, m); + exp = exp >> 1; + } + result + } + + #[test] + fn matches_malachite_muhash_prime() { + // Lehmer must agree bit-for-bit with malachite (the reference) on the MuHash prime. + let mut rng = ChaCha8Rng::seed_from_u64(42); + let mut buf = [0u8; Uint3072::BYTES]; + for _ in 0..2000 { + rng.fill_bytes(&mut buf); + let v = Uint3072::from_le_bytes(buf) % MUHASH_PRIME; + if v.is_zero() { + continue; + } + let expected = malachite_inv(v.0, MUHASH_PRIME.0).unwrap(); + let got = lehmer_inv(v.0, MUHASH_PRIME.0).unwrap(); + assert_eq!(got.as_slice(), expected.as_slice(), "v={v}"); + } + } + + #[test] + fn product_is_one_muhash_prime() { + // v * inv(v) == 1 (mod prime), independent of any oracle. + let mut rng = ChaCha8Rng::seed_from_u64(99); + let mut buf = [0u8; Uint3072::BYTES]; + for _ in 0..500 { + rng.fill_bytes(&mut buf); + let v = Uint3072::from_le_bytes(buf) % MUHASH_PRIME; + if v.is_zero() { + continue; + } + let inv = lehmer_inv(v.0, MUHASH_PRIME.0).unwrap(); + assert_eq!(mulmod(v, Uint3072(inv), MUHASH_PRIME), Uint3072::from_u64(1), "v={v}"); + } + } + + #[test] + fn small_values_force_multilimb_quotient() { + // Tiny `value` makes the first quotient ~48 limbs, exercising the + // multi-limb `div_rem` fallback and the wide `t_add_qmul` path. + for v_small in [1u64, 2, 3, 5, 7, 1103717, u64::MAX] { + let v = Uint3072::from_u64(v_small); + let inv = lehmer_inv(v.0, MUHASH_PRIME.0).unwrap(); + assert_eq!(mulmod(v, Uint3072(inv), MUHASH_PRIME), Uint3072::from_u64(1), "v={v_small}"); + } + } + + #[test] + fn edge_cases_muhash_prime() { + let one = Uint3072::from_u64(1); + let p_minus_1 = MUHASH_PRIME.overflowing_sub(one).0; + assert_eq!(lehmer_inv(one.0, MUHASH_PRIME.0).unwrap().as_slice(), one.0.as_slice()); + assert_eq!(lehmer_inv(p_minus_1.0, MUHASH_PRIME.0).unwrap().as_slice(), p_minus_1.0.as_slice()); + assert!(lehmer_inv(Uint3072::ZERO.0, MUHASH_PRIME.0).is_none()); + } + + #[test] + fn div_n_by_1_matches_uint_div() { + // The reciprocal n-by-1 division must agree with Uint3072::div_rem across + // tiny, normalized (top bit set), and arbitrary unnormalized divisors. + let mut rng = ChaCha8Rng::seed_from_u64(7); + let mut buf = [0u8; Uint3072::BYTES]; + for _ in 0..20000 { + rng.fill_bytes(&mut buf); + let x = Uint3072::from_le_bytes(buf); + let d = match rng.next_u64() % 4 { + 0 => 1 + (rng.next_u64() % 1000), // tiny (large shift) + 1 => rng.next_u64() | (1 << 63), // already normalized + 2 => (rng.next_u64() >> 7).max(1), // unnormalized + _ => rng.next_u64().max(1), // arbitrary + }; + let (q, r) = div_n_by_1(&x.0, d); + let mut dl = [0u64; N]; + dl[0] = d; + let (eq, er) = Uint3072(x.0).div_rem(Uint3072(dl)); + assert_eq!(q.as_slice(), eq.0.as_slice(), "quotient mismatch x={x} d={d}"); + assert_eq!(r, er.0[0], "remainder mismatch x={x} d={d}"); + } + // Boundary divisors: 1, 2, MAX, 2^63 (min normalized), 2^63-1 (max + // unnormalized), and the MuHash prime difference constant. + for &d in &[1u64, 2, 3, u64::MAX, 1 << 63, (1 << 63) - 1, 1103717] { + rng.fill_bytes(&mut buf); + let x = Uint3072::from_le_bytes(buf); + let (q, r) = div_n_by_1(&x.0, d); + let mut dl = [0u64; N]; + dl[0] = d; + let (eq, er) = Uint3072(x.0).div_rem(Uint3072(dl)); + assert_eq!(q.as_slice(), eq.0.as_slice(), "quotient mismatch x={x} d={d}"); + assert_eq!(r, er.0[0], "remainder mismatch x={x} d={d}"); + } + } + + #[test] + fn gcd_div_matches_u128_both_branches() { + // `gcd_div` computes `(n/d, n%d)` for 128-bit `n, d` with `d >= 2^64`. + // Cover both arms against the native `u128` oracle: the normal arm + // (`q = n1/d1 <= d1`) and the rare normalization arm (`q > d1`), which + // fires when the divisor's top word is small. The normalization arm is + // hit only occasionally through `invert`, so pin it directly here. + let oracle = |n: u128, d: u128| ((n / d) as u64, n % d); + + // Explicit cases that force `q > d1` (small divisor top word). + let forcing: &[(u128, u128)] = &[ + (u128::MAX, (1u128 << 64) | 1), // d1 = 1 + (u128::MAX, 1u128 << 64), // d = 2^64 + ((0xDEAD_BEEFu128 << 96) | 0x1234, (1u128 << 64) | u64::MAX as u128), + (u128::MAX, (2u128 << 64) | 7), // d1 = 2 + (u128::MAX - 12345, (0x1_0000u128 << 64) | 0xABCD), // d1 = 2^16 + ]; + for &(n, d) in forcing { + assert!((n >> 64) as u64 / (d >> 64) as u64 > (d >> 64) as u64, "case does not force q>d1"); + assert_eq!(gcd_div(n, d), oracle(n, d), "q>d1 arm mismatch n={n} d={d}"); + } + // Normal arm (`q <= d1`): divisor with a large top word. + for &(n, d) in &[(u128::MAX, (u64::MAX as u128) << 64 | 3), (1u128 << 127, (1u128 << 127) | 1)] { + assert!((n >> 64) as u64 / (d >> 64) as u64 <= (d >> 64) as u64, "case is not the normal arm"); + assert_eq!(gcd_div(n, d), oracle(n, d), "normal arm mismatch n={n} d={d}"); + } + + // Randomized, biased to small divisor top words so most iterations take + // the `q > d1` arm; every result checked bit-for-bit against `u128`. + let mut rng = ChaCha8Rng::seed_from_u64(2024); + let mut branch_hits = 0u64; + for _ in 0..200_000 { + let n = ((rng.next_u64() as u128) << 64) | rng.next_u64() as u128; + let d1 = (rng.next_u64() % (1 << 20)) + 1; // 1 ..= 2^20 + let d = ((d1 as u128) << 64) | rng.next_u64() as u128; + assert_eq!(gcd_div(n, d), oracle(n, d), "random gcd_div mismatch n={n} d={d}"); + if (n >> 64) as u64 / d1 > d1 { + branch_hits += 1; + } + } + assert!(branch_hits > 1000, "q>d1 arm under-exercised: only {branch_hits} hits"); + } + + #[test] + fn q_gt_d1_vector_matches_malachite() { + // A full-inverse input that drives `gcd_div` through its `q > d1` + // normalization arm, pinned as an explicit regression case (that arm is + // otherwise only hit occasionally by the random-vector tests). + let v: [u64; N] = [ + 0xbaac78a3d8d04a44, + 0x454126b8efd12383, + 0xa93bb055d701be60, + 0xe940f5627944ba89, + 0x68842c54edf3df88, + 0xf7c6332a4e2be869, + 0x396e8533e978070a, + 0x2703f08794977aad, + 0x373dc910525ad335, + 0x52520bee468a9073, + 0xf4580a27f3ca91ee, + 0x5de1060b3d65f732, + 0x4b1072d2e2cc0da1, + 0x4e2a032ba51f609a, + 0xeae5995402410005, + 0x885f9eb04a59ab3c, + 0x165446e913daa18e, + 0xf2b1311c3eb843c0, + 0x630b8232162a2fca, + 0xe4224b2c61bafd4e, + 0x25397d500e52b519, + 0x7cd55e35a4ac6022, + 0x8629892065e194c0, + 0x713ea3ded7bcd68c, + 0xf1f138b0052d1fdc, + 0x7974414d82958c31, + 0xc29243783567cf96, + 0xf9214af62a72ec1e, + 0xa1e4feb97fdbfe07, + 0x5c92a43f23a3989e, + 0x1ad829d92571a26a, + 0xd09dd31bf57bf618, + 0x00772f32162d2fd3, + 0x1c0dcedad1142715, + 0x677b857b2c9c2713, + 0xe64b41b0e187f9e7, + 0x4062f89c0309cd59, + 0x3edc673c30e10664, + 0x330b6c680b777302, + 0x423676df7dfd8ccf, + 0x347fe8b7dee8ca42, + 0x0fb4eb8629ef71c1, + 0x5e898abada4103ac, + 0x39820234bad59fb2, + 0xd780c6aeaaa89812, + 0x61bdb22ece416ba0, + 0x79f1821384ef3f44, + 0x803912aebf95eede, + ]; + let got = lehmer_inv(v, MUHASH_PRIME.0).unwrap(); + assert_eq!(got.as_slice(), malachite_inv(v, MUHASH_PRIME.0).unwrap().as_slice()); + assert_eq!(mulmod(Uint3072(v), Uint3072(got), MUHASH_PRIME), Uint3072::from_u64(1)); + } + + #[test] + fn hgcd2_none_when_top_word_below_two() { + // The guess returns `None` (caller falls back to one exact division step) + // exactly when a top word is `< 2`; otherwise it confirms a matrix. + assert!(hgcd2(1, 0, 5, 0).is_none()); // ah < 2 + assert!(hgcd2(5, 0, 0, 9).is_none()); // bh < 2 + assert!(hgcd2(u64::MAX, 0, u64::MAX >> 1, 0).is_some()); // well-separated -> confirms a quotient + } +} diff --git a/math/src/lib.rs b/math/src/lib.rs index a1a450c479..efcaa99503 100644 --- a/math/src/lib.rs +++ b/math/src/lib.rs @@ -3,6 +3,7 @@ use wasm_bindgen::JsValue; use workflow_core::sendable::Sendable; pub mod int; +pub mod lehmer; pub mod uint; pub mod wasm; From cd57641e2583c9d8f1749bca0edfa61e1ea813b2 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 22:39:28 +0400 Subject: [PATCH 04/13] bench(math): 3072-bit modular-inverse candidate benchmarks Criterion comparison of the in-repo Lehmer inverse against malachite (the current implementation), dashu, and crypto-bigint's safegcd, alongside the other production big-uint ops (blue-work sum, difficulty window, calc_work). Every candidate is asserted equal to the current malachite result before it is timed. dashu-int and crypto-bigint are bench-only dev-dependencies. --- Cargo.lock | 76 +++++++ math/Cargo.toml | 12 ++ math/benches/candidates.rs | 427 +++++++++++++++++++++++++++++++++++++ 3 files changed, 515 insertions(+) create mode 100644 math/benches/candidates.rs diff --git a/Cargo.lock b/Cargo.lock index 4b2a316269..e8deb78770 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1102,6 +1102,12 @@ version = "0.7.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a1d728cc89cf3aee9ff92b05e62b19ee65a02b5702cff7d5a377e32c6ae29d8d" +[[package]] +name = "cmov" +version = "0.5.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0c9ea0ac24bc397ab3c98583a3c9ba74fa56b09a4449bbe172b9b1ddb016027a" + [[package]] name = "cobs" version = "0.3.0" @@ -1205,6 +1211,12 @@ dependencies = [ "libc", ] +[[package]] +name = "cpubits" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "15b85f9c39137c3a891689859392b1bd49812121d0d61c9caf00d46ed5ce06ae" + [[package]] name = "cpufeatures" version = "0.2.17" @@ -1331,6 +1343,18 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7a81dae078cea95a014a339291cec439d2f232ebe854a9d672b796c6afafa9b7" +[[package]] +name = "crypto-bigint" +version = "0.7.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42a0d26b245348befa0c121944541476763dcc46ede886c88f9d12e1697d27c3" +dependencies = [ + "cpubits", + "ctutils", + "num-traits", + "rand_core 0.10.1", +] + [[package]] name = "crypto-common" version = "0.1.6" @@ -1383,6 +1407,15 @@ dependencies = [ "windows-sys 0.59.0", ] +[[package]] +name = "ctutils" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7d5515a3834141de9eafb9717ad39eea8247b5674e6066c404e8c4b365d2a29e" +dependencies = [ + "cmov", +] + [[package]] name = "curve25519-dalek" version = "4.1.3" @@ -1492,6 +1525,26 @@ dependencies = [ "parking_lot_core", ] +[[package]] +name = "dashu-base" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fab3f0756c8585395280bd81b384cd28bbd66d6b00e66124ecfd1f644938b38c" + +[[package]] +name = "dashu-int" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d6a93d93dc2aca9a071e6ecb2af883441e8b34357a1037bc536e1c52b0d3c713" +dependencies = [ + "cfg-if 1.0.0", + "dashu-base", + "num-modular", + "num-order", + "rustversion", + "static_assertions", +] + [[package]] name = "data-encoding" version = "2.6.0" @@ -3346,6 +3399,8 @@ version = "2.0.0" dependencies = [ "borsh", "criterion", + "crypto-bigint", + "dashu-int", "faster-hex 0.9.0", "js-sys", "kaspa-utils", @@ -5015,6 +5070,21 @@ dependencies = [ "num-traits", ] +[[package]] +name = "num-modular" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc41a1374056e9672221567958a66c16be12d0e2c1b408761e14d901c237d5e0" + +[[package]] +name = "num-order" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "537b596b97c40fcf8056d153049eb22f481c17ebce72a513ec9286e4986d1bb6" +dependencies = [ + "num-modular", +] + [[package]] name = "num-rational" version = "0.4.2" @@ -6634,6 +6704,12 @@ version = "1.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" +[[package]] +name = "static_assertions" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2eb9349b6444b326872e140eb1cf5e7c522154d69e7a0ffb0fb81c06b37543f" + [[package]] name = "strsim" version = "0.11.1" diff --git a/math/Cargo.toml b/math/Cargo.toml index c19ec7ea9e..59136cff3b 100644 --- a/math/Cargo.toml +++ b/math/Cargo.toml @@ -27,6 +27,14 @@ workflow-wasm.workspace = true [dev-dependencies] rand_chacha.workspace = true criterion.workspace = true +# candidate for replacing the malachite-* deps (bench-only, see benches/candidates.rs). +# the mul tuning from https://github.com/cmpute/dashu/pull/74 was also benched (via a local +# clone, the branch has a broken gitlink) and showed no effect at 48 limbs: 22.8 vs 22.4 us/inv +dashu-int = "0.4" +crypto-bigint = "0.7" +# ruled out by earlier runs (4.5x-31x slower than dashu at 3072-bit modinv): +# num-bigint = "0.4" +# ruint = "1" [features] # Bench/diagnostics only: accumulate per-inverse operation counts in `lehmer` @@ -38,5 +46,9 @@ lehmer-instrument = [] name = "bench" harness = false +[[bench]] +name = "candidates" +harness = false + [lints] workspace = true diff --git a/math/benches/candidates.rs b/math/benches/candidates.rs new file mode 100644 index 0000000000..f55927a7ac --- /dev/null +++ b/math/benches/candidates.rs @@ -0,0 +1,427 @@ +//! Evaluates replacing the malachite-* deps, based on all production big-uint usage: +//! +//! 1. `mod_inverse` (the only malachite-backed op, called from muhash u3072.rs): inverse of a +//! 3072-bit value modulo the MuHash prime 2^3072 - 1103717, with `[u64; 48]` limbs in/out. +//! Limb conversions are inside the measured loop, modulus conversion is precomputed (it is +//! a const at the call site). +//! 2. `ceiling_log_base_2` on u64 (consensus sync locator) vs plain std. +//! 3. The hand-rolled (repo-owned, license-clean) Uint ops actually used in production, each +//! against the dashu equivalent, to decide whether wholesale migration makes sense or the +//! macro should stay for everything but `mod_inverse`: +//! - blue work accumulation, Uint192 add/sum (ghostdag protocol.rs) +//! - difficulty window average, Uint320 sum / div u64 / mul u64 / min (difficulty.rs) +//! - calc_work, Uint256 `(!t / (t + 1)) + 1` (difficulty.rs) +//! - compact target bits codec round-trip (no candidate: repo-owned either way) + +use criterion::measurement::WallTime; +use criterion::{BenchmarkGroup, Criterion, black_box, criterion_group, criterion_main}; +use dashu_int::UBig; +use rand_chacha::{ + ChaCha8Rng, + rand_core::{RngCore, SeedableRng}, +}; + +use kaspa_math::{Uint192, Uint256, Uint320, Uint3072}; + +const N: usize = 32; +const LIMBS: usize = 48; + +// Same value as muhash UINT_PRIME: 2^3072 - 1103717 +const PRIME: Uint3072 = { + let mut max = Uint3072::MAX; + max.0[0] -= 1103717 - 1; + max +}; + +fn limbs_to_le_bytes(limbs: &[u64; LIMBS]) -> [u8; LIMBS * 8] { + let mut bytes = [0u8; LIMBS * 8]; + for (chunk, limb) in bytes.chunks_exact_mut(8).zip(limbs) { + chunk.copy_from_slice(&limb.to_le_bytes()); + } + bytes +} + +fn le_bytes_to_limbs(bytes: &[u8]) -> [u64; LIMBS] { + let mut limbs = [0u64; LIMBS]; + for (limb, chunk) in limbs.iter_mut().zip(bytes.chunks(8)) { + let mut buf = [0u8; 8]; + buf[..chunk.len()].copy_from_slice(chunk); + *limb = u64::from_le_bytes(buf); + } + limbs +} + +fn inputs() -> Vec<[u64; LIMBS]> { + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let mut buf = [0u8; LIMBS * 8]; + (0..N) + .map(|_| { + rng.fill_bytes(&mut buf); + (Uint3072::from_le_bytes(buf) % PRIME).0 + }) + .collect() +} + +mod current_malachite { + use super::*; + + pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { + Uint3072(*limbs).mod_inverse(PRIME).expect("0 < a < prime").0 + } +} + +// Ruled out by earlier runs (per-inverse: num-bigint ~513us, ruint ~72us vs dashu ~22us and +// malachite ~16us); kept for reference. (crypto-bigint is no longer commented out: it is wired +// into both modinv_3072 groups below to settle the "safegcd must be fastest" question directly.) +// +// mod cand_num_bigint { +// use super::*; +// use num_bigint::BigUint; +// +// pub fn prime() -> BigUint { +// BigUint::from_bytes_le(&limbs_to_le_bytes(&PRIME.0)) +// } +// +// pub fn modinv(limbs: &[u64; LIMBS], prime: &BigUint) -> [u64; LIMBS] { +// let a = BigUint::from_bytes_le(&limbs_to_le_bytes(limbs)); +// let inv = a.modinv(prime).expect("0 < a < prime"); +// le_bytes_to_limbs(&inv.to_bytes_le()) +// } +// } +// +mod cand_crypto_bigint { + use super::*; + pub use crypto_bigint::{Odd, U3072}; + + pub fn prime() -> Odd { + Odd::new(U3072::from_words(PRIME.0)).expect("prime is odd") + } + + // `from_words`/`to_words` are const `[u64; 48]` copies (uint.rs:132/148), i.e. ~free; the + // measured cost here is the safegcd inversion itself. The exact-calc group below pre-builds + // the `U3072` inputs to prove the conversion is not what makes crypto-bigint slow. + pub fn modinv_ct(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { + let a = U3072::from_words(*limbs); + let inv = Option::from(a.invert_odd_mod(prime)).expect("0 < a < prime"); + let inv: U3072 = inv; + inv.to_words() + } + + pub fn modinv_vartime(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { + let a = U3072::from_words(*limbs); + let inv = Option::from(a.invert_odd_mod_vartime(prime)).expect("0 < a < prime"); + let inv: U3072 = inv; + inv.to_words() + } +} + +// mod cand_ruint { +// use super::*; +// use ruint::Uint; +// +// pub type U3072 = Uint<3072, LIMBS>; +// +// pub fn prime() -> U3072 { +// U3072::from_limbs(PRIME.0) +// } +// +// pub fn modinv(limbs: &[u64; LIMBS], prime: &U3072) -> [u64; LIMBS] { +// let a = U3072::from_limbs(*limbs); +// let inv = a.inv_mod(*prime).expect("0 < a < prime"); +// *inv.as_limbs() +// } +// } + +// Euclidean Lehmer ext-gcd (in-repo, fixed 48-limb), the recommended permissive replacement. +// Avoids safegcd's full-width modular cofactor updates: clean integer cofactors that grow +// gradually, one sign fixup at the end. +mod cand_lehmer { + use super::*; + + pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { + let mut out = [0u64; LIMBS]; + assert!(kaspa_math::lehmer::invert(*limbs, PRIME.0, &mut out), "0 < a < prime"); + out + } +} + +mod cand_dashu { + use super::*; + use dashu_int::{IBig, UBig, ops::ExtendedGcd}; + + pub fn prime() -> UBig { + UBig::from_le_bytes(&limbs_to_le_bytes(&PRIME.0)) + } + + // dashu has no built-in modular inverse; derive it from the extended gcd + pub fn modinv(limbs: &[u64; LIMBS], prime: &UBig) -> [u64; LIMBS] { + let a = UBig::from_le_bytes(&limbs_to_le_bytes(limbs)); + let (g, x, _) = a.gcd_ext(prime.clone()); + assert_eq!(g, UBig::ONE); + let m = IBig::from(prime.clone()); + let mut x = x % &m; + if x.sign() == dashu_int::Sign::Negative { + x += &m; + } + let inv = UBig::try_from(x).unwrap(); + le_bytes_to_limbs(&inv.to_le_bytes()) + } +} + +fn bench_candidate(group: &mut BenchmarkGroup, name: &str, inputs: &[[u64; LIMBS]], f: F) +where + F: Fn(&[u64; LIMBS]) -> [u64; LIMBS], +{ + group.bench_function(name, |b| { + b.iter(|| { + for limbs in inputs { + black_box(f(black_box(limbs))); + } + }); + }); +} + +fn bench_modinv_3072(c: &mut Criterion) { + let inputs = inputs(); + + let da_prime = cand_dashu::prime(); + let cb_prime = cand_crypto_bigint::prime(); + + // all candidates must agree with the current implementation + for limbs in &inputs { + let expected = current_malachite::modinv(limbs); + assert_eq!(cand_dashu::modinv(limbs, &da_prime), expected); + assert_eq!(cand_lehmer::modinv(limbs), expected); + assert_eq!(cand_crypto_bigint::modinv_vartime(limbs, &cb_prime), expected); + assert_eq!(cand_crypto_bigint::modinv_ct(limbs, &cb_prime), expected); + } + + let mut group = c.benchmark_group("modinv_3072_muhash_prime"); + bench_candidate(&mut group, "malachite (current)", &inputs, current_malachite::modinv); + bench_candidate(&mut group, "dashu (gcd_ext)", &inputs, |l| cand_dashu::modinv(l, &da_prime)); + bench_candidate(&mut group, "lehmer (Euclidean ext-gcd)", &inputs, cand_lehmer::modinv); + bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd vartime)", &inputs, |l| cand_crypto_bigint::modinv_vartime(l, &cb_prime)); + bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd const-time)", &inputs, |l| cand_crypto_bigint::modinv_ct(l, &cb_prime)); + group.finish(); +} + +// The exact-calculation comparison the malachite-vs-crypto-bigint question hinges on: every +// candidate's inputs are converted to its NATIVE type up front (outside the timed loop), so the +// loop times only the modular inversion, never the marshalling from our `[u64; 48]` structure. +// +// This isolates the algorithm. For crypto-bigint the conversion (`U3072::from_words`) is a free +// const limb copy anyway; for malachite the `Natural` (heap bignum) is pre-built, and even its +// `&Natural::mod_inverse(&Natural)` clones the operands internally (mod_inverse.rs:210) -- that +// allocation is intrinsic to malachite's algorithm, not our boundary, so it stays in the timing. +fn bench_modinv_3072_exact(c: &mut Criterion) { + use kaspa_math::uint::malachite_base::num::arithmetic::traits::ModInverse; + use kaspa_math::uint::malachite_nz::natural::Natural; + + let inputs = inputs(); + + // malachite native: pre-built Naturals + modulus, no limb<->Natural marshalling in the loop. + let mala_prime = Natural::from_limbs_asc(&PRIME.0); + let mala_inputs: Vec = inputs.iter().map(|l| Natural::from_limbs_asc(l)).collect(); + + // crypto-bigint native: pre-built stack U3072 inputs + Odd modulus. + let cb_prime = cand_crypto_bigint::prime(); + let cb_inputs: Vec = inputs.iter().map(|l| cand_crypto_bigint::U3072::from_words(*l)).collect(); + + // sanity: the two native paths must agree before we time them. + for (m, cb) in mala_inputs.iter().zip(&cb_inputs) { + let m_inv = m.mod_inverse(&mala_prime).expect("0 < a < prime"); + let cb_inv: cand_crypto_bigint::U3072 = Option::from(cb.invert_odd_mod_vartime(&cb_prime)).expect("0 < a < prime"); + assert_eq!(Natural::from_limbs_asc(&cb_inv.to_words()), m_inv); + } + + let mut group = c.benchmark_group("modinv_3072_exact_calc"); + group.bench_function("malachite (current)", |b| { + b.iter(|| { + for v in black_box(&mala_inputs) { + black_box(black_box(v).mod_inverse(black_box(&mala_prime))); + } + }); + }); + group.bench_function("crypto-bigint 0.7 (safegcd vartime)", |b| { + b.iter(|| { + for v in black_box(&cb_inputs) { + black_box(black_box(v).invert_odd_mod_vartime(black_box(&cb_prime))); + } + }); + }); + group.bench_function("crypto-bigint 0.7 (safegcd const-time)", |b| { + b.iter(|| { + for v in black_box(&cb_inputs) { + black_box(black_box(v).invert_odd_mod(black_box(&cb_prime))); + } + }); + }); + group.finish(); +} + +// The only other malachite usage (consensus/src/processes/sync/mod.rs) is CeilingLogBase2 on a +// u64; std covers it with no crate at all. +fn bench_ceil_log2(c: &mut Criterion) { + use kaspa_math::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; + + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let values: Vec = (0..1024).map(|_| rng.next_u64() | 1).collect(); + + for v in &values { + assert_eq!(v.ceiling_log_base_2(), 64 - (v - 1).leading_zeros() as u64); + } + + let mut group = c.benchmark_group("ceil_log2_u64"); + group.bench_function("malachite (current)", |b| { + b.iter(|| { + for &v in &values { + black_box(black_box(v).ceiling_log_base_2()); + } + }); + }); + group.bench_function("std", |b| { + b.iter(|| { + for &v in &values { + black_box(64 - (black_box(v) - 1).leading_zeros() as u64); + } + }); + }); + group.finish(); +} + +// Blue work accumulation as in ghostdag protocol.rs: summing per-block work values. +// Values are kept under 2^176 so 256 of them cannot overflow Uint192. +fn bench_blue_work_sum(c: &mut Criterion) { + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let works: Vec = (0..256).map(|_| Uint192([rng.next_u64(), rng.next_u64(), rng.next_u64() >> 16])).collect(); + let works_ubig: Vec = works.iter().map(|w| UBig::from_words(&w.0)).collect(); + + let expected: Uint192 = works.iter().copied().sum(); + assert_eq!(UBig::from_words(&expected.0), works_ubig.iter().fold(UBig::ZERO, |a, b| a + b)); + + let mut group = c.benchmark_group("blue_work_sum_uint192"); + group.bench_function("kaspa-math (current)", |b| { + b.iter(|| black_box(black_box(&works).iter().copied().sum::())); + }); + group.bench_function("dashu", |b| { + b.iter(|| black_box(black_box(&works_ubig).iter().fold(UBig::ZERO, |a, b| a + b))); + }); + group.finish(); +} + +// Difficulty window average as in difficulty.rs:190-197: Uint320 sum of promoted Uint256 +// targets, divide by window size, scale by duration ratio, clamp to max target. +fn bench_difficulty_window(c: &mut Criterion) { + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let mut buf = [0u8; 32]; + // realistic mainnet-scale targets (~2^192) + let targets: Vec = (0..256) + .map(|_| { + rng.fill_bytes(&mut buf); + buf[24..].fill(0); + Uint256::from_le_bytes(buf) + }) + .collect(); + let targets_ubig: Vec = targets.iter().map(|t| UBig::from_words(&t.0)).collect(); + + let len = targets.len() as u64; + let measured_duration: u64 = 263_000; + let expected_duration: u64 = 256_000; + let max_target = Uint320::from(Uint256::from_u64(1).wrapping_shl(255) - Uint256::from_u64(1)); + let max_target_ubig = UBig::from_words(&max_target.0); + + let current = |targets: &[Uint256]| -> Uint320 { + let sum: Uint320 = targets.iter().map(|t| Uint320::from(*t)).sum(); + (sum / len * measured_duration / expected_duration).min(max_target) + }; + let dashu = |targets: &[UBig]| -> UBig { + let sum = targets.iter().fold(UBig::ZERO, |a, b| a + b); + (sum / len * measured_duration / expected_duration).min(max_target_ubig.clone()) + }; + + assert_eq!(UBig::from_words(¤t(&targets).0), dashu(&targets_ubig)); + + let mut group = c.benchmark_group("difficulty_window_uint320"); + group.bench_function("kaspa-math (current)", |b| { + b.iter(|| black_box(current(black_box(&targets)))); + }); + group.bench_function("dashu", |b| { + b.iter(|| black_box(dashu(black_box(&targets_ubig)))); + }); + group.finish(); +} + +// calc_work as in difficulty.rs:211-221: work = (!target / (target + 1)) + 1, narrowed to Uint192. +fn bench_calc_work(c: &mut Criterion) { + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let mut buf = [0u8; 32]; + let targets: Vec = (0..256) + .map(|_| { + rng.fill_bytes(&mut buf); + buf[24..].fill(0); + Uint256::from_le_bytes(buf) + }) + .collect(); + let targets_ubig: Vec = targets.iter().map(|t| UBig::from_words(&t.0)).collect(); + let max_ubig = UBig::from_words(&Uint256::MAX.0); + + let current = |t: Uint256| -> Uint192 { ((!t / (t + 1u64)) + 1u64).try_into().expect("work < 2^192") }; + let dashu = |t: &UBig| -> UBig { (&max_ubig - t) / (t + UBig::ONE) + UBig::ONE }; + + for (t, tu) in targets.iter().zip(&targets_ubig) { + assert_eq!(UBig::from_words(¤t(*t).0), dashu(tu)); + } + + let mut group = c.benchmark_group("calc_work_uint256"); + group.bench_function("kaspa-math (current)", |b| { + b.iter(|| { + for &t in black_box(&targets) { + black_box(current(black_box(t))); + } + }); + }); + group.bench_function("dashu", |b| { + b.iter(|| { + for t in black_box(&targets_ubig) { + black_box(dashu(black_box(t))); + } + }); + }); + group.finish(); +} + +// Compact target bits codec (math/src/lib.rs:64-95). Repo-owned code with no malachite +// involvement and no candidate equivalent; benched as a regression baseline only. +fn bench_compact_bits(c: &mut Criterion) { + let mut rng = ChaCha8Rng::from_seed([42u8; 32]); + let mut buf = [0u8; 32]; + let bits: Vec = (0..256) + .map(|_| { + rng.fill_bytes(&mut buf); + buf[24..].fill(0); + Uint256::from_le_bytes(buf).compact_target_bits() + }) + .collect(); + + let mut group = c.benchmark_group("compact_bits_uint256"); + group.bench_function("decode+encode (current)", |b| { + b.iter(|| { + for &v in black_box(&bits) { + black_box(Uint256::from_compact_target_bits(black_box(v)).compact_target_bits()); + } + }); + }); + group.finish(); +} + +criterion_group!( + benches, + bench_modinv_3072, + bench_modinv_3072_exact, + bench_ceil_log2, + bench_blue_work_sum, + bench_difficulty_window, + bench_calc_work, + bench_compact_bits +); +criterion_main!(benches); From a26a815f9defec40a028580f751b448e707d0b63 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 22:39:46 +0400 Subject: [PATCH 05/13] chore(math): lehmer per-inverse op-count example examples/lehmer_stats.rs reports the operation breakdown (matrices, divides, apply-limb counts) for lehmer::invert under the lehmer-instrument feature, to show where the per-inverse work goes. Inert when the feature is off. --- math/examples/lehmer_stats.rs | 65 +++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 math/examples/lehmer_stats.rs diff --git a/math/examples/lehmer_stats.rs b/math/examples/lehmer_stats.rs new file mode 100644 index 0000000000..fd9cfcaade --- /dev/null +++ b/math/examples/lehmer_stats.rs @@ -0,0 +1,65 @@ +//! Per-inverse operation-count breakdown for `lehmer::invert`, to find where the +//! work goes. Run with: +//! +//! cargo run -q --release --example lehmer_stats --features lehmer-instrument +//! +//! Build-independent (the counts come from the algorithm, not codegen). + +#[cfg(not(feature = "lehmer-instrument"))] +fn main() { + eprintln!("re-run with --features lehmer-instrument"); +} + +#[cfg(feature = "lehmer-instrument")] +fn main() { + use core::sync::atomic::Ordering::Relaxed; + use kaspa_math::Uint3072; + use kaspa_math::lehmer::{instrument as ins, invert}; + use rand_chacha::{ + ChaCha8Rng, + rand_core::{RngCore, SeedableRng}, + }; + + const PRIME: Uint3072 = { + let mut max = Uint3072::MAX; + max.0[0] -= 1103717 - 1; + max + }; + + let n: u64 = std::env::args().nth(1).and_then(|s| s.parse().ok()).unwrap_or(20000); + let mut rng = ChaCha8Rng::seed_from_u64(42); + let mut buf = [0u8; Uint3072::BYTES]; + let mut out = [0u64; 48]; + + let mut done = 0u64; + while done < n { + rng.fill_bytes(&mut buf); + let v = Uint3072::from_le_bytes(buf) % PRIME; + if v.is_zero() { + continue; + } + assert!(invert(v.0, PRIME.0, &mut out)); + done += 1; + } + + let inv = ins::INVERSES.load(Relaxed) as f64; + let m = ins::MATRICES.load(Relaxed) as f64; + let fb = ins::FALLBACK_STEPS.load(Relaxed) as f64; + let div = ins::DIVIDES.load(Relaxed) as f64; + let axy = ins::APPLY_XY_LIMBS.load(Relaxed) as f64; + let at = ins::APPLY_T_LIMBS.load(Relaxed) as f64; + let norm = ins::GCD_DIV_NORM.load(Relaxed) as f64; + + println!("inverses measured : {inv:.0}"); + println!("per inverse:"); + println!(" matrices (apply) : {:.1}", m / inv); + println!(" fallback steps : {:.2}", fb / inv); + println!(" divides (gcd_div): {:.1}", div / inv); + println!(" apply_xy limbs : {:.0} (avg len {:.1}/matrix)", axy / inv, axy / m); + println!(" apply_t limbs : {:.0} (avg len {:.1}/matrix)", at / inv, at / m); + println!(" divides/matrix : {:.1}", div / m); + // Raw count + a high-precision rate: this is ~0.0019/inverse, which a `{:.1}` + // format would misleadingly round to 0.0 (it is rare, NOT unreachable). + println!(" gcd_div q>d1 : {norm:.0} total, {:.5}/inverse ({:.2e} of divides)", norm / inv, norm / div); + println!(" bits/matrix : {:.1} (3072 / matrices)", 3072.0 / (m / inv)); +} From 0d88a59eff6a8a265763410f08088c540c9afc22 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 23:00:29 +0400 Subject: [PATCH 06/13] ci: bump CI and docker toolchains to 1.95.0 Match the workspace MSRV: the Lints job pinned rustc 1.93.0 and the docker chef images pinned 1.91, both below 1.95.0, so cargo aborted with "requires rustc 1.95.0". rust:1.95-alpine3.23 confirmed on Docker Hub. --- .github/workflows/ci.yaml | 2 +- docker/Dockerfile.kaspa-wallet | 2 +- docker/Dockerfile.kaspad | 2 +- docker/Dockerfile.rothschild | 2 +- docker/Dockerfile.simpa | 2 +- docker/Dockerfile.stratum-bridge | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 4fcb5b25df..afc48948aa 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -181,7 +181,7 @@ jobs: - name: Install stable toolchain uses: dtolnay/rust-toolchain@stable with: - toolchain: 1.93.0 + toolchain: 1.95.0 components: rustfmt, clippy - name: Set up cache diff --git a/docker/Dockerfile.kaspa-wallet b/docker/Dockerfile.kaspa-wallet index dd2dc4bc31..61fc7a1c34 100644 --- a/docker/Dockerfile.kaspa-wallet +++ b/docker/Dockerfile.kaspa-wallet @@ -1,5 +1,5 @@ # ---------------------------------------- Chef image ------------------------------------------- -FROM rust:1.91-alpine3.23 AS chef +FROM rust:1.95-alpine3.23 AS chef RUN apk --no-cache add \ musl-dev \ protobuf-dev \ diff --git a/docker/Dockerfile.kaspad b/docker/Dockerfile.kaspad index 42d8e9f134..0bb2bb9352 100644 --- a/docker/Dockerfile.kaspad +++ b/docker/Dockerfile.kaspad @@ -1,5 +1,5 @@ # ---------------------------------------- Chef image ------------------------------------------- -FROM rust:1.91-alpine3.23 AS chef +FROM rust:1.95-alpine3.23 AS chef RUN apk --no-cache add \ musl-dev \ protobuf-dev \ diff --git a/docker/Dockerfile.rothschild b/docker/Dockerfile.rothschild index 5a1c5ab86f..00001f8c6f 100644 --- a/docker/Dockerfile.rothschild +++ b/docker/Dockerfile.rothschild @@ -1,5 +1,5 @@ # ---------------------------------------- Chef image ------------------------------------------- -FROM rust:1.91-alpine3.23 AS chef +FROM rust:1.95-alpine3.23 AS chef RUN apk --no-cache add \ musl-dev \ protobuf-dev \ diff --git a/docker/Dockerfile.simpa b/docker/Dockerfile.simpa index 8d5d7528f2..9a18d02f18 100644 --- a/docker/Dockerfile.simpa +++ b/docker/Dockerfile.simpa @@ -1,5 +1,5 @@ # ---------------------------------------- Chef image ------------------------------------------- -FROM rust:1.91-alpine3.23 AS chef +FROM rust:1.95-alpine3.23 AS chef RUN apk --no-cache add \ musl-dev \ protobuf-dev \ diff --git a/docker/Dockerfile.stratum-bridge b/docker/Dockerfile.stratum-bridge index d95ff49dee..7d8b20f166 100644 --- a/docker/Dockerfile.stratum-bridge +++ b/docker/Dockerfile.stratum-bridge @@ -1,5 +1,5 @@ # ---------------------------------------- Chef image ------------------------------------------- -FROM rust:1.91-alpine3.23 AS chef +FROM rust:1.95-alpine3.23 AS chef RUN apk --no-cache add \ musl-dev \ protobuf-dev \ From 74372d82c877c5be4187dee60ecc3e8f849ada62 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sat, 13 Jun 2026 23:39:15 +0400 Subject: [PATCH 07/13] bench(math): trim modinv bench to lehmer vs malachite, drop candidate deps Comment out the dashu/crypto-bigint/num-bigint/ruint candidate arms and remove the dashu-int + crypto-bigint dev-deps; drop the dashu-comparison production-pattern groups. Bench now covers only the two malachite ops replaced on this branch (mod_inverse vs lehmer, ceil_log2 vs std). --- Cargo.lock | 76 ------- math/Cargo.toml | 17 +- math/benches/candidates.rs | 409 ++++++++++--------------------------- 3 files changed, 120 insertions(+), 382 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index e8deb78770..4b2a316269 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1102,12 +1102,6 @@ version = "0.7.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a1d728cc89cf3aee9ff92b05e62b19ee65a02b5702cff7d5a377e32c6ae29d8d" -[[package]] -name = "cmov" -version = "0.5.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0c9ea0ac24bc397ab3c98583a3c9ba74fa56b09a4449bbe172b9b1ddb016027a" - [[package]] name = "cobs" version = "0.3.0" @@ -1211,12 +1205,6 @@ dependencies = [ "libc", ] -[[package]] -name = "cpubits" -version = "0.1.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "15b85f9c39137c3a891689859392b1bd49812121d0d61c9caf00d46ed5ce06ae" - [[package]] name = "cpufeatures" version = "0.2.17" @@ -1343,18 +1331,6 @@ version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7a81dae078cea95a014a339291cec439d2f232ebe854a9d672b796c6afafa9b7" -[[package]] -name = "crypto-bigint" -version = "0.7.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "42a0d26b245348befa0c121944541476763dcc46ede886c88f9d12e1697d27c3" -dependencies = [ - "cpubits", - "ctutils", - "num-traits", - "rand_core 0.10.1", -] - [[package]] name = "crypto-common" version = "0.1.6" @@ -1407,15 +1383,6 @@ dependencies = [ "windows-sys 0.59.0", ] -[[package]] -name = "ctutils" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7d5515a3834141de9eafb9717ad39eea8247b5674e6066c404e8c4b365d2a29e" -dependencies = [ - "cmov", -] - [[package]] name = "curve25519-dalek" version = "4.1.3" @@ -1525,26 +1492,6 @@ dependencies = [ "parking_lot_core", ] -[[package]] -name = "dashu-base" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fab3f0756c8585395280bd81b384cd28bbd66d6b00e66124ecfd1f644938b38c" - -[[package]] -name = "dashu-int" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d6a93d93dc2aca9a071e6ecb2af883441e8b34357a1037bc536e1c52b0d3c713" -dependencies = [ - "cfg-if 1.0.0", - "dashu-base", - "num-modular", - "num-order", - "rustversion", - "static_assertions", -] - [[package]] name = "data-encoding" version = "2.6.0" @@ -3399,8 +3346,6 @@ version = "2.0.0" dependencies = [ "borsh", "criterion", - "crypto-bigint", - "dashu-int", "faster-hex 0.9.0", "js-sys", "kaspa-utils", @@ -5070,21 +5015,6 @@ dependencies = [ "num-traits", ] -[[package]] -name = "num-modular" -version = "0.6.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "fc41a1374056e9672221567958a66c16be12d0e2c1b408761e14d901c237d5e0" - -[[package]] -name = "num-order" -version = "1.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "537b596b97c40fcf8056d153049eb22f481c17ebce72a513ec9286e4986d1bb6" -dependencies = [ - "num-modular", -] - [[package]] name = "num-rational" version = "0.4.2" @@ -6704,12 +6634,6 @@ version = "1.2.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3" -[[package]] -name = "static_assertions" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a2eb9349b6444b326872e140eb1cf5e7c522154d69e7a0ffb0fb81c06b37543f" - [[package]] name = "strsim" version = "0.11.1" diff --git a/math/Cargo.toml b/math/Cargo.toml index 59136cff3b..a0d651e598 100644 --- a/math/Cargo.toml +++ b/math/Cargo.toml @@ -27,14 +27,15 @@ workflow-wasm.workspace = true [dev-dependencies] rand_chacha.workspace = true criterion.workspace = true -# candidate for replacing the malachite-* deps (bench-only, see benches/candidates.rs). -# the mul tuning from https://github.com/cmpute/dashu/pull/74 was also benched (via a local -# clone, the branch has a broken gitlink) and showed no effect at 48 limbs: 22.8 vs 22.4 us/inv -dashu-int = "0.4" -crypto-bigint = "0.7" -# ruled out by earlier runs (4.5x-31x slower than dashu at 3072-bit modinv): -# num-bigint = "0.4" -# ruint = "1" +# The modinv-candidate comparison is settled (Lehmer ext-gcd won), so the third-party candidate +# deps are removed. Their bench arms in benches/candidates.rs are commented out and can be revived +# by re-adding the deps: +# dashu-int = "0.4" # gcd_ext modinv (~23us) + the Uint192/256/320 production-pattern comparator. +# # the mul tuning from https://github.com/cmpute/dashu/pull/74 was also +# # benched (local clone; broken gitlink) and showed no effect at 48 limbs. +# crypto-bigint = "0.7" # safegcd vartime/const-time modinv (~114/126us, 7-8x slower) +# num-bigint = "0.4" # ~513us, ruled out +# ruint = "1" # ~72us, ruled out [features] # Bench/diagnostics only: accumulate per-inverse operation counts in `lehmer` diff --git a/math/benches/candidates.rs b/math/benches/candidates.rs index f55927a7ac..e3fdb1c5df 100644 --- a/math/benches/candidates.rs +++ b/math/benches/candidates.rs @@ -1,27 +1,24 @@ -//! Evaluates replacing the malachite-* deps, based on all production big-uint usage: +//! Benchmarks the two malachite-backed operations this branch replaces, each against its +//! replacement, to confirm the swaps are not regressions: //! -//! 1. `mod_inverse` (the only malachite-backed op, called from muhash u3072.rs): inverse of a -//! 3072-bit value modulo the MuHash prime 2^3072 - 1103717, with `[u64; 48]` limbs in/out. -//! Limb conversions are inside the measured loop, modulus conversion is precomputed (it is -//! a const at the call site). -//! 2. `ceiling_log_base_2` on u64 (consensus sync locator) vs plain std. -//! 3. The hand-rolled (repo-owned, license-clean) Uint ops actually used in production, each -//! against the dashu equivalent, to decide whether wholesale migration makes sense or the -//! macro should stay for everything but `mod_inverse`: -//! - blue work accumulation, Uint192 add/sum (ghostdag protocol.rs) -//! - difficulty window average, Uint320 sum / div u64 / mul u64 / min (difficulty.rs) -//! - calc_work, Uint256 `(!t / (t + 1)) + 1` (difficulty.rs) -//! - compact target bits codec round-trip (no candidate: repo-owned either way) +//! 1. `mod_inverse` (muhash u3072.rs): inverse of a 3072-bit value modulo the MuHash prime +//! 2^3072 - 1103717, `[u64; 48]` limbs in/out (conversions inside the measured loop, the const +//! modulus precomputed). malachite (current) vs the in-repo Lehmer ext-gcd (`lehmer.rs`), the +//! chosen replacement. The rejected third-party candidates (dashu, crypto-bigint, num-bigint, +//! ruint) are commented out and their dev-deps removed (re-add the deps and uncomment the +//! marked modules/arms below to reproduce the full sweep). +//! 2. `ceiling_log_base_2` on u64 (consensus sync locator, sync/mod.rs): malachite vs plain std. +//! +//! Per inverse = group time / N (=32). Run: `cargo bench -p kaspa-math --bench candidates`. use criterion::measurement::WallTime; use criterion::{BenchmarkGroup, Criterion, black_box, criterion_group, criterion_main}; -use dashu_int::UBig; use rand_chacha::{ ChaCha8Rng, rand_core::{RngCore, SeedableRng}, }; -use kaspa_math::{Uint192, Uint256, Uint320, Uint3072}; +use kaspa_math::Uint3072; const N: usize = 32; const LIMBS: usize = 48; @@ -33,23 +30,24 @@ const PRIME: Uint3072 = { max }; -fn limbs_to_le_bytes(limbs: &[u64; LIMBS]) -> [u8; LIMBS * 8] { - let mut bytes = [0u8; LIMBS * 8]; - for (chunk, limb) in bytes.chunks_exact_mut(8).zip(limbs) { - chunk.copy_from_slice(&limb.to_le_bytes()); - } - bytes -} - -fn le_bytes_to_limbs(bytes: &[u8]) -> [u64; LIMBS] { - let mut limbs = [0u64; LIMBS]; - for (limb, chunk) in limbs.iter_mut().zip(bytes.chunks(8)) { - let mut buf = [0u8; 8]; - buf[..chunk.len()].copy_from_slice(chunk); - *limb = u64::from_le_bytes(buf); - } - limbs -} +// Only used by the commented-out byte-marshalling candidates (dashu / num-bigint). +// fn limbs_to_le_bytes(limbs: &[u64; LIMBS]) -> [u8; LIMBS * 8] { +// let mut bytes = [0u8; LIMBS * 8]; +// for (chunk, limb) in bytes.chunks_exact_mut(8).zip(limbs) { +// chunk.copy_from_slice(&limb.to_le_bytes()); +// } +// bytes +// } +// +// fn le_bytes_to_limbs(bytes: &[u8]) -> [u64; LIMBS] { +// let mut limbs = [0u64; LIMBS]; +// for (limb, chunk) in limbs.iter_mut().zip(bytes.chunks(8)) { +// let mut buf = [0u8; 8]; +// buf[..chunk.len()].copy_from_slice(chunk); +// *limb = u64::from_le_bytes(buf); +// } +// limbs +// } fn inputs() -> Vec<[u64; LIMBS]> { let mut rng = ChaCha8Rng::from_seed([42u8; 32]); @@ -70,9 +68,24 @@ mod current_malachite { } } -// Ruled out by earlier runs (per-inverse: num-bigint ~513us, ruint ~72us vs dashu ~22us and -// malachite ~16us); kept for reference. (crypto-bigint is no longer commented out: it is wired -// into both modinv_3072 groups below to settle the "safegcd must be fastest" question directly.) +// Euclidean Lehmer ext-gcd (in-repo, fixed 48-limb), the chosen permissive replacement. +// Avoids safegcd's full-width modular cofactor updates: clean integer cofactors that grow +// gradually, one sign fixup at the end. +mod cand_lehmer { + use super::*; + + pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { + let mut out = [0u64; LIMBS]; + assert!(kaspa_math::lehmer::invert(*limbs, PRIME.0, &mut out), "0 < a < prime"); + out + } +} + +// Ruled-out third-party candidates, kept for reference (per-inverse, latest run: dashu ~23.4us, +// crypto-bigint vartime ~114us / const-time ~126us, vs lehmer ~13.6us and malachite ~15.7us; +// earlier runs: num-bigint ~513us, ruint ~72us). The `dashu-int` / `crypto-bigint` dev-deps were +// removed from Cargo.toml once Lehmer won, so these modules and their bench arms are commented out. +// Re-add the deps and uncomment to reproduce. // // mod cand_num_bigint { // use super::*; @@ -89,32 +102,31 @@ mod current_malachite { // } // } // -mod cand_crypto_bigint { - use super::*; - pub use crypto_bigint::{Odd, U3072}; - - pub fn prime() -> Odd { - Odd::new(U3072::from_words(PRIME.0)).expect("prime is odd") - } - - // `from_words`/`to_words` are const `[u64; 48]` copies (uint.rs:132/148), i.e. ~free; the - // measured cost here is the safegcd inversion itself. The exact-calc group below pre-builds - // the `U3072` inputs to prove the conversion is not what makes crypto-bigint slow. - pub fn modinv_ct(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { - let a = U3072::from_words(*limbs); - let inv = Option::from(a.invert_odd_mod(prime)).expect("0 < a < prime"); - let inv: U3072 = inv; - inv.to_words() - } - - pub fn modinv_vartime(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { - let a = U3072::from_words(*limbs); - let inv = Option::from(a.invert_odd_mod_vartime(prime)).expect("0 < a < prime"); - let inv: U3072 = inv; - inv.to_words() - } -} - +// mod cand_crypto_bigint { +// use super::*; +// pub use crypto_bigint::{Odd, U3072}; +// +// pub fn prime() -> Odd { +// Odd::new(U3072::from_words(PRIME.0)).expect("prime is odd") +// } +// +// // `from_words`/`to_words` are const `[u64; 48]` copies (uint.rs:132/148), i.e. ~free; the +// // measured cost here is the safegcd inversion itself. +// pub fn modinv_ct(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { +// let a = U3072::from_words(*limbs); +// let inv = Option::from(a.invert_odd_mod(prime)).expect("0 < a < prime"); +// let inv: U3072 = inv; +// inv.to_words() +// } +// +// pub fn modinv_vartime(limbs: &[u64; LIMBS], prime: &Odd) -> [u64; LIMBS] { +// let a = U3072::from_words(*limbs); +// let inv = Option::from(a.invert_odd_mod_vartime(prime)).expect("0 < a < prime"); +// let inv: U3072 = inv; +// inv.to_words() +// } +// } +// // mod cand_ruint { // use super::*; // use ruint::Uint; @@ -131,42 +143,29 @@ mod cand_crypto_bigint { // *inv.as_limbs() // } // } - -// Euclidean Lehmer ext-gcd (in-repo, fixed 48-limb), the recommended permissive replacement. -// Avoids safegcd's full-width modular cofactor updates: clean integer cofactors that grow -// gradually, one sign fixup at the end. -mod cand_lehmer { - use super::*; - - pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { - let mut out = [0u64; LIMBS]; - assert!(kaspa_math::lehmer::invert(*limbs, PRIME.0, &mut out), "0 < a < prime"); - out - } -} - -mod cand_dashu { - use super::*; - use dashu_int::{IBig, UBig, ops::ExtendedGcd}; - - pub fn prime() -> UBig { - UBig::from_le_bytes(&limbs_to_le_bytes(&PRIME.0)) - } - - // dashu has no built-in modular inverse; derive it from the extended gcd - pub fn modinv(limbs: &[u64; LIMBS], prime: &UBig) -> [u64; LIMBS] { - let a = UBig::from_le_bytes(&limbs_to_le_bytes(limbs)); - let (g, x, _) = a.gcd_ext(prime.clone()); - assert_eq!(g, UBig::ONE); - let m = IBig::from(prime.clone()); - let mut x = x % &m; - if x.sign() == dashu_int::Sign::Negative { - x += &m; - } - let inv = UBig::try_from(x).unwrap(); - le_bytes_to_limbs(&inv.to_le_bytes()) - } -} +// +// mod cand_dashu { +// use super::*; +// use dashu_int::{IBig, UBig, ops::ExtendedGcd}; +// +// pub fn prime() -> UBig { +// UBig::from_le_bytes(&limbs_to_le_bytes(&PRIME.0)) +// } +// +// // dashu has no built-in modular inverse; derive it from the extended gcd +// pub fn modinv(limbs: &[u64; LIMBS], prime: &UBig) -> [u64; LIMBS] { +// let a = UBig::from_le_bytes(&limbs_to_le_bytes(limbs)); +// let (g, x, _) = a.gcd_ext(prime.clone()); +// assert_eq!(g, UBig::ONE); +// let m = IBig::from(prime.clone()); +// let mut x = x % &m; +// if x.sign() == dashu_int::Sign::Negative { +// x += &m; +// } +// let inv = UBig::try_from(x).unwrap(); +// le_bytes_to_limbs(&inv.to_le_bytes()) +// } +// } fn bench_candidate(group: &mut BenchmarkGroup, name: &str, inputs: &[[u64; LIMBS]], f: F) where @@ -181,86 +180,34 @@ where }); } +// mod_inverse: the malachite incumbent vs the in-repo Lehmer ext-gcd replacement. fn bench_modinv_3072(c: &mut Criterion) { let inputs = inputs(); - let da_prime = cand_dashu::prime(); - let cb_prime = cand_crypto_bigint::prime(); - - // all candidates must agree with the current implementation + // lehmer (the chosen replacement) must agree with the current implementation for limbs in &inputs { let expected = current_malachite::modinv(limbs); - assert_eq!(cand_dashu::modinv(limbs, &da_prime), expected); assert_eq!(cand_lehmer::modinv(limbs), expected); - assert_eq!(cand_crypto_bigint::modinv_vartime(limbs, &cb_prime), expected); - assert_eq!(cand_crypto_bigint::modinv_ct(limbs, &cb_prime), expected); + // ruled-out candidates (deps removed): + // assert_eq!(cand_dashu::modinv(limbs, &da_prime), expected); + // assert_eq!(cand_crypto_bigint::modinv_vartime(limbs, &cb_prime), expected); + // assert_eq!(cand_crypto_bigint::modinv_ct(limbs, &cb_prime), expected); } let mut group = c.benchmark_group("modinv_3072_muhash_prime"); bench_candidate(&mut group, "malachite (current)", &inputs, current_malachite::modinv); - bench_candidate(&mut group, "dashu (gcd_ext)", &inputs, |l| cand_dashu::modinv(l, &da_prime)); bench_candidate(&mut group, "lehmer (Euclidean ext-gcd)", &inputs, cand_lehmer::modinv); - bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd vartime)", &inputs, |l| cand_crypto_bigint::modinv_vartime(l, &cb_prime)); - bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd const-time)", &inputs, |l| cand_crypto_bigint::modinv_ct(l, &cb_prime)); - group.finish(); -} - -// The exact-calculation comparison the malachite-vs-crypto-bigint question hinges on: every -// candidate's inputs are converted to its NATIVE type up front (outside the timed loop), so the -// loop times only the modular inversion, never the marshalling from our `[u64; 48]` structure. -// -// This isolates the algorithm. For crypto-bigint the conversion (`U3072::from_words`) is a free -// const limb copy anyway; for malachite the `Natural` (heap bignum) is pre-built, and even its -// `&Natural::mod_inverse(&Natural)` clones the operands internally (mod_inverse.rs:210) -- that -// allocation is intrinsic to malachite's algorithm, not our boundary, so it stays in the timing. -fn bench_modinv_3072_exact(c: &mut Criterion) { - use kaspa_math::uint::malachite_base::num::arithmetic::traits::ModInverse; - use kaspa_math::uint::malachite_nz::natural::Natural; - - let inputs = inputs(); - - // malachite native: pre-built Naturals + modulus, no limb<->Natural marshalling in the loop. - let mala_prime = Natural::from_limbs_asc(&PRIME.0); - let mala_inputs: Vec = inputs.iter().map(|l| Natural::from_limbs_asc(l)).collect(); - - // crypto-bigint native: pre-built stack U3072 inputs + Odd modulus. - let cb_prime = cand_crypto_bigint::prime(); - let cb_inputs: Vec = inputs.iter().map(|l| cand_crypto_bigint::U3072::from_words(*l)).collect(); - - // sanity: the two native paths must agree before we time them. - for (m, cb) in mala_inputs.iter().zip(&cb_inputs) { - let m_inv = m.mod_inverse(&mala_prime).expect("0 < a < prime"); - let cb_inv: cand_crypto_bigint::U3072 = Option::from(cb.invert_odd_mod_vartime(&cb_prime)).expect("0 < a < prime"); - assert_eq!(Natural::from_limbs_asc(&cb_inv.to_words()), m_inv); - } - - let mut group = c.benchmark_group("modinv_3072_exact_calc"); - group.bench_function("malachite (current)", |b| { - b.iter(|| { - for v in black_box(&mala_inputs) { - black_box(black_box(v).mod_inverse(black_box(&mala_prime))); - } - }); - }); - group.bench_function("crypto-bigint 0.7 (safegcd vartime)", |b| { - b.iter(|| { - for v in black_box(&cb_inputs) { - black_box(black_box(v).invert_odd_mod_vartime(black_box(&cb_prime))); - } - }); - }); - group.bench_function("crypto-bigint 0.7 (safegcd const-time)", |b| { - b.iter(|| { - for v in black_box(&cb_inputs) { - black_box(black_box(v).invert_odd_mod(black_box(&cb_prime))); - } - }); - }); + // ruled-out candidates (re-add the dev-deps + uncomment the modules above to reproduce): + // let da_prime = cand_dashu::prime(); + // let cb_prime = cand_crypto_bigint::prime(); + // bench_candidate(&mut group, "dashu (gcd_ext)", &inputs, |l| cand_dashu::modinv(l, &da_prime)); + // bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd vartime)", &inputs, |l| cand_crypto_bigint::modinv_vartime(l, &cb_prime)); + // bench_candidate(&mut group, "crypto-bigint 0.7 (safegcd const-time)", &inputs, |l| cand_crypto_bigint::modinv_ct(l, &cb_prime)); group.finish(); } -// The only other malachite usage (consensus/src/processes/sync/mod.rs) is CeilingLogBase2 on a -// u64; std covers it with no crate at all. +// The other malachite usage (consensus/src/processes/sync/mod.rs) is CeilingLogBase2 on a u64, +// replaced on this branch with std; benched here to confirm std is no slower. fn bench_ceil_log2(c: &mut Criterion) { use kaspa_math::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; @@ -289,139 +236,5 @@ fn bench_ceil_log2(c: &mut Criterion) { group.finish(); } -// Blue work accumulation as in ghostdag protocol.rs: summing per-block work values. -// Values are kept under 2^176 so 256 of them cannot overflow Uint192. -fn bench_blue_work_sum(c: &mut Criterion) { - let mut rng = ChaCha8Rng::from_seed([42u8; 32]); - let works: Vec = (0..256).map(|_| Uint192([rng.next_u64(), rng.next_u64(), rng.next_u64() >> 16])).collect(); - let works_ubig: Vec = works.iter().map(|w| UBig::from_words(&w.0)).collect(); - - let expected: Uint192 = works.iter().copied().sum(); - assert_eq!(UBig::from_words(&expected.0), works_ubig.iter().fold(UBig::ZERO, |a, b| a + b)); - - let mut group = c.benchmark_group("blue_work_sum_uint192"); - group.bench_function("kaspa-math (current)", |b| { - b.iter(|| black_box(black_box(&works).iter().copied().sum::())); - }); - group.bench_function("dashu", |b| { - b.iter(|| black_box(black_box(&works_ubig).iter().fold(UBig::ZERO, |a, b| a + b))); - }); - group.finish(); -} - -// Difficulty window average as in difficulty.rs:190-197: Uint320 sum of promoted Uint256 -// targets, divide by window size, scale by duration ratio, clamp to max target. -fn bench_difficulty_window(c: &mut Criterion) { - let mut rng = ChaCha8Rng::from_seed([42u8; 32]); - let mut buf = [0u8; 32]; - // realistic mainnet-scale targets (~2^192) - let targets: Vec = (0..256) - .map(|_| { - rng.fill_bytes(&mut buf); - buf[24..].fill(0); - Uint256::from_le_bytes(buf) - }) - .collect(); - let targets_ubig: Vec = targets.iter().map(|t| UBig::from_words(&t.0)).collect(); - - let len = targets.len() as u64; - let measured_duration: u64 = 263_000; - let expected_duration: u64 = 256_000; - let max_target = Uint320::from(Uint256::from_u64(1).wrapping_shl(255) - Uint256::from_u64(1)); - let max_target_ubig = UBig::from_words(&max_target.0); - - let current = |targets: &[Uint256]| -> Uint320 { - let sum: Uint320 = targets.iter().map(|t| Uint320::from(*t)).sum(); - (sum / len * measured_duration / expected_duration).min(max_target) - }; - let dashu = |targets: &[UBig]| -> UBig { - let sum = targets.iter().fold(UBig::ZERO, |a, b| a + b); - (sum / len * measured_duration / expected_duration).min(max_target_ubig.clone()) - }; - - assert_eq!(UBig::from_words(¤t(&targets).0), dashu(&targets_ubig)); - - let mut group = c.benchmark_group("difficulty_window_uint320"); - group.bench_function("kaspa-math (current)", |b| { - b.iter(|| black_box(current(black_box(&targets)))); - }); - group.bench_function("dashu", |b| { - b.iter(|| black_box(dashu(black_box(&targets_ubig)))); - }); - group.finish(); -} - -// calc_work as in difficulty.rs:211-221: work = (!target / (target + 1)) + 1, narrowed to Uint192. -fn bench_calc_work(c: &mut Criterion) { - let mut rng = ChaCha8Rng::from_seed([42u8; 32]); - let mut buf = [0u8; 32]; - let targets: Vec = (0..256) - .map(|_| { - rng.fill_bytes(&mut buf); - buf[24..].fill(0); - Uint256::from_le_bytes(buf) - }) - .collect(); - let targets_ubig: Vec = targets.iter().map(|t| UBig::from_words(&t.0)).collect(); - let max_ubig = UBig::from_words(&Uint256::MAX.0); - - let current = |t: Uint256| -> Uint192 { ((!t / (t + 1u64)) + 1u64).try_into().expect("work < 2^192") }; - let dashu = |t: &UBig| -> UBig { (&max_ubig - t) / (t + UBig::ONE) + UBig::ONE }; - - for (t, tu) in targets.iter().zip(&targets_ubig) { - assert_eq!(UBig::from_words(¤t(*t).0), dashu(tu)); - } - - let mut group = c.benchmark_group("calc_work_uint256"); - group.bench_function("kaspa-math (current)", |b| { - b.iter(|| { - for &t in black_box(&targets) { - black_box(current(black_box(t))); - } - }); - }); - group.bench_function("dashu", |b| { - b.iter(|| { - for t in black_box(&targets_ubig) { - black_box(dashu(black_box(t))); - } - }); - }); - group.finish(); -} - -// Compact target bits codec (math/src/lib.rs:64-95). Repo-owned code with no malachite -// involvement and no candidate equivalent; benched as a regression baseline only. -fn bench_compact_bits(c: &mut Criterion) { - let mut rng = ChaCha8Rng::from_seed([42u8; 32]); - let mut buf = [0u8; 32]; - let bits: Vec = (0..256) - .map(|_| { - rng.fill_bytes(&mut buf); - buf[24..].fill(0); - Uint256::from_le_bytes(buf).compact_target_bits() - }) - .collect(); - - let mut group = c.benchmark_group("compact_bits_uint256"); - group.bench_function("decode+encode (current)", |b| { - b.iter(|| { - for &v in black_box(&bits) { - black_box(Uint256::from_compact_target_bits(black_box(v)).compact_target_bits()); - } - }); - }); - group.finish(); -} - -criterion_group!( - benches, - bench_modinv_3072, - bench_modinv_3072_exact, - bench_ceil_log2, - bench_blue_work_sum, - bench_difficulty_window, - bench_calc_work, - bench_compact_bits -); +criterion_group!(benches, bench_modinv_3072, bench_ceil_log2); criterion_main!(benches); From a4be7a526f4804cd177eb03082c24b03798f3be3 Mon Sep 17 00:00:00 2001 From: 143672 Date: Sun, 14 Jun 2026 00:14:54 +0400 Subject: [PATCH 08/13] math: drop malachite mod_inverse for in-repo lehmer; malachite -> dev-dep - muhash 3072-bit inverse() now calls lehmer::invert; generic Uint::mod_inverse removed - malachite dropped from production deps and the uint.rs re-export; kept only as a dev-dep cross-check oracle for the lehmer/lib tests and the candidates bench - fuzz: add math/fuzz u3072 target (lehmer::invert vs num-bigint, ran clean 61k execs); remove the mod_inv arms from u128/u192/u256 - bench.rs mod_inv arm and test oracles repointed off the removed method --- Cargo.lock | 4 +-- crypto/muhash/src/u3072.rs | 7 +++- math/Cargo.toml | 7 ++-- math/benches/bench.rs | 4 ++- math/benches/candidates.rs | 14 ++++++-- math/fuzz/Cargo.toml | 6 ++++ math/fuzz/fuzz_targets/u128.rs | 48 ++------------------------ math/fuzz/fuzz_targets/u192.rs | 40 +++------------------ math/fuzz/fuzz_targets/u256.rs | 39 +++------------------ math/fuzz/fuzz_targets/u3072.rs | 61 +++++++++++++++++++++++++++++++++ math/src/lehmer.rs | 13 ++++++- math/src/lib.rs | 7 ++-- math/src/uint.rs | 28 ++++++--------- 13 files changed, 130 insertions(+), 148 deletions(-) create mode 100644 math/fuzz/fuzz_targets/u3072.rs diff --git a/Cargo.lock b/Cargo.lock index 4b2a316269..ab9b4d3236 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -7750,9 +7750,9 @@ dependencies = [ [[package]] name = "wide" -version = "1.3.0" +version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c9479f84a757f819cfab37295955906479181395de83add28f74975fde083141" +checksum = "dfdfe6a32973f2d1b268b8895845a8a96cac2f0191e72c27cc929036060dbf89" dependencies = [ "bytemuck", "safe_arch", diff --git a/crypto/muhash/src/u3072.rs b/crypto/muhash/src/u3072.rs index 3b8116beb4..b142b301b3 100644 --- a/crypto/muhash/src/u3072.rs +++ b/crypto/muhash/src/u3072.rs @@ -163,7 +163,12 @@ impl U3072 { if a == Self::zero() { return a; } - let inv = Self { limbs: Uint3072(a.limbs).mod_inverse(Self::UINT_PRIME).expect("Cannot fail, 0 < a < prime").0 }; + let inv = { + let mut out = [0u64; LIMBS]; + let invertible = kaspa_math::lehmer::invert(a.limbs, Self::UINT_PRIME.0, &mut out); + assert!(invertible, "Cannot fail, 0 < a < prime"); + Self { limbs: out } + }; if cfg!(debug_assertions) { let mut one = inv; one *= a; diff --git a/math/Cargo.toml b/math/Cargo.toml index a0d651e598..38135b53a9 100644 --- a/math/Cargo.toml +++ b/math/Cargo.toml @@ -14,8 +14,6 @@ borsh.workspace = true faster-hex.workspace = true js-sys.workspace = true kaspa-utils = { workspace = true, features = ["mem_size"] } -malachite-base.workspace = true -malachite-nz.workspace = true serde-wasm-bindgen.workspace = true serde.workspace = true thiserror.workspace = true @@ -27,6 +25,11 @@ workflow-wasm.workspace = true [dev-dependencies] rand_chacha.workspace = true criterion.workspace = true +# malachite (LGPL) is no longer a production dependency: it backed only `mod_inverse` (now the +# in-repo MIT `lehmer::invert`) and `CeilingLogBase2` (now std). It remains a dev-dependency purely +# as the cross-check oracle for `lehmer.rs` / `lib.rs` tests and the `candidates` bench baseline. +malachite-base.workspace = true +malachite-nz.workspace = true # The modinv-candidate comparison is settled (Lehmer ext-gcd won), so the third-party candidate # deps are removed. Their bench arms in benches/candidates.rs are commented out and can be revived # by re-adding the deps: diff --git a/math/benches/bench.rs b/math/benches/bench.rs index 5909463b0d..86c0387163 100644 --- a/math/benches/bench.rs +++ b/math/benches/bench.rs @@ -130,7 +130,9 @@ fn bench_uint3072(c: &mut Criterion) { uint3072_c.bench_function("mod_inv Muhash prime", |b| { b.iter(|| { for &a in &uint3072_one[..uint3072_one.len() / 4] { - black_box(a.mod_inverse(PRIME)); + // generic mod_inverse was removed; inverse the reduced value via lehmer::invert + let mut out = [0u64; 48]; + black_box(kaspa_math::lehmer::invert((a % PRIME).0, PRIME.0, &mut out)); } }); }); diff --git a/math/benches/candidates.rs b/math/benches/candidates.rs index e3fdb1c5df..b295228a23 100644 --- a/math/benches/candidates.rs +++ b/math/benches/candidates.rs @@ -62,9 +62,19 @@ fn inputs() -> Vec<[u64; LIMBS]> { mod current_malachite { use super::*; + use malachite_base::num::arithmetic::traits::ModInverse; + use malachite_nz::natural::Natural; + // The generic `Uint::mod_inverse` was removed; the malachite baseline now calls + // `Natural::mod_inverse` directly (this is exactly what the removed method did internally). pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { - Uint3072(*limbs).mod_inverse(PRIME).expect("0 < a < prime").0 + let x = Natural::from_limbs_asc(limbs); + let p = Natural::from_limbs_asc(&PRIME.0); + let inv = x.mod_inverse(p).expect("0 < a < prime"); + let mut out = [0u64; LIMBS]; + let limbs = inv.into_limbs_asc(); + out[..limbs.len()].copy_from_slice(&limbs); + out } } @@ -209,7 +219,7 @@ fn bench_modinv_3072(c: &mut Criterion) { // The other malachite usage (consensus/src/processes/sync/mod.rs) is CeilingLogBase2 on a u64, // replaced on this branch with std; benched here to confirm std is no slower. fn bench_ceil_log2(c: &mut Criterion) { - use kaspa_math::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; + use malachite_base::num::arithmetic::traits::CeilingLogBase2; let mut rng = ChaCha8Rng::from_seed([42u8; 32]); let values: Vec = (0..1024).map(|_| rng.next_u64() | 1).collect(); diff --git a/math/fuzz/Cargo.toml b/math/fuzz/Cargo.toml index 3b4743834c..7f21aa537a 100644 --- a/math/fuzz/Cargo.toml +++ b/math/fuzz/Cargo.toml @@ -38,3 +38,9 @@ name = "u192" path = "fuzz_targets/u192.rs" test = false doc = false + +[[bin]] +name = "u3072" +path = "fuzz_targets/u3072.rs" +test = false +doc = false diff --git a/math/fuzz/fuzz_targets/u128.rs b/math/fuzz/fuzz_targets/u128.rs index bdf77a2a40..4ba019bd4f 100644 --- a/math/fuzz/fuzz_targets/u128.rs +++ b/math/fuzz/fuzz_targets/u128.rs @@ -4,7 +4,6 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; use libfuzzer_sys::fuzz_target; -use std::convert::TryInto; use utils::{consume, try_opt}; construct_uint!(Uint128, 2); @@ -15,13 +14,6 @@ fn generate_ints(data: &mut &[u8]) -> Option<(Uint128, u128)> { Some((Uint128::from_le_bytes(buf), u128::from_le_bytes(buf))) } -// Consumes 16 bytes -fn generate_ints_top_bit_cleared(data: &mut &[u8]) -> Option<(Uint128, u128)> { - let mut buf = consume(data)?; - buf[15] &= 0b01111111; // clear the top/sign bit. - Some((Uint128::from_le_bytes(buf), u128::from_le_bytes(buf))) -} - fn assert_op(data: &mut &[u8], op_lib: T, op_native: U, ok_by_zero: bool) -> Option<()> where T: Fn(Uint128, Uint128) -> Uint128, @@ -131,42 +123,6 @@ fuzz_target!(|data: &[u8]| { assert_eq!(lib_bit, native_bit, "native: {native}"); } } - // mod_inv - { - // the modular inverse of 1 in Z/1Z is weird, should it be 1 or 0? - // Also, 0 never has a mod_inverse - let ((lib1, native1), (lib2, native2)) = loop { - let (lib1, native1) = try_opt!(generate_ints_top_bit_cleared(&mut data)); - let (lib2, native2) = try_opt!(generate_ints_top_bit_cleared(&mut data)); - if lib1 < lib2 && lib1 != 0u64 { - break ((lib1, native1), (lib2, native2)); - } - }; - let lib_inv = lib1.mod_inverse(lib2); - let native_inv = naive_mod_inv(native1, native2); - assert_eq!(lib_inv.is_some(), native_inv.is_some()); - if let Some(lib_inv) = lib_inv { - assert_eq!(lib_inv, native_inv.unwrap(), "lib1: {lib1}, lib2: {lib2}"); - } - } + // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse + // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. }); - -fn naive_mod_inv(x: u128, p: u128) -> Option { - let mut t = 0; - let mut newt = 1; - let p: i128 = p.try_into().unwrap(); - let mut r: i128 = p; - let mut newr: i128 = x.try_into().unwrap(); - - while newr != 0 { - let quotient = r / newr; - (t, newt) = (newt, t - quotient * newt); - (r, newr) = (newr, r - quotient * newr); - } - if r > 1 { - return None; - } else if t < 0 { - t += p; - } - Some(t.try_into().unwrap()) -} diff --git a/math/fuzz/fuzz_targets/u192.rs b/math/fuzz/fuzz_targets/u192.rs index b77d39fad7..7d7489df62 100644 --- a/math/fuzz/fuzz_targets/u192.rs +++ b/math/fuzz/fuzz_targets/u192.rs @@ -4,10 +4,8 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; use libfuzzer_sys::fuzz_target; -use num_bigint::{BigInt, BigUint}; -use num_integer::Integer; -use num_traits::{Signed, Zero}; -use std::convert::TryInto; +use num_bigint::BigUint; +use num_traits::Zero; use utils::{assert_same, consume, try_opt}; // This is important as it's a non power of two. @@ -130,36 +128,6 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv - { - // the modular inverse of 1 in Z/1Z is weird, should it be 1 or 0? - // Also, 0 never has a mod_inverse - let ((lib1, native1), (lib2, native2)) = loop { - let (lib1, native1) = try_opt!(generate_ints(&mut data)); - let (lib2, native2) = try_opt!(generate_ints(&mut data)); - if lib1 < lib2 && lib1 != 0u64 { - break ((lib1, native1), (lib2, native2)); - } - }; - let lib_inv = lib1.mod_inverse(lib2); - let native_inv = bigint_mod_inv(native1, native2); - assert_eq!(lib_inv.is_some(), native_inv.is_some()); - if let Some(lib_inv) = lib_inv { - assert_same!(lib_inv, native_inv.unwrap(), "lib1: {lib1}, lib2: {lib2}"); - } - } + // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse + // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. }); - -fn bigint_mod_inv(a: BigUint, n: BigUint) -> Option { - let a = BigInt::from(a); - let n = BigInt::from(n); - let e_gcd = a.extended_gcd(&n); - // An inverse exists iff gcd(a, n) == 1 - if e_gcd.gcd != 1u64.into() { - None - } else if e_gcd.x.is_negative() { - (e_gcd.x + n).try_into().ok() - } else { - e_gcd.x.try_into().ok() - } -} diff --git a/math/fuzz/fuzz_targets/u256.rs b/math/fuzz/fuzz_targets/u256.rs index f40146ed99..8821a86e79 100644 --- a/math/fuzz/fuzz_targets/u256.rs +++ b/math/fuzz/fuzz_targets/u256.rs @@ -4,10 +4,8 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; use libfuzzer_sys::fuzz_target; -use num_bigint::{BigInt, BigUint}; -use num_integer::Integer; -use num_traits::{Signed, Zero}; -use std::convert::TryInto; +use num_bigint::BigUint; +use num_traits::Zero; use utils::{assert_same, consume, try_opt}; construct_uint!(Uint256, 4); @@ -129,35 +127,6 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv - { - // the modular inverse of 1 in Z/1Z is weird, should it be 1 or 0? - let ((lib1, native1), (lib2, native2)) = loop { - let (lib1, native1) = try_opt!(generate_ints(&mut data)); - let (lib2, native2) = try_opt!(generate_ints(&mut data)); - if lib1 < lib2 && lib1 != 0u64 { - break ((lib1, native1), (lib2, native2)); - } - }; - let lib_inv = lib1.mod_inverse(lib2); - let native_inv = bigint_mod_inv(native1, native2); - assert_eq!(lib_inv.is_some(), native_inv.is_some()); - if let Some(lib_inv) = lib_inv { - assert_same!(lib_inv, native_inv.unwrap(), "lib1: {lib1}, lib2: {lib2}"); - } - } + // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse + // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. }); - -fn bigint_mod_inv(a: BigUint, n: BigUint) -> Option { - let a = BigInt::from(a); - let n = BigInt::from(n); - let e_gcd = a.extended_gcd(&n); - // An inverse exists iff gcd(a, n) == 1 - if e_gcd.gcd != 1u64.into() { - None - } else if e_gcd.x.is_negative() { - (e_gcd.x + n).try_into().ok() - } else { - e_gcd.x.try_into().ok() - } -} diff --git a/math/fuzz/fuzz_targets/u3072.rs b/math/fuzz/fuzz_targets/u3072.rs new file mode 100644 index 0000000000..0e13583eec --- /dev/null +++ b/math/fuzz/fuzz_targets/u3072.rs @@ -0,0 +1,61 @@ +#![no_main] +// Fuzzes the in-repo 3072-bit modular inverse (`kaspa_math::lehmer::invert`) that replaced the +// malachite-backed generic `mod_inverse`. This is the only production modular inverse (muhash +// finalize, mod the MuHash prime 2^3072 - 1103717). Each nonzero reduced value is inverted and +// cross-checked bit-for-bit against a num-bigint extended-gcd reference. +use kaspa_math::{lehmer, Uint3072}; +use libfuzzer_sys::fuzz_target; +use num_bigint::{BigInt, Sign}; +use num_integer::Integer; +use num_traits::{One, Signed}; + +const LIMBS: usize = 48; +const PRIME_DIFF: u64 = 1103717; + +// MuHash prime: 2^3072 - 1103717 +fn prime_uint() -> Uint3072 { + let mut max = Uint3072::MAX; + max.0[0] -= PRIME_DIFF - 1; + max +} + +fn prime_num() -> BigInt { + let mut prime = BigInt::one(); + prime <<= 3072; + prime -= PRIME_DIFF; + prime +} + +fuzz_target!(|data: &[u8]| { + if data.len() < Uint3072::BYTES { + return; + } + let prime = prime_uint(); + let buf: [u8; Uint3072::BYTES] = data[..Uint3072::BYTES].try_into().unwrap(); + let value = Uint3072::from_le_bytes(buf) % prime; + if value.is_zero() { + return; // 0 has no modular inverse + } + + // The production inverse under test. + let mut out = [0u64; LIMBS]; + let invertible = lehmer::invert(value.0, prime.0, &mut out); + // The MuHash prime is prime, so every nonzero reduced value is invertible. + assert!(invertible, "value={value}"); + + // num-bigint extended-gcd reference, compared bit-for-bit. + let value_num = BigInt::from_bytes_le(Sign::Plus, &value.to_le_bytes()); + let expected = inverse_num(&value_num, &prime_num()); + let got = BigInt::from_bytes_le(Sign::Plus, &Uint3072(out).to_le_bytes()); + assert_eq!(got, expected, "value={value}"); +}); + +fn inverse_num(n: &BigInt, prime: &BigInt) -> BigInt { + let e_gcd = n.extended_gcd(prime); + assert!(e_gcd.gcd.is_one()); + if e_gcd.x.is_negative() { + e_gcd.x + prime + } else { + e_gcd.x + } +} diff --git a/math/src/lehmer.rs b/math/src/lehmer.rs index 13ee7e5c20..3aed8bc247 100644 --- a/math/src/lehmer.rs +++ b/math/src/lehmer.rs @@ -889,8 +889,19 @@ mod tests { invert(value, modulus, &mut out).then_some(out) } + // Malachite reference oracle. The generic `Uint::mod_inverse` was removed, so this calls + // malachite's `Natural::mod_inverse` directly (malachite is a dev-dependency, available here). fn malachite_inv(value: [u64; N], modulus: [u64; N]) -> Option<[u64; N]> { - Uint3072(value).mod_inverse(Uint3072(modulus)).map(|inv| inv.0) + use malachite_base::num::arithmetic::traits::ModInverse; + use malachite_nz::natural::Natural; + let x = Natural::from_limbs_asc(&value); + let m = Natural::from_limbs_asc(&modulus); + x.mod_inverse(m).map(|inv| { + let mut out = [0u64; N]; + let limbs = inv.into_limbs_asc(); + out[..limbs.len()].copy_from_slice(&limbs); + out + }) } fn addmod(a: Uint3072, b: Uint3072, m: Uint3072) -> Uint3072 { diff --git a/math/src/lib.rs b/math/src/lib.rs index efcaa99503..c85b49efae 100644 --- a/math/src/lib.rs +++ b/math/src/lib.rs @@ -256,10 +256,9 @@ mod ceil_log_2_tests { #[test] fn compatibility_matches_malachite() { - // Direct equivalence with the malachite `CeilingLogBase2` this replaced. Tied to the - // malachite dependency; remove alongside it once `mod_inverse` is migrated off malachite. - // Permanent correctness coverage is provided by `correctness_matches_oracle`. - use crate::uint::malachite_base::num::arithmetic::traits::CeilingLogBase2; + // Direct equivalence with the malachite `CeilingLogBase2` this replaced. malachite is a + // dev-dependency only (the oracle); permanent coverage is `correctness_matches_oracle`. + use malachite_base::num::arithmetic::traits::CeilingLogBase2; for x in sample_inputs() { assert_eq!(ceil_log_2(x), x.ceiling_log_base_2(), "x={x}"); } diff --git a/math/src/uint.rs b/math/src/uint.rs index 0a07ef3ef5..c59f25db70 100644 --- a/math/src/uint.rs +++ b/math/src/uint.rs @@ -1,5 +1,5 @@ #[doc(hidden)] -pub use {faster_hex, malachite_base, malachite_nz, serde}; +pub use {faster_hex, serde}; // TODO: Add u32 support for optimization on 32 bit machines. @@ -342,23 +342,9 @@ macro_rules! construct_uint { (Self(ret), sub_copy) } - /// Assumes self < prime - #[inline] - pub fn mod_inverse(self, prime: Self) -> Option { - use $crate::uint::malachite_nz::natural::Natural; - use $crate::uint::malachite_base::num::arithmetic::traits::ModInverse; - - let x = Natural::from_limbs_asc(&self.0); - let p = Natural::from_limbs_asc(&prime.0); - let mod_inv = x.mod_inverse(p); - - mod_inv.map(|n| { - let mut res = [0u64; Self::LIMBS]; - let limbs = n.into_limbs_asc(); - res[..limbs.len()].copy_from_slice(&limbs); - Self(res) - }) - } + // The general malachite-backed `mod_inverse` was removed: the only production caller + // is muhash's 3072-bit inverse, now served by the in-repo `lehmer::invert` + // (`math/src/lehmer.rs`, MIT, faster than malachite). No other Uint size needs it. #[inline] pub fn iter_be_bits(self) -> impl ExactSizeIterator + core::iter::FusedIterator { @@ -1175,6 +1161,11 @@ mod tests { assert_eq!(u2.saturating_add(Uint128::from_u64(1)), Uint128::from_u128(u64::MAX as u128 + 1)); } + // `test_mod_inv` is commented out: it exercised the generic `mod_inverse` (Uint128), which was + // removed along with the malachite-backed implementation. The only production inverse is the + // 3072-bit MuHash case, now covered by the `lehmer.rs` unit tests (bit-for-bit vs malachite + + // `v * inv == 1` self-check) and the `math/fuzz` `u3072` target. Kept here for history. + /* #[test] fn test_mod_inv() { use core::cmp::Ordering; @@ -1214,4 +1205,5 @@ mod tests { res } } + */ } From dadfea234080f255c095657e2cbd16ab66cdbbb7 Mon Sep 17 00:00:00 2001 From: 143672 Date: Fri, 10 Jul 2026 18:52:11 +0400 Subject: [PATCH 09/13] math: generalize lehmer modinv over uint widths via traits - LimbArray/LehmerOps/LehmerInvert in lehmer.rs; construct_uint! impls LehmerOps, plus primitive u64/u128 impls (native-division oracles) - construct_uint! made $crate-hygienic (js_sys/wasm_bindgen/kaspa_utils/serde) so external crates can use it standalone - pre-trait fixed-48 snapshot kept in benches/support as a candidates bench arm for A/B and asm diffing; bench parity within noise - width parity tests: exhaustive u64 sweep, u128/Uint192/256/320 vs malachite --- crypto/muhash/src/u3072.rs | 2 +- math/benches/bench.rs | 2 +- math/benches/candidates.rs | 35 +- math/benches/support/lehmer_fixed48.rs | 822 +++++++++++++++++++++++++ math/examples/lehmer_stats.rs | 2 +- math/src/lehmer.rs | 494 +++++++++++---- math/src/uint.rs | 52 +- 7 files changed, 1246 insertions(+), 163 deletions(-) create mode 100644 math/benches/support/lehmer_fixed48.rs diff --git a/crypto/muhash/src/u3072.rs b/crypto/muhash/src/u3072.rs index b142b301b3..fbc40a2db4 100644 --- a/crypto/muhash/src/u3072.rs +++ b/crypto/muhash/src/u3072.rs @@ -165,7 +165,7 @@ impl U3072 { } let inv = { let mut out = [0u64; LIMBS]; - let invertible = kaspa_math::lehmer::invert(a.limbs, Self::UINT_PRIME.0, &mut out); + let invertible = kaspa_math::lehmer::invert::(a.limbs, Self::UINT_PRIME.0, &mut out); assert!(invertible, "Cannot fail, 0 < a < prime"); Self { limbs: out } }; diff --git a/math/benches/bench.rs b/math/benches/bench.rs index 86c0387163..63c1b2cad4 100644 --- a/math/benches/bench.rs +++ b/math/benches/bench.rs @@ -132,7 +132,7 @@ fn bench_uint3072(c: &mut Criterion) { for &a in &uint3072_one[..uint3072_one.len() / 4] { // generic mod_inverse was removed; inverse the reduced value via lehmer::invert let mut out = [0u64; 48]; - black_box(kaspa_math::lehmer::invert((a % PRIME).0, PRIME.0, &mut out)); + black_box(kaspa_math::lehmer::invert::((a % PRIME).0, PRIME.0, &mut out)); } }); }); diff --git a/math/benches/candidates.rs b/math/benches/candidates.rs index b295228a23..5f723bda2c 100644 --- a/math/benches/candidates.rs +++ b/math/benches/candidates.rs @@ -4,9 +4,11 @@ //! 1. `mod_inverse` (muhash u3072.rs): inverse of a 3072-bit value modulo the MuHash prime //! 2^3072 - 1103717, `[u64; 48]` limbs in/out (conversions inside the measured loop, the const //! modulus precomputed). malachite (current) vs the in-repo Lehmer ext-gcd (`lehmer.rs`), the -//! chosen replacement. The rejected third-party candidates (dashu, crypto-bigint, num-bigint, -//! ruint) are commented out and their dev-deps removed (re-add the deps and uncomment the -//! marked modules/arms below to reproduce the full sweep). +//! chosen replacement, plus the frozen pre-generalization fixed 48-limb snapshot +//! (`support/lehmer_fixed48.rs`) that isolates the cost of the trait/generic plumbing. +//! The rejected third-party candidates (dashu, crypto-bigint, num-bigint, ruint) are +//! commented out and their dev-deps removed (re-add the deps and uncomment the marked +//! modules/arms below to reproduce the full sweep). //! 2. `ceiling_log_base_2` on u64 (consensus sync locator, sync/mod.rs): malachite vs plain std. //! //! Per inverse = group time / N (=32). Run: `cargo bench -p kaspa-math --bench candidates`. @@ -20,6 +22,11 @@ use rand_chacha::{ use kaspa_math::Uint3072; +// Frozen pre-generalization 48-limb Lehmer implementation, an A/B baseline for +// the generic trait-based `src/lehmer.rs` (bench arm below + asm diffing). +#[path = "support/lehmer_fixed48.rs"] +mod lehmer_fixed48; + const N: usize = 32; const LIMBS: usize = 48; @@ -78,15 +85,27 @@ mod current_malachite { } } -// Euclidean Lehmer ext-gcd (in-repo, fixed 48-limb), the chosen permissive replacement. -// Avoids safegcd's full-width modular cofactor updates: clean integer cofactors that grow -// gradually, one sign fixup at the end. +// Euclidean Lehmer ext-gcd (in-repo, generic over the limb width), the chosen permissive +// replacement. Avoids safegcd's full-width modular cofactor updates: clean integer cofactors +// that grow gradually, one sign fixup at the end. mod cand_lehmer { use super::*; pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { let mut out = [0u64; LIMBS]; - assert!(kaspa_math::lehmer::invert(*limbs, PRIME.0, &mut out), "0 < a < prime"); + assert!(kaspa_math::lehmer::invert::(*limbs, PRIME.0, &mut out), "0 < a < prime"); + out + } +} + +// The pre-generalization fixed 48-limb snapshot (benches/support/lehmer_fixed48.rs): any gap +// between this arm and `cand_lehmer` is the cost of the trait/generic plumbing. +mod cand_lehmer_fixed { + use super::*; + + pub fn modinv(limbs: &[u64; LIMBS]) -> [u64; LIMBS] { + let mut out = [0u64; LIMBS]; + assert!(lehmer_fixed48::invert(*limbs, PRIME.0, &mut out), "0 < a < prime"); out } } @@ -198,6 +217,7 @@ fn bench_modinv_3072(c: &mut Criterion) { for limbs in &inputs { let expected = current_malachite::modinv(limbs); assert_eq!(cand_lehmer::modinv(limbs), expected); + assert_eq!(cand_lehmer_fixed::modinv(limbs), expected); // ruled-out candidates (deps removed): // assert_eq!(cand_dashu::modinv(limbs, &da_prime), expected); // assert_eq!(cand_crypto_bigint::modinv_vartime(limbs, &cb_prime), expected); @@ -207,6 +227,7 @@ fn bench_modinv_3072(c: &mut Criterion) { let mut group = c.benchmark_group("modinv_3072_muhash_prime"); bench_candidate(&mut group, "malachite (current)", &inputs, current_malachite::modinv); bench_candidate(&mut group, "lehmer (Euclidean ext-gcd)", &inputs, cand_lehmer::modinv); + bench_candidate(&mut group, "lehmer fixed-48 (pre-trait snapshot)", &inputs, cand_lehmer_fixed::modinv); // ruled-out candidates (re-add the dev-deps + uncomment the modules above to reproduce): // let da_prime = cand_dashu::prime(); // let cb_prime = cand_crypto_bigint::prime(); diff --git a/math/benches/support/lehmer_fixed48.rs b/math/benches/support/lehmer_fixed48.rs new file mode 100644 index 0000000000..9b451b1c05 --- /dev/null +++ b/math/benches/support/lehmer_fixed48.rs @@ -0,0 +1,822 @@ +//! Frozen pre-generalization snapshot of `src/lehmer.rs`: the fixed 48-limb +//! (Uint3072-only) implementation, kept verbatim so the candidates bench can +//! A/B it against the generic trait-based version in `src/lehmer.rs` and the +//! two binaries' asm can be diffed. Not production code; instrumentation +//! counters are compiled out. See `src/lehmer.rs` for the algorithm notes. + +use core::cmp::Ordering; +use kaspa_math::Uint3072; + +macro_rules! count { + ($name:ident += $v:expr) => {{}}; +} + +/// Limb count of the operands (3072 bits in base 2^64). +const N: usize = 48; +/// Limb count for the cofactor buffers. The cofactor of `value` is bounded in +/// magnitude by the modulus (< 2^3072, i.e. <= `N` limbs); the extra headroom +/// absorbs carry-out during a matrix application and the multi-limb fallback. +const T: usize = N + 8; +// Matrix entries from `hgcd2` are below `2^63` by construction (the half-limb +// refinement discards the least significant half limb), so the `i128`/`u128` +// accumulators in the matrix applications cannot overflow +// (`a*x + b*y + carry < 2^128` with `a, b < 2^63` and `x, y < 2^64`). + +/// Compute the multiplicative inverse of `value` modulo `modulus` (an odd +/// modulus, e.g. the MuHash prime), via Lehmer's extended GCD. `value`, +/// `modulus` and `out` are little-endian base-2^64 limbs (`N = 48` words), and +/// `value` must already be reduced (`value < modulus`). +/// +/// Returns `true` and writes the inverse into `out` when +/// `gcd(value, modulus) == 1`; returns `false` (leaving `out` unspecified) +/// otherwise. A zero `value` yields `false`. +pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { + count!(INVERSES += 1); + if is_zero(&value) { + core::hint::cold_path(); + return false; + } + + // Euclidean state on (x, y) with x >= y, plus the cofactor of `value`. + // Invariant: x ≡ -t0·value (mod m) and y ≡ t1·value (mod m) when + // `swapped == false`, and the roles flip with each swap. `x`, `y` are owned + // working buffers (they are reduced in place). + let mut x: [u64; N] = modulus; + let mut y: [u64; N] = value; + let mut x_len = top_len(&x); + let mut y_len = top_len(&y); + let mut t0 = [0u64; T]; + let mut t1 = [0u64; T]; + t1[0] = 1; + let mut t0_len = 1usize; // 0 has conventional length 1 + let mut t1_len = 1usize; + let mut swapped = false; + + // Multi-limb phase: run until y collapses to a single limb, then finish + // with a u64 extended-Euclidean tail. + while y_len > 1 { + let (x_hi, y_hi) = highest_two_words_normalized(&x, &y, x_len); + let guess = hgcd2((x_hi >> 64) as u64, x_hi as u64, (y_hi >> 64) as u64, y_hi as u64); + + if let Some((a, b, c, d)) = guess { + count!(MATRICES += 1); + (x_len, y_len) = apply_matrix_xy(&mut x, &mut y, x_len, a, b, c, d); + apply_matrix_t(&mut t0, &mut t1, &mut t0_len, &mut t1_len, a, b, c, d); + + // A partial guess can leave x < y; restore x >= y, toggling the sign. + if x_len < y_len || (x_len == y_len && cmp_prefix(&x, &y, x_len).is_lt()) { + core::mem::swap(&mut x, &mut y); + core::mem::swap(&mut x_len, &mut y_len); + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + } else { + // The guess almost always confirms a quotient, so this exact-division + // path is cold; keep it out of the hot loop body. + core::hint::cold_path(); + count!(FALLBACK_STEPS += 1); + // Guess could not confirm a quotient: one exact Euclidean step. + let (q, r) = div_rem_step(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); + x = y; + x_len = y_len; + y = r; + y_len = top_len(&y); + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + } + + // Single-limb tail. `y` is now <= 1 limb; one u64 extended Euclidean step + // replaces the remaining ~tens of `t_add_qmul` iterations. + if y_len == 1 && y[0] != 0 { + if x[1..N].iter().any(|&w| w != 0) { + let (q, r) = div_rem_step(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); + let old_y = y[0]; + x = [0u64; N]; + x[0] = old_y; + y = r; + core::mem::swap(&mut t0, &mut t1); + core::mem::swap(&mut t0_len, &mut t1_len); + swapped = !swapped; + } + + let (g, cx, cy) = gcd_ext_u64(x[0], y[0]); + + // Fold the sign of the chosen cofactor into `swapped`, then combine + // magnitudes: new_t0 = |cx|·t0 + |cy|·t1. + if cx < 0 || (cx == 0 && cy > 0) { + swapped = !swapped; + } + let mut new_t0 = [0u64; T]; + let mut new_t0_len = 1usize; + add_mul_word(&mut new_t0, &mut new_t0_len, cx.unsigned_abs(), &t0, t0_len); + add_mul_word(&mut new_t0, &mut new_t0_len, cy.unsigned_abs(), &t1, t1_len); + t0 = new_t0; + + x = [0u64; N]; + x[0] = g; + x_len = if g == 0 { 0 } else { 1 }; + } + + // For an odd prime modulus and `value` in [1, m), the gcd is 1. + if !(x_len == 1 && x[0] == 1) { + return false; + } + + // x = 1 = ±t0·value (mod m). With `swapped`, the inverse is +t0; otherwise + // it is -t0 ≡ m - t0. Reduce t0 into [0, m) first. + let mut inv = reduce_mod(&t0, &modulus); + if !swapped && !is_zero(&inv) { + inv = sub_n(&modulus, &inv); + } + *out = inv; + true +} + +/// Two-word subtract `(h1·2^64 + l1) - (h2·2^64 + l2)`, returning `(hi, lo)`. +/// Callers guarantee the minuend is the larger, so no underflow off the top. +#[inline(always)] +fn sub2(h1: u64, l1: u64, h2: u64, l2: u64) -> (u64, u64) { + let (lo, borrow) = l1.overflowing_sub(l2); + (h1.wrapping_sub(h2).wrapping_sub(borrow as u64), lo) +} + +/// Fold the next half limb up into a fresh full word, `(hi << 32) | (lo >> 32)`, +/// for the handoff from the main hgcd2 loop (working on top *words*) to the +/// half-limb refinement loop. `hi < 2^32` here, so the result fits a `u64`. +#[inline(always)] +fn fold_half(hi: u64, lo: u64) -> u64 { + (hi << 32).wrapping_add(lo >> 32) +} + +/// The half-GCD guess: from the aligned top *two* words of the remainders +/// (`(ah, al) >=~ (bh, bl)`), run Euclidean steps on the 128-bit approximation +/// and return the 2x2 reduction matrix in this module's `(a, b, c, d)` +/// convention, `(x', y') = (a·x - b·y, d·y - c·x)`. Returns `None` when no +/// quotient can be confirmed (the caller then takes one exact division step). +/// +/// A quotient of 1 is resolved by one two-word subtraction and a top-word +/// compare, so [`gcd_div`] (the hardware divide) is reached only for `q >= 2`. +/// The trailing half-limb refinement (the `HALF_LIMIT_2` loop) keeps every +/// matrix entry below `2^63`, which is why the `apply_matrix_*` accumulators +/// need no per-step overflow guard. All matrix arithmetic is `wrapping_*` +/// (entries are `< 2^63`, so it never actually wraps) to avoid the workspace +/// `overflow-checks = true` panic pads. +#[inline] +fn hgcd2(mut ah: u64, mut al: u64, mut bh: u64, mut bl: u64) -> Option<(u64, u64, u64, u64)> { + // Need at least two leading bits of headroom in both top words. The guess + // almost always confirms a quotient, so every `return None` here is cold. + if ah < 2 || bh < 2 { + core::hint::cold_path(); + return None; + } + + // Reduction matrix M = (m00 m01; m10 m11), accumulated by the same column + // updates as the full-width apply; returned remapped to (a,b,c,d). + let (mut m00, mut m01, mut m10, mut m11): (u64, u64, u64, u64); + if ah > bh || (ah == bh && al > bl) { + (ah, al) = sub2(ah, al, bh, bl); // a -= b + if ah < 2 { + return None; + } + (m01, m10) = (1, 0); + } else { + (bh, bl) = sub2(bh, bl, ah, al); // b -= a + if bh < 2 { + return None; + } + (m01, m10) = (0, 1); + } + (m00, m11) = (1, 1); + + const HALF: u32 = 32; + const HALF_LIMIT_1: u64 = 1 << HALF; + + let mut subtract_a = ah < bh; + let mut subtract_a1 = false; + let mut done = false; + loop { + if subtract_a { + subtract_a = false; + } else { + if ah == bh { + done = true; + break; + } + if ah < HALF_LIMIT_1 { + // Top words collapsed to a half limb: fold in the next half from + // the low words and finish in the refinement loop below. + ah = fold_half(ah, al); + bh = fold_half(bh, bl); + break; + } + // Subtract a -= q·b (affects the second column of M). + (ah, al) = sub2(ah, al, bh, bl); + if ah < 2 { + done = true; + break; + } + if ah <= bh { + m01 = m01.wrapping_add(m00); // q = 1: no divide + m11 = m11.wrapping_add(m10); + } else { + let n = ((ah as u128) << 64) | al as u128; + let d = ((bh as u128) << 64) | bl as u128; + let (mut q, r) = gcd_div(n, d); + ah = (r >> 64) as u64; + al = r as u64; + if ah < 2 { + m01 = m01.wrapping_add(q.wrapping_mul(m00)); // q correct, a too small + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + done = true; + break; + } + q = q.wrapping_add(1); // one subtraction already taken above + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + } + } + if ah == bh { + done = true; + break; + } + if bh < HALF_LIMIT_1 { + ah = fold_half(ah, al); + bh = fold_half(bh, bl); + subtract_a1 = true; + break; + } + // Subtract b -= q·a (affects the first column of M). + (bh, bl) = sub2(bh, bl, ah, al); + if bh < 2 { + done = true; + break; + } + if bh <= ah { + m00 = m00.wrapping_add(m01); // q = 1: no divide + m10 = m10.wrapping_add(m11); + } else { + let n = ((bh as u128) << 64) | bl as u128; + let d = ((ah as u128) << 64) | al as u128; + let (mut q, r) = gcd_div(n, d); + bh = (r >> 64) as u64; + bl = r as u64; + if bh < 2 { + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + done = true; + break; + } + q = q.wrapping_add(1); + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + } + } + + // Half-limb refinement: peel a bit more (single-word divides on the folded + // top words) until |a - b| fits in one limb + 1 bit. This is what keeps the + // matrix entries below 2^63, discarding the least significant half limb. + if !done { + const HALF_LIMIT_2: u64 = 1 << (HALF + 1); + loop { + if subtract_a1 { + subtract_a1 = false; + } else { + ah = ah.wrapping_sub(bh); + if ah < HALF_LIMIT_2 { + break; + } + if ah <= bh { + m01 = m01.wrapping_add(m00); + m11 = m11.wrapping_add(m10); + } else { + let q = ah / bh; + ah %= bh; + if ah < HALF_LIMIT_2 { + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + break; + } + let q = q.wrapping_add(1); + m01 = m01.wrapping_add(q.wrapping_mul(m00)); + m11 = m11.wrapping_add(q.wrapping_mul(m10)); + } + } + bh = bh.wrapping_sub(ah); + if bh < HALF_LIMIT_2 { + break; + } + if ah >= bh { + m00 = m00.wrapping_add(m01); + m10 = m10.wrapping_add(m11); + } else { + let q = bh / ah; + bh %= ah; + if bh < HALF_LIMIT_2 { + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + break; + } + let q = q.wrapping_add(1); + m00 = m00.wrapping_add(q.wrapping_mul(m01)); + m10 = m10.wrapping_add(q.wrapping_mul(m11)); + } + } + } + + // Remap M to this module's apply convention: + // x' = m11·x - m01·y, y' = m00·y - m10·x. + Some((m11, m01, m10, m00)) +} + +/// `(q, r) = (n / d, n % d)` for 128-bit `n, d` with `d >= 2^64`, computed from +/// hardware `u64` divisions (via [`div_2by1`]) instead of the `u128 / u128` +/// software routine. The quotient is `< 2^64` since `n < 2^128 <= d * 2^64`. +#[inline(always)] +fn gcd_div(n: u128, d: u128) -> (u64, u128) { + count!(DIVIDES += 1); + let (n1, n0) = ((n >> 64) as u64, n as u64); + let (mut d1, mut d0) = ((d >> 64) as u64, d as u64); + debug_assert!(d1 != 0); + let (q, r) = (n1 / d1, n1 % d1); + if q > d1 { + // Non-normalized divisor: when `d1`'s top word is small, the single-word + // estimate `q = n1/d1` can exceed `d1`. Rare but genuinely reachable (a + // b-step can drop the divisor below 2^32 before an a-step divides by it), + // so this is required for correctness, not merely defensive; marked cold + // because it is rare. Normalize `d1` to have its top bit set, then do one + // 128/64 division. + core::hint::cold_path(); + count!(GCD_DIV_NORM += 1); + let c = d1.leading_zeros(); + let wc = 64 - c; + let n2 = n1 >> wc; + let n1n = (n1 << c) | (n0 >> wc); + let n0n = n0 << c; + d1 = (d1 << c) | (d0 >> wc); + d0 <<= c; + let (mut q, rem1) = div_2by1(n2, n1n, d1); + let prod = (q as u128) * (d0 as u128); + let (mut t1, mut t0) = ((prod >> 64) as u64, prod as u64); + if t1 > rem1 || (t1 == rem1 && t0 > n0n) { + q -= 1; + let (nt0, br) = t0.overflowing_sub(d0); + t0 = nt0; + t1 = t1.wrapping_sub(d1).wrapping_sub(br as u64); + } + let (rr0, br) = n0n.overflowing_sub(t0); + let rr1 = rem1.wrapping_sub(t1).wrapping_sub(br as u64); + let rr = ((rr1 as u128) << 64) | rr0 as u128; + (q, rr >> c) // undo normalization + } else { + let mut q = q; + let prod = (q as u128).wrapping_mul(d0 as u128); + let (mut t1, mut t0) = ((prod >> 64) as u64, prod as u64); + if t1 > r || (t1 == r && t0 > n0) { + q -= 1; + let (nt0, br) = t0.overflowing_sub(d0); + t0 = nt0; + t1 = t1.wrapping_sub(d1).wrapping_sub(br as u64); + } + let (rr0, br) = n0.overflowing_sub(t0); + let rr1 = r.wrapping_sub(t1).wrapping_sub(br as u64); + (q, ((rr1 as u128) << 64) | rr0 as u128) + } +} + +/// Multiply-accumulate: `acc + m·w + carry`, returning `(low 64 bits, high 64 +/// bits)`. The high word is the carry into the next limb. A single widening +/// `64x64->128` multiply plus an add-with-carry; the carry is one word. +#[inline(always)] +fn mac(acc: u64, m: u64, w: u64, carry: u64) -> (u64, u64) { + let t = (m as u128).wrapping_mul(w as u128).wrapping_add(acc as u128).wrapping_add(carry as u128); + (t as u64, (t >> 64) as u64) +} + +/// Multiply-subtract-borrow: `acc - m·w - borrow`, returning `(low 64 bits, +/// borrow)`. The borrow is the magnitude subtracted off the next limb. It fits +/// in one word because callers pass `m < 2^63` (the [`hgcd2`] entry bound), so +/// `m·w + borrow < 2^127 + 2^64` and its high word stays below `2^63`. +#[inline(always)] +fn msb(acc: u64, m: u64, w: u64, borrow: u64) -> (u64, u64) { + let sub = (m as u128).wrapping_mul(w as u128).wrapping_add(borrow as u128); + let (lo, b) = acc.overflowing_sub(sub as u64); + (lo, ((sub >> 64) as u64).wrapping_add(b as u64)) +} + +/// `(x, y) <- (a·x - b·y, d·y - c·x)`, over the `len` live limbs only (the +/// current larger length), returning the new `(x_len, y_len)`. Both results are +/// non-negative and fit in `len` limbs for matrices from [`hgcd2`], so +/// the limbs at `[len..N]` (already zero) stay zero. Tracking lengths here lets +/// passes shrink with the remainders instead of always touching all `N` limbs. +/// +/// Each output uses a multiply-accumulate / multiply-subtract structure: a `mac` +/// carry chain for the additive term and an independent `msb` borrow chain for +/// the subtractive one, both single-word. Keeping the two chains independent +/// shortens the loop-carried dependency versus a fused signed-`i128` accumulator, +/// which would serialize a two-register carry plus a per-limb sign-extension +/// every step. +#[inline] +fn apply_matrix_xy(x: &mut [u64; N], y: &mut [u64; N], len: usize, a: u64, b: u64, c: u64, d: u64) -> (usize, usize) { + // A single up-front bound: once the compiler knows `len <= N`, every `[i]` + // with `i < len` below is provably in-bounds, so the per-limb bounds checks + // collapse into this one branch and no `unsafe` is needed. `len <= N` always + // holds here, so this never actually panics. + assert!(len <= N, "apply_matrix_xy: len exceeds N"); + count!(APPLY_XY_LIMBS += len as u64); + // x' = a·x - b·y: `carry_ax` accumulates a·x, `borrow_x` peels off b·y. + let mut carry_ax: u64 = 0; + let mut borrow_x: u64 = 0; + // y' = d·y - c·x. + let mut carry_dy: u64 = 0; + let mut borrow_y: u64 = 0; + + for i in 0..len { + let xi = x[i]; + let yi = y[i]; + + let (px, ncx) = mac(0, a, xi, carry_ax); + let (nx, nbx) = msb(px, b, yi, borrow_x); + carry_ax = ncx; + borrow_x = nbx; + + let (py, ncy) = mac(0, d, yi, carry_dy); + let (ny, nby) = msb(py, c, xi, borrow_y); + carry_dy = ncy; + borrow_y = nby; + + x[i] = nx; + y[i] = ny; + } + // Non-negative result fitting in `len` limbs: the positive carry-out exactly + // cancels the negative borrow-out (both name the same high limb, which is 0). + debug_assert_eq!(carry_ax, borrow_x, "apply_matrix_xy: x went negative or overflowed len limbs"); + debug_assert_eq!(carry_dy, borrow_y, "apply_matrix_xy: y went negative or overflowed len limbs"); + + let nlx = x[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + let nly = y[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + (nlx, nly) +} + +/// `(t0, t1) <- (a·t0 + b·t1, c·t0 + d·t1)`, growing by at most one limb. +/// All-positive (the cofactor matrix uses `+` signs), so `u128` accumulators +/// suffice. Lengths are tracked inline; the cell at `max_len` is always written +/// (carry or zero) to keep `t[len..] == 0` without a separate clearing pass. +#[allow(clippy::too_many_arguments)] +#[inline] +fn apply_matrix_t(t0: &mut [u64; T], t1: &mut [u64; T], t0_len: &mut usize, t1_len: &mut usize, a: u64, b: u64, c: u64, d: u64) { + let max_len = (*t0_len).max(*t1_len); + count!(APPLY_T_LIMBS += max_len as u64); + // Single up-front bound; see `apply_matrix_xy`. With `max_len < T` known, the + // body `[i]` (i < max_len), the carry write `[max_len]`, and the length scan + // (`max_len + 1 <= T`) are all provably in-bounds, collapsing the per-limb + // checks into this one branch. The cofactor magnitude is bounded by the + // modulus with `T = N + 8` limbs of headroom, so `max_len < T` always holds. + assert!(max_len < T, "cofactor exceeded T limbs"); + let (a, b, c, d) = (a as u128, b as u128, c as u128, d as u128); + + let mut c0: u128 = 0; + let mut c1: u128 = 0; + + for i in 0..max_len { + let ti0 = t0[i] as u128; + let ti1 = t1[i] as u128; + + let p_at0 = a.wrapping_mul(ti0); + let p_bt1 = b.wrapping_mul(ti1); + let p_ct0 = c.wrapping_mul(ti0); + let p_dt1 = d.wrapping_mul(ti1); + + let n0 = c0.wrapping_add(p_at0).wrapping_add(p_bt1); + let n1 = c1.wrapping_add(p_ct0).wrapping_add(p_dt1); + + t0[i] = n0 as u64; + t1[i] = n1 as u64; + c0 = n0 >> 64; + c1 = n1 >> 64; + } + // Carry-out (provably < 2^64); also clears any stale cell at `max_len`. + t0[max_len] = c0 as u64; + t1[max_len] = c1 as u64; + + let nl0 = t0[..max_len + 1].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + let nl1 = t1[..max_len + 1].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); + *t0_len = nl0.max(1); + *t1_len = nl1.max(1); +} + +/// `t0 += q · t1`, where `q` is up to `N` limbs (the exact Euclidean quotient). +fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T], t1_len: usize) { + for (i, &qi) in q.iter().enumerate() { + if qi == 0 { + continue; + } + let qi = qi as u128; + let mut carry: u64 = 0; + for (j, &t1j) in t1.iter().enumerate().take(t1_len) { + let k = i + j; + if k >= T { + debug_assert_eq!(carry, 0, "t_add_qmul: truncated nonzero limb"); + break; + } + let prod = qi.wrapping_mul(t1j as u128).wrapping_add(t0[k] as u128).wrapping_add(carry as u128); + t0[k] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = i + t1_len; + while carry > 0 && k < T { + let (s, c) = t0[k].overflowing_add(carry); + t0[k] = s; + carry = c as u64; + k += 1; + } + debug_assert!(carry == 0, "t_add_qmul: carry lost off the top"); + } + *t0_len = top_len_t(t0); +} + +/// `out += w · t` (multi-limb unsigned), used by the single-limb tail. +fn add_mul_word(out: &mut [u64; T], out_len: &mut usize, w: u64, t: &[u64; T], t_len: usize) { + if w != 0 { + let w = w as u128; + let mut carry: u64 = 0; + for j in 0..t_len { + // `j < t_len <= T` always, but the bound also lets the compiler clamp + // the trip count and skip the per-limb bounds check. + if j >= T { + break; + } + let prod = w.wrapping_mul(t[j] as u128).wrapping_add(out[j] as u128).wrapping_add(carry as u128); + out[j] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = t_len; + while carry > 0 && k < T { + let (s, c) = out[k].overflowing_add(carry); + out[k] = s; + carry = c as u64; + k += 1; + } + debug_assert!(carry == 0, "add_mul_word: carry lost off the top"); + } + *out_len = top_len_t(out); +} + +/// `(q, r) = (x / y, x % y)`. Fast path for a single-limb divisor (the common +/// case once Lehmer has shrunk `y`); otherwise [`Uint3072::div_rem`]. +fn div_rem_step(x: &[u64; N], y: &[u64; N]) -> ([u64; N], [u64; N]) { + if y[1..N].iter().all(|&w| w == 0) { + let (q, r0) = div_n_by_1(x, y[0]); + let mut r_arr = [0u64; N]; + r_arr[0] = r0; + return (q, r_arr); + } + let (q, r) = Uint3072(*x).div_rem(Uint3072(*y)); + (q.0, r.0) +} + +/// `(q, r) = (x / d, x % d)` for a nonzero single-limb divisor `d` (the common +/// case once Lehmer has collapsed `y`). The divisor is constant across all `N` +/// limbs, so a reciprocal is computed once and each limb costs two multiplies +/// instead of a hardware divide. Precomputed-reciprocal division +/// (Moller-Granlund, 2011). +fn div_n_by_1(x: &[u64; N], d: u64) -> ([u64; N], u64) { + debug_assert!(d != 0); + let mut q = [0u64; N]; + let s = d.leading_zeros(); + if s == 0 { + // `d` already has its top bit set: no normalization shift needed. + let v = invert_limb(d); + let mut r: u64 = 0; + for i in (0..N).rev() { + (q[i], r) = div_2by1_preinv(r, x[i], d, v); + } + (q, r) + } else { + // Normalize: divide `x << s` by `d << s` (same quotient); the bits shifted + // off the top of `x` seed the running remainder, and `x % d = r >> s`. + let d_norm = d << s; + let v = invert_limb(d_norm); + let mut r = x[N - 1] >> (64 - s); + for i in (0..N).rev() { + let lo = if i == 0 { 0 } else { x[i - 1] }; + let cur = (x[i] << s) | (lo >> (64 - s)); + (q[i], r) = div_2by1_preinv(r, cur, d_norm, v); + } + (q, r >> s) + } +} + +/// Reciprocal of a normalized 64-bit divisor `d` (top bit set): +/// `v = floor((2^128 - 1) / d) - 2^64`, the "3/2" reciprocal consumed by +/// [`div_2by1_preinv`] (Moller-Granlund, "Improved division by invariant +/// integers", 2011). The lone wide divide here runs once per [`div_n_by_1`] +/// call (amortized over all `N` limbs), so it is off the per-limb hot path. +#[inline] +const fn invert_limb(d: u64) -> u64 { + debug_assert!(d >> 63 == 1, "invert_limb requires a normalized divisor"); + // v = floor((2^128 - 1) / d) - 2^64 = floor(((!d)·2^64 + !0) / d): subtracting + // 2^64·d from the numerator drops the quotient by exactly 2^64, and that + // numerator's high word `!d` is `< d` for a normalized `d`, so the quotient + // fits a u64 and the software 128/64 [`div_2by1`] computes it -- avoiding the + // slow `u128 / u128` (`__udivti3`). + div_2by1(!d, u64::MAX, d).0 +} + +/// `(nh·2^64 + nl) / d -> (q, r)` via the precomputed reciprocal `di` from +/// [`invert_limb`]. Requires `d` normalized (top bit set) and `nh < d` (so the +/// quotient fits a `u64`); the returned `r` satisfies `r < d`, sustaining that +/// precondition limb to limb. One widening multiply, a branchless first +/// correction, and a rarely-taken second. All `wrapping_*` (the workspace builds +/// with `overflow-checks = true`). +#[inline(always)] +fn div_2by1_preinv(nh: u64, nl: u64, d: u64, di: u64) -> (u64, u64) { + debug_assert!(d >> 63 == 1 && nh < d); + // (qh:ql) = nh·di + (nh + 1)·2^64 + nl + let prod = (nh as u128).wrapping_mul(di as u128); + let (mut qh, ql) = ((prod >> 64) as u64, prod as u64); + let (ql, carry) = ql.overflowing_add(nl); + qh = qh.wrapping_add(nh.wrapping_add(1)).wrapping_add(carry as u64); + + let mut r = nl.wrapping_sub(qh.wrapping_mul(d)); + // First correction: if the estimate `qh` is one too high, `r` wraps above + // `ql`; `mask` is all-ones then, folding the `-1`/`+d` fixup in branch-free. + let mask = 0u64.wrapping_sub((r > ql) as u64); + qh = qh.wrapping_add(mask); + r = r.wrapping_add(mask & d); + // Second correction (rare): a still-too-small quotient leaves `r >= d`. + if r >= d { + r = r.wrapping_sub(d); + qh = qh.wrapping_add(1); + } + (qh, r) +} + +/// `(hi·2^64 + lo) / d -> (q, r)`. Knuth Algorithm D specialized to 128/64. +/// Precondition: `d != 0` and `hi < d` (so the quotient fits in a `u64`). +/// +/// Not on the per-limb hot path (that goes through [`div_n_by_1`]'s reciprocal +/// step, [`div_2by1_preinv`]): this generic divide is reached only once per +/// [`div_n_by_1`] (to seed the reciprocal, via [`invert_limb`]) and from +/// [`gcd_div`]'s rare normalization branch, so a plain software long division is +/// fine. +#[inline] +const fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) { + debug_assert!(d != 0); + debug_assert!(hi < d, "div_2by1 quotient would overflow u64"); + + let s = d.leading_zeros(); + let d = if s == 0 { d } else { d << s }; + let un32 = if s == 0 { hi } else { (hi << s) | (lo >> (64 - s)) }; + let un10 = if s == 0 { lo } else { lo << s }; + + let vn1 = d >> 32; + let vn0 = d & 0xFFFF_FFFF; + let un1 = un10 >> 32; + let un0 = un10 & 0xFFFF_FFFF; + + let mut q1 = un32 / vn1; + let mut rhat = un32 - q1 * vn1; + while q1 >= (1u64 << 32) || q1 * vn0 > (rhat << 32) | un1 { + q1 -= 1; + rhat += vn1; + if rhat >= (1u64 << 32) { + break; + } + } + + let un21 = (un32 << 32).wrapping_add(un1).wrapping_sub(q1.wrapping_mul(d)); + let mut q0 = un21 / vn1; + let mut rhat = un21 - q0 * vn1; + while q0 >= (1u64 << 32) || q0 * vn0 > (rhat << 32) | un0 { + q0 -= 1; + rhat += vn1; + if rhat >= (1u64 << 32) { + break; + } + } + + let r = (un21 << 32).wrapping_add(un0).wrapping_sub(q0.wrapping_mul(d)) >> s; + ((q1 << 32) | q0, r) +} + +/// Extended Euclidean GCD of two `u64`s: `(gcd, cx, cy)` with `cx·x + cy·y = g`. +/// Called once per inversion, so `i128` intermediates are fine. +fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { + let (mut a, mut b, mut c, mut d) = (1i128, 0i128, 0i128, 1i128); + while y != 0 { + let q = (x / y) as i128; + let r = x.wrapping_sub((q as u64).wrapping_mul(y)); + let (nc, nd) = (a.wrapping_sub(q.wrapping_mul(c)), b.wrapping_sub(q.wrapping_mul(d))); + a = c; + b = d; + c = nc; + d = nd; + x = y; + y = r; + } + (x, a as i64, b as i64) +} + +/// Reduce a cofactor (`<= N` significant limbs) into `[0, m)`. The cofactor is +/// bounded by `m` at loop exit, so a few subtractions suffice; the `div_rem` +/// branch is a defensive fallback. +fn reduce_mod(value: &[u64; T], m: &[u64; N]) -> [u64; N] { + debug_assert!(value[N..T].iter().all(|&w| w == 0), "cofactor exceeded N limbs"); + let mut out = [0u64; N]; + out.copy_from_slice(&value[..N]); + for _ in 0..8 { + if cmp_n(&out, m).is_lt() { + return out; + } + out = sub_n(&out, m); + } + (Uint3072(out) % Uint3072(*m)).0 +} + +/// Top 128-bit prefix of `x` and of `y`, both shifted left by the same amount so +/// the prefix of the larger operand `x` has its top bit set (the normalization +/// [`hgcd2`] expects). `y` is read at the same limb positions as `x`, so its +/// prefix is naturally the smaller. Requires `x` the larger operand and +/// `x_len >= 2`. +#[inline] +fn highest_two_words_normalized(x: &[u64; N], y: &[u64; N], x_len: usize) -> (u128, u128) { + debug_assert!(x_len >= 2); + let i = x_len - 1; + let lz = x[i].leading_zeros(); + let lo_idx_ok = i >= 2; + // Combine the limbs at positions (i, i-1, i-2) into the top 128 bits after a + // left shift by `lz`. `arr[i]`'s `lz` leading zeros guarantee no overflow. + let combine = |hi: u64, mid: u64, lo: u64| -> u128 { + let base = ((hi as u128) << 64) | (mid as u128); + if lz == 0 { base } else { (base << lz) | ((lo >> (64 - lz)) as u128) } + }; + let x_lo = if lo_idx_ok { x[i - 2] } else { 0 }; + let y_lo = if lo_idx_ok { y[i - 2] } else { 0 }; + (combine(x[i], x[i - 1], x_lo), combine(y[i], y[i - 1], y_lo)) +} + +#[inline] +fn is_zero(x: &[u64; N]) -> bool { + x.iter().all(|&w| w == 0) +} + +#[inline] +fn top_len(x: &[u64; N]) -> usize { + for i in (0..N).rev() { + if x[i] != 0 { + return i + 1; + } + } + 0 +} + +/// Significant length of a cofactor buffer; 0 maps to 1 by convention. +#[inline] +fn top_len_t(x: &[u64; T]) -> usize { + for i in (0..T).rev() { + if x[i] != 0 { + return i + 1; + } + } + 1 +} + +#[inline] +fn sub_n(a: &[u64; N], b: &[u64; N]) -> [u64; N] { + let mut out = [0u64; N]; + let mut borrow = 0u64; + for i in 0..N { + let (d1, b1) = a[i].overflowing_sub(b[i]); + let (d2, b2) = d1.overflowing_sub(borrow); + out[i] = d2; + borrow = (b1 | b2) as u64; + } + out +} + +#[inline] +fn cmp_n(a: &[u64; N], b: &[u64; N]) -> Ordering { + for i in (0..N).rev() { + match a[i].cmp(&b[i]) { + Ordering::Equal => continue, + other => return other, + } + } + Ordering::Equal +} + +#[inline] +fn cmp_prefix(a: &[u64; N], b: &[u64; N], n: usize) -> Ordering { + for i in (0..n).rev() { + match a[i].cmp(&b[i]) { + Ordering::Equal => continue, + other => return other, + } + } + Ordering::Equal +} diff --git a/math/examples/lehmer_stats.rs b/math/examples/lehmer_stats.rs index fd9cfcaade..a775eb9a19 100644 --- a/math/examples/lehmer_stats.rs +++ b/math/examples/lehmer_stats.rs @@ -38,7 +38,7 @@ fn main() { if v.is_zero() { continue; } - assert!(invert(v.0, PRIME.0, &mut out)); + assert!(invert::(v.0, PRIME.0, &mut out)); done += 1; } diff --git a/math/src/lehmer.rs b/math/src/lehmer.rs index 3aed8bc247..79921a29ab 100644 --- a/math/src/lehmer.rs +++ b/math/src/lehmer.rs @@ -1,5 +1,9 @@ -//! Variable-time modular inversion via Lehmer's extended-GCD algorithm, -//! specialized for the 3072-bit MuHash element and fully stack-allocated. +//! Variable-time modular inversion via Lehmer's extended-GCD algorithm, fully +//! stack-allocated and generic over the operand width (little-endian base-2^64 +//! limbs). [`LehmerInvert::lehmer_invert`] is the entry point, available on +//! every `construct_uint!` type and on primitive `u64`/`u128` (whose native +//! division doubles as a reference oracle in tests and fuzzing); the production +//! user is the 3072-bit MuHash element. //! //! This is a *Euclidean* extended GCD (true quotients), not the binary //! Bernstein-Yang "safegcd" divstep. The distinction matters: safegcd divides by @@ -28,11 +32,17 @@ //! Algorithm: Knuth, TAOCP Vol. 2, sec. 4.5.2, Algorithm L. Adapted from the //! MIT-licensed 256-bit implementation by Dean Little (Copyright (c) 2026 Dean //! Little, , `src/lehmer.rs`), -//! generalized to 48 limbs with a two-word half-GCD guess. Uses native -//! `u128`/`i128` accumulators plus this crate's `Uint3072` for the rare -//! multi-limb-divisor fallback; no external dependency. +//! generalized to arbitrary limb counts with a two-word half-GCD guess. Uses +//! native `u128`/`i128` accumulators plus the uint type's own division +//! ([`LehmerOps::div_rem`]) for the rare multi-limb-divisor fallback; no +//! external dependency. +//! +//! The limb loops are generic over the concrete `[u64; N]` buffer types +//! ([`LimbArray`]) rather than slices, so each width monomorphizes with +//! compile-time trip counts and elided bounds checks, keeping the 3072-bit +//! codegen equivalent to the pre-generalization fixed-width version (kept +//! frozen in `benches/support/lehmer_fixed48.rs` for A/B comparison). -use crate::Uint3072; use core::cmp::Ordering; /// Diagnostics-only per-inverse operation counters (enabled with the @@ -58,12 +68,117 @@ macro_rules! count { }}; } -/// Limb count of the operands (3072 bits in base 2^64). -const N: usize = 48; -/// Limb count for the cofactor buffers. The cofactor of `value` is bounded in -/// magnitude by the modulus (< 2^3072, i.e. <= `N` limbs); the extra headroom -/// absorbs carry-out during a matrix application and the multi-limb fallback. -const T: usize = N + 8; +/// Fixed-size little-endian limb buffer, blanket-implemented for every +/// `[u64; N]`. The inversion's working storage; keeping the concrete array +/// type generic (rather than passing slices around) monomorphizes the limb +/// loops per width, so trip counts and bounds stay compile-time constants. +pub trait LimbArray: Copy { + const LEN: usize; + const ZERO: Self; + fn as_slice(&self) -> &[u64]; + fn as_mut_slice(&mut self) -> &mut [u64]; +} + +impl LimbArray for [u64; N] { + const LEN: usize = N; + const ZERO: Self = [0; N]; + #[inline(always)] + fn as_slice(&self) -> &[u64] { + self + } + #[inline(always)] + fn as_mut_slice(&mut self) -> &mut [u64] { + self + } +} + +/// Extra limbs every cofactor buffer ([`LehmerOps::Cofactor`]) carries beyond +/// the operand width, absorbing carry-out during a matrix application and the +/// multi-limb fallback. Compile-time checked per width in [`invert`]. +pub const COFACTOR_HEADROOM: usize = 8; + +/// What Lehmer inversion needs from a uint type beyond plain limb arithmetic: +/// its limb representation and an exact multi-limb division (the rare fallback +/// when the two-word guess cannot confirm a quotient, and the defensive branch +/// of the final reduction). +/// +/// Implemented by `construct_uint!` for every generated type, and manually for +/// primitive `u64`/`u128`. +pub trait LehmerOps: Copy { + /// The operand representation, `[u64; N]`. + type Limbs: LimbArray; + /// The cofactor working buffer, `[u64; N + COFACTOR_HEADROOM]`. The + /// cofactor of `value` is bounded in magnitude by the modulus (<= `N` + /// limbs); the extra headroom absorbs carry-out during a matrix + /// application and the multi-limb fallback. Supplied per-impl because + /// stable Rust cannot spell the sum generically. + type Cofactor: LimbArray; + fn from_limbs(limbs: Self::Limbs) -> Self; + fn into_limbs(self) -> Self::Limbs; + /// `(self / other, self % other)`. + fn div_rem(self, other: Self) -> (Self, Self); +} + +impl LehmerOps for u64 { + type Limbs = [u64; 1]; + type Cofactor = [u64; 1 + COFACTOR_HEADROOM]; + #[inline] + fn from_limbs(limbs: Self::Limbs) -> Self { + limbs[0] + } + #[inline] + fn into_limbs(self) -> Self::Limbs { + [self] + } + #[inline] + fn div_rem(self, other: Self) -> (Self, Self) { + (self / other, self % other) + } +} + +impl LehmerOps for u128 { + type Limbs = [u64; 2]; + type Cofactor = [u64; 2 + COFACTOR_HEADROOM]; + #[inline] + fn from_limbs(limbs: Self::Limbs) -> Self { + ((limbs[1] as u128) << 64) | (limbs[0] as u128) + } + #[inline] + fn into_limbs(self) -> Self::Limbs { + [self as u64, (self >> 64) as u64] + } + #[inline] + fn div_rem(self, other: Self) -> (Self, Self) { + (self / other, self % other) + } +} + +/// Modular inversion as a method on any [`LehmerOps`] type. +pub trait LehmerInvert: Sized { + /// The multiplicative inverse of `self` modulo `modulus`, in `[0, modulus)`. + /// + /// Total over all inputs: `self` is reduced first when `self >= modulus`, + /// and `None` is returned when no inverse exists (`gcd != 1`, which covers + /// a zero residue, or `modulus < 2`). + fn lehmer_invert(self, modulus: Self) -> Option; +} + +impl LehmerInvert for U { + fn lehmer_invert(self, modulus: Self) -> Option { + let m = modulus.into_limbs(); + let ms = m.as_slice(); + if ms[0] < 2 && ms[1..].iter().all(|&w| w == 0) { + return None; + } + let mut value = self.into_limbs(); + if cmp_n(value.as_slice(), ms) != Ordering::Less { + value = self.div_rem(modulus).1.into_limbs(); + } + let mut out = U::Limbs::ZERO; + invert::(value, m, &mut out).then(|| U::from_limbs(out)) + } +} + // Matrix entries from `hgcd2` are below `2^63` by construction (the half-limb // refinement discards the least significant half limb), so the `i128`/`u128` // accumulators in the matrix applications cannot overflow @@ -71,15 +186,20 @@ const T: usize = N + 8; /// Compute the multiplicative inverse of `value` modulo `modulus` (an odd /// modulus, e.g. the MuHash prime), via Lehmer's extended GCD. `value`, -/// `modulus` and `out` are little-endian base-2^64 limbs (`N = 48` words), and -/// `value` must already be reduced (`value < modulus`). +/// `modulus` and `out` are little-endian base-2^64 limb arrays, and `value` +/// must already be reduced (`value < modulus`). +/// [`LehmerInvert::lehmer_invert`] is the totalized wrapper without the +/// reduction precondition. /// /// Returns `true` and writes the inverse into `out` when /// `gcd(value, modulus) == 1`; returns `false` (leaving `out` unspecified) /// otherwise. A zero `value` yields `false`. -pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { +pub fn invert(value: U::Limbs, modulus: U::Limbs, out: &mut U::Limbs) -> bool { + const { + assert!(U::Cofactor::LEN >= U::Limbs::LEN + COFACTOR_HEADROOM, "cofactor buffer lacks the required headroom"); + } count!(INVERSES += 1); - if is_zero(&value) { + if is_zero(value.as_slice()) { core::hint::cold_path(); return false; } @@ -88,13 +208,12 @@ pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { // Invariant: x ≡ -t0·value (mod m) and y ≡ t1·value (mod m) when // `swapped == false`, and the roles flip with each swap. `x`, `y` are owned // working buffers (they are reduced in place). - let mut x: [u64; N] = modulus; - let mut y: [u64; N] = value; - let mut x_len = top_len(&x); - let mut y_len = top_len(&y); - let mut t0 = [0u64; T]; - let mut t1 = [0u64; T]; - t1[0] = 1; + let mut x: U::Limbs = modulus; + let mut y: U::Limbs = value; + let mut x_len = top_len(x.as_slice()); + let mut y_len = top_len(y.as_slice()); + let mut t0 = U::Cofactor::ZERO; + let mut t1 = from_word::(1); let mut t0_len = 1usize; // 0 has conventional length 1 let mut t1_len = 1usize; let mut swapped = false; @@ -102,7 +221,7 @@ pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { // Multi-limb phase: run until y collapses to a single limb, then finish // with a u64 extended-Euclidean tail. while y_len > 1 { - let (x_hi, y_hi) = highest_two_words_normalized(&x, &y, x_len); + let (x_hi, y_hi) = highest_two_words_normalized(x.as_slice(), y.as_slice(), x_len); let guess = hgcd2((x_hi >> 64) as u64, x_hi as u64, (y_hi >> 64) as u64, y_hi as u64); if let Some((a, b, c, d)) = guess { @@ -111,7 +230,7 @@ pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { apply_matrix_t(&mut t0, &mut t1, &mut t0_len, &mut t1_len, a, b, c, d); // A partial guess can leave x < y; restore x >= y, toggling the sign. - if x_len < y_len || (x_len == y_len && cmp_prefix(&x, &y, x_len).is_lt()) { + if x_len < y_len || (x_len == y_len && cmp_prefix(x.as_slice(), y.as_slice(), x_len).is_lt()) { core::mem::swap(&mut x, &mut y); core::mem::swap(&mut x_len, &mut y_len); core::mem::swap(&mut t0, &mut t1); @@ -124,12 +243,12 @@ pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { core::hint::cold_path(); count!(FALLBACK_STEPS += 1); // Guess could not confirm a quotient: one exact Euclidean step. - let (q, r) = div_rem_step(&x, &y); - t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); + let (q, r) = div_rem_step::(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, q.as_slice(), &t1, t1_len); x = y; x_len = y_len; y = r; - y_len = top_len(&y); + y_len = top_len(y.as_slice()); core::mem::swap(&mut t0, &mut t1); core::mem::swap(&mut t0_len, &mut t1_len); swapped = !swapped; @@ -138,46 +257,43 @@ pub fn invert(value: [u64; N], modulus: [u64; N], out: &mut [u64; N]) -> bool { // Single-limb tail. `y` is now <= 1 limb; one u64 extended Euclidean step // replaces the remaining ~tens of `t_add_qmul` iterations. - if y_len == 1 && y[0] != 0 { - if x[1..N].iter().any(|&w| w != 0) { - let (q, r) = div_rem_step(&x, &y); - t_add_qmul(&mut t0, &mut t0_len, &q, &t1, t1_len); - let old_y = y[0]; - x = [0u64; N]; - x[0] = old_y; + if y_len == 1 && y.as_slice()[0] != 0 { + if x.as_slice()[1..].iter().any(|&w| w != 0) { + let (q, r) = div_rem_step::(&x, &y); + t_add_qmul(&mut t0, &mut t0_len, q.as_slice(), &t1, t1_len); + x = from_word(y.as_slice()[0]); y = r; core::mem::swap(&mut t0, &mut t1); core::mem::swap(&mut t0_len, &mut t1_len); swapped = !swapped; } - let (g, cx, cy) = gcd_ext_u64(x[0], y[0]); + let (g, cx, cy) = gcd_ext_u64(x.as_slice()[0], y.as_slice()[0]); // Fold the sign of the chosen cofactor into `swapped`, then combine // magnitudes: new_t0 = |cx|·t0 + |cy|·t1. if cx < 0 || (cx == 0 && cy > 0) { swapped = !swapped; } - let mut new_t0 = [0u64; T]; + let mut new_t0 = U::Cofactor::ZERO; let mut new_t0_len = 1usize; add_mul_word(&mut new_t0, &mut new_t0_len, cx.unsigned_abs(), &t0, t0_len); add_mul_word(&mut new_t0, &mut new_t0_len, cy.unsigned_abs(), &t1, t1_len); t0 = new_t0; - x = [0u64; N]; - x[0] = g; + x = from_word(g); x_len = if g == 0 { 0 } else { 1 }; } // For an odd prime modulus and `value` in [1, m), the gcd is 1. - if !(x_len == 1 && x[0] == 1) { + if !(x_len == 1 && x.as_slice()[0] == 1) { return false; } // x = 1 = ±t0·value (mod m). With `swapped`, the inverse is +t0; otherwise // it is -t0 ≡ m - t0. Reduce t0 into [0, m) first. - let mut inv = reduce_mod(&t0, &modulus); - if !swapped && !is_zero(&inv) { + let mut inv = reduce_mod::(&t0, &modulus); + if !swapped && !is_zero(inv.as_slice()) { inv = sub_n(&modulus, &inv); } *out = inv; @@ -458,8 +574,8 @@ fn msb(acc: u64, m: u64, w: u64, borrow: u64) -> (u64, u64) { /// `(x, y) <- (a·x - b·y, d·y - c·x)`, over the `len` live limbs only (the /// current larger length), returning the new `(x_len, y_len)`. Both results are /// non-negative and fit in `len` limbs for matrices from [`hgcd2`], so -/// the limbs at `[len..N]` (already zero) stay zero. Tracking lengths here lets -/// passes shrink with the remainders instead of always touching all `N` limbs. +/// the limbs above `len` (already zero) stay zero. Tracking lengths here lets +/// passes shrink with the remainders instead of always touching all limbs. /// /// Each output uses a multiply-accumulate / multiply-subtract structure: a `mac` /// carry chain for the additive term and an independent `msb` borrow chain for @@ -468,12 +584,14 @@ fn msb(acc: u64, m: u64, w: u64, borrow: u64) -> (u64, u64) { /// which would serialize a two-register carry plus a per-limb sign-extension /// every step. #[inline] -fn apply_matrix_xy(x: &mut [u64; N], y: &mut [u64; N], len: usize, a: u64, b: u64, c: u64, d: u64) -> (usize, usize) { - // A single up-front bound: once the compiler knows `len <= N`, every `[i]` - // with `i < len` below is provably in-bounds, so the per-limb bounds checks - // collapse into this one branch and no `unsafe` is needed. `len <= N` always - // holds here, so this never actually panics. - assert!(len <= N, "apply_matrix_xy: len exceeds N"); +fn apply_matrix_xy(x: &mut A, y: &mut A, len: usize, a: u64, b: u64, c: u64, d: u64) -> (usize, usize) { + let (x, y) = (x.as_mut_slice(), y.as_mut_slice()); + // A single up-front bound: once the compiler knows `len` is within both + // buffers, every `[i]` with `i < len` below is provably in-bounds, so the + // per-limb bounds checks collapse into this one branch and no `unsafe` is + // needed. `len` never exceeds the limb count here, so this never actually + // panics. + assert!(len <= x.len() && x.len() == y.len(), "apply_matrix_xy: len exceeds the limb count"); count!(APPLY_XY_LIMBS += len as u64); // x' = a·x - b·y: `carry_ax` accumulates a·x, `borrow_x` peels off b·y. let mut carry_ax: u64 = 0; @@ -515,15 +633,17 @@ fn apply_matrix_xy(x: &mut [u64; N], y: &mut [u64; N], len: usize, a: u64, b: u6 /// (carry or zero) to keep `t[len..] == 0` without a separate clearing pass. #[allow(clippy::too_many_arguments)] #[inline] -fn apply_matrix_t(t0: &mut [u64; T], t1: &mut [u64; T], t0_len: &mut usize, t1_len: &mut usize, a: u64, b: u64, c: u64, d: u64) { +fn apply_matrix_t(t0: &mut C, t1: &mut C, t0_len: &mut usize, t1_len: &mut usize, a: u64, b: u64, c: u64, d: u64) { + let (t0, t1) = (t0.as_mut_slice(), t1.as_mut_slice()); let max_len = (*t0_len).max(*t1_len); count!(APPLY_T_LIMBS += max_len as u64); - // Single up-front bound; see `apply_matrix_xy`. With `max_len < T` known, the - // body `[i]` (i < max_len), the carry write `[max_len]`, and the length scan - // (`max_len + 1 <= T`) are all provably in-bounds, collapsing the per-limb - // checks into this one branch. The cofactor magnitude is bounded by the - // modulus with `T = N + 8` limbs of headroom, so `max_len < T` always holds. - assert!(max_len < T, "cofactor exceeded T limbs"); + // Single up-front bound; see `apply_matrix_xy`. With `max_len` known to be + // within both buffers, the body `[i]` (i < max_len), the carry write + // `[max_len]`, and the length scan (`max_len + 1` limbs) are all provably + // in-bounds, collapsing the per-limb checks into this one branch. The + // cofactor magnitude is bounded by the modulus with 8 limbs of headroom, + // so the bound always holds. + assert!(max_len < t0.len() && t0.len() == t1.len(), "cofactor exceeded the cofactor buffer"); let (a, b, c, d) = (a as u128, b as u128, c as u128, d as u128); let mut c0: u128 = 0; @@ -557,7 +677,9 @@ fn apply_matrix_t(t0: &mut [u64; T], t1: &mut [u64; T], t0_len: &mut usize, t1_l } /// `t0 += q · t1`, where `q` is up to `N` limbs (the exact Euclidean quotient). -fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T], t1_len: usize) { +fn t_add_qmul(t0: &mut C, t0_len: &mut usize, q: &[u64], t1: &C, t1_len: usize) { + let (t0, t1) = (t0.as_mut_slice(), t1.as_slice()); + let cap = t0.len(); for (i, &qi) in q.iter().enumerate() { if qi == 0 { continue; @@ -566,7 +688,7 @@ fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T] let mut carry: u64 = 0; for (j, &t1j) in t1.iter().enumerate().take(t1_len) { let k = i + j; - if k >= T { + if k >= cap { debug_assert_eq!(carry, 0, "t_add_qmul: truncated nonzero limb"); break; } @@ -575,7 +697,7 @@ fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T] carry = (prod >> 64) as u64; } let mut k = i + t1_len; - while carry > 0 && k < T { + while carry > 0 && k < cap { let (s, c) = t0[k].overflowing_add(carry); t0[k] = s; carry = c as u64; @@ -587,14 +709,16 @@ fn t_add_qmul(t0: &mut [u64; T], t0_len: &mut usize, q: &[u64; N], t1: &[u64; T] } /// `out += w · t` (multi-limb unsigned), used by the single-limb tail. -fn add_mul_word(out: &mut [u64; T], out_len: &mut usize, w: u64, t: &[u64; T], t_len: usize) { +fn add_mul_word(out: &mut C, out_len: &mut usize, w: u64, t: &C, t_len: usize) { + let (out, t) = (out.as_mut_slice(), t.as_slice()); + let cap = out.len(); if w != 0 { let w = w as u128; let mut carry: u64 = 0; for j in 0..t_len { - // `j < t_len <= T` always, but the bound also lets the compiler clamp - // the trip count and skip the per-limb bounds check. - if j >= T { + // `j < t_len <= cap` always, but the bound also lets the compiler + // clamp the trip count and skip the per-limb bounds check. + if j >= cap { break; } let prod = w.wrapping_mul(t[j] as u128).wrapping_add(out[j] as u128).wrapping_add(carry as u128); @@ -602,7 +726,7 @@ fn add_mul_word(out: &mut [u64; T], out_len: &mut usize, w: u64, t: &[u64; T], t carry = (prod >> 64) as u64; } let mut k = t_len; - while carry > 0 && k < T { + while carry > 0 && k < cap { let (s, c) = out[k].overflowing_add(carry); out[k] = s; carry = c as u64; @@ -614,48 +738,51 @@ fn add_mul_word(out: &mut [u64; T], out_len: &mut usize, w: u64, t: &[u64; T], t } /// `(q, r) = (x / y, x % y)`. Fast path for a single-limb divisor (the common -/// case once Lehmer has shrunk `y`); otherwise [`Uint3072::div_rem`]. -fn div_rem_step(x: &[u64; N], y: &[u64; N]) -> ([u64; N], [u64; N]) { - if y[1..N].iter().all(|&w| w == 0) { - let (q, r0) = div_n_by_1(x, y[0]); - let mut r_arr = [0u64; N]; - r_arr[0] = r0; - return (q, r_arr); +/// case once Lehmer has shrunk `y`); otherwise the uint type's own +/// [`LehmerOps::div_rem`]. +fn div_rem_step(x: &U::Limbs, y: &U::Limbs) -> (U::Limbs, U::Limbs) { + if y.as_slice()[1..].iter().all(|&w| w == 0) { + let (q, r0) = div_n_by_1(x, y.as_slice()[0]); + return (q, from_word(r0)); } - let (q, r) = Uint3072(*x).div_rem(Uint3072(*y)); - (q.0, r.0) + let (q, r) = U::from_limbs(*x).div_rem(U::from_limbs(*y)); + (q.into_limbs(), r.into_limbs()) } /// `(q, r) = (x / d, x % d)` for a nonzero single-limb divisor `d` (the common -/// case once Lehmer has collapsed `y`). The divisor is constant across all `N` +/// case once Lehmer has collapsed `y`). The divisor is constant across all /// limbs, so a reciprocal is computed once and each limb costs two multiplies /// instead of a hardware divide. Precomputed-reciprocal division /// (Moller-Granlund, 2011). -fn div_n_by_1(x: &[u64; N], d: u64) -> ([u64; N], u64) { +fn div_n_by_1(x: &A, d: u64) -> (A, u64) { debug_assert!(d != 0); - let mut q = [0u64; N]; + let x = x.as_slice(); + let n = x.len(); + let mut q = A::ZERO; + let qs = q.as_mut_slice(); let s = d.leading_zeros(); - if s == 0 { + let r = if s == 0 { // `d` already has its top bit set: no normalization shift needed. let v = invert_limb(d); let mut r: u64 = 0; - for i in (0..N).rev() { - (q[i], r) = div_2by1_preinv(r, x[i], d, v); + for i in (0..n).rev() { + (qs[i], r) = div_2by1_preinv(r, x[i], d, v); } - (q, r) + r } else { // Normalize: divide `x << s` by `d << s` (same quotient); the bits shifted // off the top of `x` seed the running remainder, and `x % d = r >> s`. let d_norm = d << s; let v = invert_limb(d_norm); - let mut r = x[N - 1] >> (64 - s); - for i in (0..N).rev() { + let mut r = x[n - 1] >> (64 - s); + for i in (0..n).rev() { let lo = if i == 0 { 0 } else { x[i - 1] }; let cur = (x[i] << s) | (lo >> (64 - s)); - (q[i], r) = div_2by1_preinv(r, cur, d_norm, v); + (qs[i], r) = div_2by1_preinv(r, cur, d_norm, v); } - (q, r >> s) - } + r >> s + }; + (q, r) } /// Reciprocal of a normalized 64-bit divisor `d` (top bit set): @@ -769,20 +896,21 @@ fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { (x, a as i64, b as i64) } -/// Reduce a cofactor (`<= N` significant limbs) into `[0, m)`. The cofactor is -/// bounded by `m` at loop exit, so a few subtractions suffice; the `div_rem` -/// branch is a defensive fallback. -fn reduce_mod(value: &[u64; T], m: &[u64; N]) -> [u64; N] { - debug_assert!(value[N..T].iter().all(|&w| w == 0), "cofactor exceeded N limbs"); - let mut out = [0u64; N]; - out.copy_from_slice(&value[..N]); +/// Reduce a cofactor (with at most the modulus' significant limbs) into +/// `[0, m)`. The cofactor is bounded by `m` at loop exit, so a few +/// subtractions suffice; the `div_rem` branch is a defensive fallback. +fn reduce_mod(value: &U::Cofactor, m: &U::Limbs) -> U::Limbs { + let n = U::Limbs::LEN; + debug_assert!(value.as_slice()[n..].iter().all(|&w| w == 0), "cofactor exceeded the modulus width"); + let mut out = U::Limbs::ZERO; + out.as_mut_slice().copy_from_slice(&value.as_slice()[..n]); for _ in 0..8 { - if cmp_n(&out, m).is_lt() { + if cmp_n(out.as_slice(), m.as_slice()).is_lt() { return out; } out = sub_n(&out, m); } - (Uint3072(out) % Uint3072(*m)).0 + U::from_limbs(out).div_rem(U::from_limbs(*m)).1.into_limbs() } /// Top 128-bit prefix of `x` and of `y`, both shifted left by the same amount so @@ -791,7 +919,7 @@ fn reduce_mod(value: &[u64; T], m: &[u64; N]) -> [u64; N] { /// prefix is naturally the smaller. Requires `x` the larger operand and /// `x_len >= 2`. #[inline] -fn highest_two_words_normalized(x: &[u64; N], y: &[u64; N], x_len: usize) -> (u128, u128) { +fn highest_two_words_normalized(x: &[u64], y: &[u64], x_len: usize) -> (u128, u128) { debug_assert!(x_len >= 2); let i = x_len - 1; let lz = x[i].leading_zeros(); @@ -807,58 +935,53 @@ fn highest_two_words_normalized(x: &[u64; N], y: &[u64; N], x_len: usize) -> (u1 (combine(x[i], x[i - 1], x_lo), combine(y[i], y[i - 1], y_lo)) } +/// A limb buffer holding the single word `w`. #[inline] -fn is_zero(x: &[u64; N]) -> bool { +fn from_word(w: u64) -> A { + let mut a = A::ZERO; + a.as_mut_slice()[0] = w; + a +} + +#[inline] +fn is_zero(x: &[u64]) -> bool { x.iter().all(|&w| w == 0) } #[inline] -fn top_len(x: &[u64; N]) -> usize { - for i in (0..N).rev() { - if x[i] != 0 { - return i + 1; - } - } - 0 +fn top_len(x: &[u64]) -> usize { + x.iter().rposition(|&w| w != 0).map_or(0, |i| i + 1) } /// Significant length of a cofactor buffer; 0 maps to 1 by convention. #[inline] -fn top_len_t(x: &[u64; T]) -> usize { - for i in (0..T).rev() { - if x[i] != 0 { - return i + 1; - } - } - 1 +fn top_len_t(x: &[u64]) -> usize { + x.iter().rposition(|&w| w != 0).map_or(1, |i| i + 1) } #[inline] -fn sub_n(a: &[u64; N], b: &[u64; N]) -> [u64; N] { - let mut out = [0u64; N]; +fn sub_n(a: &A, b: &A) -> A { + let mut out = A::ZERO; + let (a, b, o) = (a.as_slice(), b.as_slice(), out.as_mut_slice()); let mut borrow = 0u64; - for i in 0..N { + for i in 0..o.len() { let (d1, b1) = a[i].overflowing_sub(b[i]); let (d2, b2) = d1.overflowing_sub(borrow); - out[i] = d2; + o[i] = d2; borrow = (b1 | b2) as u64; } out } +/// Full-width limb compare of two equal-length buffers. #[inline] -fn cmp_n(a: &[u64; N], b: &[u64; N]) -> Ordering { - for i in (0..N).rev() { - match a[i].cmp(&b[i]) { - Ordering::Equal => continue, - other => return other, - } - } - Ordering::Equal +fn cmp_n(a: &[u64], b: &[u64]) -> Ordering { + debug_assert_eq!(a.len(), b.len()); + cmp_prefix(a, b, a.len()) } #[inline] -fn cmp_prefix(a: &[u64; N], b: &[u64; N], n: usize) -> Ordering { +fn cmp_prefix(a: &[u64], b: &[u64], n: usize) -> Ordering { for i in (0..n).rev() { match a[i].cmp(&b[i]) { Ordering::Equal => continue, @@ -877,6 +1000,8 @@ mod tests { rand_core::{RngCore, SeedableRng}, }; + const N: usize = 48; + // MuHash prime: 2^3072 - 1103717. const MUHASH_PRIME: Uint3072 = { let mut max = Uint3072::MAX; @@ -886,21 +1011,22 @@ mod tests { fn lehmer_inv(value: [u64; N], modulus: [u64; N]) -> Option<[u64; N]> { let mut out = [0u64; N]; - invert(value, modulus, &mut out).then_some(out) + invert::(value, modulus, &mut out).then_some(out) } - // Malachite reference oracle. The generic `Uint::mod_inverse` was removed, so this calls - // malachite's `Natural::mod_inverse` directly (malachite is a dev-dependency, available here). - fn malachite_inv(value: [u64; N], modulus: [u64; N]) -> Option<[u64; N]> { + // Malachite reference oracle over any LehmerOps type. The generic `Uint::mod_inverse` was + // removed, so this calls malachite's `Natural::mod_inverse` directly (malachite is a + // dev-dependency, available here). Requires `value < modulus`, like `invert`. + fn malachite_inv(value: U, modulus: U) -> Option { use malachite_base::num::arithmetic::traits::ModInverse; use malachite_nz::natural::Natural; - let x = Natural::from_limbs_asc(&value); - let m = Natural::from_limbs_asc(&modulus); + let x = Natural::from_limbs_asc(value.into_limbs().as_slice()); + let m = Natural::from_limbs_asc(modulus.into_limbs().as_slice()); x.mod_inverse(m).map(|inv| { - let mut out = [0u64; N]; + let mut out = U::Limbs::ZERO; let limbs = inv.into_limbs_asc(); - out[..limbs.len()].copy_from_slice(&limbs); - out + out.as_mut_slice()[..limbs.len()].copy_from_slice(&limbs); + U::from_limbs(out) }) } @@ -935,9 +1061,9 @@ mod tests { if v.is_zero() { continue; } - let expected = malachite_inv(v.0, MUHASH_PRIME.0).unwrap(); + let expected = malachite_inv(v, MUHASH_PRIME).unwrap(); let got = lehmer_inv(v.0, MUHASH_PRIME.0).unwrap(); - assert_eq!(got.as_slice(), expected.as_slice(), "v={v}"); + assert_eq!(got.as_slice(), expected.0.as_slice(), "v={v}"); } } @@ -1112,7 +1238,7 @@ mod tests { 0x803912aebf95eede, ]; let got = lehmer_inv(v, MUHASH_PRIME.0).unwrap(); - assert_eq!(got.as_slice(), malachite_inv(v, MUHASH_PRIME.0).unwrap().as_slice()); + assert_eq!(got.as_slice(), malachite_inv(Uint3072(v), MUHASH_PRIME).unwrap().0.as_slice()); assert_eq!(mulmod(Uint3072(v), Uint3072(got), MUHASH_PRIME), Uint3072::from_u64(1)); } @@ -1124,4 +1250,100 @@ mod tests { assert!(hgcd2(5, 0, 0, 9).is_none()); // bh < 2 assert!(hgcd2(u64::MAX, 0, u64::MAX >> 1, 0).is_some()); // well-separated -> confirms a quotient } + + /// Extended-Euclid reference for u64 moduli, in i128 (all magnitudes fit). + fn egcd_inv_u64(v: u64, m: u64) -> Option { + let (mut old_r, mut r) = (m as i128, v as i128); + let (mut old_t, mut t) = (0i128, 1i128); + while r != 0 { + let q = old_r / r; + (old_r, r) = (r, old_r - q * r); + (old_t, t) = (t, old_t - q * t); + } + (old_r == 1).then(|| old_t.rem_euclid(m as i128) as u64) + } + + #[test] + fn u64_exhaustive_small_range() { + // Every (value, modulus) with modulus in [2, 512) and value in [0, 4m): + // agreement with the extended-Euclid oracle, covering gcd != 1 -> None + // and the value-reduction path. + for m in 2u64..512 { + for v in 0..m * 4 { + let got = v.lehmer_invert(m); + assert_eq!(got, egcd_inv_u64(v % m, m), "v={v} m={m}"); + if let Some(inv) = got { + assert_eq!((inv as u128 * (v % m) as u128 % m as u128) as u64, 1, "v={v} m={m}"); + } + } + } + } + + #[test] + fn u64_random_matches_oracle() { + let mut rng = ChaCha8Rng::seed_from_u64(17); + for _ in 0..200_000 { + let m = rng.next_u64().max(2); + let v = rng.next_u64(); + assert_eq!(v.lehmer_invert(m), egcd_inv_u64(v % m, m), "v={v} m={m}"); + } + } + + #[test] + fn u128_random_matches_malachite() { + // Primitive u128 impl vs malachite, mixing full-width and single-limb + // operands to cross the multi-limb/single-limb phase boundary. + let mut rng = ChaCha8Rng::seed_from_u64(23); + let mut sample = |wide: bool| { + if wide { ((rng.next_u64() as u128) << 64) | rng.next_u64() as u128 } else { rng.next_u64() as u128 } + }; + for i in 0..50_000u32 { + let m = sample(i % 4 != 0).max(2); + let v = sample(i % 3 != 0); + let got = v.lehmer_invert(m); + assert_eq!(got, malachite_inv(v % m, m), "v={v} m={m}"); + } + } + + #[test] + fn uint_sizes_match_malachite() { + // The construct_uint! impls at the other production widths, vs malachite. + use crate::{Uint192, Uint256, Uint320}; + let mut rng = ChaCha8Rng::seed_from_u64(29); + macro_rules! check { + ($ty:ty, $iters:expr) => {{ + let mut buf = [0u8; <$ty>::BYTES]; + for _ in 0..$iters { + rng.fill_bytes(&mut buf); + let m = <$ty>::from_le_bytes(buf); + if m < 2u64 { + continue; + } + rng.fill_bytes(&mut buf); + let v = <$ty>::from_le_bytes(buf); + assert_eq!(v.lehmer_invert(m), malachite_inv(v % m, m), "v={v} m={m}"); + } + }}; + } + check!(Uint192, 3000); + check!(Uint256, 3000); + check!(Uint320, 3000); + } + + #[test] + fn lehmer_invert_edge_semantics() { + // Totality of the trait method: degenerate moduli, unreduced values, + // non-coprime pairs. + assert_eq!(0u64.lehmer_invert(0), None); + assert_eq!(5u64.lehmer_invert(0), None); + assert_eq!(5u64.lehmer_invert(1), None); + assert_eq!(1u64.lehmer_invert(2), Some(1)); + assert_eq!(7u64.lehmer_invert(7), None); // value == modulus reduces to 0 + assert_eq!(9u64.lehmer_invert(7), 2u64.lehmer_invert(7)); // reduction path + assert_eq!(4u64.lehmer_invert(6), None); // gcd = 2 + assert_eq!(3u64.lehmer_invert(8), Some(3)); // even modulus, coprime value + assert_eq!(3u128.lehmer_invert(8), Some(3)); + assert_eq!(Uint3072::from_u64(5).lehmer_invert(Uint3072::from_u64(1)), None); + assert_eq!(Uint3072::from_u64(3).lehmer_invert(Uint3072::from_u64(8)), Some(Uint3072::from_u64(3))); + } } diff --git a/math/src/uint.rs b/math/src/uint.rs index c59f25db70..0f95c14d9f 100644 --- a/math/src/uint.rs +++ b/math/src/uint.rs @@ -1,5 +1,5 @@ #[doc(hidden)] -pub use {faster_hex, serde}; +pub use {faster_hex, js_sys, kaspa_utils, serde, wasm_bindgen}; // TODO: Add u32 support for optimization on 32 bit machines. @@ -342,9 +342,10 @@ macro_rules! construct_uint { (Self(ret), sub_copy) } - // The general malachite-backed `mod_inverse` was removed: the only production caller - // is muhash's 3072-bit inverse, now served by the in-repo `lehmer::invert` - // (`math/src/lehmer.rs`, MIT, faster than malachite). No other Uint size needs it. + // The general malachite-backed `mod_inverse` was removed; modular inversion is + // provided for every Uint by the in-repo `lehmer` module (`math/src/lehmer.rs`, MIT, + // faster than malachite) through the `LehmerInvert` extension trait, backed by the + // `LehmerOps` impl below. #[inline] pub fn iter_be_bits(self) -> impl ExactSizeIterator + core::iter::FusedIterator { @@ -418,37 +419,54 @@ macro_rules! construct_uint { } #[inline] - pub fn as_bigint(&self) -> Result { + pub fn as_bigint(&self) -> Result<$crate::uint::js_sys::BigInt, $crate::Error> { self.try_into() } #[inline] - pub fn to_bigint(self) -> Result { + pub fn to_bigint(self) -> Result<$crate::uint::js_sys::BigInt, $crate::Error> { self.try_into() } } - impl kaspa_utils::mem_size::MemSizeEstimator for $name { + impl $crate::lehmer::LehmerOps for $name { + type Limbs = [u64; $n_words]; + type Cofactor = [u64; $n_words + $crate::lehmer::COFACTOR_HEADROOM]; + #[inline] + fn from_limbs(limbs: Self::Limbs) -> Self { + Self(limbs) + } + #[inline] + fn into_limbs(self) -> Self::Limbs { + self.0 + } + #[inline] + fn div_rem(self, other: Self) -> (Self, Self) { + $name::div_rem(self, other) + } + } + + impl $crate::uint::kaspa_utils::mem_size::MemSizeEstimator for $name { fn estimate_mem_units(&self) -> usize { 1 } } - impl kaspa_utils::hex::ToHex for $name { + impl $crate::uint::kaspa_utils::hex::ToHex for $name { fn to_hex(&self) -> String { self.to_be_bytes().as_slice().to_hex() } } - impl kaspa_utils::hex::ToHex for &$name { + impl $crate::uint::kaspa_utils::hex::ToHex for &$name { fn to_hex(&self) -> String { self.to_be_bytes().as_slice().to_hex() } } - impl kaspa_utils::hex::FromHex for $name { + impl $crate::uint::kaspa_utils::hex::FromHex for $name { type Error = $crate::Error; fn from_hex(hex: &str) -> Result<$name, Self::Error> { Ok($name::from_hex(hex)?) @@ -848,7 +866,7 @@ macro_rules! construct_uint { #[inline] fn deserialize>(deserializer: D) -> Result { if deserializer.is_human_readable() { - let hex = ::deserialize(deserializer)?; + let hex = ::deserialize(deserializer)?; Ok(Self::from_hex(&hex).map_err($crate::uint::serde::de::Error::custom)?) } else { use core::{fmt, marker::PhantomData}; @@ -946,27 +964,27 @@ macro_rules! construct_uint { } - impl TryFrom<&$name> for js_sys::BigInt { + impl TryFrom<&$name> for $crate::uint::js_sys::BigInt { type Error = $crate::Error; #[inline] - fn try_from(value: &$name) -> Result { + fn try_from(value: &$name) -> Result<$crate::uint::js_sys::BigInt, Self::Error> { use $crate::wasm::*; BigInt::new(&JsValue::from_str(&format!("0x{value:x}"))).map_err(|err|$crate::Error::JsSys(Sendable(err))) } } - impl TryFrom<$name> for js_sys::BigInt { + impl TryFrom<$name> for $crate::uint::js_sys::BigInt { type Error = $crate::Error; #[inline] - fn try_from(value: $name) -> Result { + fn try_from(value: $name) -> Result<$crate::uint::js_sys::BigInt, Self::Error> { use $crate::wasm::*; BigInt::try_from(&value) } } - impl TryFrom for $name { + impl TryFrom<$crate::uint::wasm_bindgen::JsValue> for $name { type Error = $crate::Error; - fn try_from(js_value: wasm_bindgen::JsValue) -> Result { + fn try_from(js_value: $crate::uint::wasm_bindgen::JsValue) -> Result { use $crate::wasm::*; if js_value.is_string() || js_value.is_array() { From 56ec8f5365f7b05ac83ffc57cc645bc54bdbb117 Mon Sep 17 00:00:00 2001 From: 143672 Date: Fri, 10 Jul 2026 18:52:18 +0400 Subject: [PATCH 10/13] fuzz(math): recover mod_inv arms via lehmer_invert - u128/u192/u256 fuzz both value and modulus full-width against a shared num-bigint oracle in utils.rs - u128 also cross-checks the primitive u128 impl --- math/fuzz/fuzz_targets/u128.rs | 22 +++++++++++++++++++--- math/fuzz/fuzz_targets/u192.rs | 24 +++++++++++++++++++++--- math/fuzz/fuzz_targets/u256.rs | 24 +++++++++++++++++++++--- math/fuzz/fuzz_targets/u3072.rs | 2 +- math/fuzz/fuzz_targets/utils.rs | 18 ++++++++++++++++++ 5 files changed, 80 insertions(+), 10 deletions(-) diff --git a/math/fuzz/fuzz_targets/u128.rs b/math/fuzz/fuzz_targets/u128.rs index 4ba019bd4f..8e4a193e07 100644 --- a/math/fuzz/fuzz_targets/u128.rs +++ b/math/fuzz/fuzz_targets/u128.rs @@ -3,8 +3,10 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; +use kaspa_math::lehmer::LehmerInvert; use libfuzzer_sys::fuzz_target; -use utils::{consume, try_opt}; +use num_bigint::BigUint; +use utils::{bigint_mod_inv, consume, try_opt}; construct_uint!(Uint128, 2); @@ -123,6 +125,20 @@ fuzz_target!(|data: &[u8]| { assert_eq!(lib_bit, native_bit, "native: {native}"); } } - // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse - // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. + // lehmer_invert (the mod_inv arm, recovered): both the Uint128 and the primitive u128 + // impls against the shared num-bigint extended-gcd oracle. `lehmer_invert` is total, so + // both value and modulus are fuzzed freely; only the trivial-ring modulus < 2 is filtered + // (the inverse in Z/1Z is convention: the oracle says 0, lehmer says None). + { + let ((lib1, native1), (lib2, native2)) = loop { + let (lib1, native1) = try_opt!(generate_ints(&mut data)); + let (lib2, native2) = try_opt!(generate_ints(&mut data)); + if native2 >= 2 { + break ((lib1, native1), (lib2, native2)); + } + }; + let expected = bigint_mod_inv(BigUint::from(native1), BigUint::from(native2)).map(|x| u128::try_from(x).unwrap()); + assert_eq!(native1.lehmer_invert(native2), expected, "native1: {native1}, native2: {native2}"); + assert_eq!(lib1.lehmer_invert(lib2).map(|x| x.as_u128()), expected, "lib1: {lib1}, lib2: {lib2}"); + } }); diff --git a/math/fuzz/fuzz_targets/u192.rs b/math/fuzz/fuzz_targets/u192.rs index 7d7489df62..79872a05e8 100644 --- a/math/fuzz/fuzz_targets/u192.rs +++ b/math/fuzz/fuzz_targets/u192.rs @@ -3,10 +3,11 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; +use kaspa_math::lehmer::LehmerInvert; use libfuzzer_sys::fuzz_target; use num_bigint::BigUint; use num_traits::Zero; -use utils::{assert_same, consume, try_opt}; +use utils::{assert_same, bigint_mod_inv, consume, try_opt}; // This is important as it's a non power of two. construct_uint!(Uint192, 3); @@ -128,6 +129,23 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse - // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. + // lehmer_invert (the mod_inv arm, recovered) against a num-bigint extended-gcd + // reference. `lehmer_invert` is total, so both value and modulus are fuzzed freely; + // only the trivial-ring modulus < 2 is filtered (the inverse in Z/1Z is convention: + // the oracle says 0, lehmer says None). + { + let ((lib1, native1), (lib2, native2)) = loop { + let (lib1, native1) = try_opt!(generate_ints(&mut data)); + let (lib2, native2) = try_opt!(generate_ints(&mut data)); + if lib2 > 1u64 { + break ((lib1, native1), (lib2, native2)); + } + }; + let lib_inv = lib1.lehmer_invert(lib2); + let native_inv = bigint_mod_inv(native1, native2); + assert_eq!(lib_inv.is_some(), native_inv.is_some(), "lib1: {lib1}, lib2: {lib2}"); + if let Some(lib_inv) = lib_inv { + assert_same!(lib_inv, native_inv.unwrap(), "lib1: {lib1}, lib2: {lib2}"); + } + } }); diff --git a/math/fuzz/fuzz_targets/u256.rs b/math/fuzz/fuzz_targets/u256.rs index 8821a86e79..743d53addd 100644 --- a/math/fuzz/fuzz_targets/u256.rs +++ b/math/fuzz/fuzz_targets/u256.rs @@ -3,10 +3,11 @@ mod utils; use core::ops::{Add, BitAnd, BitOr, BitXor, Div, Mul, Rem}; use kaspa_math::construct_uint; +use kaspa_math::lehmer::LehmerInvert; use libfuzzer_sys::fuzz_target; use num_bigint::BigUint; use num_traits::Zero; -use utils::{assert_same, consume, try_opt}; +use utils::{assert_same, bigint_mod_inv, consume, try_opt}; construct_uint!(Uint256, 4); @@ -127,6 +128,23 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv was removed with the generic `Uint::mod_inverse`; the 3072-bit MuHash inverse - // (the only production case) is fuzzed in `u3072.rs` against a num-bigint reference. + // lehmer_invert (the mod_inv arm, recovered) against a num-bigint extended-gcd + // reference. `lehmer_invert` is total, so both value and modulus are fuzzed freely; + // only the trivial-ring modulus < 2 is filtered (the inverse in Z/1Z is convention: + // the oracle says 0, lehmer says None). + { + let ((lib1, native1), (lib2, native2)) = loop { + let (lib1, native1) = try_opt!(generate_ints(&mut data)); + let (lib2, native2) = try_opt!(generate_ints(&mut data)); + if lib2 > 1u64 { + break ((lib1, native1), (lib2, native2)); + } + }; + let lib_inv = lib1.lehmer_invert(lib2); + let native_inv = bigint_mod_inv(native1, native2); + assert_eq!(lib_inv.is_some(), native_inv.is_some(), "lib1: {lib1}, lib2: {lib2}"); + if let Some(lib_inv) = lib_inv { + assert_same!(lib_inv, native_inv.unwrap(), "lib1: {lib1}, lib2: {lib2}"); + } + } }); diff --git a/math/fuzz/fuzz_targets/u3072.rs b/math/fuzz/fuzz_targets/u3072.rs index 0e13583eec..988edcc175 100644 --- a/math/fuzz/fuzz_targets/u3072.rs +++ b/math/fuzz/fuzz_targets/u3072.rs @@ -39,7 +39,7 @@ fuzz_target!(|data: &[u8]| { // The production inverse under test. let mut out = [0u64; LIMBS]; - let invertible = lehmer::invert(value.0, prime.0, &mut out); + let invertible = lehmer::invert::(value.0, prime.0, &mut out); // The MuHash prime is prime, so every nonzero reduced value is invertible. assert!(invertible, "value={value}"); diff --git a/math/fuzz/fuzz_targets/utils.rs b/math/fuzz/fuzz_targets/utils.rs index 309b9bbe2c..d3efb52a01 100644 --- a/math/fuzz/fuzz_targets/utils.rs +++ b/math/fuzz/fuzz_targets/utils.rs @@ -26,4 +26,22 @@ pub fn consume(data: &mut &[u8]) -> Option<[u8; N]> { } } +/// num-bigint extended-gcd modular-inverse reference: `a^-1 (mod n)` when `gcd(a, n) == 1`. +pub fn bigint_mod_inv(a: num_bigint::BigUint, n: num_bigint::BigUint) -> Option { + use num_bigint::BigInt; + use num_integer::Integer; + use num_traits::Signed; + let a = BigInt::from(a); + let n = BigInt::from(n); + let e_gcd = a.extended_gcd(&n); + // An inverse exists iff gcd(a, n) == 1 + if e_gcd.gcd != 1u64.into() { + None + } else if e_gcd.x.is_negative() { + (e_gcd.x + n).try_into().ok() + } else { + e_gcd.x.try_into().ok() + } +} + pub(crate) use {assert_same, try_opt}; From c910b6a619151a7b9697b7c6967773b8fd21a6bd Mon Sep 17 00:00:00 2001 From: 143672 Date: Fri, 10 Jul 2026 19:22:21 +0400 Subject: [PATCH 11/13] math: ceil_log_2 returns 0 for x=0 via leading_zeros Replace ilog2() (which panics on 0) with leading_zeros() arithmetic, guarded to return 0 when x=0. Eliminates panic case and maintains performance. --- math/src/lib.rs | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/math/src/lib.rs b/math/src/lib.rs index c85b49efae..04b5f68c32 100644 --- a/math/src/lib.rs +++ b/math/src/lib.rs @@ -13,13 +13,10 @@ construct_uint!(Uint320, 5); construct_uint!(Uint3072, 48); /// Returns the ceiling of the base-2 logarithm of `x`, i.e. the smallest `k` such that `2^k >= x`. -/// -/// # Panics -/// Panics if `x` is 0 (the base-2 logarithm of 0 is undefined). +/// If `x` is 0, returns 0. #[inline] pub const fn ceil_log_2(x: u64) -> u64 { - // power of two -> floor; not a power of two -> floor + 1 - x.ilog2() as u64 + (!x.is_power_of_two()) as u64 + (u64::BITS - x.saturating_sub(1).leading_zeros()) as u64 } #[derive(thiserror::Error, Debug)] @@ -204,9 +201,11 @@ mod ceil_log_2_tests { use crate::ceil_log_2; /// Independent reference: the smallest `k` such that `2^k >= x`. Computed in `u128` so it - /// stays correct for `x` near `u64::MAX` (where the answer is 64). + /// stays correct for `x` near `u64::MAX` (where the answer is 64). Returns 0 if `x` is 0. fn oracle(x: u64) -> u64 { - assert!(x != 0); + if x == 0 { + return 0; + } let mut k = 0u64; while (1u128 << k) < x as u128 { k += 1; @@ -216,11 +215,11 @@ mod ceil_log_2_tests { /// A spread of inputs exercising the dense low range, every power-of-2 boundary, and the top. fn sample_inputs() -> impl Iterator { - let dense = 1..=8192u64; + let dense = 0..=8192u64; let boundaries = (0..64u32).flat_map(|k| { let p = 1u64 << k; - // p-1 (clamped away from 0), p, p+1 - [p.saturating_sub(1).max(1), p, p + 1] + // p-1 (clamped away from negative), p, p+1 + [p.saturating_sub(1), p, p + 1] }); dense.chain(boundaries).chain([u64::MAX]) } @@ -228,7 +227,8 @@ mod ceil_log_2_tests { #[test] fn known_values() { for (x, expected) in [ - (1u64, 0u64), + (0u64, 0u64), + (1, 0), (2, 1), (3, 2), (4, 2), @@ -258,16 +258,12 @@ mod ceil_log_2_tests { fn compatibility_matches_malachite() { // Direct equivalence with the malachite `CeilingLogBase2` this replaced. malachite is a // dev-dependency only (the oracle); permanent coverage is `correctness_matches_oracle`. + // Skip x=0 since malachite panics on it; our implementation intentionally returns 0. use malachite_base::num::arithmetic::traits::CeilingLogBase2; for x in sample_inputs() { - assert_eq!(ceil_log_2(x), x.ceiling_log_base_2(), "x={x}"); + if x != 0 { + assert_eq!(ceil_log_2(x), x.ceiling_log_base_2(), "x={x}"); + } } } - - #[test] - #[should_panic] - fn panics_on_zero() { - // Matches malachite's "Cannot take the base-2 logarithm of 0." panic. - let _ = ceil_log_2(0); - } } From 5160f8dd5d9c095d40092c1c4a8fe20a6cedc15a Mon Sep 17 00:00:00 2001 From: Maxim <59533214+biryukovmaxim@users.noreply.github.com> Date: Fri, 10 Jul 2026 19:35:16 +0400 Subject: [PATCH 12/13] apply black box Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com> --- math/benches/bench.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/math/benches/bench.rs b/math/benches/bench.rs index 63c1b2cad4..e5389dcc8d 100644 --- a/math/benches/bench.rs +++ b/math/benches/bench.rs @@ -132,7 +132,9 @@ fn bench_uint3072(c: &mut Criterion) { for &a in &uint3072_one[..uint3072_one.len() / 4] { // generic mod_inverse was removed; inverse the reduced value via lehmer::invert let mut out = [0u64; 48]; - black_box(kaspa_math::lehmer::invert::((a % PRIME).0, PRIME.0, &mut out)); + let invertible = kaspa_math::lehmer::invert::((a % PRIME).0, PRIME.0, &mut out); + black_box(invertible); + black_box(&out); } }); }); From 46b0e2cc5956d9adde44391485eec3bb7c45ce9b Mon Sep 17 00:00:00 2001 From: 143672 Date: Fri, 10 Jul 2026 19:57:29 +0400 Subject: [PATCH 13/13] math: strict-asserts feature; promote lehmer debug_asserts, add invariant checks --- math/Cargo.toml | 3 ++ math/fuzz/Cargo.toml | 1 + math/fuzz/fuzz_targets/u128.rs | 2 +- math/src/lehmer.rs | 95 ++++++++++++++++++++++++++-------- 4 files changed, 77 insertions(+), 24 deletions(-) diff --git a/math/Cargo.toml b/math/Cargo.toml index 38135b53a9..83acf545b7 100644 --- a/math/Cargo.toml +++ b/math/Cargo.toml @@ -45,6 +45,9 @@ malachite-nz.workspace = true # (matrices, divides, apply-limb iterations) into atomics. Off by default; has # zero effect on the production path. lehmer-instrument = [] +# Promote the `lehmer` invariant checks from `debug_assert!` to always-on +# `assert!`, independent of `debug-assertions`. +strict-asserts = [] [[bench]] name = "bench" diff --git a/math/fuzz/Cargo.toml b/math/fuzz/Cargo.toml index 7f21aa537a..93fcc6902c 100644 --- a/math/fuzz/Cargo.toml +++ b/math/fuzz/Cargo.toml @@ -16,6 +16,7 @@ num-integer = "0.1" [dependencies.kaspa-math] path = ".." +features = ["strict-asserts"] # Prevent this from interfering with workspaces [workspace] diff --git a/math/fuzz/fuzz_targets/u128.rs b/math/fuzz/fuzz_targets/u128.rs index 8e4a193e07..ba2c811cf5 100644 --- a/math/fuzz/fuzz_targets/u128.rs +++ b/math/fuzz/fuzz_targets/u128.rs @@ -98,7 +98,7 @@ fuzz_target!(|data: &[u8]| { // as u128 { let (lib, native) = try_opt!(generate_ints(&mut data)); - assert_eq!(lib.as_u128(), native as u128, "native: {native}"); + assert_eq!(lib.as_u128(), native, "native: {native}"); } // as f64 { diff --git a/math/src/lehmer.rs b/math/src/lehmer.rs index 79921a29ab..b8f0108bb9 100644 --- a/math/src/lehmer.rs +++ b/math/src/lehmer.rs @@ -68,6 +68,26 @@ macro_rules! count { }}; } +// Invariant checks: `debug_assert!` by default (free in release), promoted to +// always-on `assert!` under the `strict-asserts` feature so optimized fuzz and +// test builds verify every invariant regardless of `debug-assertions`. +#[cfg(not(feature = "strict-asserts"))] +macro_rules! strict_assert { + ($($arg:tt)*) => { debug_assert!($($arg)*) }; +} +#[cfg(feature = "strict-asserts")] +macro_rules! strict_assert { + ($($arg:tt)*) => { assert!($($arg)*) }; +} +#[cfg(not(feature = "strict-asserts"))] +macro_rules! strict_assert_eq { + ($($arg:tt)*) => { debug_assert_eq!($($arg)*) }; +} +#[cfg(feature = "strict-asserts")] +macro_rules! strict_assert_eq { + ($($arg:tt)*) => { assert_eq!($($arg)*) }; +} + /// Fixed-size little-endian limb buffer, blanket-implemented for every /// `[u64; N]`. The inversion's working storage; keeping the concrete array /// type generic (rather than passing slices around) monomorphizes the limb @@ -203,6 +223,7 @@ pub fn invert(value: U::Limbs, modulus: U::Limbs, out: &mut U::Lim core::hint::cold_path(); return false; } + strict_assert!(cmp_n(value.as_slice(), modulus.as_slice()).is_lt(), "invert requires value < modulus"); // Euclidean state on (x, y) with x >= y, plus the cofactor of `value`. // Invariant: x ≡ -t0·value (mod m) and y ≡ t1·value (mod m) when @@ -296,6 +317,10 @@ pub fn invert(value: U::Limbs, modulus: U::Limbs, out: &mut U::Lim if !swapped && !is_zero(inv.as_slice()) { inv = sub_n(&modulus, &inv); } + strict_assert!( + !is_zero(inv.as_slice()) && cmp_n(inv.as_slice(), modulus.as_slice()).is_lt(), + "invert: result outside (0, modulus)" + ); *out = inv; true } @@ -491,6 +516,12 @@ fn hgcd2(mut ah: u64, mut al: u64, mut bh: u64, mut bl: u64) -> Option<(u64, u64 } } + // Every entry stays below 2^63 (the half-limb refinement's exit bound, which + // the `apply_matrix_*` overflow analyses rely on), and the column updates + // preserve the unit determinant of the elementary step matrices. + strict_assert!(m00 < 1 << 63 && m01 < 1 << 63 && m10 < 1 << 63 && m11 < 1 << 63, "hgcd2: matrix entry reached 2^63"); + strict_assert_eq!((m00 as i128) * (m11 as i128) - (m01 as i128) * (m10 as i128), 1, "hgcd2: determinant is not 1"); + // Remap M to this module's apply convention: // x' = m11·x - m01·y, y' = m00·y - m10·x. Some((m11, m01, m10, m00)) @@ -504,9 +535,9 @@ fn gcd_div(n: u128, d: u128) -> (u64, u128) { count!(DIVIDES += 1); let (n1, n0) = ((n >> 64) as u64, n as u64); let (mut d1, mut d0) = ((d >> 64) as u64, d as u64); - debug_assert!(d1 != 0); + strict_assert!(d1 != 0); let (q, r) = (n1 / d1, n1 % d1); - if q > d1 { + let (q, rem) = if q > d1 { // Non-normalized divisor: when `d1`'s top word is small, the single-word // estimate `q = n1/d1` can exceed `d1`. Rare but genuinely reachable (a // b-step can drop the divisor below 2^32 before an a-step divides by it), @@ -548,7 +579,10 @@ fn gcd_div(n: u128, d: u128) -> (u64, u128) { let (rr0, br) = n0.overflowing_sub(t0); let rr1 = r.wrapping_sub(t1).wrapping_sub(br as u64); (q, ((rr1 as u128) << 64) | rr0 as u128) - } + }; + // Exact reconstruction: cannot wrap for a correct result since q·d + r = n < 2^128. + strict_assert!(rem < d && (q as u128).wrapping_mul(d).wrapping_add(rem) == n, "gcd_div: bad quotient or remainder"); + (q, rem) } /// Multiply-accumulate: `acc + m·w + carry`, returning `(low 64 bits, high 64 @@ -619,8 +653,8 @@ fn apply_matrix_xy(x: &mut A, y: &mut A, len: usize, a: u64, b: u6 } // Non-negative result fitting in `len` limbs: the positive carry-out exactly // cancels the negative borrow-out (both name the same high limb, which is 0). - debug_assert_eq!(carry_ax, borrow_x, "apply_matrix_xy: x went negative or overflowed len limbs"); - debug_assert_eq!(carry_dy, borrow_y, "apply_matrix_xy: y went negative or overflowed len limbs"); + strict_assert_eq!(carry_ax, borrow_x, "apply_matrix_xy: x went negative or overflowed len limbs"); + strict_assert_eq!(carry_dy, borrow_y, "apply_matrix_xy: y went negative or overflowed len limbs"); let nlx = x[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); let nly = y[..len].iter().rposition(|&w| w != 0).map_or(0, |i| i + 1); @@ -667,6 +701,7 @@ fn apply_matrix_t(t0: &mut C, t1: &mut C, t0_len: &mut usize, t1_l c1 = n1 >> 64; } // Carry-out (provably < 2^64); also clears any stale cell at `max_len`. + strict_assert!(c0 >> 64 == 0 && c1 >> 64 == 0, "apply_matrix_t: carry-out exceeded one limb"); t0[max_len] = c0 as u64; t1[max_len] = c1 as u64; @@ -689,7 +724,7 @@ fn t_add_qmul(t0: &mut C, t0_len: &mut usize, q: &[u64], t1: &C, t for (j, &t1j) in t1.iter().enumerate().take(t1_len) { let k = i + j; if k >= cap { - debug_assert_eq!(carry, 0, "t_add_qmul: truncated nonzero limb"); + strict_assert_eq!(carry, 0, "t_add_qmul: truncated nonzero limb"); break; } let prod = qi.wrapping_mul(t1j as u128).wrapping_add(t0[k] as u128).wrapping_add(carry as u128); @@ -703,7 +738,7 @@ fn t_add_qmul(t0: &mut C, t0_len: &mut usize, q: &[u64], t1: &C, t carry = c as u64; k += 1; } - debug_assert!(carry == 0, "t_add_qmul: carry lost off the top"); + strict_assert!(carry == 0, "t_add_qmul: carry lost off the top"); } *t0_len = top_len_t(t0); } @@ -732,7 +767,7 @@ fn add_mul_word(out: &mut C, out_len: &mut usize, w: u64, t: &C, t carry = c as u64; k += 1; } - debug_assert!(carry == 0, "add_mul_word: carry lost off the top"); + strict_assert!(carry == 0, "add_mul_word: carry lost off the top"); } *out_len = top_len_t(out); } @@ -755,7 +790,7 @@ fn div_rem_step(x: &U::Limbs, y: &U::Limbs) -> (U::Limbs, U::Limbs /// instead of a hardware divide. Precomputed-reciprocal division /// (Moller-Granlund, 2011). fn div_n_by_1(x: &A, d: u64) -> (A, u64) { - debug_assert!(d != 0); + strict_assert!(d != 0); let x = x.as_slice(); let n = x.len(); let mut q = A::ZERO; @@ -782,6 +817,7 @@ fn div_n_by_1(x: &A, d: u64) -> (A, u64) { } r >> s }; + strict_assert!(r < d, "div_n_by_1: remainder not below the divisor"); (q, r) } @@ -792,7 +828,7 @@ fn div_n_by_1(x: &A, d: u64) -> (A, u64) { /// call (amortized over all `N` limbs), so it is off the per-limb hot path. #[inline] const fn invert_limb(d: u64) -> u64 { - debug_assert!(d >> 63 == 1, "invert_limb requires a normalized divisor"); + strict_assert!(d >> 63 == 1, "invert_limb requires a normalized divisor"); // v = floor((2^128 - 1) / d) - 2^64 = floor(((!d)·2^64 + !0) / d): subtracting // 2^64·d from the numerator drops the quotient by exactly 2^64, and that // numerator's high word `!d` is `< d` for a normalized `d`, so the quotient @@ -809,7 +845,7 @@ const fn invert_limb(d: u64) -> u64 { /// with `overflow-checks = true`). #[inline(always)] fn div_2by1_preinv(nh: u64, nl: u64, d: u64, di: u64) -> (u64, u64) { - debug_assert!(d >> 63 == 1 && nh < d); + strict_assert!(d >> 63 == 1 && nh < d); // (qh:ql) = nh·di + (nh + 1)·2^64 + nl let prod = (nh as u128).wrapping_mul(di as u128); let (mut qh, ql) = ((prod >> 64) as u64, prod as u64); @@ -827,6 +863,7 @@ fn div_2by1_preinv(nh: u64, nl: u64, d: u64, di: u64) -> (u64, u64) { r = r.wrapping_sub(d); qh = qh.wrapping_add(1); } + strict_assert!(r < d, "div_2by1_preinv: remainder not below the divisor"); (qh, r) } @@ -840,16 +877,16 @@ fn div_2by1_preinv(nh: u64, nl: u64, d: u64, di: u64) -> (u64, u64) { /// fine. #[inline] const fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) { - debug_assert!(d != 0); - debug_assert!(hi < d, "div_2by1 quotient would overflow u64"); + strict_assert!(d != 0); + strict_assert!(hi < d, "div_2by1 quotient would overflow u64"); let s = d.leading_zeros(); - let d = if s == 0 { d } else { d << s }; + let dn = if s == 0 { d } else { d << s }; let un32 = if s == 0 { hi } else { (hi << s) | (lo >> (64 - s)) }; let un10 = if s == 0 { lo } else { lo << s }; - let vn1 = d >> 32; - let vn0 = d & 0xFFFF_FFFF; + let vn1 = dn >> 32; + let vn0 = dn & 0xFFFF_FFFF; let un1 = un10 >> 32; let un0 = un10 & 0xFFFF_FFFF; @@ -863,7 +900,7 @@ const fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) { } } - let un21 = (un32 << 32).wrapping_add(un1).wrapping_sub(q1.wrapping_mul(d)); + let un21 = (un32 << 32).wrapping_add(un1).wrapping_sub(q1.wrapping_mul(dn)); let mut q0 = un21 / vn1; let mut rhat = un21 - q0 * vn1; while q0 >= (1u64 << 32) || q0 * vn0 > (rhat << 32) | un0 { @@ -874,13 +911,19 @@ const fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) { } } - let r = (un21 << 32).wrapping_add(un0).wrapping_sub(q0.wrapping_mul(d)) >> s; - ((q1 << 32) | q0, r) + let r = (un21 << 32).wrapping_add(un0).wrapping_sub(q0.wrapping_mul(dn)) >> s; + let q = (q1 << 32) | q0; + strict_assert!( + r < d && (q as u128).wrapping_mul(d as u128).wrapping_add(r as u128) == ((hi as u128) << 64) | lo as u128, + "div_2by1: bad quotient or remainder" + ); + (q, r) } /// Extended Euclidean GCD of two `u64`s: `(gcd, cx, cy)` with `cx·x + cy·y = g`. /// Called once per inversion, so `i128` intermediates are fine. fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { + let (x0, y0) = (x as i128, y as i128); let (mut a, mut b, mut c, mut d) = (1i128, 0i128, 0i128, 1i128); while y != 0 { let q = (x / y) as i128; @@ -893,6 +936,10 @@ fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { x = y; y = r; } + // Cofactors are bounded by the inputs (so the i64 narrowing is lossless) and + // satisfy the Bezout identity for the returned gcd. + strict_assert!(i64::try_from(a).is_ok() && i64::try_from(b).is_ok(), "gcd_ext_u64: cofactor overflowed i64"); + strict_assert_eq!(a.wrapping_mul(x0).wrapping_add(b.wrapping_mul(y0)), x as i128, "gcd_ext_u64: Bezout identity failed"); (x, a as i64, b as i64) } @@ -901,7 +948,7 @@ fn gcd_ext_u64(mut x: u64, mut y: u64) -> (u64, i64, i64) { /// subtractions suffice; the `div_rem` branch is a defensive fallback. fn reduce_mod(value: &U::Cofactor, m: &U::Limbs) -> U::Limbs { let n = U::Limbs::LEN; - debug_assert!(value.as_slice()[n..].iter().all(|&w| w == 0), "cofactor exceeded the modulus width"); + strict_assert!(value.as_slice()[n..].iter().all(|&w| w == 0), "cofactor exceeded the modulus width"); let mut out = U::Limbs::ZERO; out.as_mut_slice().copy_from_slice(&value.as_slice()[..n]); for _ in 0..8 { @@ -920,7 +967,7 @@ fn reduce_mod(value: &U::Cofactor, m: &U::Limbs) -> U::Limbs { /// `x_len >= 2`. #[inline] fn highest_two_words_normalized(x: &[u64], y: &[u64], x_len: usize) -> (u128, u128) { - debug_assert!(x_len >= 2); + strict_assert!(x_len >= 2); let i = x_len - 1; let lz = x[i].leading_zeros(); let lo_idx_ok = i >= 2; @@ -932,7 +979,9 @@ fn highest_two_words_normalized(x: &[u64], y: &[u64], x_len: usize) -> (u128, u1 }; let x_lo = if lo_idx_ok { x[i - 2] } else { 0 }; let y_lo = if lo_idx_ok { y[i - 2] } else { 0 }; - (combine(x[i], x[i - 1], x_lo), combine(y[i], y[i - 1], y_lo)) + let (x_hi, y_hi) = (combine(x[i], x[i - 1], x_lo), combine(y[i], y[i - 1], y_lo)); + strict_assert!(x_hi >> 127 == 1 && y_hi <= x_hi, "highest_two_words_normalized: prefixes unnormalized or misordered"); + (x_hi, y_hi) } /// A limb buffer holding the single word `w`. @@ -976,7 +1025,7 @@ fn sub_n(a: &A, b: &A) -> A { /// Full-width limb compare of two equal-length buffers. #[inline] fn cmp_n(a: &[u64], b: &[u64]) -> Ordering { - debug_assert_eq!(a.len(), b.len()); + strict_assert_eq!(a.len(), b.len()); cmp_prefix(a, b, a.len()) }