/root/kaspa/kaspa-fuzz/rusty-kaspa/math/src/lehmer.rs

Source
//! 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, <https://github.com/blueshift-gg/solana-secp256k1>, `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<const N: usize> 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<Self>;
}

impl<U: LehmerOps> LehmerInvert for U {
    fn lehmer_invert(self, modulus: Self) -> Option<Self> {
        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::<U>(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<U: LehmerOps>(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::<U::Cofactor>(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::<U>(&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::<U>(&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::<U>(&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<A: LimbArray>(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<C: LimbArray>(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<C: LimbArray>(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<C: LimbArray>(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<U: LehmerOps>(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<A: LimbArray>(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<U: LehmerOps>(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<A: LimbArray>(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: LimbArray>(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::<Uint3072>(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<U: LehmerOps>(value: U, modulus: U) -> Option<U> {
        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<u64> {
        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)));
    }
}

/root/kaspa/kaspa-fuzz/rusty-kaspa/math/src/uint.rs

Source
#[doc(hidden)]
pub use {faster_hex, js_sys, kaspa_utils, serde, wasm_bindgen};

// TODO: Add u32 support for optimization on 32 bit machines.

#[macro_export]
macro_rules! construct_uint {
    ($name:ident, $n_words:literal $(, $derive_trait:ty)*) => {
        /// Little-endian large integer type
        #[derive(Copy, Clone, PartialEq, Eq, Hash, Debug$(, $derive_trait )*)]
        pub struct $name(pub [u64; $n_words]);
        #[allow(unused)]
        impl $name {
            pub const ZERO: Self = $name([0; $n_words]);
            pub const MIN: Self = Self::ZERO;
            pub const MAX: Self = $name([u64::MAX; $n_words]);
            pub const BITS: u32 = $n_words * u64::BITS;
            pub const BYTES: usize = $n_words * size_of::<u64>();
            pub const LIMBS: usize = $n_words;

            #[inline]
            pub fn from_u64(n: u64) -> Self {
                let mut ret = Self::ZERO;
                ret.0[0] = n;
                ret
            }
            #[inline]
            pub fn from_u128(n: u128) -> Self {
                let mut ret = Self::ZERO;
                ret.0[0] = n as u64;
                ret.0[1] = (n >> 64) as u64;
                ret
            }

            #[inline]
            pub fn as_u128(self) -> u128 {
                self.0[0] as u128 | ((self.0[1] as u128) << 64)
            }

            #[inline]
            pub fn as_u64(self) -> u64 {
                self.0[0] as u64
            }

            #[inline(always)]
            pub fn is_zero(self) -> bool {
                self.0.iter().all(|&a| a == 0)
            }

            /// Return the least number of bits needed to represent the number
            #[inline(always)]
            pub fn bits(&self) -> u32 {
                for (i, &word) in self.0.iter().enumerate().rev() {
                    if word != 0 {
                        return u64::BITS * (i as u32 + 1) - word.leading_zeros();
                    }
                }
                0
            }

            #[inline(always)]
            pub fn leading_zeros(&self) -> u32 {
                return Self::BITS - self.bits();
            }

            #[inline]
            pub fn overflowing_shl(self, mut s: u32) -> (Self, bool) {
                let overflows = s >= Self::BITS;
                s %= Self::BITS;
                let mut ret = [0u64; $n_words];
                let left_words = (s / 64) as usize;
                let left_shifts = s % 64;

                for i in left_words..$n_words {
                    ret[i] = self.0[i - left_words] << left_shifts;
                }
                if left_shifts > 0 {
                    let left_over = 64 - left_shifts;
                    for i in left_words + 1..$n_words {
                        ret[i] |= self.0[i - 1 - left_words] >> left_over;
                    }
                }
                (Self(ret), overflows)
            }

            #[inline]
            pub fn wrapping_shl(self, s: u32) -> Self {
                self.overflowing_shl(s).0
            }

            #[inline]
            pub fn overflowing_shr(self, mut s: u32) -> (Self, bool) {
                let overflows = s >= Self::BITS;
                s %= Self::BITS;
                let mut ret = [0u64; Self::LIMBS];
                let left_words = (s / 64) as usize;
                let left_shifts = s % 64;

                for i in left_words..Self::LIMBS {
                    ret[i - left_words] = self.0[i] >> left_shifts;
                }
                if left_shifts > 0 {
                    let left_over = 64 - left_shifts;
                    for i in left_words + 1..Self::LIMBS {
                        ret[i - left_words - 1] |= self.0[i] << left_over;
                    }
                }
                (Self(ret), overflows)
            }

            #[inline]
            pub fn overflowing_add(mut self, other: Self) -> (Self, bool) {
                // Replace with std once stabilized:https://github.com/rust-lang/rust/issues/85532
                #[inline(always)]
                pub const fn carrying_add_u64(lhs: u64, rhs: u64, carry: bool) -> (u64, bool) {
                    let (a, b) = lhs.overflowing_add(rhs);
                    let (c, d) = a.overflowing_add(carry as u64);
                    (c, b != d)
                }
                let mut carry = false;
                let mut carry_out;
                for i in 0..Self::LIMBS {
                    (self.0[i], carry_out) = carrying_add_u64(self.0[i], other.0[i], carry);
                    carry = carry_out;
                }
                (self, carry)
            }

            #[inline]
            pub fn overflowing_add_u64(mut self, other: u64) -> (Self, bool) {
                let mut carry: bool;
                (self.0[0], carry) = self.0[0].overflowing_add(other);
                for i in 1..Self::LIMBS {
                    if !carry {
                        break;
                    }
                    (self.0[i], carry) = self.0[i].overflowing_add(1);
                }
                (self, carry)
            }

            #[inline]
            pub fn overflowing_sub(mut self, other: Self) -> (Self, bool) {
                // Replace with std once stabilized:https://github.com/rust-lang/rust/issues/85532
                #[inline(always)]
                pub const fn borrowing_sub_u64(lhs: u64, rhs: u64, borrow: bool) -> (u64, bool) {
                    let (a, b) = lhs.overflowing_sub(rhs);
                    let (c, d) = a.overflowing_sub(borrow as u64);
                    (c, b != d)
                }

                let mut carry = false;
                let mut carry_out;
                for i in 0..Self::LIMBS {
                    (self.0[i], carry_out) = borrowing_sub_u64(self.0[i], other.0[i], carry);
                    carry = carry_out;
                }
                (self, carry)
            }

            #[inline]
            pub fn saturating_sub(self, other: Self) -> Self {
                let (sum, carry) = self.overflowing_sub(other);
                if carry { Self::ZERO } else { sum }
            }

            #[inline]
            pub fn saturating_add(self, other: Self) -> Self {
                let (sum, carry) = self.overflowing_add(other);
                if carry { Self::MAX } else { sum }
            }

            /// Multiplication by u64
            #[inline]
            pub fn overflowing_mul_u64(self, other: u64) -> (Self, bool) {
                let (this, carry) = self.carrying_mul_u64(other);
                (this, carry != 0)
            }

            #[inline]
            pub fn carrying_mul_u64(mut self, other: u64) -> (Self, u64) {
                let mut carry: u128 = 0;
                for i in 0..Self::LIMBS {
                    // TODO: Use `carrying_mul` when stabilized: https://github.com/rust-lang/rust/issues/85532
                    let n = carry + (other as u128) * (self.0[i] as u128);
                    self.0[i] = n as u64;
                    carry = (n >> 64) & u64::MAX as u128;
                }
                (self, carry as u64)
            }

            #[inline]
            pub fn overflowing_mul(self, other: Self) -> (Self, bool) {
                // We should probably replace this with a Montgomery multiplication algorithm
                let mut result = Self::ZERO;
                let mut carry_out = false;
                for j in 0..Self::LIMBS {
                    let mut carry = 0;
                    let mut i = 0;
                    while i + j < Self::LIMBS {
                        let n = (self.0[i] as u128) * (other.0[j] as u128) + (result.0[i + j] as u128) + (carry as u128);
                        result.0[i + j] = n as u64;
                        carry = (n >> 64) as u64;
                        i += 1;
                    }
                    carry_out |= carry != 0;
                }
                (result, carry_out)
            }
            /// Creates big integer value from a byte slice using
            /// little-endian encoding
            #[inline(always)]
            pub fn from_le_bytes(bytes: [u8; Self::BYTES]) -> Self {
                let mut out = [0u64; Self::LIMBS];
                // This should optimize to basically a transmute.
                out.iter_mut()
                    .zip(bytes.chunks_exact(8))
                    .for_each(|(word, bytes)| *word = u64::from_le_bytes(bytes.try_into().unwrap()));
                Self(out)
            }

            /// Creates big integer value from a byte slice using
            /// big-endian encoding
            #[inline(always)]
            pub fn from_be_bytes(bytes: [u8; Self::BYTES]) -> Self {
                let mut out = [0u64; Self::LIMBS];
                out.iter_mut()
                    .rev()
                    .zip(bytes.chunks_exact(8))
                    .for_each(|(word, bytes)| *word = u64::from_be_bytes(bytes.try_into().unwrap()));
                Self(out)
            }

            /// Convert's the Uint into little endian byte array
            #[inline(always)]
            pub fn to_le_bytes(self) -> [u8; Self::BYTES] {
                let mut out = [0u8; Self::BYTES];
                // This should optimize to basically a transmute.
                out.chunks_exact_mut(8).zip(self.0).for_each(|(bytes, word)| bytes.copy_from_slice(&word.to_le_bytes()));
                out
            }

            /// Convert's the Uint into big endian byte array
            #[inline(always)]
            pub fn to_be_bytes(self) -> [u8; Self::BYTES] {
                let mut out = [0u8; Self::BYTES];
                // This should optimize to basically a transmute.
                out.chunks_exact_mut(8)
                    .zip(self.0.into_iter().rev())
                    .for_each(|(bytes, word)| bytes.copy_from_slice(&word.to_be_bytes()));
                out
            }

            #[inline(always)]
            pub fn to_be_bytes_var(self) -> Vec<u8> {
                let bytes = self.to_be_bytes();
                let start = bytes.iter().copied().position(|b| b != 0).unwrap_or(bytes.len());
                Vec::from(&bytes[start..])
            }

            #[inline]
            pub fn div_rem_u64(mut self, other: u64) -> (Self, u64) {
                let mut rem = 0u64;
                self.0.iter_mut().rev().for_each(|d| {
                    let n = (rem as u128) << 64 | (*d as u128);
                    *d = (n / other as u128) as u64;
                    rem = (n % other as u128) as u64;
                });
                (self, rem)
            }

            #[inline]
            pub fn as_f64(&self) -> f64 {
                // Reference: https://blog.m-ou.se/floats/
                // Step 1: Get leading zeroes
                let leading_zeroes =  self.leading_zeros();
                // Step 2: Align the bits to the left, so the highest bit will be 1.
                let left_aligned = self.wrapping_shl(leading_zeroes);
                // Step 3: Take the highest 53 bits as the mantissa (equivalent to shifting by (Self::BITS - 53))
                let mut mantissa = left_aligned.0[Self::LIMBS - 1] >> 11;
                // Step 4: Get the dropped bits, which are the bits that are not part of the mantissa
                // The dropped bits are left_aligned << 53 (everything except the highest 53 bits).
                // Unlike the blog here we split the highest bit and the rest of the bits into 2 variables.
                // We first take the highest 11 bits that were dropped.
                let highest_dropped_bits = left_aligned.0[Self::LIMBS - 1] << 53;
                let highest_dropped_bit = highest_dropped_bits >> 63 != 0;
                // Now we OR together the rest of the bits.
                let mut rest_dropped_bits = highest_dropped_bits << 1; // Remove the highest.
                for &word in &left_aligned.0[..Self::LIMBS - 1] {
                    rest_dropped_bits |= word;
                }
                // This is true if the dropped bits are higher than half the int.
                let higher_than_half = highest_dropped_bit & (rest_dropped_bits != 0);
                let exactly_half_but_mantissa_odd = highest_dropped_bit & (rest_dropped_bits == 0) & (mantissa & 1 == 1);
                // Step 5: if the dropped bits are higher than half the int, we add 1 to the mantissa.
                // If the dropped bits are exactly half the int, we add 1 to the mantissa only if the mantissa is odd. (IEEE-754)
                mantissa += (higher_than_half | exactly_half_but_mantissa_odd) as u64;
                // Step 6: Calculate the exponent
                // If self is 0, exponent should be 0 (special meaning) and mantissa will end up 0 too
                // Otherwise, (Self::BITS - 1 - leading_zeros) + 1022 so it simplifies to Self::BITS + 1021 - leading_zeroes
                // 1023 and 1022 are the cutoffs for the exponent having the msb next to the decimal point
                let exponent = if self.is_zero() { 0 } else { u64::from(Self::BITS) + 1021 - u64::from(leading_zeroes) };
                // Step 7: sign bit is always 0, exponent is shifted into place
                // Use addition instead of bitwise OR to saturate the exponent if mantissa overflows
                f64::from_bits((exponent << 52) + mantissa)
            }

            // divmod like operation, returns (quotient, remainder)
            #[inline]
            pub fn div_rem(self, other: Self) -> (Self, Self) {
                let mut sub_copy = self;
                let mut shift_copy = other;
                let mut ret = [0u64; Self::LIMBS];

                let my_bits = self.bits();
                let your_bits = other.bits();

                // Check for division by 0
                assert_ne!(your_bits, 0, "attempted to divide {} by zero", self);

                // Early return in case we are dividing by a larger number than us
                if my_bits < your_bits {
                    return (Self(ret), sub_copy);
                }

                // Bitwise long division
                let mut shift = my_bits - your_bits;
                shift_copy = shift_copy << shift;
                loop {
                    if sub_copy >= shift_copy {
                        let (shift_index, shift_val) = ((shift / 64) as usize, shift % 64);
                        ret[shift_index] |= 1 << shift_val;
                        sub_copy = sub_copy - shift_copy;
                    }
                    shift_copy = shift_copy >> 1;
                    if shift == 0 {
                        break;
                    }
                    shift -= 1;
                }

                (Self(ret), sub_copy)
            }

            // 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<Item = bool> + core::iter::FusedIterator {
                struct BinaryIterator {
                    array: [u64; $n_words],
                    bit: usize,
                }

                impl Iterator for BinaryIterator {
                    type Item = bool;

                    #[inline]
                    fn next(&mut self) -> Option<Self::Item> {
                        if self.bit >= 64 * $n_words {
                            return None;
                        }
                        let (word, subbit) = (self.bit / 64, self.bit % 64);
                        let current_bit = self.array[$n_words - word - 1] & (1 << 64 - subbit - 1);
                        self.bit += 1;
                        Some(current_bit != 0)
                    }

                    #[inline]
                    fn nth(&mut self, n: usize) -> Option<Self::Item> {
                        match self.bit.checked_add(n) {
                            Some(bit) => {
                                self.bit = bit;
                                self.next()
                            }
                            None => {
                                self.bit = usize::MAX;
                                None
                            }
                        }
                    }
                    #[inline]
                    fn size_hint(&self) -> (usize, Option<usize>) {
                        let remaining_bits = $n_words * (u64::BITS as usize) - self.bit;
                        (remaining_bits, Some(remaining_bits))
                    }
                }
                impl ExactSizeIterator for BinaryIterator {}
                impl core::iter::FusedIterator for BinaryIterator {}

                BinaryIterator { array: self.0, bit: 0 }
            }

            /// Converts a Self::BYTES*2 hex string interpreted as big endian, into a Uint
            #[inline]
            pub fn from_hex(hex: &str) -> Result<Self, $crate::uint::faster_hex::Error> {
                if hex.len() > Self::BYTES * 2 {
                    return Err($crate::uint::faster_hex::Error::InvalidLength(hex.len()));
                }
                let mut out = [0u8; Self::BYTES];
                let mut input = [b'0'; Self::BYTES * 2];
                let start = input.len() - hex.len();
                input[start..].copy_from_slice(hex.as_bytes());
                $crate::uint::faster_hex::hex_decode(&input, &mut out)?;
                Ok(Self::from_be_bytes(out))
            }

            #[inline]
            pub fn from_be_bytes_var(bytes: &[u8]) -> Result<Self, $crate::uint::TryFromSliceError> {
                if bytes.len() > Self::BYTES {
                    return Err($crate::uint::TryFromSliceError);
                }
                let mut out = [0u8; Self::BYTES];
                let start = Self::BYTES - bytes.len();
                out[start..].copy_from_slice(bytes);
                Ok(Self::from_be_bytes(out))
            }

            #[inline]
            pub fn as_bigint(&self) -> Result<$crate::uint::js_sys::BigInt, $crate::Error> {
                self.try_into()
            }

            #[inline]
            pub fn to_bigint(self) -> Result<$crate::uint::js_sys::BigInt, $crate::Error> {
                self.try_into()
            }

        }

        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 $crate::uint::kaspa_utils::hex::ToHex for $name {
            fn to_hex(&self) -> String {
                self.to_be_bytes().as_slice().to_hex()
            }
        }

        impl $crate::uint::kaspa_utils::hex::ToHex for &$name {
            fn to_hex(&self) -> String {
                self.to_be_bytes().as_slice().to_hex()
            }
        }

        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)?)
            }
        }

        impl PartialEq<u64> for $name {
            #[inline]
            fn eq(&self, other: &u64) -> bool {
                let bigger = self.0[1..].iter().any(|&x| x != 0);
                !bigger && self.0[0] == *other
            }
        }
        impl PartialOrd<u64> for $name {
            #[inline]
            fn partial_cmp(&self, other: &u64) -> Option<core::cmp::Ordering> {
                let bigger = self.0[1..].iter().any(|&x| x != 0);
                if bigger {
                    Some(core::cmp::Ordering::Greater)
                } else {
                    self.0[0].partial_cmp(other)
                }
            }
        }

        impl PartialEq<u128> for $name {
            #[inline]
            fn eq(&self, other: &u128) -> bool {
                let bigger = self.0[2..].iter().any(|&x| x != 0);
                !bigger && self.0[0] == (*other as u64) && self.0[1] == ((*other >> 64) as u64)
            }
        }
        impl PartialOrd<u128> for $name {
            #[inline]
            fn partial_cmp(&self, other: &u128) -> Option<core::cmp::Ordering> {
                let bigger = self.0[2..].iter().any(|&x| x != 0);
                if bigger {
                    Some(core::cmp::Ordering::Greater)
                } else {
                    self.as_u128().partial_cmp(other)
                }
            }
        }

        impl PartialOrd for $name {
            #[inline]
            fn partial_cmp(&self, other: &$name) -> Option<core::cmp::Ordering> {
                Some(self.cmp(&other))
            }
        }

        impl Ord for $name {
            #[inline]
            fn cmp(&self, other: &$name) -> core::cmp::Ordering {
                // We need to manually implement ordering because we use little-endian
                // and the auto derive is a lexicographic ordering(i.e. memcmp)
                // which with numbers is equivalent to big-endian
                Iterator::cmp(self.0.iter().rev(), other.0.iter().rev())
            }
        }

        impl core::ops::Add<$name> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn add(self, other: $name) -> $name {
                let (sum, carry) = self.overflowing_add(other);
                debug_assert!(!carry, "attempt to add with overflow"); // Check in debug that it didn't overflow
                sum
            }
        }

        impl core::ops::Add<u64> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn add(self, other: u64) -> $name {
                let (sum, carry) = self.overflowing_add_u64(other);
                debug_assert!(!carry, "attempt to add with overflow"); // Check in debug that it didn't overflow
                sum
            }
        }

        impl core::ops::Sub<$name> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn sub(self, other: $name) -> $name {
                let (sum, carry) = self.overflowing_sub(other);
                debug_assert!(!carry, "attempt to subtract with overflow"); // Check in debug that it didn't overflow
                sum
            }
        }

        impl core::ops::Mul<$name> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn mul(self, other: $name) -> $name {
                let (product, carry) = self.overflowing_mul(other);
                debug_assert!(!carry, "attempt to multiply with overflow"); // Check in debug that it didn't overflow
                product
            }
        }

        impl core::ops::Mul<u64> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn mul(self, other: u64) -> $name {
                let (product, carry) = self.overflowing_mul_u64(other);
                debug_assert!(!carry, "attempt to multiply with overflow"); // Check in debug that it didn't overflow
                product
            }
        }

        impl core::ops::Div<$name> for $name {
            type Output = $name;

            #[inline]
            fn div(self, other: $name) -> $name {
                self.div_rem(other).0
            }
        }

        impl core::ops::Rem<$name> for $name {
            type Output = $name;

            #[inline]
            fn rem(self, other: $name) -> $name {
                self.div_rem(other).1
            }
        }

        impl core::ops::Div<u64> for $name {
            type Output = $name;

            #[inline]
            fn div(self, other: u64) -> $name {
                self.div_rem_u64(other).0
            }
        }

        impl core::ops::Rem<u64> for $name {
            type Output = u64;

            fn rem(self, other: u64) -> u64 {
                self.div_rem_u64(other).1
            }
        }

        impl core::ops::BitAnd<$name> for $name {
            type Output = $name;

            #[inline]
            fn bitand(mut self, other: $name) -> $name {
                self.0.iter_mut().zip(other.0.iter()).for_each(|(a, b)| *a &= *b);
                self
            }
        }

        impl core::ops::BitXor<$name> for $name {
            type Output = $name;

            #[inline]
            fn bitxor(mut self, other: $name) -> $name {
                self.0.iter_mut().zip(other.0.iter()).for_each(|(a, b)| *a ^= *b);
                self
            }
        }

        impl core::ops::BitOr<$name> for $name {
            type Output = $name;

            #[inline]
            fn bitor(mut self, other: $name) -> $name {
                self.0.iter_mut().zip(other.0.iter()).for_each(|(a, b)| *a |= *b);
                self
            }
        }

        impl core::ops::Not for $name {
            type Output = $name;

            #[inline]
            fn not(mut self) -> $name {
                self.0.iter_mut().for_each(|a| *a = !*a);
                self
            }
        }

        impl core::ops::Shl<u32> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn shl(self, shift: u32) -> $name {
                let (res, carry) = self.overflowing_shl(shift);
                debug_assert!(!carry, "attempt to shift left with overflow"); // Check in debug that it didn't overflow
                res
            }
        }

        impl core::ops::Shr<u32> for $name {
            type Output = $name;

            #[inline]
            #[track_caller]
            fn shr(self, shift: u32) -> $name {
                let (res, carry) = self.overflowing_shr(shift);
                debug_assert!(!carry, "attempt to shift left with overflow"); // Check in debug that it didn't overflow
                res
            }
        }

        impl core::iter::Sum for $name {
            #[inline]
            #[track_caller]
            fn sum<I: Iterator<Item = Self>>(mut iter: I) -> Self {
                let first = iter.next().unwrap_or_else(|| Self::ZERO);
                iter.fold(first, |a, b| a + b)
            }
        }

        impl core::iter::Product for $name {
            #[inline]
            #[track_caller]
            fn product<I: Iterator<Item = Self>>(mut iter: I) -> Self {
                let first = iter.next().unwrap_or_else(|| Self::from_u64(1));
                iter.fold(first, |a, b| a * b)
            }
        }

        impl<'a> core::iter::Sum<&'a $name> for $name {
            #[inline]
            #[track_caller]
            fn sum<I: Iterator<Item = &'a Self>>(mut iter: I) -> Self {
                let first = iter.next().copied().unwrap_or_else(|| Self::ZERO);
                iter.fold(first, |a, &b| a + b)
            }
        }

        impl<'a> core::iter::Product<&'a $name> for $name {
            #[inline]
            #[track_caller]
            fn product<I: Iterator<Item = &'a Self>>(mut iter: I) -> Self {
                let first = iter.next().copied().unwrap_or_else(|| Self::from_u64(1));
                iter.fold(first, |a, &b| a * b)
            }
        }

        impl Default for $name {
            #[inline]
            fn default() -> Self {
                Self::ZERO
            }
        }

        impl From<u64> for $name {
            #[inline]
            fn from(x: u64) -> Self {
                Self::from_u64(x)
            }
        }

        impl core::convert::TryFrom<$name> for u128 {
            type Error = $crate::uint::TryFromIntError;

            #[inline]
            fn try_from(value: $name) -> Result<Self, Self::Error> {
                if value.0[2..].iter().any(|&x| x != 0) {
                    Err($crate::uint::TryFromIntError)
                } else {
                    Ok(value.as_u128())
                }
            }
        }

        impl core::fmt::LowerHex for $name {
            #[inline]
            fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
                let mut hex = [0u8; Self::BYTES * 2];
                let bytes = self.to_be_bytes();
                $crate::uint::faster_hex::hex_encode(&bytes, &mut hex).expect("The output is exactly twice the size of the input");
                let first_non_zero = hex.iter().position(|&x| x != b'0').unwrap_or(hex.len() - 1);
                // The string is hex encoded so must be valid UTF8.
                let str = unsafe { core::str::from_utf8_unchecked(&hex[first_non_zero..]) };
                f.pad_integral(true, "0x", str)
            }
        }

        // Based on https://github.com/rust-lang/rust/blob/2e44c17c12cec45b6a682b1e53a04ac5b5fcc9d2/library/core/src/fmt/num.rs#L209
        impl core::fmt::Display for $name {
            #[inline]
            fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
                // 2 digit decimal look up table
                static DEC_DIGITS_LUT: &[u8; 200] = b"0001020304050607080910111213141516171819\
            2021222324252627282930313233343536373839\
            4041424344454647484950515253545556575859\
            6061626364656667686970717273747576777879\
            8081828384858687888990919293949596979899";

                let mut buf = [0u8; $name::LIMBS * 20]; // 2**64-1 takes 20 digits to represent.
                let mut n = *self;
                let mut curr = buf.len();

                // eagerly decode 4 characters at a time
                const STEP: u64 = 10_000;
                while n >= STEP {
                    let rem: u64;
                    (n, rem) = n.div_rem_u64(STEP);
                    let rem = rem as usize;
                    let d1 = (rem / 100) << 1;
                    let d2 = (rem % 100) << 1;
                    curr -= 4;

                    buf[curr] = DEC_DIGITS_LUT[d1];
                    buf[curr + 1] = DEC_DIGITS_LUT[d1 + 1];
                    buf[curr + 2] = DEC_DIGITS_LUT[d2];
                    buf[curr + 3] = DEC_DIGITS_LUT[d2 + 1];
                }
                // if we reach here numbers are <= 9999, so at most 4 chars long
                let mut n = n.as_u64() as usize; // possibly reduce 64bit math

                // decode 2 more chars, if > 2 chars
                if n >= 100 {
                    let d1 = (n % 100) << 1;
                    n /= 100;
                    curr -= 2;
                    buf[curr] = DEC_DIGITS_LUT[d1 as usize];
                    buf[curr + 1] = DEC_DIGITS_LUT[d1 + 1 as usize];
                }

                // decode last 1 or 2 chars
                if n < 10 {
                    curr -= 1;
                    buf[curr] = (n as u8) + b'0'
                } else {
                    let d1 = n << 1;
                    curr -= 2;
                    buf[curr] = DEC_DIGITS_LUT[d1];
                    buf[curr + 1] = DEC_DIGITS_LUT[d1 + 1];
                }

                // SAFETY: everything up to `curr` is valid UTF8 because `DEC_DIGITS_LUT` is.
                let buf_str = unsafe { std::str::from_utf8_unchecked(&buf[curr..]) };
                f.pad_integral(true, "", buf_str)
            }
        }

        impl core::fmt::Binary for $name {
            #[inline]
            fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
                const BIN_LEN: usize = $name::BITS as usize;
                let mut buf = [0u8; BIN_LEN];
                let mut first_one = BIN_LEN - 1;
                for (index, (bit, char)) in self.iter_be_bits().zip(buf.iter_mut()).enumerate() {
                    *char = bit as u8 + b'0';
                    if first_one == BIN_LEN - 1 && bit {
                        first_one = index;
                    }
                }
                // We only wrote '0' and '1' so this is always valid UTF-8
                let buf_str = unsafe { std::str::from_utf8_unchecked(&buf[first_one..]) };
                f.pad_integral(true, "0b", buf_str)
            }
        }

