fin_decimal 0.4.0

//! Multi-word integer primitives backing the wide decimal types
//! ([`Decimal128`](crate::Decimal128) and [`Decimal256`](crate::Decimal256)).
//!
//! Values are held as little-endian arrays of 64-bit limbs. Two division
//! flavors exist:
//!
//! * [`div_words_by_pow10`] / [`div_rem_u128_pow10`] divide by the
//!   compile-time constant 10^DIGITS via Möller–Granlund reciprocal
//!   multiplication, with the reciprocal computed at compile time: no
//!   division instructions and no `__udivti3` on any architecture. (LLVM
//!   does NOT strength-reduce 128-bit division by a constant, so this must
//!   be done by hand; `scripts/check_asm.sh` verifies it stays that way.)
//!   All decimal re-scaling — the hot path — goes through these, and they
//!   are `const fn`, so whole expressions can be evaluated at compile time.
//! * For runtime divisors the key primitive is [`div_2by1`]: a 128-bit-by-
//!   64-bit division with a 64-bit quotient. On `x86_64` with the `asm`
//!   feature it is the native `div` instruction; the portable fallback is
//!   Knuth's long division on 32-bit half-digits, substantially cheaper than
//!   the compiler's full 128-by-128 `__udivti3`. Word-based long division and
//!   Knuth's algorithm D are built on top of it.

use crate::AmountErrorKind;
use crate::AmountSign;
use crate::Rounding;
use core::cmp::Ordering;

/// Computes 10^pow as u64. Valid for `pow <= 19` (10^19 < 2^64).
///
/// Call sites that need the value folded into generated code must wrap the
/// call in an inline `const { .. }` block: without it, a cross-crate call
/// can survive into codegen and defeat constant propagation entirely.
#[inline]
pub(crate) const fn upow10(pow: u32) -> u64 {
    const P10: [u64; 20] = [
        1,
        10,
        100,
        1_000,
        10_000,
        100_000,
        1_000_000,
        10_000_000,
        100_000_000,
        1_000_000_000,
        10_000_000_000,
        100_000_000_000,
        1_000_000_000_000,
        10_000_000_000_000,
        100_000_000_000_000,
        1_000_000_000_000_000,
        10_000_000_000_000_000,
        100_000_000_000_000_000,
        1_000_000_000_000_000_000,
        10_000_000_000_000_000_000,
    ];
    P10[pow as usize]
}

/// Divide `(hi << 64) | lo` by `d`. Requires `hi < d` so the quotient fits in
/// a single limb. Returns `(quotient, remainder)`.
#[cfg(all(target_arch = "x86_64", feature = "asm"))]
#[inline]
pub(crate) fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) {
    debug_assert!(hi < d);
    let q: u64;
    let r: u64;
    unsafe {
        core::arch::asm!(
            "div {d}",
            d = in(reg) d,
            inout("rax") lo => q,
            inout("rdx") hi => r,
            options(pure, nomem, nostack)
        );
    }
    (q, r)
}

/// Divide the 128-bit value `(n1 << 64) | n0` by `d`, which must be
/// normalized (top bit set) with `n1 < d`. Returns `(quotient, remainder)`.
///
/// Knuth long division on 32-bit half-digits (the classic `udiv_qrnnd`
/// construction); avoids the compiler's 128-by-128 division builtin.
#[cfg(any(not(target_arch = "x86_64"), not(feature = "asm")))]
fn udiv_qrnnd(n1: u64, n0: u64, d: u64) -> (u64, u64) {
    debug_assert!(d >= 1 << 63 && n1 < d);
    let d_high = d >> 32;
    let d_low = d & 0xffff_ffff;

    let mut q_high = n1 / d_high;
    let mut r_high = n1 - q_high * d_high;
    let m = q_high * d_low;
    r_high = (r_high << 32) | (n0 >> 32);
    // Fine tune the estimate. `wrapping_add(d)` may overflow; when it does the
    // logical value of r_high exceeds any m, detected via `r_high >= d`.
    if r_high < m {
        q_high -= 1;
        r_high = r_high.wrapping_add(d);
        if r_high >= d && r_high < m {
            q_high -= 1;
            r_high = r_high.wrapping_add(d);
        }
    }
    r_high = r_high.wrapping_sub(m);

    let mut q_low = r_high / d_high;
    let mut r_low = r_high - q_low * d_high;
    let m = q_low * d_low;
    r_low = (r_low << 32) | (n0 & 0xffff_ffff);
    if r_low < m {
        q_low -= 1;
        r_low = r_low.wrapping_add(d);
        if r_low >= d && r_low < m {
            q_low -= 1;
            r_low = r_low.wrapping_add(d);
        }
    }
    r_low = r_low.wrapping_sub(m);

    ((q_high << 32) | q_low, r_low)
}

/// Divide `(hi << 64) | lo` by `d`. Requires `hi < d` so the quotient fits in
/// a single limb. Returns `(quotient, remainder)`.
#[cfg(any(not(target_arch = "x86_64"), not(feature = "asm")))]
#[inline]
pub(crate) fn div_2by1(hi: u64, lo: u64, d: u64) -> (u64, u64) {
    debug_assert!(hi < d);
    if hi == 0 {
        return (lo / d, lo % d);
    }
    let s = d.leading_zeros();
    if s == 0 {
        udiv_qrnnd(hi, lo, d)
    } else {
        let (q, r) = udiv_qrnnd((hi << s) | (lo >> (64 - s)), lo << s, d << s);
        (q, r >> s)
    }
}

/// In-place long division of `u` by the single word `d`; `u` becomes the
/// quotient. Returns the remainder.
#[inline]
pub(crate) fn div_words_by_word(u: &mut [u64], d: u64) -> u64 {
    debug_assert!(d != 0);
    let mut rem: u64 = 0;
    for x in u.iter_mut().rev() {
        // Leading limbs smaller than the divisor produce a zero quotient limb
        // and pass straight into the remainder; real values rarely fill all
        // limbs, so this skips most divisions.
        if rem == 0 && *x < d {
            rem = *x;
            *x = 0;
            continue;
        }
        let (q, r) = div_2by1(rem, *x, d);
        *x = q;
        rem = r;
    }
    rem
}

/// Per-scale constants for reciprocal division by 10^DIGITS: the normalized
/// divisor (top bit set), its normalization shift, and the Möller–Granlund
/// reciprocal `floor((2^128 - 1) / dn) - 2^64`. Evaluated at compile time
/// (the 128-bit division here is const-folded, it never reaches codegen).
const fn pow10_norm(digits: u32) -> (u64, u32, u64) {
    let d = upow10(digits);
    let s = d.leading_zeros();
    let dn = d << s;
    let v = (u128::MAX / (dn as u128) - (1u128 << 64)) as u64;
    (dn, s, v)
}

/// Möller–Granlund 2-by-1 division by a normalized divisor `dn` (top bit
/// set) with the precomputed reciprocal `v`: one 64x64->128 multiply plus
/// fixups, no division instruction. Requires `u1 < dn`.
#[inline]
const fn div_2by1_recip(u1: u64, u0: u64, dn: u64, v: u64) -> (u64, u64) {
    // (q1, q0) estimate = v * u1 + (u1:u0); cannot overflow u128 since
    // u1 <= dn - 1 and v + 2^64 = floor((2^128 - 1) / dn).
    let q = (v as u128) * (u1 as u128) + (((u1 as u128) << 64) | u0 as u128);
    let mut q1 = ((q >> 64) as u64).wrapping_add(1);
    let q0 = q as u64;
    let mut r = u0.wrapping_sub(q1.wrapping_mul(dn));
    // First fixup, branchless: a data-dependent branch here would mispredict.
    let over = (r > q0) as u64;
    q1 = q1.wrapping_sub(over);
    r = r.wrapping_add(dn & over.wrapping_neg());
    // Second fixup is almost never taken; keep it a predictable branch.
    if r >= dn {
        q1 += 1;
        r -= dn;
    }
    (q1, r)
}

/// One step of long division by 10^DIGITS: divides `(rem << 64) | limb` for
/// `rem < 10^DIGITS`, normalizing into the reciprocal division above. All
/// per-scale constants fold at monomorphization time.
#[inline]
const fn div_2by1_pow10<const DIGITS: u8>(rem: u64, limb: u64) -> (u64, u64) {
    let (dn, s, v) = const { pow10_norm(DIGITS as u32) };
    // (rem:limb) << s; for s == 0 (DIGITS = 19) the shifts fold away.
    let u1 = if s == 0 {
        rem
    } else {
        (rem << s) | (limb >> (64 - s))
    };
    let (q, r) = div_2by1_recip(u1, limb << s, dn, v);
    (q, r >> s)
}

/// In-place long division of `u` by the compile-time constant 10^DIGITS;
/// `u` becomes the quotient. Returns the remainder.
///
/// Every step is a reciprocal multiplication with compile-time constants —
/// no division instructions and no `__udivti3` on any architecture. (LLVM
/// does NOT strength-reduce 128-bit division by a constant, so writing
/// `u128 / const` would lower to a libcall; `scripts/check_asm.sh` guards
/// against regressing this.) Leading limbs smaller than the divisor are
/// skipped entirely, which covers typical financial magnitudes.
#[inline]
pub(crate) const fn div_words_by_pow10<const DIGITS: u8>(u: &mut [u64]) -> u64 {
    if DIGITS == 0 {
        return 0; // dividing by 1
    }
    let d = const { upow10(DIGITS as u32) };
    let mut rem: u64 = 0;
    let mut i = u.len();
    while i > 0 {
        i -= 1;
        // Leading limbs smaller than the divisor produce a zero quotient limb
        // and pass straight into the remainder.
        if rem == 0 && u[i] < d {
            rem = u[i];
            u[i] = 0;
            continue;
        }
        let (q, r) = div_2by1_pow10::<DIGITS>(rem, u[i]);
        u[i] = q;
        rem = r;
    }
    rem
}

/// Divides a full-range `u128` by the compile-time constant 10^DIGITS,
/// returning `(quotient, remainder)`, using the same reciprocal
/// multiplication as [`div_words_by_pow10`] — division-instruction-free.
#[inline]
pub(crate) const fn div_rem_u128_pow10<const DIGITS: u8>(n: u128) -> (u128, u64) {
    if DIGITS == 0 {
        return (n, 0);
    }
    let d = const { upow10(DIGITS as u32) };
    let hi = (n >> 64) as u64;
    let lo = n as u64;
    let (q1, r1) = if hi < d {
        (0, hi) // covers the common case of values below 2^64
    } else {
        div_2by1_pow10::<DIGITS>(0, hi)
    };
    let (q0, r) = div_2by1_pow10::<DIGITS>(r1, lo);
    (((q1 as u128) << 64) | q0 as u128, r)
}

/// `acc = acc * mul + add`, in place. Returns `true` if a carry overflowed
/// out of the limb array.
#[inline]
pub(crate) const fn mul_add_word(acc: &mut [u64], mul: u64, add: u64) -> bool {
    let mut carry: u128 = add as u128;
    let mut i = 0;
    while i < acc.len() {
        let v = (acc[i] as u128) * (mul as u128) + carry;
        acc[i] = v as u64;
        carry = v >> 64;
        i += 1;
    }
    carry != 0
}

/// `c = a * b` (schoolbook). Requires `c.len() >= a.len() + b.len()`; `c` is
/// fully overwritten (excess top limbs are zeroed).
#[inline]
pub(crate) const fn mul_words(c: &mut [u64], a: &[u64], b: &[u64]) {
    debug_assert!(c.len() >= a.len() + b.len());
    let mut i = 0;
    while i < c.len() {
        c[i] = 0;
        i += 1;
    }
    let mut i = 0;
    while i < a.len() {
        // c[i..i+b.len()] += a[i] * b, tracking the carry into c[i+b.len()].
        let mut carry: u128 = 0;
        let mut j = 0;
        while j < b.len() {
            let v = (a[i] as u128) * (b[j] as u128) + (c[i + j] as u128) + carry;
            c[i + j] = v as u64;
            carry = v >> 64;
            j += 1;
        }
        c[i + b.len()] = carry as u64;
        i += 1;
    }
}

/// Number of significant limbs (at least 1, even for zero).
#[inline]
pub(crate) fn sig_limbs(l: &[u64]) -> usize {
    let mut n = l.len();
    while n > 1 && l[n - 1] == 0 {
        n -= 1;
    }
    n
}

/// True if every limb is zero.
#[inline]
pub(crate) const fn is_zero(l: &[u64]) -> bool {
    let mut i = 0;
    while i < l.len() {
        if l[i] != 0 {
            return false;
        }
        i += 1;
    }
    true
}

/// Lexicographic comparison of equal-length limb arrays as unsigned integers.
#[inline]
pub(crate) fn cmp_words(a: &[u64], b: &[u64]) -> Ordering {
    debug_assert_eq!(a.len(), b.len());
    for i in (0..a.len()).rev() {
        match a[i].cmp(&b[i]) {
            Ordering::Equal => {}
            o => return o,
        }
    }
    Ordering::Equal
}

/// Shifts left by one bit in place; returns the bit shifted out of the top.
#[inline]
pub(crate) fn shl1(l: &mut [u64]) -> u64 {
    let mut carry = 0u64;
    for x in l.iter_mut() {
        let next = *x >> 63;
        *x = (*x << 1) | carry;
        carry = next;
    }
    carry
}

/// Maximum dividend size (in limbs) supported by [`div_knuth`].
pub(crate) const KNUTH_MAX_M: usize = 8;

