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/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/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" 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/crypto/muhash/src/u3072.rs b/crypto/muhash/src/u3072.rs index 3b8116beb4..fbc40a2db4 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/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 \ diff --git a/math/Cargo.toml b/math/Cargo.toml index 3c40dbe756..83acf545b7 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,10 +25,37 @@ 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: +# 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` +# (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" harness = false +[[bench]] +name = "candidates" +harness = false + [lints] workspace = true diff --git a/math/benches/bench.rs b/math/benches/bench.rs index 5909463b0d..e5389dcc8d 100644 --- a/math/benches/bench.rs +++ b/math/benches/bench.rs @@ -130,7 +130,11 @@ 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]; + let invertible = kaspa_math::lehmer::invert::((a % PRIME).0, PRIME.0, &mut out); + black_box(invertible); + black_box(&out); } }); }); diff --git a/math/benches/candidates.rs b/math/benches/candidates.rs new file mode 100644 index 0000000000..5f723bda2c --- /dev/null +++ b/math/benches/candidates.rs @@ -0,0 +1,271 @@ +//! Benchmarks the two malachite-backed operations this branch replaces, each against its +//! replacement, to confirm the swaps are not regressions: +//! +//! 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, 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`. + +use criterion::measurement::WallTime; +use criterion::{BenchmarkGroup, Criterion, black_box, criterion_group, criterion_main}; +use rand_chacha::{ + ChaCha8Rng, + rand_core::{RngCore, SeedableRng}, +}; + +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; + +// Same value as muhash UINT_PRIME: 2^3072 - 1103717 +const PRIME: Uint3072 = { + let mut max = Uint3072::MAX; + max.0[0] -= 1103717 - 1; + max +}; + +// 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]); + 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::*; + 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] { + 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 + } +} + +// 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"); + 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 + } +} + +// 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::*; +// 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. +// 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() +// } +// } +// +// 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))); + } + }); + }); +} + +// mod_inverse: the malachite incumbent vs the in-repo Lehmer ext-gcd replacement. +fn bench_modinv_3072(c: &mut Criterion) { + let inputs = inputs(); + + // lehmer (the chosen replacement) must agree with the current implementation + 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); + // 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, "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(); + // 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 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 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(); +} + +criterion_group!(benches, bench_modinv_3072, bench_ceil_log2); +criterion_main!(benches); 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 new file mode 100644 index 0000000000..a775eb9a19 --- /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)); +} diff --git a/math/fuzz/Cargo.toml b/math/fuzz/Cargo.toml index 3b4743834c..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] @@ -38,3 +39,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..ba2c811cf5 100644 --- a/math/fuzz/fuzz_targets/u128.rs +++ b/math/fuzz/fuzz_targets/u128.rs @@ -3,9 +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 std::convert::TryInto; -use utils::{consume, try_opt}; +use num_bigint::BigUint; +use utils::{bigint_mod_inv, consume, try_opt}; construct_uint!(Uint128, 2); @@ -15,13 +16,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, @@ -104,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 { @@ -131,42 +125,20 @@ fuzz_target!(|data: &[u8]| { assert_eq!(lib_bit, native_bit, "native: {native}"); } } - // mod_inv + // 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). { - // 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 { + 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 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}"); - } + 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}"); } }); - -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..79872a05e8 100644 --- a/math/fuzz/fuzz_targets/u192.rs +++ b/math/fuzz/fuzz_targets/u192.rs @@ -3,12 +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::{BigInt, BigUint}; -use num_integer::Integer; -use num_traits::{Signed, Zero}; -use std::convert::TryInto; -use utils::{assert_same, consume, try_opt}; +use num_bigint::BigUint; +use num_traits::Zero; +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); @@ -130,36 +129,23 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv + // 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). { - // 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 { + if lib2 > 1u64 { break ((lib1, native1), (lib2, native2)); } }; - let lib_inv = lib1.mod_inverse(lib2); + 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()); + 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}"); } } }); - -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..743d53addd 100644 --- a/math/fuzz/fuzz_targets/u256.rs +++ b/math/fuzz/fuzz_targets/u256.rs @@ -3,12 +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::{BigInt, BigUint}; -use num_integer::Integer; -use num_traits::{Signed, Zero}; -use std::convert::TryInto; -use utils::{assert_same, consume, try_opt}; +use num_bigint::BigUint; +use num_traits::Zero; +use utils::{assert_same, bigint_mod_inv, consume, try_opt}; construct_uint!(Uint256, 4); @@ -129,35 +128,23 @@ fuzz_target!(|data: &[u8]| { } } - // mod_inv + // 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). { - // 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 { + if lib2 > 1u64 { break ((lib1, native1), (lib2, native2)); } }; - let lib_inv = lib1.mod_inverse(lib2); + 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()); + 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}"); } } }); - -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..988edcc175 --- /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/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}; diff --git a/math/src/lehmer.rs b/math/src/lehmer.rs new file mode 100644 index 0000000000..b8f0108bb9 --- /dev/null +++ b/math/src/lehmer.rs @@ -0,0 +1,1398 @@ +//! 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 +//! 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 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 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); + }}; +} + +// 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 +/// 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 +// (`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 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: 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.as_slice()) { + 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 + // `swapped == false`, and the roles flip with each swap. `x`, `y` are owned + // working buffers (they are reduced in place). + 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; + + // 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.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 { + 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.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); + 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.as_slice(), &t1, t1_len); + x = y; + x_len = y_len; + y = r; + 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; + } + } + + // 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.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.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 = 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 = 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.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.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 +} + +/// 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)); + } + } + } + + // 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)) +} + +/// `(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); + strict_assert!(d1 != 0); + let (q, r) = (n1 / d1, n1 % 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), + // 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) + }; + // 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 +/// 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 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 +/// 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 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; + 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). + 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); + (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 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` 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; + 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`. + 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; + + 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 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; + } + 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 >= cap { + 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); + t0[k] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = i + t1_len; + while carry > 0 && k < cap { + let (s, c) = t0[k].overflowing_add(carry); + t0[k] = s; + carry = c as u64; + k += 1; + } + strict_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 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 <= 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); + out[j] = prod as u64; + carry = (prod >> 64) as u64; + } + let mut k = t_len; + while carry > 0 && k < cap { + let (s, c) = out[k].overflowing_add(carry); + out[k] = s; + carry = c as u64; + k += 1; + } + strict_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 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) = 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 +/// 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: &A, d: u64) -> (A, u64) { + strict_assert!(d != 0); + 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(); + 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() { + (qs[i], r) = div_2by1_preinv(r, x[i], d, v); + } + 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)); + (qs[i], r) = div_2by1_preinv(r, cur, d_norm, v); + } + r >> s + }; + strict_assert!(r < d, "div_n_by_1: remainder not below the divisor"); + (q, r) +} + +/// 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 { + 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 + // 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) { + 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); + 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); + } + strict_assert!(r < d, "div_2by1_preinv: remainder not below the divisor"); + (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) { + strict_assert!(d != 0); + strict_assert!(hi < d, "div_2by1 quotient would overflow u64"); + + let s = d.leading_zeros(); + 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 = dn >> 32; + let vn0 = dn & 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(dn)); + 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(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; + 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; + } + // 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) +} + +/// 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; + 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 { + if cmp_n(out.as_slice(), m.as_slice()).is_lt() { + return out; + } + out = sub_n(&out, m); + } + 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 +/// 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], y: &[u64], x_len: usize) -> (u128, u128) { + strict_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 }; + 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`. +#[inline] +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]) -> 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]) -> usize { + x.iter().rposition(|&w| w != 0).map_or(1, |i| i + 1) +} + +#[inline] +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..o.len() { + let (d1, b1) = a[i].overflowing_sub(b[i]); + let (d2, b2) = d1.overflowing_sub(borrow); + o[i] = d2; + borrow = (b1 | b2) as u64; + } + out +} + +/// Full-width limb compare of two equal-length buffers. +#[inline] +fn cmp_n(a: &[u64], b: &[u64]) -> Ordering { + strict_assert_eq!(a.len(), b.len()); + cmp_prefix(a, b, a.len()) +} + +#[inline] +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, + other => return other, + } + } + Ordering::Equal +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::Uint3072; + use rand_chacha::{ + ChaCha8Rng, + rand_core::{RngCore, SeedableRng}, + }; + + const N: usize = 48; + + // 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) + } + + // 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.into_limbs().as_slice()); + let m = Natural::from_limbs_asc(modulus.into_limbs().as_slice()); + x.mod_inverse(m).map(|inv| { + let mut out = U::Limbs::ZERO; + let limbs = inv.into_limbs_asc(); + out.as_mut_slice()[..limbs.len()].copy_from_slice(&limbs); + U::from_limbs(out) + }) + } + + 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, MUHASH_PRIME).unwrap(); + let got = lehmer_inv(v.0, MUHASH_PRIME.0).unwrap(); + assert_eq!(got.as_slice(), expected.0.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(Uint3072(v), MUHASH_PRIME).unwrap().0.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 + } + + /// 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/lib.rs b/math/src/lib.rs index e8c319eea0..04b5f68c32 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; @@ -11,6 +12,13 @@ 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`. +/// If `x` is 0, returns 0. +#[inline] +pub const fn ceil_log_2(x: u64) -> u64 { + (u64::BITS - x.saturating_sub(1).leading_zeros()) as u64 +} + #[derive(thiserror::Error, Debug)] pub enum Error { #[error("{0:?}")] @@ -187,3 +195,75 @@ 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). Returns 0 if `x` is 0. + fn oracle(x: u64) -> u64 { + if x == 0 { + return 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 = 0..=8192u64; + let boundaries = (0..64u32).flat_map(|k| { + let p = 1u64 << k; + // p-1 (clamped away from negative), p, p+1 + [p.saturating_sub(1), p, p + 1] + }); + dense.chain(boundaries).chain([u64::MAX]) + } + + #[test] + fn known_values() { + for (x, expected) in [ + (0u64, 0u64), + (1, 0), + (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. 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() { + if x != 0 { + 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..0f95c14d9f 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, js_sys, kaspa_utils, serde, wasm_bindgen}; // TODO: Add u32 support for optimization on 32 bit machines. @@ -342,23 +342,10 @@ 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; 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 { @@ -432,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)?) @@ -862,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}; @@ -960,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() { @@ -1175,6 +1179,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 +1223,5 @@ mod tests { res } } + */ }