        // We can't derive because the array might be bigger than 32,
        // so we just implement it the same as arrays.
        impl $crate::uint::serde::Serialize for $name {
            #[inline]
            fn serialize<S: $crate::uint::serde::Serializer>(&self, serializer: S) -> Result<S::Ok, S::Error> {
                if serializer.is_human_readable() {
                    let mut hex = [0u8; Self::BYTES * 2];
                    let bytes = self.to_be_bytes();
                    $crate::uint::faster_hex::hex_encode(&bytes, &mut hex).expect("The output is exactly twice the size of the input");
                    let hex_str = unsafe { std::str::from_utf8_unchecked(&hex) };
                    serializer.serialize_str(hex_str)
                } else {
                    use $crate::uint::serde::ser::SerializeTuple;
                    let mut seq = serializer.serialize_tuple(Self::LIMBS)?;
                    for limb in &self.0 {
                        seq.serialize_element(limb)?;
                    }
                    seq.end()
                }
            }
        }

        impl<'de> $crate::uint::serde::Deserialize<'de> for $name {
            #[inline]
            fn deserialize<D: $crate::uint::serde::Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error> {
                if deserializer.is_human_readable() {
                    let hex = <std::string::String as $crate::uint::serde::Deserialize>::deserialize(deserializer)?;
                    Ok(Self::from_hex(&hex).map_err($crate::uint::serde::de::Error::custom)?)
                } else {
                    use core::{fmt, marker::PhantomData};
                    use $crate::uint::serde::de::{Error, SeqAccess, Visitor};
                    struct EmptyVisitor(PhantomData<$name>);
                    impl<'de> Visitor<'de> for EmptyVisitor {
                        type Value = $name;
                        #[inline]

                        fn expecting(&self, formatter: &mut fmt::Formatter) -> fmt::Result {
                            formatter.write_str(concat!("an integer with ", $n_words, " limbs"))
                        }

                        #[inline]
                        fn visit_seq<A: SeqAccess<'de>>(self, mut seq: A) -> Result<Self::Value, A::Error> {
                            let mut ret = $name::ZERO;
                            for (i, limb) in ret.0.iter_mut().enumerate() {
                                *limb = seq.next_element()?.ok_or_else(|| Error::invalid_length(i, &self))?;
                            }
                            Ok(ret)
                        }
                    }
                    deserializer.deserialize_tuple(Self::LIMBS, EmptyVisitor(PhantomData))
                }
            }