/// Knuth's algorithm D: divides `u` (m limbs) by `v` (n limbs, `2 <= n <= m`,
/// top limb of `v` non-zero). Writes the quotient into `q` (needs `m - n + 1`
/// limbs; excess limbs are zeroed) and the remainder into `r` (needs `n`
/// limbs; excess limbs are zeroed).
///
/// Direct port of the classic 32-bit `bn_div_ex` construction, widened to
/// 64-bit limbs with `u128`/`i128` intermediates.
#[inline]
pub(crate) fn div_knuth(q: &mut [u64], r: &mut [u64], u: &[u64], v: &[u64]) {
    let m = u.len();
    let n = v.len();
    debug_assert!((2..=m).contains(&n), "use div_words_by_word for n == 1");
    debug_assert!(v[n - 1] != 0);
    debug_assert!(m <= KNUTH_MAX_M);
    debug_assert!(q.len() > m - n && r.len() >= n);

    let mut un = [0u64; KNUTH_MAX_M + 1]; // normalized u, gains one limb
    let mut vn = [0u64; KNUTH_MAX_M]; // normalized v

    // Normalize so the divisor's top bit is set (required by div_2by1's
    // estimate quality guarantees).
    let s = v[n - 1].leading_zeros();
    if s == 0 {
        vn[..n].copy_from_slice(v);
        un[..m].copy_from_slice(u);
        un[m] = 0;
    } else {
        for i in (1..n).rev() {
            vn[i] = (v[i] << s) | (v[i - 1] >> (64 - s));
        }
        vn[0] = v[0] << s;
        un[m] = u[m - 1] >> (64 - s);
        for i in (1..m).rev() {
            un[i] = (u[i] << s) | (u[i - 1] >> (64 - s));
        }
        un[0] = u[0] << s;
    }

    for x in q.iter_mut() {
        *x = 0;
    }

    const B: u128 = 1u128 << 64;
    let vtop = vn[n - 1];
    let vnext = vn[n - 2];

    // Main loop, reducing un digit by digit from the top.
    let mut j = m - n;
    loop {
        debug_assert!(un[j + n] <= vtop);
        // Estimate the quotient digit from the top two dividend limbs.
        let (mut qhat, mut rhat): (u128, u128) = if un[j + n] == vtop {
            // The generic path needs un[j+n] < vtop; here the division splits
            // exactly: (vtop*B + lo) / vtop == B + lo/vtop, rem lo%vtop.
            (
                B + (un[j + n - 1] / vtop) as u128,
                (un[j + n - 1] % vtop) as u128,
            )
        } else {
            let (q0, r0) = div_2by1(un[j + n], un[j + n - 1], vtop);
            (q0 as u128, r0 as u128)
        };
        // Fine tune: while the estimate is too large by the next-limb test.
        // The `qhat >= B` arm must be checked first so the multiplication
        // below never overflows u128.
        loop {
            if qhat >= B || qhat * (vnext as u128) > ((rhat << 64) | un[j + n - 2] as u128) {
                qhat -= 1;
                rhat += vtop as u128;
                if rhat < B {
                    continue;
                }
            }
            break;
        }

        // Multiply and subtract: un[j..j+n+1] -= qhat * vn.
        let mut qd = qhat as u64;
        let mut borrow: i128 = 0;
        for i in 0..n {
            let p = (qd as u128) * (vn[i] as u128);
            let t = (un[i + j] as i128) - borrow - ((p as u64) as i128);
            un[i + j] = t as u64;
            borrow = ((p >> 64) as i128) - (t >> 64);
        }
        let t = (un[j + n] as i128) - borrow;
        un[j + n] = t as u64;

        // If we borrowed past the top, the estimate was one too large:
        // add the divisor back and adjust.
        if t < 0 {
            qd -= 1;
            let mut carry: u128 = 0;
            for i in 0..n {
                let sum = un[i + j] as u128 + vn[i] as u128 + carry;
                un[i + j] = sum as u64;
                carry = sum >> 64;
            }
            un[j + n] = un[j + n].wrapping_add(carry as u64);
        }
        q[j] = qd;

        if j == 0 {
            break;
        }
        j -= 1;
    }

    // Denormalize the remainder (shift right by s bits).
    if s == 0 {
        r[..n].copy_from_slice(&un[..n]);
    } else {
        for i in 0..n - 1 {
            r[i] = (un[i] >> s) | (un[i + 1] << (64 - s));
        }
        r[n - 1] = un[n - 1] >> s;
    }
    for x in r.iter_mut().skip(n) {
        *x = 0;
    }
}

/// Multiplies two sign-magnitude decimals and re-scales by 10^DIGITS with
/// the given rounding mode. This is THE multiplication core shared by every
/// backing width: `W` is the magnitude size in limbs (1 for `Decimal`, 2 for
/// `Decimal128`, 4 for `Decimal256`) and `W2` must equal `2 * W` (passed
/// separately because stable Rust cannot compute const-generic expressions).
///
/// Returns `None` on overflow of the symmetric range: thanks to `MIN ==
/// -MAX` on every backing, "fits" is uniformly "top bit of the top limb is
/// clear", regardless of sign.
#[inline]
pub(crate) const fn dec_mul<const DIGITS: u8, const W: usize, const W2: usize>(
    a_neg: bool,
    a_mag: &[u64; W],
    b_neg: bool,
    b_mag: &[u64; W],
    mode: Rounding,
) -> Option<(bool, [u64; W])> {
    debug_assert!(W2 == 2 * W);
    // Single-limb operands take the two-limb core: this replaces the
    // WxW-limb schoolbook product with a 2x2 one for typical magnitudes.
    // Only 1-limb operands qualify because there the narrow core provably
    // never overflows (product < 2^128, re-scaled result < 2^115 for
    // DIGITS >= 1), so no work is ever thrown away; the rare DIGITS = 0
    // overflow falls through and recomputes wide.
    if W > 2 && is_zero(a_mag.split_at(1).1) && is_zero(b_mag.split_at(1).1) {
        let a2 = [a_mag[0], 0];
        let b2 = [b_mag[0], 0];
        if let Some((n2, m2)) = dec_mul::<DIGITS, 2, 4>(a_neg, &a2, b_neg, &b2, mode) {
            let mut mag = [0u64; W];
            mag[0] = m2[0];
            mag[1] = m2[1];
            return Some((n2, mag));
        }
    }
    let neg = a_neg != b_neg;
    let mut prod = [0u64; W2];
    mul_words(&mut prod, a_mag, b_mag);
    // HalfUp (the operators' default) folds its rounding into the division:
    // adding half the divisor before truncating rounds half away from zero
    // with no data-dependent branch. `mode` is statically known at almost
    // every call site, so this `if` disappears after inlining.
    let half_up = matches!(mode, Rounding::HalfUp);
    if half_up {
        // The product is at most (2^(64W-1) - 1)^2: no carry out of W2 limbs.
        mul_add_word(&mut prod, 1, const { upow10(DIGITS as u32) / 2 });
    }
    let rem = div_words_by_pow10::<DIGITS>(&mut prod);
    if !is_zero(prod.split_at(W).1) {
        return None;
    }
    let mut mag = [0u64; W];
    let mut i = 0;
    while i < W {
        mag[i] = prod[i];
        i += 1;
    }
    if !half_up
        && round_up_by_cmp(
            cmp_twice_rem_u64(rem, const { upow10(DIGITS as u32) }),
            rem == 0,
            mag[0] & 1 != 0,
            neg,
            mode,
        )
        && mul_add_word(&mut mag, 1, 1)
    {
        return None;
    }
    if mag[W - 1] >> 63 != 0 {
        return None;
    }
    Some((neg, mag))
}

/// Divides two sign-magnitude decimals (scaling the dividend by 10^DIGITS
/// first) with the given rounding mode. The division core shared by every
/// backing width; `WP1` must equal `W + 1` (the scaled dividend gains one
/// limb). Single-limb divisors take word division (through [`div_2by1`], so
/// the `asm` feature applies); wider divisors take Knuth's algorithm D over
/// significant limbs only.
///
/// Returns `None` if `b` is zero or on overflow of the symmetric range.
#[inline]
pub(crate) fn dec_div<const DIGITS: u8, const W: usize, const WP1: usize>(
    a_neg: bool,
    a_mag: &[u64; W],
    b_neg: bool,
    b_mag: &[u64; W],
    mode: Rounding,
) -> Option<(bool, [u64; W])> {
    debug_assert!(WP1 == W + 1);
    if is_zero(b_mag) {
        return None;
    }
    // Same narrow-operand tiering as dec_mul; `b` is known non-zero here, so
    // `None` from the narrow core only means quotient overflow: fall through.
    if W > 2 && is_zero(a_mag.split_at(2).1) && is_zero(b_mag.split_at(2).1) {
        let a2 = [a_mag[0], a_mag[1]];
        let b2 = [b_mag[0], b_mag[1]];
        if let Some((n2, m2)) = dec_div::<DIGITS, 2, 3>(a_neg, &a2, b_neg, &b2, mode) {
            let mut mag = [0u64; W];
            mag[0] = m2[0];
            mag[1] = m2[1];
            return Some((n2, mag));
        }
    }
    let neg = a_neg != b_neg;
    let mut num = [0u64; WP1];
    mul_words(&mut num, a_mag, &[const { upow10(DIGITS as u32) }]);
    let (rem_cmp, rem_zero) = div_mag_in_place(&mut num, b_mag);
    round_quotient::<W>(&num, rem_cmp, rem_zero, neg, mode)
}

/// Divides `num` in place by the significant limbs of `div` (`num` becomes
/// the quotient), returning the ordering of `2 * remainder` against the
/// divisor and whether the remainder was zero — exactly what the rounding
/// decision needs. The divide-and-classify tail shared by [`dec_div`] and
/// [`dec_mul_div`]. `div` must be non-zero; single-limb divisors take word
/// division (through [`div_2by1`], so the `asm` feature applies), wider
/// ones Knuth's algorithm D over significant limbs only.
///
/// `#[inline(always)]`: each caller invokes this exactly once per
/// monomorphization, and inlining is what lets the slice lengths become
/// compile-time constants again (without it, every width shares one
/// dynamic-length body and the division paths measurably regress).
#[inline(always)]
fn div_mag_in_place(num: &mut [u64], div: &[u64]) -> (Ordering, bool) {
    let n = sig_limbs(div);
    if n == 1 {
        // Single-limb divisor: word-based long division. For 1-limb callers
        // this is the only possible path and the rest folds away.
        let d = div[0];
        let rem = div_words_by_word(num, d);
        return (cmp_twice_rem_u64(rem, d), rem == 0);
    }
    let mut q = [0u64; KNUTH_MAX_M];
    let mut r = [0u64; KNUTH_MAX_M];
    let m = sig_limbs(num);
    if m < n {
        // Fewer dividend limbs than divisor limbs: quotient is zero and the
        // whole dividend is the remainder.
        r[..m].copy_from_slice(&num[..m]);
        num[..m].fill(0);
    } else {
        div_knuth(&mut q[..m - n + 1], &mut r[..n], &num[..m], &div[..n]);
        num.copy_from_slice(&q[..num.len()]);
    }
    // Order 2 * rem against the divisor: doubling in place is safe because
    // the zero test happens first.
    let rem_zero = is_zero(&r[..n]);
    let carry = shl1(&mut r[..n]);
    let cmp = if carry != 0 {
        Ordering::Greater
    } else {
        cmp_words(&r[..n], &div[..n])
    };
    (cmp, rem_zero)
}

/// Applies the rounding decision to a quotient's low `W` limbs and runs the
/// symmetric-range checks shared by [`dec_div`] and [`dec_mul_div`]: `None`
/// if the quotient needs more than `W` limbs, if rounding up carries out of
/// them, or if the result's top bit is set.
#[inline]
fn round_quotient<const W: usize>(
    q: &[u64],
    rem_cmp: Ordering,
    rem_zero: bool,
    neg: bool,
    mode: Rounding,
) -> Option<(bool, [u64; W])> {
    if !is_zero(q.split_at(W).1) {
        return None;
    }
    let mut mag = [0u64; W];
    mag.copy_from_slice(&q[..W]);
    if round_up_by_cmp(rem_cmp, rem_zero, mag[0] & 1 != 0, neg, mode)
        && mul_add_word(&mut mag, 1, 1)
    {
        return None;
    }
    if mag[W - 1] >> 63 != 0 {
        return None;
    }
    Some((neg, mag))
}

/// `a * b / c` on the exact double-width product with a single rounding at
/// the end: the width-generic muldiv core shared by [`Decimal128`](crate::
/// Decimal128) and [`Decimal256`](crate::Decimal256). `W2` must equal
/// `2 * W`. The `b / c` scale cancels at the call site, so unlike
/// [`dec_div`] there is no re-scaling by 10^DIGITS: the dividend is the
/// plain 2W-limb product (which is why `KNUTH_MAX_M` is 2W for the widest
/// backing). Single-limb divisors take word division (through [`div_2by1`],
/// so the `asm` feature applies); wider divisors take Knuth's algorithm D
/// over significant limbs only.
///
/// Returns `None` if `c` is zero or the result overflows the symmetric
/// range.
#[inline]
pub(crate) fn dec_mul_div<const W: usize, const W2: usize>(
    a_neg: bool,
    a_mag: &[u64; W],
    b_neg: bool,
    b_mag: &[u64; W],
    c_neg: bool,
    c_mag: &[u64; W],
    mode: Rounding,
) -> Option<(bool, [u64; W])> {
    debug_assert!(W2 == 2 * W);
    if is_zero(c_mag) {
        return None;
    }
    // Same narrow-operand tiering as dec_mul/dec_div; `c` is known non-zero
    // here, so `None` from the narrow core only means quotient overflow:
    // fall through and recompute wide.
    if W > 2
        && is_zero(a_mag.split_at(2).1)
        && is_zero(b_mag.split_at(2).1)
        && is_zero(c_mag.split_at(2).1)
    {
        let a2 = [a_mag[0], a_mag[1]];
        let b2 = [b_mag[0], b_mag[1]];
        let c2 = [c_mag[0], c_mag[1]];
        if let Some((n2, m2)) = dec_mul_div::<2, 4>(a_neg, &a2, b_neg, &b2, c_neg, &c2, mode) {
            let mut mag = [0u64; W];
            mag[0] = m2[0];
            mag[1] = m2[1];
            return Some((n2, mag));
        }
    }
    let neg = (a_neg != b_neg) != c_neg;
    // Full-width fixed-length product on purpose: bounding the slices by
    // sig_limbs was measured slower (dynamic lengths defeat unrolling, and
    // the division below already skips leading zero limbs via sig_limbs).
    let mut prod = [0u64; W2];
    mul_words(&mut prod, a_mag, b_mag);
    let (rem_cmp, rem_zero) = div_mag_in_place(&mut prod, c_mag);
    round_quotient::<W>(&prod, rem_cmp, rem_zero, neg, mode)
}