            #[inline]
            fn deserialize_in_place<D: $crate::uint::serde::Deserializer<'de>>(
                deserializer: D,
                place: &mut Self,
            ) -> Result<(), D::Error> {
                if deserializer.is_human_readable() {

                    use core::fmt;
                    use $crate::uint::serde::de::{Error, Visitor};
                    struct InPlaceVisitor<'a>(&'a mut $name);

                    impl<'de, 'a> Visitor<'de> for InPlaceVisitor<'a> {
                        type Value = ();
                        #[inline]
                        fn expecting(&self, formatter: &mut fmt::Formatter) -> fmt::Result {
                            formatter.write_str("a hex string")
                        }
                        #[inline]
                        fn visit_str<E>(self, hex: &str) -> Result<Self::Value, E>
                        where
                            E: Error,
                        {
                            if hex.len() > $name::BYTES * 2 {
                                return Err($crate::uint::serde::de::Error::custom("invalid hex string length"));
                            }
                            let mut bytes = [0u8; $name::BYTES];
                            let mut input = [b'0'; $name::BYTES * 2];
                            let start = input.len() - hex.len();
                            input[start..].copy_from_slice(hex.as_bytes());
                            $crate::uint::faster_hex::hex_decode(&input, &mut bytes).map_err(Error::custom)?;

                            self.0.0.iter_mut()
                                .rev()
                                .zip(bytes.chunks_exact(8))
                                .for_each(|(word, bytes)| *word = u64::from_be_bytes(bytes.try_into().unwrap()));

                            Ok(())
                        }
                    }
                    deserializer.deserialize_str(InPlaceVisitor(place))

                } else {
                    use core::fmt;
                    use $crate::uint::serde::de::{Error, SeqAccess, Visitor};
                    struct InPlaceVisitor<'a>(&'a mut $name);

                    impl<'de, 'a> Visitor<'de> for InPlaceVisitor<'a> {
                        type Value = ();
                        #[inline]
                        fn expecting(&self, formatter: &mut fmt::Formatter) -> fmt::Result {
                            formatter.write_str(concat!("an integer with ", $n_words, " limbs"))
                        }
                        #[inline]
                        fn visit_seq<A: SeqAccess<'de>>(self, mut seq: A) -> Result<Self::Value, A::Error> {
                            for (idx, dest) in self.0 .0[..].iter_mut().enumerate() {
                                match seq.next_element()? {
                                    Some(elem) => *dest = elem,
                                    None => {
                                        return Err(Error::invalid_length(idx, &self));
                                    }
                                }
                            }
                            Ok(())
                        }
                    }
                    deserializer.deserialize_tuple(Self::LIMBS, InPlaceVisitor(place))
                }
            }

        }

        impl TryFrom<&$name> for $crate::uint::js_sys::BigInt {
            type Error = $crate::Error;
            #[inline]
            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 $crate::uint::js_sys::BigInt {
            type Error = $crate::Error;
            #[inline]
            fn try_from(value: $name) -> Result<$crate::uint::js_sys::BigInt, Self::Error> {
                use $crate::wasm::*;
                BigInt::try_from(&value)
            }
        }

        impl TryFrom<$crate::uint::wasm_bindgen::JsValue> for $name {
            type Error = $crate::Error;
            fn try_from(js_value: $crate::uint::wasm_bindgen::JsValue) -> Result<Self, Self::Error> {
                use $crate::wasm::*;

                if js_value.is_string() || js_value.is_array() {
                    let bytes = js_value.try_as_vec_u8()?;
                    Ok(Self::from_be_bytes_var(&bytes)?)
                } else if js_value.is_bigint() {
                    let v: &BigInt = js_value.dyn_ref().unwrap();
                    let hex = String::from(v.to_string(16)?);
                    Ok(Self::from_hex(hex.as_str())?)
                } else {
                    return Err(Self::Error::NotCompatible);
                }
            }
        }

    };
}