/// Remainder of `a` by non-zero `b` on `W`-limb magnitudes (`|a| mod |b|`):
/// the shared core of the wide backings' `%` operator. Single-limb divisors
/// take word division; a dividend narrower than the divisor passes through;
/// wider divisors take Knuth's algorithm D over significant limbs only.
#[inline]
pub(crate) fn rem_mag<const W: usize>(a: &[u64; W], b: &[u64; W]) -> [u64; W] {
    let n = sig_limbs(b);
    let mut r = [0u64; W];
    if n == 1 {
        let mut q = *a;
        r[0] = div_words_by_word(&mut q, b[0]);
        return r;
    }
    // |a| < |b|: the dividend is the remainder. A value compare (not just a
    // limb-count one) — this fast path is measurably load-bearing for
    // remainders of same-width operands.
    if cmp_words(a, b) == Ordering::Less {
        return *a;
    }
    let m = sig_limbs(a);
    let mut q = [0u64; W];
    div_knuth(&mut q[..m - n + 1], &mut r[..n], &a[..m], &b[..n]);
    r
}

/// Divides a `W`-limb magnitude by a single non-zero word, applying the
/// shared rounding decision: the wide backings' integer-divisor core
/// (`div_int_rounded`, and the division half of `round_dp`). The result
/// magnitude never exceeds the input's: rounding up requires `d >= 2`.
#[inline]
pub(crate) fn div_mag_by_word<const W: usize>(
    neg: bool,
    mag: &[u64; W],
    d: u64,
    mode: Rounding,
) -> [u64; W] {
    let mut q = *mag;
    let r = div_words_by_word(&mut q, d);
    if round_up_by_cmp(cmp_twice_rem_u64(r, d), r == 0, q[0] & 1 != 0, neg, mode) {
        mul_add_word(&mut q, 1, 1);
    }
    q
}

/// Rounds a `W`-limb magnitude to a multiple of `step` (a single-limb
/// runtime divisor, `10^(DIGITS - dp)` at the call sites): divide with the
/// shared rounding decision, then rebuild `(q + up) * step`. The caller
/// applies its own top-of-range check; the multiply-back itself cannot
/// carry out of `W` limbs (`(q + up) * step <= mag + step`, and `mag`'s top
/// bit is clear on every backing).
#[inline]
pub(crate) fn round_mag_to_step<const W: usize>(
    neg: bool,
    mag: &[u64; W],
    step: u64,
    mode: Rounding,
) -> [u64; W] {
    let mut q = div_mag_by_word(neg, mag, step, mode);
    let carry = mul_add_word(&mut q, step, 0);
    debug_assert!(!carry);
    q
}

/// Largest-remainder even split of a `W`-limb magnitude into `n > 0`
/// parts: returns `(low, high, high_count)` where `high_count` parts of
/// `high` plus `n - high_count` parts of `low` sum back to `mag` exactly.
/// When the split is inexact, `high == low + 1` ulp; when exact,
/// `high == low` and `high_count == 0`. The invariant-bearing core shared
/// by the wide backings (the i64 backing keeps a const copy).
#[inline]
pub(crate) fn split_mag_even<const W: usize>(mag: &[u64; W], n: u32) -> ([u64; W], [u64; W], u32) {
    let mut low = *mag;
    let r = div_words_by_word(&mut low, n as u64) as u32;
    let mut high = low;
    if r != 0 {
        // An inexact split implies n >= 2, hence low <= MAX / 2: the bump
        // cannot carry out of W limbs.
        let carry = mul_add_word(&mut high, 1, 1);
        debug_assert!(!carry);
    }
    (low, high, r)
}

/// The square-root rounding decision shared by every backing: whether to
/// bump the floor root, given the remainder facts from the isqrt cores.
/// Half-mode ties cannot occur (`(s + 1/2)^2` is never an integer), so the
/// three half modes agree.
#[inline]
pub(crate) const fn sqrt_round_up(mode: Rounding, rem_zero: bool, rem_gt_root: bool) -> bool {
    // Exhaustive on purpose, like round_up_by_cmp: a new Rounding variant
    // must force a decision here, not silently become a half mode.
    match mode {
        Rounding::Down => false,
        Rounding::Up => !rem_zero,
        Rounding::HalfUp | Rounding::HalfDown | Rounding::HalfEven => rem_gt_root,
    }
}

/// Decides whether a truncated magnitude should be rounded up (away from the
/// integer's current value on the magnitude scale), given the ordering of
/// `2 * remainder` against the divisor.
///
/// The `Down`/`Up` modes are directional (floor/ceil), matching
/// `mul_rounded`/`div_rounded` on [`Decimal`](crate::Decimal): flooring a
/// negative result means rounding its magnitude up.
#[inline]
pub(crate) const fn round_up_by_cmp(
    twice_rem_vs_div: Ordering,
    rem_is_zero: bool,
    quo_is_odd: bool,
    is_neg: bool,
    mode: Rounding,
) -> bool {
    if rem_is_zero {
        return false;
    }
    match mode {
        Rounding::HalfEven => {
            matches!(twice_rem_vs_div, Ordering::Greater)
                || (matches!(twice_rem_vs_div, Ordering::Equal) && quo_is_odd)
        }
        Rounding::HalfUp => !matches!(twice_rem_vs_div, Ordering::Less),
        Rounding::HalfDown => matches!(twice_rem_vs_div, Ordering::Greater),
        Rounding::Down => is_neg,
        Rounding::Up => !is_neg,
    }
}

/// Orders `2 * rem` against `d` without overflow, given `rem < d`.
#[inline]
pub(crate) const fn cmp_twice_rem_u64(rem: u64, d: u64) -> Ordering {
    debug_assert!(rem < d);
    let other = d - rem;
    if rem < other {
        Ordering::Less
    } else if rem > other {
        Ordering::Greater
    } else {
        Ordering::Equal
    }
}

/// Floor square root of a `u128`, seeded from an `f64` reciprocal-sqrt
/// estimate and refined with multiply-only Newton steps — no division
/// instructions and no libcalls anywhere (`core` has no `f64::sqrt`, and a
/// `u128` division here would be a `__udivti3` call), so it is usable in
/// const and stays fast at runtime.
///
/// Steps: the classic exponent-halving bit trick gives `y ~ 1/sqrt(n)` to
/// ~3%; four Newton iterations for the *reciprocal* square root
/// (`y *= 1.5 - 0.5*n*y*y`, multiply-only) converge quadratically to f64
/// precision; `n * y` is then `sqrt(n)` within a few thousand ULPs (f64
/// carries 53 bits, the result up to 64). One further Newton step on the
/// exact integer residual — `s += (n - s^2) / (2s)`, with the division again
/// a multiplication by `y` — lands within ~1, and a short walk pins the
/// exact floor.
///
/// Requires `n < 2^127` (all callers pass at most `i64::MAX * 10^19 <
/// 2^127`), which keeps every intermediate square inside `u128`.
///
/// Benchmarks (`cargo bench`, "Square root") put this at ~31 ns/op flat
/// across magnitudes — beating the pure integer Newton on wide inputs.
/// Core's table-seeded `u128::isqrt` is faster still on x86, but its
/// codegen contains a `__udivti3` libcall, which the crate bans; the
/// public `sqrt` API therefore uses `u64::isqrt` (native word divisions
/// only) for single-limb inputs and this core for wider ones. Pinned to
/// `u128::isqrt` by tests along with [`isqrt_newton`].
pub const fn isqrt_f64(n: u128) -> u64 {
    if n == 0 {
        return 0;
    }
    debug_assert!(n >> 127 == 0);
    let nf = n as f64;
    let mut y = f64::from_bits(0x5FE6EB50C7B537A9u64.wrapping_sub(nf.to_bits() >> 1));
    let mut i = 0;
    while i < 4 {
        y = y * (1.5 - 0.5 * nf * y * y);
        i += 1;
    }
    // `nf * y` is within 2^-49 relative of sqrt(n) < 2^63.5, so `s` stays
    // well below 2^64 and `s * s` cannot overflow (or wrap) below.
    let s0 = (nf * y) as u128;
    let s2 = s0 * s0;
    let mut s = if n >= s2 {
        s0 + (((n - s2) as f64) * (0.5 * y)) as u128
    } else {
        s0.saturating_sub((((s2 - n) as f64) * (0.5 * y)) as u128)
    };
    while s > 0 && s * s > n {
        s -= 1;
    }
    while (s + 1) * (s + 1) <= n {
        s += 1;
    }
    s as u64
}

/// Floor square root of any `u128` (unlike [`isqrt_f64`], the top bit may
/// be set): single-limb inputs take core's `u64::isqrt` (native word
/// divisions only), wider ones the division-free f64-seeded core, with
/// inputs at or above `2^127` reduced by `isqrt(n) = 2 * isqrt(n / 4)`
/// plus a one-step walk. Const-evaluable and division-instruction-free
/// above one limb, like [`isqrt_f64`].
pub(crate) const fn isqrt_u128_any(n: u128) -> u64 {
    if n >> 64 == 0 {
        return (n as u64).isqrt();
    }
    let mut s = if n >> 127 == 0 {
        isqrt_f64(n) as u128
    } else {
        (isqrt_f64(n >> 2) as u128) << 1
    };
    // The reduction loses at most one: walk to the exact floor. `s` stays
    // below 2^64, so `s * s` and `(s + 1)^2` fit `u128`.
    while s > 0 && s * s > n {
        s -= 1;
    }
    while s < u64::MAX as u128 && (s + 1) * (s + 1) <= n {
        s += 1;
    }
    s as u64
}

/// Floor square root of any `u128` plus the remainder facts the decimal
/// rounding needs: `(root, rem_zero, rem_gt_root)` as in [`isqrt_words`].
/// The shared single/two-limb fast tier of every backing's `sqrt`.
#[inline]
pub(crate) const fn isqrt_u128_with_rem(n: u128) -> (u64, bool, bool) {
    let s = isqrt_u128_any(n);
    let rem = n - (s as u128) * (s as u128);
    (s, rem == 0, rem > s as u128)
}

/// Floor square root of an `N`-limb magnitude, plus the remainder facts the
/// decimal rounding needs: returns `(root, rem_zero, rem_gt_root)` with
/// `root^2 <= n < (root + 1)^2`, `rem = n - root^2`, and `rem_gt_root` the
/// tie-free half-mode test (`rem > root` means `n` is past
/// `(root + 1/2)^2`; a half tie cannot occur since `(root + 1/2)^2` is
/// never an integer).
///
/// `R` is the root width and `R2 == 2 * R >= N` (the root's square is
/// formed in `R2` limbs). Two-limb inputs take the `u128` cores; wider
/// ones an f64-seeded Newton iteration whose division steps run through
/// [`div_knuth`] (so the `asm` feature applies and no 128-bit builtins are
/// emitted).
pub(crate) fn isqrt_words<const N: usize, const R: usize, const R2: usize>(
    n: &[u64; N],
) -> ([u64; R], bool, bool) {
    const {
        assert!(R2 == 2 * R && R2 >= N && N <= KNUTH_MAX_M && N >= 3);
    }
    let m = sig_limbs(n);
    let mut root = [0u64; R];
    if m <= 2 {
        let v = (n[0] as u128) | ((n[1] as u128) << 64);
        let (s, rem_zero, rem_gt_root) = isqrt_u128_with_rem(v);
        root[0] = s;
        return (root, rem_zero, rem_gt_root);
    }

    // Seed from the top 127..=128 bits: r0 = isqrt(n >> e) << (e / 2) for an
    // even `e` is a floor-accurate estimate with ~2^-63 relative error.
    let bits = 64 * m - n[m - 1].leading_zeros() as usize;
    let e = (bits - 128).next_multiple_of(2);
    let (wo, bo) = (e / 64, e % 64);
    let w = |i: usize| -> u128 { if i < N { n[i] as u128 } else { 0 } };
    let t: u128 = if bo == 0 {
        w(wo) | (w(wo + 1) << 64)
    } else {
        (w(wo) >> bo) | (w(wo + 1) << (64 - bo)) | (w(wo + 2) << (128 - bo))
    };
    let s0 = isqrt_u128_any(t);
    let (hw, hb) = (e / 2 / 64, e / 2 % 64);
    root[hw] = s0 << hb;
    if hb != 0 {
        root[hw + 1] = s0 >> (64 - hb);
    }

    // Newton (Heron) steps: root = (root + n / root) / 2. The seed carries
    // ~63 correct bits and each step doubles them, so one step suffices for
    // 3-limb inputs (root < 2^96) and two pin 5-limb ones (root < 2^160) —
    // and the second step is a whole div_knuth, so skipping it on the
    // narrow shape matters. `n >= 2^128` makes `root >= 2^64`, so the
    // divisor always has at least two significant limbs.
    let steps = if m <= 3 { 1 } else { 2 };
    let mut i = 0;
    while i < steps {
        // The true root is at least 2^64 (n >= 2^128 here); clamp an iterate
        // that undershot below one limb so div_knuth always sees a two-limb
        // divisor.
        if sig_limbs(&root) < 2 {
            root[0] = 0;
            root[1] = 1;
        }
        let rl = sig_limbs(&root);
        let mut q = [0u64; N];
        let mut r = [0u64; R];
        div_knuth(&mut q[..m - rl + 1], &mut r[..rl], &n[..m], &root[..rl]);
        // sum = root + q, then halve. The sum fits N limbs: both terms are
        // within a hair of the true root (< 2^(64R - 1)).
        let mut carry: u128 = 0;
        let mut sum = [0u64; N];
        let mut j = 0;
        while j < N {
            let s = q[j] as u128 + if j < R { root[j] as u128 } else { 0 } + carry;
            sum[j] = s as u64;
            carry = s >> 64;
            j += 1;
        }
        debug_assert!(carry == 0);
        let mut j = 0;
        while j < N {
            let lo = sum[j] >> 1;
            let hi = if j + 1 < N { sum[j + 1] << 63 } else { 0 };
            sum[j] = lo | hi;
            j += 1;
        }
        if !is_zero(&sum[R..]) {
            // The iterate is at most ceil(sqrt(n)) + O(1), so exceeding R
            // limbs only happens when the true root sits at the very top of
            // the range (e.g. all-ones input, root 2^(64R) - 1): saturate
            // and let the fixup walk land on the exact floor.
            root = [u64::MAX; R];
        } else {
            root.copy_from_slice(&sum[..R]);
        }
        i += 1;
    }

    // Walk to the exact floor (the Newton iterate is within a couple of
    // ulps), then extract the remainder facts.
    let mut n_pad = [0u64; R2];
    n_pad[..N].copy_from_slice(n);
    let mut sq = [0u64; R2];
    loop {
        mul_words(&mut sq, &root, &root);
        if cmp_words(&sq, &n_pad) == Ordering::Greater {
            // root -= 1
            let mut j = 0;
            while j < R {
                let (v, borrow) = root[j].overflowing_sub(1);
                root[j] = v;
                if !borrow {
                    break;
                }
                j += 1;
            }
        } else {
            break;
        }
    }
    loop {
        let mut r1 = root;
        if mul_add_word(&mut r1, 1, 1) {
            break; // root + 1 does not fit R limbs, so its square exceeds n
        }
        mul_words(&mut sq, &r1, &r1);
        if cmp_words(&sq, &n_pad) == Ordering::Greater {
            break;
        }
        root = r1;
    }

    // rem = n - root^2 (fits R2 limbs, non-negative by construction).
    mul_words(&mut sq, &root, &root);
    let mut rem = [0u64; R2];
    let mut borrow: u128 = 0;
    let mut j = 0;
    while j < R2 {
        let d = n_pad[j] as u128;
        let s = sq[j] as u128 + borrow;
        if d >= s {
            rem[j] = (d - s) as u64;
            borrow = 0;
        } else {
            rem[j] = ((d + (1u128 << 64)) - s) as u64;
            borrow = 1;
        }
        j += 1;
    }
    debug_assert!(borrow == 0);
    let rem_zero = is_zero(&rem);
    let mut root_pad = [0u64; R2];
    root_pad[..R].copy_from_slice(&root);
    let rem_gt_root = cmp_words(&rem, &root_pad) == Ordering::Greater;
    (root, rem_zero, rem_gt_root)
}

/// Floor square root of a `u128` by pure integer Newton (Heron) iteration:
/// the seed `2^ceil(bits/2)` overshoots `sqrt(n)`, and each step
/// `x = (x + n/x) / 2` runs `n / x` through [`div_words_by_word`] (so the
/// `asm` feature's `div_2by1` applies and no `__udivti3` is ever emitted).
/// The strictly decreasing sequence stops exactly at the floor. ~6-7
/// division steps against [`isqrt_f64`]'s zero: ~20 ns/op on money-sized
/// inputs but ~60 ns at 2^120 (the iteration count grows with magnitude).
/// Kept as the measured comparison implementation, pinned to
/// `u128::isqrt` by tests.
pub fn isqrt_newton(n: u128) -> u64 {
    if n == 0 {
        return 0;
    }
    let bits = 128 - n.leading_zeros();
    let shift = bits.div_ceil(2);
    let div = |x: u128| {
        let mut q = [n as u64, (n >> 64) as u64];
        // The seed can be exactly 2^64 (one bit past a word); its quotient
        // is a plain shift. Every later iterate fits a word.
        if x >> 64 != 0 {
            return n >> 64;
        }
        div_words_by_word(&mut q, x as u64);
        (q[0] as u128) | ((q[1] as u128) << 64)
    };
    let mut x0: u128 = 1 << shift;
    let mut x1 = (x0 + div(x0)) >> 1;
    while x1 < x0 {
        x0 = x1;
        x1 = (x0 + div(x0)) >> 1;
    }
    x0 as u64
}

/// Converts a decimal string into a sign and a 4-limb magnitude scaled by
/// 10^`scale`, rounding excess fractional digits with `mode`.
///
/// The limb-based counterpart of
/// [`parse_decimal_i64_rounded`](crate::parse_decimal_i64_rounded), with the
/// same grammar, rounding semantics (first dropped digit + sticky flag +
/// parity for `HalfEven`) and error kinds. Range checks against a concrete
/// backing type are the caller's job; `Overflow` is returned only when the
/// scaled magnitude exceeds 256 bits.
pub(crate) const fn parse_decimal_mag_rounded<const W: usize>(
    src: &str,
    scale: u8,
    mode: Rounding,
) -> Result<(bool, [u64; W]), AmountErrorKind> {
    let bytes: &[u8] = src.as_bytes();
    let scale = scale as i64;
    debug_assert!(scale <= 19);
    // Kept digits are buffered in one word and folded into the wide `acc` once
    // per `CHUNK_DIGITS`. That flush is where `Overflow` is *detected*, so 19
    // is not adjustable: a shorter budget changes which error a string that
    // both overflows and carries a stray byte returns ("35 digits then 'x'" is
    // `InvalidDigit` at 19, `Overflow` at 16).
    //
    // 19 is also the **overflow licence** for the wrapping accumulation: the
    // `chunk_len == CHUNK_DIGITS` branch below empties `chunk` on the 19th
    // digit, so every `chunk * 10 + d` runs with `chunk_len < 19`, hence
    // `chunk <= 10^18 - 1` and `chunk * 10 + 9 < 10^19 < u64::MAX`. The
    // `debug_assert`s at the two accumulation sites restate that as code.
    // Wrapping matters because the crate builds with `overflow-checks = true`:
    // a checked mul+add puts two branches on the loop-carried `chunk` chain
    // the next digit waits on, measured at ~15% of the whole parse.
    const CHUNK_DIGITS: u32 = 19;
    let mut acc = [0u64; W];
    let mut chunk: u64 = 0;
    let mut chunk_len: u32 = 0;
    // Set at the flush site itself, never derived: the inline finish below is
    // valid only while `acc` is still zero.
    let mut flushed = false;

    if bytes.is_empty() {
        return Err(AmountErrorKind::Empty);
    }
    let (sign, digits): (bool, &[u8]) = match bytes[0] {
        b'+' => (false, bytes.split_at(1).1),
        b'-' => (true, bytes.split_at(1).1),
        _ => (false, bytes),
    };
    let n = digits.len();
    let mut di = 0;

    // Integer phase: always kept, never a round digit, never moves the
    // fractional counter, so the body is pure accumulation. `wrapping_sub`
    // classifies and converts in one compare and doubles as the terminator —
    // the loop stops on the '.', on a stray byte and on end-of-input alike.
    while di < n {
        let d = digits[di].wrapping_sub(b'0');
        if d > 9 {
            break;
        }
        debug_assert!(chunk_len < CHUNK_DIGITS);
        chunk = chunk.wrapping_mul(10).wrapping_add(d as u64);
        chunk_len += 1;
        di += 1;
        if chunk_len == CHUNK_DIGITS {
            if mul_add_word(&mut acc, upow10(CHUNK_DIGITS), chunk) {
                return Err(AmountErrorKind::Overflow);
            }
            chunk = 0;
            chunk_len = 0;
            flushed = true;
        }
    }
    let int_digits = di;

    // Terminator, then the fraction phase. `nfrac < scale` is the one test the
    // seed's three (`point <= scale`, `point == scale + 1`, `point += 1`)
    // collapse to: digits past it are dropped, and dropping them is the cold
    // tail's business, which is what keeps the over-scale shape off a re-parse.
    // The loop condition is tested *after* the flush, so a chunk that fills on
    // the very digit that exhausts `scale` is flushed on that digit as the seed
    // does — the flush reports `Overflow` ahead of any stray byte to come.
    let mut nfrac: i64 = 0;
    if di < n {
        if digits[di] != b'.' {
            return Err(AmountErrorKind::InvalidDigit);
        }
        di += 1;
        while di < n && nfrac < scale {
            let d = digits[di].wrapping_sub(b'0');
            if d > 9 {
                break;
            }
            debug_assert!(chunk_len < CHUNK_DIGITS);
            chunk = chunk.wrapping_mul(10).wrapping_add(d as u64);
            chunk_len += 1;
            nfrac += 1;
            di += 1;
            if chunk_len == CHUNK_DIGITS {
                if mul_add_word(&mut acc, upow10(CHUNK_DIGITS), chunk) {
                    return Err(AmountErrorKind::Overflow);
                }
                chunk = 0;
                chunk_len = 0;
                flushed = true;
            }
        }
    }

    if di == n {
        // Only at end of input, so `Empty` cannot pre-empt an `InvalidDigit`:
        // a remaining byte is either a digit or an error already returned.
        if int_digits == 0 && nfrac == 0 {
            return Err(AmountErrorKind::Empty);
        }
        if !flushed {
            // Nothing was dropped, so `nfrac <= scale`, the round digit is 0
            // and every mode yields `round_up == false` — `mode` is dead here.
            // `acc` is still zero, so the partial-chunk flush and the
            // trailing-zero scaling collapse into one exact 128-bit multiply:
            // `chunk < 10^19`, exponent at most 19, product below 10^38.
            let wide = (chunk as u128) * (upow10((scale - nfrac) as u32) as u128);
            let lo = wide as u64;
            let hi = (wide >> 64) as u64;
            if W == 1 && hi != 0 {
                // Exactly where `mul_add_word` would have carried out.
                return Err(AmountErrorKind::Overflow);
            }
            let out: &mut [u64] = &mut acc;
            out[0] = lo;
            if W > 1 {
                out[1] = hi;
            }
            return Ok((sign, acc));
        }
    }
    // `nfrac + 1` is the seed's `point`, including the value its
    // `if point == 0 { point = 1 }` fixup lands on when no '.' was seen.
    parse_tail(
        sign,
        digits,
        di,
        scale,
        mode,
        acc,
        chunk,
        chunk_len,
        nfrac + 1,
        int_digits != 0 || nfrac != 0,
    )
}

/// Everything after the digit scan that a money-sized string never executes:
/// the dropped-precision remainder, the partial-chunk flush, the trailing-zero
/// scaling, and the one and only five-mode rounding decision.
///
/// `#[inline(never)]` is the point of the split — code that never runs still
/// costs, by setting the hot path's register allocation, spill slots and branch
/// layout. The line is drawn at the *end* of the scan rather than the front, so
/// the scan's state is handed forward instead of thrown away and the input is
/// never walked twice. Entered on two cold exits: `flushed`, or over-scale
/// (`di < n`, fractional budget spent), where `point == scale + 1` makes the
/// first remaining digit the round digit exactly as the seed's counter would.
#[inline(never)]
#[allow(clippy::too_many_arguments)] // the scan state the split hands forward
const fn parse_tail<const W: usize>(
    sign: bool,
    digits: &[u8],
    mut di: usize,
    scale: i64,
    mode: Rounding,
    mut acc: [u64; W],
    chunk: u64,
    chunk_len: u32,
    mut point: i64,
    mut digit: bool,
) -> Result<(bool, [u64; W]), AmountErrorKind> {
    let mut round_digit: u64 = 0; // first fractional digit dropped beyond `scale`
    let mut sticky = false; // any non-zero digit beyond `round_digit`
    while di < digits.len() {
        let x = digits[di].wrapping_sub(b'0') as u64;
        if x > 9 {
            // A second '.' included: `point > 0` here, so the seed rejects it.
            return Err(AmountErrorKind::InvalidDigit);
        }
        digit = true;
        if point == scale + 1 {
            round_digit = x;
        } else if x != 0 {
            sticky = true;
        }
        point += 1;
        di += 1;
    }
    if !digit {
        return Err(AmountErrorKind::Empty);
    }
    // Flush the partial chunk before scaling and rounding (the HalfEven parity
    // below inspects the fully-accumulated magnitude), then append zeros so the
    // magnitude carries exactly `scale` fractional digits.
    if chunk_len > 0 && mul_add_word(&mut acc, upow10(chunk_len), chunk) {
        return Err(AmountErrorKind::Overflow);
    }
    if point <= scale && mul_add_word(&mut acc, upow10((scale + 1 - point) as u32), 0) {
        return Err(AmountErrorKind::Overflow);
    }
    let round_up = match mode {
        Rounding::HalfUp => round_digit >= 5,
        Rounding::HalfDown => round_digit > 5 || (round_digit == 5 && sticky),
        Rounding::HalfEven => round_digit > 5 || (round_digit == 5 && (sticky || acc[0] & 1 != 0)),
        Rounding::Down => false,
        Rounding::Up => round_digit != 0 || sticky,
    };
    if round_up && mul_add_word(&mut acc, 1, 1) {
        return Err(AmountErrorKind::Overflow);
    }
    Ok((sign, acc))
}