/// The error type returned when a checked integral type conversion fails.
#[derive(Debug, Copy, Clone, PartialEq, Eq)]
pub struct TryFromIntError;

impl std::error::Error for TryFromIntError {}

impl core::fmt::Display for TryFromIntError {
    fn fmt(&self, fmt: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        "out of range integral type conversion attempted".fmt(fmt)
    }
}

impl From<core::convert::Infallible> for TryFromIntError {
    fn from(x: core::convert::Infallible) -> TryFromIntError {
        match x {}
    }
}

/// The error type returned when a slice conversion fails.
#[derive(Debug, Copy, Clone, PartialEq, Eq)]
pub struct TryFromSliceError;

impl std::error::Error for TryFromSliceError {}

impl core::fmt::Display for TryFromSliceError {
    fn fmt(&self, fmt: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        "conversion attempted from a slice too large".fmt(fmt)
    }
}

impl From<core::convert::Infallible> for TryFromSliceError {
    fn from(x: core::convert::Infallible) -> TryFromSliceError {
        match x {}
    }
}

#[cfg(test)]
mod tests {
    use rand_chacha::{
        ChaCha8Rng,
        rand_core::{RngCore, SeedableRng},
    };
    use std::fmt::Write;
    construct_uint!(Uint128, 2);

    #[test]
    fn test_u128() {
        use core::fmt::Arguments;
        let mut fmt_buf = String::with_capacity(256);
        let mut fmt_buf2 = String::with_capacity(256);
        let mut assert_equal_args = |arg1: Arguments, arg2: Arguments| {
            fmt_buf.clear();
            fmt_buf2.clear();
            fmt_buf.write_fmt(arg1).unwrap();
            fmt_buf2.write_fmt(arg2).unwrap();
            assert_eq!(fmt_buf, fmt_buf2);
        };
        let mut assert_equal = |a: Uint128, b: u128, check_fmt: bool| {
            assert_eq!(a, b);
            assert_eq!(a.to_le_bytes(), b.to_le_bytes());
            if !check_fmt {
                return;
            }

            assert_equal_args(format_args!("{a:}"), format_args!("{b:}"));
            assert_equal_args(format_args!("{a:b}"), format_args!("{b:b}")); // Test Binary
            assert_equal_args(format_args!("{a:#b}"), format_args!("{b:#b}")); // Test Binary with prefix
            assert_equal_args(format_args!("{a:0128b}"), format_args!("{b:0128b}")); // Test binary with length
            assert_equal_args(format_args!("{a:x}"), format_args!("{b:x}")); // Test LowerHex
            assert_equal_args(format_args!("{a:#x}"), format_args!("{b:#x}")); // Test LowerHex with prefix
            // Test LowerHex with padding
            assert_equal_args(format_args!("{a:032x}"), format_args!("{b:032x}"));
        };
        let mut rng = ChaCha8Rng::from_seed([0; 32]);
        let mut buf = [0u8; 16];
        let mut str_buf = String::with_capacity(32);
        for i in 0..80_000 {
            // Checking all the fmt's is quite expensive.
            let check_fmt = i % 8 == 1;
            rng.fill_bytes(&mut buf);
            let mine = Uint128::from_le_bytes(buf);
            let default = u128::from_le_bytes(buf);
            rng.fill_bytes(&mut buf);
            let mine2 = Uint128::from_le_bytes(buf);
            let default2 = u128::from_le_bytes(buf);
            assert_equal(mine, default, check_fmt);
            assert_equal(mine2, default2, check_fmt);

            let mine = mine.overflowing_add(mine2).0.overflowing_mul(mine2).0;
            let default = default.overflowing_add(default2).0.overflowing_mul(default2).0;
            assert_equal(mine, default, check_fmt);
            let shift = rng.next_u32() % 4096;
            {
                let mine_overflow_shl = mine.overflowing_shl(shift);
                let default_overflow_shl = default.overflowing_shl(shift);
                assert_equal(mine_overflow_shl.0, default_overflow_shl.0, check_fmt);
                assert_eq!(mine_overflow_shl.1, default_overflow_shl.1);
            }
            {
                let mine_overflow_shr = mine.overflowing_shl(shift);
                let default_overflow_shr = default.overflowing_shl(shift);
                assert_equal(mine_overflow_shr.0, default_overflow_shr.0, check_fmt);
                assert_eq!(mine_overflow_shr.1, default_overflow_shr.1);
            }
            {
                let mine_divrem = mine.div_rem(mine2);
                let default_divrem = (default / default2, default % default2);
                assert_equal(mine_divrem.0, default_divrem.0, check_fmt);
                assert_equal(mine_divrem.1, default_divrem.1, check_fmt);
            }
            // Test conversion to f64
            {
                let mine_f64 = mine.as_f64();
                let default_f64 = default as f64;
                assert_eq!(mine_f64, default_f64);
            }
            // Test fast u64 division.
            {
                let rand_u64 = rng.next_u64();
                let mine_divrem = mine.div_rem_u64(rand_u64);
                let default_divrem = (default / u128::from(rand_u64), default % u128::from(rand_u64));
                assert_equal(mine_divrem.0, default_divrem.0, check_fmt);
                assert_eq!(mine_divrem.1, u64::try_from(default_divrem.1).unwrap());
            }
            // Test fast u64 multiplication
            {
                let rand_u64 = rng.next_u64();
                let mine_mult = mine.overflowing_mul_u64(rand_u64);
                let default_mult = default.overflowing_mul(rand_u64 as u128);
                assert_equal(mine_mult.0, default_mult.0, check_fmt);
                assert_eq!(mine_mult.1, default_mult.1);
            }
            // Test fast u64 addition
            {
                let rand_u64 = rng.next_u64();
                let mine_add = mine.overflowing_add_u64(rand_u64);
                let default_add = default.overflowing_add(rand_u64 as u128);
                assert_equal(mine_add.0, default_add.0, check_fmt);
                assert_eq!(mine_add.1, default_add.1);
            }
            // Roundtrip Little-Endian bytes conversion
            {
                let mine_le = mine.to_le_bytes();
                let default_le = default.to_le_bytes();
                assert_eq!(mine_le, default_le);
                assert_eq!(mine, Uint128::from_le_bytes(mine_le));
            }
            // Roundtrip Big-Endian bytes conversion
            {
                let mine_le = mine.to_be_bytes();
                let default_le = default.to_be_bytes();
                assert_eq!(mine_le, default_le);
                assert_eq!(mine, Uint128::from_be_bytes(mine_le));
            }
            // Roundtrip hex
            if check_fmt {
                str_buf.clear();
                str_buf.write_fmt(format_args!("{mine:032x}")).unwrap();
                assert_eq!(mine, Uint128::from_hex(&str_buf).unwrap());
            }
        }
    }

    #[test]
    fn test_saturating_ops() {
        let u1 = Uint128::from_u128(u128::MAX);
        let u2 = Uint128::from_u64(u64::MAX);
        // Sub
        assert_eq!(u1.saturating_sub(u2), Uint128::from_u128(u128::MAX - u64::MAX as u128));
        assert_eq!(u1.saturating_sub(u2).as_u128(), u128::MAX - u64::MAX as u128);
        assert_eq!(u2.saturating_sub(u1), Uint128::ZERO);
        // Add
        assert_eq!(u1.saturating_add(Uint128::from_u64(1)), Uint128::MAX);
        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;
        let mut rng = ChaCha8Rng::from_seed([0; 32]);
        let mut buf = [0u8; 16];
        for _ in 0..50_000 {
            rng.fill_bytes(&mut buf);
            let uint1 = Uint128::from_le_bytes(buf);
            rng.fill_bytes(&mut buf);
            let uint2 = Uint128::from_le_bytes(buf);
            let (bigger, smaller) = match uint1.cmp(&uint2) {
                Ordering::Greater => (uint1, uint2),
                Ordering::Less => (uint2, uint1),
                Ordering::Equal => continue,
            };
            let inv = smaller.mod_inverse(bigger);
            if let Some(inv) = inv {
                assert_eq!(prod_bin(inv, smaller, bigger), 1u64);
            }
        }

        fn sum(x: Uint128, y: Uint128, m: Uint128) -> Uint128 {
            let res = x.overflowing_add(y).0;
            if res < x || res >= m { res.overflowing_sub(m).0 } else { res }
        }
        fn prod_bin(x: Uint128, y: Uint128, m: Uint128) -> Uint128 {
            if y == 1u64 {
                return x;
            } else if y == 0u64 {
                return Uint128::ZERO;
            }
            let mut res = prod_bin(x, y >> 1, m);
            res = sum(res, res, m);
            if (y.as_u64() & 1) == 1 {
                res = sum(res, x, m);
            }
            res
        }
    }
    */
}