/// ASCII digit pairs `"00"`..`"99"` (most significant digit first), so digit
/// extraction below can take two digits per strength-reduced division by 100,
/// halving the serial divide chain that dominates formatting latency.
const DIGIT_PAIRS: [u8; 200] = {
    let mut t = [0u8; 200];
    let mut i = 0;
    while i < 100 {
        t[2 * i] = b'0' + (i / 10) as u8;
        t[2 * i + 1] = b'0' + (i % 10) as u8;
        i += 1;
    }
    t
};

/// Formats a scaled decimal magnitude into `buf`, mirroring the semantics of
/// [`str_i64`](crate::str_i64): trailing fractional zeros are trimmed when no
/// precision is given, an explicit precision rounds (half away from zero) or
/// zero-pads, and the string is built right-to-left in `buf`.
pub(crate) fn str_mag<'a>(
    mag: &[u64],
    neg: bool,
    frac_digits: usize,
    precision: Option<usize>,
    sign: AmountSign,
    buf: &'a mut [u8],
) -> Option<&'a str> {
    debug_assert!(frac_digits <= 19 && mag.len() <= 4);

    // Fast path: a single-limb magnitude with default precision — every
    // i64-backed value and typical wide-type magnitudes. It reads `mag`
    // directly, so the common case skips the working-copy memcpy entirely,
    // and writes digits straight into `buf` right-to-left: no intermediate
    // digit buffer and no second pass. At most 21 digits, the point and a
    // sign are written (23 bytes), so with `buf.len() >= 24` the positions
    // never underflow. (The `== 1` guard also proves `mag` is non-empty.)
    if precision.is_none() && buf.len() >= 24 && sig_limbs(mag) == 1 {
        let mut pos = buf.len();
        let mut v = mag[0];
        let is_zero_val = v == 0;
        // Trim trailing fractional zeros, then emit the remaining fraction.
        let mut fd = frac_digits;
        while fd > 0 && v.is_multiple_of(10) {
            v /= 10;
            fd -= 1;
        }
        if fd > 0 {
            while fd >= 2 {
                let p = (v % 100) as usize * 2;
                v /= 100;
                pos -= 2;
                buf[pos] = DIGIT_PAIRS[p];
                buf[pos + 1] = DIGIT_PAIRS[p + 1];
                fd -= 2;
            }
            if fd > 0 {
                pos -= 1;
                buf[pos] = (v % 10) as u8 + b'0';
                v /= 10;
            }
            pos -= 1;
            buf[pos] = b'.';
        }
        // Integer part, at least one digit (possibly zero).
        while v >= 100 {
            let p = (v % 100) as usize * 2;
            v /= 100;
            pos -= 2;
            buf[pos] = DIGIT_PAIRS[p];
            buf[pos + 1] = DIGIT_PAIRS[p + 1];
        }
        if v >= 10 {
            let p = v as usize * 2;
            pos -= 2;
            buf[pos] = DIGIT_PAIRS[p];
            buf[pos + 1] = DIGIT_PAIRS[p + 1];
        } else {
            pos -= 1;
            buf[pos] = v as u8 + b'0';
        }
        match (sign, neg && !is_zero_val) {
            (AmountSign::None, _) => {}
            (_, true) => {
                pos -= 1;
                buf[pos] = b'-';
            }
            (AmountSign::Always, false) if !is_zero_val => {
                pos -= 1;
                buf[pos] = b'+';
            }
            _ => {}
        }
        // All bytes written are ASCII digits, '.', '+' or '-'.
        return core::str::from_utf8(&buf[pos..]).ok();
    }

    // Remaining paths mutate the magnitude, so take a working copy. One spare
    // limb: rounding up at a precision boundary can carry.
    let mut work = [0u64; 5];
    work[..mag.len()].copy_from_slice(mag);

    if let Some(precision) = precision {
        // If the requested precision cannot fit, bail out immediately.
        if precision >= buf.len() {
            return None;
        }
        // Round to the requested precision first (half away from zero), then
        // scale back so digit extraction below stays uniform.
        if precision < frac_digits {
            let round_scale = upow10((frac_digits - precision) as u32);
            let dropped = div_words_by_word(&mut work, round_scale);
            if dropped.wrapping_shl(1) >= round_scale {
                mul_add_word(&mut work, 1, 1);
            }
            mul_add_word(&mut work, round_scale, 0);
        }
    }

    // Extract ASCII digits, least significant first. While the value spans
    // several limbs, full 19-digit chunks come off with one reciprocal
    // division by 10^19 per pass; the last limb's digits come straight off
    // that word. The value is < 10^78 (and each chunk emits real positional
    // digits), so `nd` stays well within the 96-digit buffer.
    let mut digits = [b'0'; 96];
    let mut nd = 0usize;
    while sig_limbs(&work) > 1 {
        let mut rem = div_words_by_pow10::<19>(&mut work);
        // Exactly 19 digits: nine pairs, then the leftover top digit
        // (rem < 10^19 / 100^9 = 10).
        for _ in 0..9 {
            let p = (rem % 100) as usize * 2;
            rem /= 100;
            digits[nd] = DIGIT_PAIRS[p + 1];
            digits[nd + 1] = DIGIT_PAIRS[p];
            nd += 2;
        }
        digits[nd] = rem as u8 + b'0';
        nd += 1;
    }
    let mut rem = work[0];
    while rem >= 100 {
        let p = (rem % 100) as usize * 2;
        rem /= 100;
        digits[nd] = DIGIT_PAIRS[p + 1];
        digits[nd + 1] = DIGIT_PAIRS[p];
        nd += 2;
    }
    if rem >= 10 {
        let p = rem as usize * 2;
        digits[nd] = DIGIT_PAIRS[p + 1];
        digits[nd + 1] = DIGIT_PAIRS[p];
        nd += 2;
    } else if rem != 0 {
        digits[nd] = rem as u8 + b'0';
        nd += 1;
    }
    // Guarantee at least one integer digit (possibly zero).
    if nd < frac_digits + 1 {
        nd = frac_digits + 1;
    }

    // Compose right-to-left: [pad zeros][fraction]['.'][integer][sign].
    let mut pos = buf.len();
    macro_rules! push {
        ($c:expr) => {{
            if pos == 0 {
                return None;
            }
            pos -= 1;
            buf[pos] = $c;
        }};
    }

    // Number of fractional digits to print.
    let print_frac = match precision {
        Some(p) => p,
        None => {
            // Trim trailing zeros of the fractional part.
            let mut skip = 0;
            while skip < frac_digits && digits[skip] == b'0' {
                skip += 1;
            }
            frac_digits - skip
        }
    };

    if print_frac > 0 {
        // Zero-pad on the right if the precision exceeds the scale.
        for _ in frac_digits..print_frac {
            push!(b'0');
        }
        let start = frac_digits - print_frac.min(frac_digits);
        for digit in digits.iter().take(frac_digits).skip(start) {
            push!(*digit);
        }
        push!(b'.');
    }

    for digit in digits.iter().take(nd).skip(frac_digits) {
        push!(*digit);
    }

    let is_zero_val = is_zero(mag);
    match (sign, neg && !is_zero_val) {
        (AmountSign::None, _) => {}
        (_, true) => push!(b'-'),
        (AmountSign::Always, false) if !is_zero_val => push!(b'+'),
        _ => {}
    }

    // All bytes written are ASCII digits, '.', '+' or '-'.
    core::str::from_utf8(&buf[pos..]).ok()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::AmountSign;

    /// 10^19, the largest power of ten that fits in a u64.
    const POW10_19: u64 = upow10(19);

    /// xorshift64* pseudo-random generator for deterministic tests.
    struct Rng(u64);
    impl Rng {
        fn next(&mut self) -> u64 {
            let mut x = self.0;
            x ^= x >> 12;
            x ^= x << 25;
            x ^= x >> 27;
            self.0 = x;
            x.wrapping_mul(0x2545F4914F6CDD1D)
        }
    }

    #[test]
    fn test_div_2by1() {
        // Edge and random cases, checked against native u128 division.
        let mut rng = Rng(0x12345678DEADBEEF);
        let check = |hi: u64, lo: u64, d: u64| {
            let n = ((hi as u128) << 64) | lo as u128;
            let (q, r) = div_2by1(hi, lo, d);
            assert_eq!(q as u128, n / d as u128, "q for {n} / {d}");
            assert_eq!(r as u128, n % d as u128, "r for {n} / {d}");
        };
        check(0, 0, 1);
        check(0, u64::MAX, 1);
        check(0, u64::MAX, u64::MAX);
        check(u64::MAX - 1, u64::MAX, u64::MAX);
        check(1, 0, 2);
        check(1, 1, 2);
        check(0x7FFF_FFFF_FFFF_FFFF, u64::MAX, 1 << 63);
        check(POW10_19 - 1, u64::MAX, POW10_19);
        for _ in 0..20000 {
            let d = rng.next() | 1;
            let hi = rng.next() % d;
            let lo = rng.next();
            check(hi, lo, d);
            // small divisors too
            let ds = (rng.next() % 1000) + 1;
            check(rng.next() % ds, rng.next(), ds);
        }
    }

    #[test]
    fn test_div_words_by_pow10() {
        // Cross-check the constant-divisor path against the generic word
        // division for several scales, sizes and magnitude shapes.
        fn check<const DIGITS: u8>(rng: &mut Rng) {
            let d = upow10(DIGITS as u32);
            for shape in 0..4 {
                // 1, 2, 3 or 4 significant limbs
                let mut u = [0u64; 4];
                for x in u.iter_mut().take(shape + 1) {
                    *x = rng.next();
                }
                let mut a = u;
                let ra = div_words_by_word(&mut a, d);
                let mut b = u;
                let rb = div_words_by_pow10::<DIGITS>(&mut b);
                assert_eq!((a, ra), (b, rb), "DIGITS={DIGITS} u={u:?}");
                // Also with an 8-limb buffer as used by dec256_mul.
                let mut w = [0u64; 8];
                w[..4].copy_from_slice(&u);
                let mut wa = w;
                let ra = div_words_by_word(&mut wa, d);
                let mut wb = w;
                let rb = div_words_by_pow10::<DIGITS>(&mut wb);
                assert_eq!((wa, ra), (wb, rb), "DIGITS={DIGITS} w={w:?}");
            }
        }
        let mut rng = Rng(0x9E3779B97F4A7C15);
        for _ in 0..2000 {
            check::<0>(&mut rng);
            check::<1>(&mut rng);
            check::<4>(&mut rng);
            check::<8>(&mut rng);
            check::<19>(&mut rng);
        }
        // Zero value stays zero on the fast path.
        let mut z = [0u64; 4];
        assert_eq!(div_words_by_pow10::<4>(&mut z), 0);
        assert_eq!(z, [0u64; 4]);
    }

    #[test]
    fn test_div_words_by_word() {
        let mut u = [u64::MAX, u64::MAX, u64::MAX, u64::MAX];
        let rem = div_words_by_word(&mut u, 10);
        // (2^256 - 1) % 10 == 5
        assert_eq!(rem, 5);
        // Rebuild: q * 10 + 5 == 2^256 - 1.
        assert!(!mul_add_word(&mut u, 10, 5));
        assert_eq!(u, [u64::MAX; 4]);
    }

    #[test]
    fn test_mul_words() {
        let mut c = [0u64; 4];
        mul_words(&mut c, &[u64::MAX, u64::MAX], &[u64::MAX, u64::MAX]);
        // (2^128 - 1)^2 = 2^256 - 2^129 + 1
        assert_eq!(c, [1, 0, u64::MAX - 1, u64::MAX]);

        let mut rng = Rng(0xC0FFEE);
        for _ in 0..5000 {
            let a = ((rng.next() as u128) << 64) | rng.next() as u128;
            let b = rng.next() as u128;
            let mut c = [0u64; 3];
            mul_words(&mut c, &[a as u64, (a >> 64) as u64], &[b as u64]);
            // Verify low 128 bits against native u128 multiplication.
            let lo = a.wrapping_mul(b);
            assert_eq!(c[0], lo as u64);
            assert_eq!(c[1], (lo >> 64) as u64);
        }
    }

    #[test]
    fn test_div_knuth_reconstruct() {
        // Validate q*v + r == u and r < v over random shapes.
        let mut rng = Rng(0xFEEDFACE);
        for iter in 0..20000 {
            let n = 2 + (rng.next() as usize) % 3; // 2..=4 divisor limbs
            let m = n + (rng.next() as usize) % (5 - n) + 1; // n+1..=5
            let mut u = [0u64; 5];
            let mut v = [0u64; 4];
            for x in u.iter_mut().take(m) {
                *x = rng.next();
            }
            for x in v.iter_mut().take(n) {
                *x = rng.next();
            }
            if v[n - 1] == 0 {
                v[n - 1] = 1 + (rng.next() >> 32);
            }
            // Sometimes force interesting top-limb relationships.
            if iter % 7 == 0 {
                u[m - 1] = v[n - 1];
            }
            if iter % 11 == 0 {
                v[n - 1] = 1 << 63;
            }

            let mut q = [0u64; 4];
            let mut r = [0u64; 4];
            div_knuth(&mut q[..m - n + 1], &mut r[..n], &u[..m], &v[..n]);

            // r < v
            assert_eq!(cmp_words(&r[..n], &v[..n]), Ordering::Less);
            // q * v + r == u  (reconstruct in 9 limbs)
            let mut prod = [0u64; 9];
            mul_words(&mut prod[..m + 1], &q[..m - n + 1], &v[..n]);
            let mut carry: u128 = 0;
            for i in 0..9 {
                let add = if i < n { r[i] as u128 } else { 0 };
                let s = prod[i] as u128 + add + carry;
                prod[i] = s as u64;
                carry = s >> 64;
            }
            assert_eq!(&prod[..m], &u[..m], "reconstruction failed");
            assert!(is_zero(&prod[m..]));
        }
    }

    /// Path flags reported by [`div_knuth_probed`], used to assert that the
    /// hard branches of algorithm D are actually exercised by the tests.
    #[derive(Default)]
    struct DivKnuthProbe {
        /// Largest number of D3 correction-loop decrements at any position.
        max_corrections: u32,
        /// The add-back step (D6, `t < 0`) fired.
        addback: bool,
        /// The `un[j+n] == vtop` estimate branch (`qhat` starts `>= B`) ran.
        eq_estimate: bool,
        /// `qhat >= B` survived the correction loop (must never happen: the
        /// multiply-subtract truncates `qhat` to one limb).
        qhat_overflow: bool,
    }

    /// Instrumented copy of [`div_knuth`] reporting which hard paths ran.
    ///
    /// Kept in lockstep with the real function by
    /// `test_div_knuth_d3_adversarial`, which asserts identical `q`/`r` on
    /// every case, so any drift between the two is caught immediately.
    fn div_knuth_probed(q: &mut [u64], r: &mut [u64], u: &[u64], v: &[u64]) -> DivKnuthProbe {
        let m = u.len();
        let n = v.len();
        let mut un = [0u64; KNUTH_MAX_M + 1];
        let mut vn = [0u64; KNUTH_MAX_M];
        let s = v[n - 1].leading_zeros();
        if s == 0 {
            vn[..n].copy_from_slice(v);
            un[..m].copy_from_slice(u);
            un[m] = 0;
        } else {
            for i in (1..n).rev() {
                vn[i] = (v[i] << s) | (v[i - 1] >> (64 - s));
            }
            vn[0] = v[0] << s;
            un[m] = u[m - 1] >> (64 - s);
            for i in (1..m).rev() {
                un[i] = (u[i] << s) | (u[i - 1] >> (64 - s));
            }
            un[0] = u[0] << s;
        }
        for x in q.iter_mut() {
            *x = 0;
        }
        const B: u128 = 1u128 << 64;
        let vtop = vn[n - 1];
        let vnext = vn[n - 2];
        let mut probe = DivKnuthProbe::default();
        let mut j = m - n;
        loop {
            let (mut qhat, mut rhat): (u128, u128) = if un[j + n] == vtop {
                probe.eq_estimate = true;
                (
                    B + (un[j + n - 1] / vtop) as u128,
                    (un[j + n - 1] % vtop) as u128,
                )
            } else {
                let (q0, r0) = div_2by1(un[j + n], un[j + n - 1], vtop);
                (q0 as u128, r0 as u128)
            };
            let mut corrections = 0u32;
            loop {
                if qhat >= B || qhat * (vnext as u128) > ((rhat << 64) | un[j + n - 2] as u128) {
                    qhat -= 1;
                    corrections += 1;
                    rhat += vtop as u128;
                    if rhat < B {
                        continue;
                    }
                }
                break;
            }
            probe.max_corrections = probe.max_corrections.max(corrections);
            if qhat >= B {
                probe.qhat_overflow = true;
            }
            let mut qd = qhat as u64;
            let mut borrow: i128 = 0;
            for i in 0..n {
                let p = (qd as u128) * (vn[i] as u128);
                let t = (un[i + j] as i128) - borrow - ((p as u64) as i128);
                un[i + j] = t as u64;
                borrow = ((p >> 64) as i128) - (t >> 64);
            }
            let t = (un[j + n] as i128) - borrow;
            un[j + n] = t as u64;
            if t < 0 {
                probe.addback = true;
                qd -= 1;
                let mut carry: u128 = 0;
                for i in 0..n {
                    let sum = un[i + j] as u128 + vn[i] as u128 + carry;
                    un[i + j] = sum as u64;
                    carry = sum >> 64;
                }
                un[j + n] = un[j + n].wrapping_add(carry as u64);
            }
            q[j] = qd;
            if j == 0 {
                break;
            }
            j -= 1;
        }
        if s == 0 {
            r[..n].copy_from_slice(&un[..n]);
        } else {
            for i in 0..n - 1 {
                r[i] = (un[i] >> s) | (un[i + 1] << (64 - s));
            }
            r[n - 1] = un[n - 1] >> s;
        }
        for x in r.iter_mut().skip(n) {
            *x = 0;
        }
        probe
    }

    /// Oracle for a single division: checks `r < v` and `q * v + r == u`,
    /// runs the probed copy, and asserts it agrees with [`div_knuth`].
    fn check_div_knuth(u: &[u64], v: &[u64]) -> DivKnuthProbe {
        let m = u.len();
        let n = v.len();
        let mut q = [0u64; KNUTH_MAX_M];
        let mut r = [0u64; KNUTH_MAX_M];
        div_knuth(&mut q[..m - n + 1], &mut r[..n], u, v);

        let mut q2 = [0u64; KNUTH_MAX_M];
        let mut r2 = [0u64; KNUTH_MAX_M];
        let probe = div_knuth_probed(&mut q2[..m - n + 1], &mut r2[..n], u, v);
        assert_eq!(q, q2, "probed copy diverged: u={u:?} v={v:?}");
        assert_eq!(r, r2, "probed copy diverged: u={u:?} v={v:?}");
        assert!(
            !probe.qhat_overflow,
            "qhat >= B survived the correction loop: u={u:?} v={v:?}"
        );

        assert_eq!(cmp_words(&r[..n], &v[..n]), Ordering::Less, "r >= v");
        let mut prod = [0u64; 2 * KNUTH_MAX_M + 1];
        mul_words(&mut prod[..m + 1], &q[..m - n + 1], &v[..n]);
        let mut carry: u128 = 0;
        for (i, p) in prod.iter_mut().enumerate() {
            let add = if i < n { r[i] as u128 } else { 0 };
            let s = *p as u128 + add + carry;
            *p = s as u64;
            carry = s >> 64;
        }
        assert_eq!(&prod[..m], u, "reconstruction failed: u={u:?} v={v:?}");
        assert!(is_zero(&prod[m..]));
        probe
    }

    /// Adversarial stress test for the D3 correction step, prompted by
    /// kolja.rs/algorithm-d: for an *odd* base, Knuth's proven bound
    /// `qhat <= q + 2` fails (three corrections can be needed) and
    /// implementations that truncate `qhat` to one digit break. Our base
    /// `2^64` is even, the correction loop is unbounded, and add-back covers
    /// the residual `q + 1` case, so the bug cannot apply; this test pins all
    /// of that empirically by hammering the article's worst-case digit
    /// pattern mapped to base 2^64 (divisor top limb at the normalization
    /// boundary, next limbs saturated, dividend top equal to the divisor
    /// top). At the end it asserts every hard path was actually exercised.
    #[test]
    fn test_div_knuth_d3_adversarial() {
        let iters: u64 = if cfg!(debug_assertions) {
            200_000
        } else {
            2_000_000
        };
        let mut rng = Rng(0xA160B17D);
        let specials: [u64; 8] = [
            1 << 63,
            (1 << 63) | 1,
            (1 << 63) - 1, // top bit clear: forces a normalization shift
            u64::MAX,
            u64::MAX - 1,
            (1 << 63) + (1 << 62),
            1, // maximal normalization shift (s == 63)
            (u64::MAX >> 1) + 1,
        ];
        let mut max_corrections = 0u32;
        let mut addback_seen = false;
        let mut eq_estimate_seen = false;
        let mut s_zero_seen = false;
        let mut s_shift_seen = false;
        let mut m_eq_n_seen = false;
        for iter in 0..iters {
            let n = 2 + (rng.next() as usize) % 3; // 2..=4 divisor limbs
            let m = n + (rng.next() as usize) % (6 - n); // n..=5 dividend limbs
            let mut u = [0u64; 5];
            let mut v = [0u64; 4];
            for x in u.iter_mut().take(m) {
                *x = match rng.next() % 4 {
                    0 => specials[(rng.next() as usize) % specials.len()],
                    1 => u64::MAX,
                    _ => rng.next(),
                };
            }
            for x in v.iter_mut().take(n) {
                *x = match rng.next() % 4 {
                    0 => specials[(rng.next() as usize) % specials.len()],
                    1 => u64::MAX,
                    _ => rng.next(),
                };
            }
            if v[n - 1] == 0 {
                v[n - 1] = specials[(rng.next() as usize) % specials.len()];
            }
            // Deterministic hammering of the article's digit pattern
            // (v_{n-1} = b/2, v_{n-2} = b-1, u_top = v_{n-1}, u_next = b-1)
            // and of the `qhat >= B` estimate path.
            match iter % 5 {
                0 => {
                    v[n - 1] = 1 << 63;
                    v[n - 2] = u64::MAX;
                    u[m - 1] = v[n - 1];
                    u[m - 2] = u64::MAX;
                }
                1 => u[m - 1] = v[n - 1],
                2 => {
                    v[n - 1] = 1 << 63;
                    u[m - 1] = (1 << 63) - 1;
                    u[m - 2] = u64::MAX;
                }
                _ => {}
            }

            let probe = check_div_knuth(&u[..m], &v[..n]);
            max_corrections = max_corrections.max(probe.max_corrections);
            addback_seen |= probe.addback;
            eq_estimate_seen |= probe.eq_estimate;
            s_zero_seen |= v[n - 1] >= 1 << 63;
            s_shift_seen |= v[n - 1] < 1 << 63;
            m_eq_n_seen |= m == n;
        }
        // Knuth's bound for an even base: at most two corrections. Both must
        // be reachable, and every hard path must have been hit at least once.
        assert_eq!(max_corrections, 2, "expected the 2-correction worst case");
        assert!(addback_seen, "add-back path never exercised");
        assert!(eq_estimate_seen, "qhat >= B estimate path never exercised");
        assert!(s_zero_seen && s_shift_seen, "both normalization paths");
        assert!(m_eq_n_seen, "m == n (single quotient limb) never exercised");
    }

    /// Fixed vectors pinning each hard path of algorithm D, so a regression
    /// in any of them fails deterministically (the adversarial test reaches
    /// the same paths, but through a PRNG). Vectors were mined from the
    /// adversarial search; each is asserted to still hit its path.
    #[test]
    fn test_isqrt_words_self_check() {
        // The floor-root property is self-verifying: root^2 <= n <
        // (root+1)^2, rem flags consistent with n - root^2. Random shapes
        // across both instantiations, biased toward boundary limb patterns.
        fn check<const N: usize, const R: usize, const R2: usize>(n: &[u64; N]) {
            let (root, rem_zero, rem_gt) = isqrt_words::<N, R, R2>(n);
            let mut n_pad = [0u64; R2];
            n_pad[..N].copy_from_slice(n);
            let mut sq = [0u64; R2];
            mul_words(&mut sq, &root, &root);
            assert_ne!(cmp_words(&sq, &n_pad), Ordering::Greater, "root^2 > n");
            // (root + 1)^2 > n, unless root + 1 overflows R limbs.
            let mut r1 = root;
            if !mul_add_word(&mut r1, 1, 1) {
                mul_words(&mut sq, &r1, &r1);
                assert_eq!(cmp_words(&sq, &n_pad), Ordering::Greater, "(root+1)^2 <= n");
            }
            // rem facts.
            mul_words(&mut sq, &root, &root);
            let mut rem = [0u64; R2];
            let mut borrow: i128 = 0;
            for j in 0..R2 {
                let d = n_pad[j] as i128 - sq[j] as i128 - borrow;
                rem[j] = d as u64;
                borrow = (d < 0) as i128;
            }
            assert_eq!(borrow, 0);
            assert_eq!(rem_zero, is_zero(&rem), "rem_zero");
            let mut rp = [0u64; R2];
            rp[..R].copy_from_slice(&root);
            assert_eq!(
                rem_gt,
                cmp_words(&rem, &rp) == Ordering::Greater,
                "rem_gt_root"
            );
        }

        // Deterministic extremes: all-ones inputs give an all-ones root, so
        // the "root + 1 does not fit R limbs" fixup break is exercised.
        check::<4, 2, 4>(&[u64::MAX; 4]);
        check::<6, 3, 6>(&[u64::MAX; 6]);
        check::<4, 2, 4>(&[0, 0, 0, u64::MAX]);
        check::<6, 3, 6>(&[0, 0, 0, 0, 0, u64::MAX]);

        let mut rng = Rng(0x15C4A57B00B5);
        for iter in 0..20000 {
            let mut n4 = [0u64; 4];
            let mut n6 = [0u64; 6];
            // Vary significant width, including 1..2-limb (u128 tier), the
            // >= 2^127 two-limb corner, and full-width values.
            let w4 = 1 + (rng.next() as usize) % 3; // callers cap at 3 limbs
            let w6 = 1 + (rng.next() as usize) % 5; // callers cap at 5 limbs
            for x in n4.iter_mut().take(w4) {
                *x = rng.next();
            }
            for x in n6.iter_mut().take(w6) {
                *x = rng.next();
            }
            match iter % 5 {
                0 => {
                    n4[w4 - 1] |= 1 << 63; // saturate the top bit
                    n6[w6 - 1] |= 1 << 63;
                }
                1 => {
                    n4[w4 - 1] = 1; // barely into the next limb
                    n6[w6 - 1] = 1;
                }
                2 => {
                    // Perfect squares (and their neighbors on iter % 8).
                    let s = rng.next() as u128;
                    let sq = s * s;
                    n4 = [sq as u64, (sq >> 64) as u64, 0, 0];
                    n6 = [sq as u64, (sq >> 64) as u64, 0, 0, 0, 0];
                    if iter % 8 == 6 {
                        n4[0] = n4[0].wrapping_add(1);
                        n6[0] = n6[0].wrapping_add(1);
                    }
                }
                3 => {
                    // Wide perfect squares and their +-1 / +2s neighbors:
                    // the exact-boundary shapes that force the post-Newton
                    // fixup walks in both directions.
                    let r4 = [rng.next(), rng.next() >> 33, 0, 0]; // ~95-bit root
                    let r6 = [rng.next(), rng.next(), rng.next() >> 34, 0, 0, 0]; // ~158-bit
                    let mut sq4 = [0u64; 4];
                    mul_words(&mut sq4, &r4[..2], &r4[..2]);
                    let mut sq6 = [0u64; 6];
                    mul_words(&mut sq6, &r6[..3], &r6[..3]);
                    n4 = sq4;
                    n6 = sq6;
                    match iter % 3 {
                        0 => {
                            // n = root^2 - 1 (previous root, maximal rem).
                            let b4 = n4[0] == 0;
                            n4[0] = n4[0].wrapping_sub(1);
                            if b4 {
                                n4[1] = n4[1].wrapping_sub(1);
                            }
                            let b6 = n6[0] == 0;
                            n6[0] = n6[0].wrapping_sub(1);
                            if b6 {
                                n6[1] = n6[1].wrapping_sub(1);
                            }
                        }
                        1 => {
                            // n = root^2 + 1 (minimal non-zero rem).
                            n4[0] = n4[0].wrapping_add(1);
                            n6[0] = n6[0].wrapping_add(1);
                        }
                        _ => {}
                    }
                }
                _ => {}
            }
            check::<4, 2, 4>(&n4);
            check::<6, 3, 6>(&n6);
        }
    }

    #[test]
    fn test_div_knuth_hard_path_fixtures() {
        // Add-back (D6) with an already-normalized divisor (s == 0).
        let u = [
            9223372036854775808,
            8792676872254834457,
            18446744073709551615,
            9223372036854775808,
        ];
        let v = [
            14320622906001072715,
            9223372036854775808,
            18446744073709551615,
            9223372036854775808,
        ];
        assert!(check_div_knuth(&u, &v).addback, "add-back, s == 0");

        // Add-back (D6) behind a normalization shift (s > 0).
        let u = [
            15833662907273503877,
            1193572382016947486,
            9223372036854775808,
            18446744073709551615,
            9223372036854775807,
        ];
        let v = [
            10296164541181878266,
            10908281225356688296,
            18446744073709551615,
            9223372036854775807,
        ];
        assert!(check_div_knuth(&u, &v).addback, "add-back, s > 0");

        // Two D3 corrections: Knuth's worst case for an even base.
        let u = [
            10216183818892624758,
            12442129881179520362,
            18446744073709551615,
            9223372036854775808,
        ];
        let v = [
            16550061096035433150,
            18446744073709551615,
            9223372036854775808,
        ];
        assert_eq!(check_div_knuth(&u, &v).max_corrections, 2, "2 corrections");

        // `un[j+n] == vtop` estimate branch: qhat starts at B or B + 1.
        let u = [
            7873722608829012771,
            1077289598778399907,
            18446744073709551615,
        ];
        let v = [18446744073709551615, 18446744073709551615];
        assert!(check_div_knuth(&u, &v).eq_estimate, "qhat >= B estimate");

        // Maximal normalization shift (s == 63, divisor top limb is 1).
        let u = [0x123456789ABCDEF0, u64::MAX, 0xFEDCBA9876543210];
        let v = [u64::MAX, 1];
        check_div_knuth(&u, &v);

        // Dividend smaller than divisor at equal width: q == 0, r == u.
        let u = [7u64, 3];
        let v = [9u64, u64::MAX];
        let mut q = [0xAAAA_AAAA_AAAA_AAAAu64; 3];
        let mut r = [0xAAAA_AAAA_AAAA_AAAAu64; 4];
        div_knuth(&mut q, &mut r, &u, &v);
        assert_eq!(q, [0, 0, 0], "quotient and excess limbs zeroed");
        assert_eq!(r, [7, 3, 0, 0], "remainder and excess limbs zeroed");
    }

    /// Rounding and overflow edges of [`dec_div`] that random inputs almost
    /// never reach: the exact tie (`2 * rem == divisor`) through the
    /// multi-limb compare, the `2 * rem` left-shift carrying out of the
    /// divisor width, round-up carrying out of an all-ones quotient, and the
    /// narrow-tier fall-through on 2-limb quotient overflow.
    #[test]
    fn test_dec_div_rounding_edges() {
        let b = 1u64 << 63;
        // Tie: 7*2^64 / 2*2^64 = 3 rem 2^64, and 2*rem == divisor exactly
        // (multi-limb cmp_words returns Equal). HalfEven sees an odd
        // quotient and rounds up; Down truncates.
        let a2 = [0u64, 7];
        let b2 = [0u64, 2];
        let r = dec_div::<0, 2, 3>(false, &a2, false, &b2, Rounding::HalfEven);
        assert_eq!(r, Some((false, [4, 0])), "tie rounds odd quotient up");
        let r = dec_div::<0, 2, 3>(false, &a2, false, &b2, Rounding::Down);
        assert_eq!(r, Some((false, [3, 0])), "tie truncates under Down");

        // 2 * rem overflows the divisor's limb count: divisor top limb has
        // its top bit set (legal when the divisor is narrower than W), and
        // rem == 2^127 so the doubling shl1 carries out => "round up".
        let a4 = [3u64, 0, 2, 0]; // 3 * b + 2^127 with b = 2^127 + 1
        let b4 = [1u64, b, 0, 0];
        let r = dec_div::<0, 4, 5>(false, &a4, false, &b4, Rounding::HalfUp);
        assert_eq!(r, Some((false, [4, 0, 0, 0])), "shl1 carry means rem > b/2");
        let r = dec_div::<0, 4, 5>(false, &a4, false, &b4, Rounding::Down);
        assert_eq!(r, Some((false, [3, 0, 0, 0])));

        // Round-up carry out of an all-ones quotient: a * 100 = 9 * (2^128 -
        // 1) + 5, so q is all ones and 2 * 5 > 9 rounds up, overflowing W
        // limbs => None (not a wrong wrapped value).
        let a2 = [8116567392432202711u64, 1660206966633859645];
        let b2 = [9u64, 0];
        let r = dec_div::<2, 2, 3>(false, &a2, false, &b2, Rounding::HalfUp);
        assert_eq!(r, None, "round-up past all-ones quotient overflows");
        // (Down also overflows: all-ones already exceeds MAX = 2^127 - 1,
        // caught by the top-bit check instead of the round-up carry.)
        let r = dec_div::<2, 2, 3>(false, &a2, false, &b2, Rounding::Down);
        assert_eq!(r, None);

        // Narrow-tier fall-through: 2-limb operands whose quotient has the
        // top bit set overflow the 2-limb core but fit the 4-limb range.
        let a4 = [u64::MAX, u64::MAX, 0, 0];
        let b4 = [1u64, 0, 0, 0];
        let r = dec_div::<0, 4, 5>(false, &a4, false, &b4, Rounding::Down);
        assert_eq!(r, Some((false, [u64::MAX, u64::MAX, 0, 0])));
    }

    /// Overflow edges of [`dec_mul`]: the re-scaled product exceeding W
    /// limbs, non-HalfUp round-up carrying out of an all-ones magnitude, and
    /// the single-limb fast path falling through to the wide core when the
    /// 2-limb result overflows (only possible for DIGITS == 0).
    #[test]
    fn test_dec_mul_overflow_edges() {
        // 2^126 * 4 == 2^128: the product needs a third limb => None.
        let a2 = [0u64, 1 << 62];
        let b2 = [4u64, 0];
        let r = dec_mul::<0, 2, 4>(false, &a2, false, &b2, Rounding::Down);
        assert_eq!(r, None, "re-scaled product exceeds W limbs");

        // 155 * b == 10 * (2^64 - 1) + 5: quotient all ones, exact tie, odd,
        // so HalfEven rounds up and the increment carries out => None.
        let b1 = [(u64::MAX as u128 * 2 + 1).div_euclid(31) as u64]; // (2^65-1)/31
        let r = dec_mul::<1, 1, 2>(false, &[155], false, &b1, Rounding::HalfEven);
        assert_eq!(r, None, "round-up past all-ones product overflows");
        // (Down also overflows: all-ones already exceeds 2^63 - 1, caught by
        // the top-bit check instead of the round-up carry.)
        let r = dec_mul::<1, 1, 2>(false, &[155], false, &b1, Rounding::Down);
        assert_eq!(r, None);

        // Single-limb fast path fall-through: (2^64-1)^2 overflows the
        // 2-limb symmetric range but fits 4 limbs.
        let a4 = [u64::MAX, 0, 0, 0];
        let r = dec_mul::<0, 4, 8>(false, &a4, false, &a4, Rounding::Down);
        assert_eq!(r, Some((false, [1, u64::MAX - 1, 0, 0])), "(2^64-1)^2");
    }

    #[test]
    fn test_div_rem_u128_pow10_edges() {
        // DIGITS == 0 is the identity.
        assert_eq!(div_rem_u128_pow10::<0>(12345), (12345, 0));
        assert_eq!(div_rem_u128_pow10::<0>(u128::MAX), (u128::MAX, 0));
        // High half >= 10^DIGITS exercises the wide first reduction.
        let n = (10_000_000u128 << 64) + 123_456_789;
        assert_eq!(
            div_rem_u128_pow10::<4>(n),
            (n / 10_000, (n % 10_000) as u64)
        );
        assert_eq!(
            div_rem_u128_pow10::<19>(u128::MAX),
            (
                u128::MAX / POW10_19 as u128,
                (u128::MAX % POW10_19 as u128) as u64,
            )
        );
    }

    #[test]
    fn test_div_knuth_exact_and_edges() {
        // 2^192 / 2^128 == 2^64
        let u = [0u64, 0, 0, 1];
        let v = [0u64, 0, 1];
        let mut q = [0u64; 2];
        let mut r = [0u64; 3];
        div_knuth(&mut q, &mut r, &u, &v);
        assert_eq!(q, [0, 1]);
        assert!(is_zero(&r));

        // Divisor with top bit set (s == 0 path).
        let u = [5u64, 7, 9];
        let v = [1u64, 1 << 63];
        let mut q = [0u64; 2];
        let mut r = [0u64; 2];
        div_knuth(&mut q, &mut r, &u, &v);
        let mut prod = [0u64; 4];
        mul_words(&mut prod, &q, &v);
        let mut carry: u128 = 0;
        for i in 0..4 {
            let add = if i < 2 { r[i] as u128 } else { 0 };
            let s = prod[i] as u128 + add + carry;
            prod[i] = s as u64;
            carry = s >> 64;
        }
        assert_eq!(&prod[..3], &u);
        assert_eq!(prod[3], 0);
    }

    #[test]
    fn test_parse_mag() {
        use crate::Rounding::*;
        assert_eq!(
            parse_decimal_mag_rounded::<4>("1.0001", 4, HalfUp),
            Ok((false, [10001, 0, 0, 0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>("-1.00005", 4, HalfUp),
            Ok((true, [10001, 0, 0, 0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>("1.00005", 4, HalfEven),
            Ok((false, [10000, 0, 0, 0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>("", 4, HalfUp),
            Err(AmountErrorKind::Empty)
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>("1.2.3", 4, HalfUp),
            Err(AmountErrorKind::InvalidDigit)
        );
        // A magnitude needing more than one limb:
        // 2^128 == 340282366920938463463374607431768211456
        let (neg, mag) =
            parse_decimal_mag_rounded::<4>("340282366920938463463374607431768211456", 0, HalfUp)
                .unwrap();
        assert!(!neg);
        assert_eq!(mag, [0, 0, 1, 0]);
        // Too many digits for 256 bits.
        let huge = "9".repeat(80);
        assert_eq!(
            parse_decimal_mag_rounded::<4>(&huge, 0, HalfUp),
            Err(AmountErrorKind::Overflow)
        );
    }

    /// The five-mode rounding decision, over every `(round digit, sticky,
    /// parity)` the parser can see.
    ///
    /// The expectation is derived from the exact remainder — `2 * rem` against
    /// the divisor, the crate's magnitude rounding rule — not from the parser's
    /// own first-dropped-digit-plus-sticky formulation, so agreeing with it is
    /// evidence rather than a restatement.
    #[test]
    fn test_parse_rounding_matrix() {
        use crate::Rounding::*;
        // `1.2{k}{r}{s}` at scale 2 keeps `12k` and drops the two-digit
        // remainder `rs` out of a divisor of 100. `k` selects the parity of
        // the kept magnitude, which is the only thing `HalfEven` adds.
        for k in *b"45" {
            let kept = 120 + (k - b'0') as u64;
            for r in 0..10u8 {
                for s in 0..2u8 {
                    let rem = (r as u64) * 10 + s as u64;
                    let src = [b'1', b'.', b'2', k, b'0' + r, b'0' + s];
                    let src = core::str::from_utf8(&src).unwrap();
                    let want = |up: bool| Ok((false, [kept + up as u64]));
                    let got = |mode| parse_decimal_mag_rounded::<1>(src, 2, mode);
                    assert_eq!(got(Down), want(false), "{src}");
                    assert_eq!(got(Up), want(rem != 0), "{src}");
                    assert_eq!(got(HalfUp), want(2 * rem >= 100), "{src}");
                    assert_eq!(got(HalfDown), want(2 * rem > 100), "{src}");
                    let even = 2 * rem > 100 || (2 * rem == 100 && kept % 2 == 1);
                    assert_eq!(got(HalfEven), want(even), "{src}");
                }
            }
        }
    }

    /// Rounding that carries out of the magnitude.
    ///
    /// The final `mul_add_word(&mut acc, 1, 1)` is the only overflow that
    /// happens *after* every digit is accumulated: the digits are in range and
    /// the increment is not, so the same string parses under `Down`.
    #[test]
    fn test_parse_round_carry_overflow() {
        use crate::Rounding::{Down, HalfUp};
        // 2^64 - 1 rounds up to 2^64: out of range at one limb, exact at two.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("18446744073709551615.5", 0, HalfUp),
            Err(AmountErrorKind::Overflow)
        );
        assert_eq!(
            parse_decimal_mag_rounded::<2>("18446744073709551615.5", 0, HalfUp),
            Ok((false, [0, 1]))
        );
        // The value itself fits — only the carry does not.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("18446744073709551615.5", 0, Down),
            Ok((false, [u64::MAX]))
        );
        // One less: the carry stops inside the limb.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("18446744073709551614.5", 0, HalfUp),
            Ok((false, [u64::MAX]))
        );
        // Same boundary two and four limbs up.
        assert_eq!(
            parse_decimal_mag_rounded::<2>("340282366920938463463374607431768211455.5", 0, HalfUp),
            Err(AmountErrorKind::Overflow)
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>("340282366920938463463374607431768211455.5", 0, HalfUp),
            Ok((false, [0, 0, 1, 0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<4>(
                "115792089237316195423570985008687907853269984665640564039457584007913129639935.5",
                0,
                HalfUp
            ),
            Err(AmountErrorKind::Overflow)
        );
    }

    /// Inputs with no integer digits at all.
    ///
    /// A leading '.' hands the cold tail a scan state in which nothing has
    /// been seen yet, so "was there a digit" and "where is the round digit"
    /// both have to survive the hand-off with the integer counter at zero.
    #[test]
    fn test_parse_no_integer_digits() {
        use crate::Rounding::*;
        let p = |s, scale, mode| parse_decimal_mag_rounded::<1>(s, scale, mode);
        // Every mode against a bare half, with a kept magnitude of zero
        // (even, so `HalfEven` goes down).
        assert_eq!(p(".5", 0, HalfUp), Ok((false, [1])));
        assert_eq!(p(".5", 0, HalfDown), Ok((false, [0])));
        assert_eq!(p(".5", 0, HalfEven), Ok((false, [0])));
        assert_eq!(p(".5", 0, Down), Ok((false, [0])));
        assert_eq!(p(".5", 0, Up), Ok((false, [1])));
        assert_eq!(p(".4", 0, HalfUp), Ok((false, [0])));
        assert_eq!(p("-.5", 0, HalfUp), Ok((true, [1])));
        // Over-scale, still with no integer digits.
        assert_eq!(p(".00005", 4, HalfUp), Ok((false, [1])));
        assert_eq!(p(".000049", 4, HalfUp), Ok((false, [0])));
        assert_eq!(p(".00005", 4, HalfEven), Ok((false, [0])));
        // Within scale: the same shape on the inline finish instead.
        assert_eq!(p(".5", 4, HalfUp), Ok((false, [5000])));
        assert_eq!(p("-.5", 4, HalfUp), Ok((true, [5000])));
        // A point and nothing else is empty, not a zero.
        assert_eq!(p(".", 4, HalfUp), Err(AmountErrorKind::Empty));
        assert_eq!(p("-.", 4, HalfUp), Err(AmountErrorKind::Empty));
        assert_eq!(p(".", 0, HalfUp), Err(AmountErrorKind::Empty));
        // Strays reached through the fraction phase and through the tail.
        assert_eq!(p(".x", 0, HalfUp), Err(AmountErrorKind::InvalidDigit));
        assert_eq!(p(".x", 4, HalfUp), Err(AmountErrorKind::InvalidDigit));
        assert_eq!(p(".5x", 0, HalfUp), Err(AmountErrorKind::InvalidDigit));
        assert_eq!(p("..5", 4, HalfUp), Err(AmountErrorKind::InvalidDigit));
        assert_eq!(p(".5.5", 4, HalfUp), Err(AmountErrorKind::InvalidDigit));
    }

    /// The 19-digit chunk boundary where it meets the scale boundary.
    ///
    /// The flush is where `Overflow` is detected, so its position relative to
    /// the fractional budget is observable: a chunk that fills on the very
    /// digit that exhausts `scale` must flush on that digit, ahead of any
    /// stray byte still to come.
    #[test]
    fn test_parse_chunk_boundary() {
        use crate::Rounding::HalfUp;
        // 15 integer + 4 fractional digits: the chunk fills on the last
        // fractional digit, which is also the one that exhausts `scale`.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("123456789012345.6789", 4, HalfUp),
            Ok((false, [1234567890123456789]))
        );
        // Same string, stray byte after it: the value is in range, so the
        // stray wins and the error is `InvalidDigit`.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("123456789012345.6789x", 4, HalfUp),
            Err(AmountErrorKind::InvalidDigit)
        );
        // Digits past the boundary round against the already-flushed value.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("123456789012345.67891", 4, HalfUp),
            Ok((false, [1234567890123456789]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("123456789012345.67895", 4, HalfUp),
            Ok((false, [1234567890123456790]))
        );
        // 19 kept digits but fewer than `scale` fractional ones: the tail has
        // to append the missing zeros to an accumulator that is not zero.
        // One digit short of `scale` is the boundary of that padding —
        // `point == scale` exactly, the only case a `<` there would drop.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("1234567890123456.789", 4, HalfUp),
            Ok((false, [12345678901234567890]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<2>("12345678901234567.89", 4, HalfUp),
            Ok((false, [12776324570088369204, 6]))
        );
        // 18 digits take the inline finish, 19 take the flush; both exact.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("999999999999999999", 0, HalfUp),
            Ok((false, [999999999999999999]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("9999999999999999999", 0, HalfUp),
            Ok((false, [9999999999999999999]))
        );
        // Leading zeros spend the digit budget without changing the value,
        // so the flush path has to stay exact across them.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("00000000000000000000000001", 0, HalfUp),
            Ok((false, [1]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("0.00000000000000000000000001", 4, HalfUp),
            Ok((false, [0]))
        );
    }

    /// The inline finish: one 128-bit product standing in for the partial-chunk
    /// flush and the trailing-zero scaling, taken when nothing was flushed.
    #[test]
    fn test_parse_inline_finish_limits() {
        use crate::Rounding::HalfUp;
        // 18 kept digits scaled by 10^2, either side of 2^64: at one limb the
        // high half of the product is the overflow test, at two it is limb 1.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("1844674407370955.16", 4, HalfUp),
            Ok((false, [18446744073709551600]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("1844674407370955.17", 4, HalfUp),
            Err(AmountErrorKind::Overflow)
        );
        assert_eq!(
            parse_decimal_mag_rounded::<2>("1844674407370955.17", 4, HalfUp),
            Ok((false, [84, 1]))
        );
        // The widest product the inline finish can form: 18 digits by 10^19,
        // just under 10^38 — the bound that keeps it inside one u128.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("999999999999999999", 19, HalfUp),
            Err(AmountErrorKind::Overflow)
        );
        assert_eq!(
            parse_decimal_mag_rounded::<2>("999999999999999999", 19, HalfUp),
            Ok((false, [8515484028849618944, 542101086242752216]))
        );
        // Zero at both ends of the scale range, and a signed zero.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("0", 0, HalfUp),
            Ok((false, [0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("0", 19, HalfUp),
            Ok((false, [0]))
        );
        assert_eq!(
            parse_decimal_mag_rounded::<1>("-0.0", 4, HalfUp),
            Ok((true, [0]))
        );
        // A trailing point is a trailing zero.
        assert_eq!(
            parse_decimal_mag_rounded::<1>("7.", 4, HalfUp),
            Ok((false, [70000]))
        );
    }

    /// Randomised cross-check against an oracle that works on the whole digit
    /// string at once: form the integer of every digit, divide by the power of
    /// ten that places the point at `scale`, and round by comparing `2 * rem`
    /// to the divisor. It shares no structure with the incremental scan, so it
    /// covers digit placement, scaling and rounding together.
    #[test]
    fn test_parse_random_vs_exact_reference() {
        use crate::Rounding::*;
        let modes = [HalfUp, HalfDown, HalfEven, Down, Up];
        let mut rng = Rng(0x0BADC0DE5EED1234);
        for _ in 0..20000 {
            // Up to 18 digits and scale up to 19 keeps the oracle's product
            // (10^18 * 10^19) inside a u128.
            let m = 1 + (rng.next() % 18) as usize;
            let ip = (rng.next() % (m as u64 + 1)) as usize;
            let scale = (rng.next() % 20) as u8;
            let neg = rng.next() & 1 != 0;
            let mode = modes[(rng.next() % 5) as usize];

            let mut buf = [0u8; 24];
            let mut len = 0;
            if neg {
                buf[len] = b'-';
                len += 1;
            }
            let mut d: u128 = 0;
            for i in 0..m {
                if i == ip {
                    buf[len] = b'.';
                    len += 1;
                }
                let digit = (rng.next() % 10) as u8;
                d = d * 10 + digit as u128;
                buf[len] = b'0' + digit;
                len += 1;
            }
            let src = core::str::from_utf8(&buf[..len]).unwrap();

            let frac = (m - ip) as u32;
            let (mag, rem, div) = if frac <= scale as u32 {
                (d * upow10(scale as u32 - frac) as u128, 0u128, 1u128)
            } else {
                let div = upow10(frac - scale as u32) as u128;
                (d / div, d % div, div)
            };
            let up = match mode {
                Down => false,
                Up => rem != 0,
                HalfUp => 2 * rem >= div,
                HalfDown => 2 * rem > div,
                HalfEven => 2 * rem > div || (2 * rem == div && mag & 1 == 1),
            };
            let want = mag + up as u128;
            assert_eq!(
                parse_decimal_mag_rounded::<2>(src, scale, mode),
                Ok((neg, [want as u64, (want >> 64) as u64])),
                "{src} at scale {scale} with {mode:?}"
            );
        }
    }

    #[test]
    fn test_str_mag() {
        extern crate std;
        let mut buf = [0u8; 128];
        let m = |v: u64| [v, 0, 0, 0];
        assert_eq!(
            str_mag(&m(10000), false, 4, None, AmountSign::Negative, &mut buf),
            Some("1")
        );
        assert_eq!(
            str_mag(&m(10001), false, 4, None, AmountSign::Negative, &mut buf),
            Some("1.0001")
        );
        assert_eq!(
            str_mag(&m(10001), true, 4, None, AmountSign::Negative, &mut buf),
            Some("-1.0001")
        );
        assert_eq!(
            str_mag(&m(10050), false, 4, Some(2), AmountSign::Negative, &mut buf),
            Some("1.01")
        );
        assert_eq!(
            str_mag(&m(10000), false, 4, Some(5), AmountSign::Negative, &mut buf),
            Some("1.00000")
        );
        assert_eq!(
            str_mag(&m(10000), false, 4, None, AmountSign::Always, &mut buf),
            Some("+1")
        );
        assert_eq!(
            str_mag(&m(0), false, 4, None, AmountSign::Always, &mut buf),
            Some("0")
        );
        // 2^128 with scale 4: 34028236692093846346337460743176821.1456
        assert_eq!(
            str_mag(
                &[0, 0, 1, 0],
                false,
                4,
                None,
                AmountSign::Negative,
                &mut buf
            ),
            Some("34028236692093846346337460743176821.1456")
        );
        // 20-digit single-limb magnitude: the direct path, no chunk division.
        assert_eq!(
            str_mag(&m(u64::MAX), false, 4, None, AmountSign::Negative, &mut buf),
            Some("1844674407370955.1615")
        );
        // 2^64, the smallest two-limb magnitude: one 19-digit chunk plus a
        // single leading digit.
        assert_eq!(
            str_mag(
                &[0, 1, 0, 0],
                false,
                0,
                None,
                AmountSign::Negative,
                &mut buf
            ),
            Some("18446744073709551616")
        );
        // 10^20 = [7766279631452241920, 5]: the low chunk is all zeros, so
        // the chunk path must still emit its 19 positional zero digits.
        assert_eq!(
            str_mag(
                &[7766279631452241920, 5, 0, 0],
                false,
                0,
                None,
                AmountSign::Negative,
                &mut buf
            ),
            Some("100000000000000000000")
        );
        // Round-trip with the parser.
        let s = "-1234567890123456789012345678901234567890123456789012345.67891";
        let (neg, mag) = parse_decimal_mag_rounded::<4>(s, 5, crate::Rounding::HalfUp).unwrap();
        assert_eq!(
            str_mag(&mag, neg, 5, None, AmountSign::Negative, &mut buf),
            Some(s)
        );
        // Buffers under the fast-path minimum fall back to the general path
        // with identical output, and still reject values that do not fit.
        let mut small = [0u8; 8];
        assert_eq!(
            str_mag(&m(10001), true, 4, None, AmountSign::Negative, &mut small),
            Some("-1.0001")
        );
        assert_eq!(
            str_mag(
                &m(u64::MAX),
                false,
                4,
                None,
                AmountSign::Negative,
                &mut small
            ),
            None
        );
    }

    #[test]
    fn test_isqrt_cores() {
        // Both custom cores must agree exactly with core's u128::isqrt.
        let check = |n: u128| {
            let want = n.isqrt() as u64;
            assert_eq!(isqrt_f64(n), want, "isqrt_f64({n})");
            assert_eq!(isqrt_newton(n), want, "isqrt_newton({n})");
        };
        for n in 0..2000u128 {
            check(n);
        }
        // Perfect squares and their neighbours across the whole width.
        for bit in 1..63u32 {
            let k = (1u128 << bit) - 1;
            for s in [k * k, k * k + 1, k * k + 2 * k] {
                check(s);
            }
        }
        // The largest value any backing feeds in: i64::MAX * 10^19.
        let top = i64::MAX as u128 * POW10_19 as u128;
        for n in [top, top - 1, 1 << 126, (1 << 126) + 1, (1 << 127) - 1] {
            check(n);
        }
        // Random values over mixed magnitudes.
        let mut rng = Rng(0x5851F42D4C957F2D);
        for i in 0..20000u32 {
            let bits = 1 + (i % 127);
            let top = 1u128 << (bits - 1);
            let n = top | (((rng.next() as u128) << 64 | rng.next() as u128) & (top - 1));
            check(n);
        }
    }
}