Line	Count	Source
1		use core::fmt::{self, Display};
2		use core::ops::{Add, Div, Mul, Sub};
3
4		#[derive(Copy, Clone, Debug)]
5		pub struct SignedInteger<T> {
6		abs: T,
7		negative: bool,
8		}
9
10		impl<T> From<T> for SignedInteger<T> {
11		#[inline]
12	0	fn from(u: T) -> Self {
13	0	Self { abs: u, negative: false }
14	0	}
15		}
16		impl<T: From<u64>> SignedInteger<T> {
17		#[inline]
18	0	pub fn positive_u64(u: u64) -> Self {
19	0	Self { abs: T::from(u), negative: false }
20	0	}
21		}
22
23		impl<T: Display> Display for SignedInteger<T> {
24		#[inline]
25	0	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
26	0	if self.negative {
27	0	write!(f, "-")?;
28	0	}
29	0	write!(f, "{}", self.abs)
30	0	}
31		}
32
33		impl<T: Copy> SignedInteger<T> {
34		#[inline]
35	0	pub const fn abs(&self) -> T {
36	0	self.abs
37	0	}
38
39		#[inline]
40	0	pub const fn negative(&self) -> bool {
41	0	self.negative
42	0	}
43		}
44
45		impl<T: Sub<Output = T> + Add<Output = T> + Ord> Sub for SignedInteger<T> {
46		type Output = Self;
47		#[inline]
48		#[track_caller]
49	0	fn sub(self, other: Self) -> Self::Output {
50	0	match (self.negative, other.negative) {
51		(false, false) \| (true, true) => {
52	0	if self.abs < other.abs {
53	0	Self { negative: !self.negative, abs: other.abs - self.abs }
54		} else {
55	0	Self { negative: self.negative, abs: self.abs - other.abs }
56		}
57		}
58	0	(false, true) \| (true, false) => Self { negative: self.negative, abs: self.abs + other.abs },
59		}
60	0	}
61		}
62
63		impl<T: Mul<Output = T>> Mul for SignedInteger<T> {
64		type Output = Self;
65		#[inline]
66		#[track_caller]
67	0	fn mul(self, rhs: Self) -> Self::Output {
68	0	Self { negative: self.negative ^ rhs.negative, abs: self.abs * rhs.abs }
69	0	}
70		}
71
72		impl<T: Div<Output = T>> Div for SignedInteger<T> {
73		type Output = Self;
74		#[inline]
75		#[track_caller]
76	0	fn div(self, rhs: Self) -> Self::Output {
77	0	Self { negative: self.negative ^ rhs.negative, abs: self.abs / rhs.abs }
78	0	}
79		}
80
81		impl<T: PartialEq + PartialEq<u64>> PartialEq for SignedInteger<T> {
82	0	fn eq(&self, other: &Self) -> bool {
83	0	if self.abs == 0 && other.abs == 0 {
84		// neg/pos zeros are considered equal
85	0	return true;
86	0	}
87	0	self.negative == other.negative && self.abs == other.abs
88	0	}
89		}
90
91		impl<T: PartialEq + PartialEq<u64>> Eq for SignedInteger<T> {}
92
93		impl<T: PartialOrd + PartialEq<u64>> PartialOrd for SignedInteger<T> {
94	0	fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
95	0	if self.abs == 0 && other.abs == 0 {
96		// neg/pos zeros are considered equal
97	0	return Some(std::cmp::Ordering::Equal);
98	0	}
99	0	match (self.negative, other.negative) {
100	0	(false, false) => self.abs.partial_cmp(&other.abs),
101	0	(true, true) => other.abs.partial_cmp(&self.abs),
102	0	(true, false) => Some(std::cmp::Ordering::Less),
103	0	(false, true) => Some(std::cmp::Ordering::Greater),
104		}
105	0	}
106		}
107
108		#[cfg(test)]
109		mod tests {
110		use crate::{Uint192, int::SignedInteger};
111
112		fn from_u64(val: u64) -> SignedInteger<Uint192> {
113		SignedInteger::from(Uint192::from_u64(val))
114		}
115
116		#[test]
117		fn test_partial_eq() {
118		assert_eq!(from_u64(0), SignedInteger::from(Uint192::ZERO));
119		assert_eq!(from_u64(0), from_u64(10) - from_u64(10));
120		assert_eq!(from_u64(0), from_u64(10) - from_u64(20) - from_u64(10) * (from_u64(0) - from_u64(1))); // 0 == 10 - 20 -(-10)
121		assert_eq!(from_u64(0) - from_u64(1000), from_u64(0) - from_u64(1000)); // -1000 = -1000
122		assert_eq!(from_u64(1000), from_u64(1000));
123		}
124
125		#[test]
126		fn test_partial_cmp() {
127		// Test cases related to 0 and equality
128		assert!(from_u64(0) >= from_u64(10) - from_u64(20) - from_u64(10) * (from_u64(0) - from_u64(1))); // pos 0 >= neg 0
129		assert!(from_u64(0) <= from_u64(10) - from_u64(20) - from_u64(10) * (from_u64(0) - from_u64(1))); // pos 0 <= neg 0
130
131		// Test all possible neg/pos combinations
132		assert!(from_u64(100) > from_u64(0) - from_u64(1000)); // pos > neg
133		assert!(from_u64(0) - from_u64(100) < from_u64(10)); // neg < pos
134		assert!(from_u64(0) - from_u64(1000) < from_u64(0) - from_u64(100)); // -1000 < -100
135		assert!(from_u64(0) - from_u64(1000) <= from_u64(0) - from_u64(1000)); // -1000 <= -1000
136		assert!(from_u64(0) - from_u64(1000) >= from_u64(0) - from_u64(1000)); // -1000 >= -1000
137		assert!(from_u64(1000) > from_u64(100));
138		assert!(from_u64(100) < from_u64(1000));
139		assert!(from_u64(1000) >= from_u64(1000));
140		assert!(from_u64(100) <= from_u64(100));
141		}
142		}