use super::{
    Rng, bpsw, carmichael, continued_fraction, crt, crt2, extended_gcd, factorize, gcd, iroot,
    is_lucas_probable_prime, is_prime, is_prime_certified, is_strong_probable_prime, isqrt, jacobi,
    kronecker, lcm, modinv, modpow, mulmod, next_prime, order, prime_certificate, primitive_root,
    random_prime, sqrtmod, totient, verify_certificate,
};

// ── Implied coverage — algebraic laws entailed by the defining invariants ──
//
// Each function's *defining* invariant is asserted over its full domain: dense
// grids + 100k–200k LCG inputs in the `*_wide` tests, unbounded coverage-guided
// input in `coprime/fuzz`, and — for the base functions — symbolically in
// `coprime/proofs`. A law that is a mathematical *consequence* of such an
// invariant cannot fail while the invariant holds, so it is recorded here as
// covered-by-implication rather than re-asserted (a test that cannot fail when
// the suite is green tests a theorem, not the code). Two families below are
// NOT consequences and therefore get their own checks — flagged inline.
//
// gcd — entailed by `gcd == the true greatest common divisor` (verified vs
// `brute_gcd` over the grid and vs the independent Stein `binary_gcd` over
// unbounded fuzz):
//   • commutativity  gcd(a,b)==gcd(b,a)  — brute_gcd is symmetric and the
//       0..200² grid covers both orders (also asserted directly in fuzz)
//   • associativity  gcd(a,gcd(b,c))==gcd(gcd(a,b),c)
//   • homogeneity    gcd(k·a,k·b)==k·gcd(a,b)
//   • Euclid step    gcd(a,b)==gcd(a, b mod a)
//   • idempotence/unit/zero  gcd(a,a)==a, gcd(a,1)==1, gcd(a,0)==a
//   each a theorem of the true gcd, which is pinned by two independent
//   algorithms over the full domain. (gcd·lcm==a·b waits on an `lcm` fn.)
//
// extended_gcd — Bézout `s·a+t·b==d` and `d∣a,b` are verified
// (extended_gcd_bezout_wide + Kani). The cofactor *magnitude* bounds
// |s|≤b/2, |t|≤a/2 are NOT entailed (Bézout does not pin the magnitude) and the
// no-overflow safety rests on them, so they are asserted directly in
// extended_gcd_cofactor_bounds, NOT covered here. Sign-opposition of s·a,t·b
// and swap-symmetry extgcd(a,b)=(d,s,t) ⇒ extgcd(b,a)=(d,t,s) follow once
// minimality is checked.
//
// modinv — entailed by `a·inv(a) ≡ 1 (mod m)` over the full domain
// (modinv_inverts_wide + fuzz) plus uniqueness of inverses:
//   • involution     inv(inv(a))==a
//   • homomorphism   inv(a·b)==inv(a)·inv(b)
//   • self-inverse units  inv(1)==1, inv(m-1)==m-1
//   • invertible ⇔ coprime — the fuzz None-arm asserts None ⇒ gcd≠1 (so
//       gcd==1 ⇒ Some), and the Some-arm's a·inv≡1 forces gcd==1
//   • Euler form inv(a)==a^(φ(m)-1) — asserted directly in euler_carmichael_laws
//
// modpow — entailed by `modpow == a^e mod m` (verified vs naive, L2R, and an
// addition-chain mulmod over grid + unbounded fuzz):
//   • additivity a^(x+y)==a^x·a^y   • exp-mult a^(x·y)==(a^x)^y
//   • base-mult (a·b)^x==a^x·b^x     • identities a^1==a, 1^e==1
//   • Fermat a^p==a incl. p∣a        • Euler a^φ(m)==1 (in euler_carmichael_laws)
//
// is_strong_probable_prime — base reduction spp(n,a)==spp(n, a mod n) holds by
// construction (the body reduces `a` first) and spp ⇒ Fermat-PRP is a theorem.
// BUT correctness over *arbitrary* (n,a) — large n with a base outside the
// fixed set {2..37} — rests only on `modpow correct + a simple wrapper`, never
// an independent SPRP oracle; that quadrant gets its own differential checks
// (spp_matches_independent_oracle for completeness, prime⇒liar; and
// spp_soundness_vs_reference_mr for soundness, witness⇒false over composite n
// × arbitrary base / the fuzz target), NOT covered here.
//
// is_prime — entailed:
//   • no false negative — a prime is never an MR witness (theorem); two-sided
//       is_prime==trial_is_prime is fuzzed for n<2³² and Feitsma rules out
//       false positives over all u64
//   • is_prime == AND of its 12 rounds — true by construction
//   • prefilter soundness — multiples of 2..37 are composite, caught by the
//       sieve (0..10⁶) and the composite-rejection fuzz
//
// crt / crt2 — entailed by `result satisfies every congruence and lies in
// [0,∏m)` (crt2_matches_brute over every residue pair of coprime moduli;
// congruence+range over unbounded fuzz) plus CRT uniqueness (coprime ⇒
// Z/∏m ≅ ∏ Z/mᵢ is a bijection):
//   • commutativity crt2(r1,m1,r2,m2)==crt2(r2,m2,r1,m1) for moduli ≥2
//   • uniqueness — exactly one solution in range, so the value found IS it
//   • surjectivity — crt2_matches_brute sweeps every residue pair
//   • input reduction crt2(r1,m1,…)==crt2(r1 mod m1,m1,…) (reduced in-body)
//   • k-ary order-independence crt(permuted)==crt(original) — entailed (the
//       result meets all congruences + range, hence is THE unique solution
//       regardless of fold), and also asserted directly over shuffled
//       coprime-prime subsets in crt_satisfies_congruences_wide
//   • fold-consistency crt==iterated crt2; solvable ⇔ pairwise coprime (fuzz
//       None-arm); identities crt(&[])==Some(0), crt(&[(r,m)])==Some(r mod m)
//
// next_prime / random_prime — entailed by `is_prime` being exhaustively correct
// (Feitsma) over a trivial scan:
//   • next_prime(prime)==Some(prime); least-prime / gap-is-composite at any
//       magnitude — the scan is trivial given a correct is_prime
//   • random_prime: p∈[2^(bits-1),2^bits) ≡ the bit-length check; determinism
//       (same seed ⇒ same prime) is shown by the rustdoc example. coprime's own
//       candidate construction `(rng & mask) | top | 1` is pinned directly in
//       random_prime_reaches_every_band_prime (every odd band prime reachable,
//       in-band, odd, exact length). Only the statistical *distribution*
//       (uniformity/bias of the word stream) is the injected RNG's
//       responsibility — out of scope by mechanism-not-policy.

/// Brute-force greatest common divisor (the largest `d` dividing both) — an
/// oracle independent of Euclid's algorithm.
fn brute_gcd(a: u64, b: u64) -> u64 {
    if a == 0 && b == 0 {
        return 0;
    }
    let mut g = 0;
    let hi = a.max(b);
    let mut d = 1;
    while d <= hi {
        if a.is_multiple_of(d) && b.is_multiple_of(d) {
            g = d;
        }
        d += 1;
    }
    g
}

/// PROOF: `gcd` equals the brute-force greatest common divisor and divides
/// both inputs, exhaustively over 0..200 × 0..200; `gcd(n,0)=n`, `gcd(0,0)=0`.
#[test]
fn gcd_is_greatest_common_divisor() {
    assert_eq!(gcd(0, 0), 0);
    for a in 0..200u64 {
        assert_eq!(gcd(a, 0), a, "gcd({a},0)");
        for b in 0..200u64 {
            let d = gcd(a, b);
            assert_eq!(d, brute_gcd(a, b), "a={a} b={b}");
            if d != 0 {
                assert_eq!(a % d, 0, "a={a} b={b}: d does not divide a");
                assert_eq!(b % d, 0, "a={a} b={b}: d does not divide b");
            }
        }
    }
}

/// PROOF of the Bézout identity: `extended_gcd(a,b)` returns `(d,s,t)` with
/// `d == gcd(a,b)` and `s*a + t*b == d`, exhaustively over 0..150 × 0..150,
/// plus large inputs (>= 2^63) that would silently wrap under an `i64`
/// accumulator. (The all-`u64` symbolic version is the Kani proof.)
#[test]
fn extended_gcd_satisfies_bezout() {
    for a in 0..150u64 {
        for b in 0..150u64 {
            let (d, s, t) = extended_gcd(a, b);
            assert_eq!(d, gcd(a, b), "a={a} b={b} gcd");
            assert_eq!(
                s * i128::from(a) + t * i128::from(b),
                i128::from(d),
                "a={a} b={b} bezout"
            );
        }
    }
    // >= 2^63: exact under i128, silently wrong under an i64 accumulator.
    assert_eq!(extended_gcd(1u64 << 63, 1), (1, 0, 1));
    let (d, s, t) = extended_gcd(u64::MAX, u64::MAX - 1);
    assert_eq!(d, 1, "gcd of consecutive integers");
    assert_eq!(s * i128::from(u64::MAX) + t * i128::from(u64::MAX - 1), 1);
}

/// Naive modular exponentiation: `exp` repeated modular multiplications, no
/// squaring — an oracle independent of the square-and-multiply bit walk.
#[allow(clippy::cast_possible_truncation)] // acc < modulus <= u64::MAX fits u64
fn naive_modpow(base: u64, exp: u64, modulus: u64) -> u64 {
    if modulus <= 1 {
        return 0;
    }
    let m = u128::from(modulus);
    let b = u128::from(base) % m;
    let mut acc: u128 = 1;
    for _ in 0..exp {
        acc = acc * b % m;
    }
    acc as u64
}

/// PROOF: `modpow` equals naive repeated-multiplication exponentiation,
/// exhaustively over 0..40 base × 0..40 exp × 1..40 modulus, plus the
/// `modulus <= 1` convention and the `0^0 = 1` / `x^0 = 1` / `0^e = 0` edges.
#[test]
fn modpow_matches_naive() {
    assert_eq!(modpow(7, 13, 0), 0, "mod 0 convention");
    assert_eq!(modpow(7, 13, 1), 0, "mod 1");
    assert_eq!(modpow(0, 0, 7), 1, "0^0 = 1 by convention");
    assert_eq!(modpow(0, 5, 7), 0, "0^5 = 0");
    assert_eq!(modpow(123_456_789, 0, 1_000_000_007), 1, "x^0 = 1");
    for base in 0..40u64 {
        for exp in 0..40u64 {
            for modulus in 1..40u64 {
                assert_eq!(
                    modpow(base, exp, modulus),
                    naive_modpow(base, exp, modulus),
                    "base={base} exp={exp} mod={modulus}",
                );
            }
        }
    }
}

/// PROOF via Fermat's little theorem: for prime `p` and `a` not divisible by
/// `p`, `a^(p-1) == 1 (mod p)`. The exponents here reach ~2^31 — far beyond
/// what the linear oracle could check — exercising the square-and-multiply
/// fast path and the wide `u128` products under a near-`2^31` modulus.
#[test]
fn modpow_fermat_little_theorem() {
    // each value is prime: 7919 is the 1000th prime, 1_000_000_007 (10^9+7)
    // is the well-known competitive-programming prime, 2_147_483_647 is the
    // Mersenne prime 2^31 - 1.
    for &p in &[2u64, 3, 5, 7, 13, 101, 7919, 1_000_000_007, 2_147_483_647] {
        for a in 1..50u64 {
            if a % p != 0 {
                assert_eq!(modpow(a, p - 1, p), 1, "fermat a={a} p={p}");
            }
        }
    }
}

/// `mulmod`'s public contract, including the degenerate moduli that its
/// internal callers (always `m >= 2`) never reach: `m <= 1` returns `0`, and the
/// `u128` widening is exact at the `u64` ceiling. The `m == 0` arm is a runtime
/// check on the total contract that `mulmod_total_spec` proves in Lean.
#[test]
fn mulmod_contract() {
    assert_eq!(mulmod(7, 8, 5), 1); // 56 mod 5
    assert_eq!(mulmod(0, 0, 7), 0);
    assert_eq!(mulmod(u64::MAX, u64::MAX, 1_000_000_007), 114_944_269); // no overflow
    assert_eq!(mulmod(u64::MAX, u64::MAX, u64::MAX), 0); // exact at the ceiling
    // degenerate moduli: total, no panic — the branch no caller exercises
    assert_eq!(mulmod(7, 8, 0), 0);
    assert_eq!(mulmod(7, 8, 1), 0);
    assert_eq!(mulmod(u64::MAX, u64::MAX, 0), 0);
}

/// Deterministic LCG (Knuth MMIX constants) — full-range pseudo-random test
/// inputs without an external rng dependency.
fn lcg(state: &mut u64) -> u64 {
    *state = state
        .wrapping_mul(6_364_136_223_846_793_005)
        .wrapping_add(1_442_695_040_888_963_407);
    *state
}

/// Left-to-right square-and-multiply: scan exponent bits high to low. An
/// oracle independent of `modpow`'s right-to-left direction.
#[allow(clippy::cast_possible_truncation)] // result < modulus <= u64::MAX fits u64
fn modpow_l2r(base: u64, exp: u64, modulus: u64) -> u64 {
    if modulus <= 1 {
        return 0;
    }
    let m = u128::from(modulus);
    let b = u128::from(base) % m;
    let mut result: u128 = 1;
    let mut bit = u64::BITS;
    while bit > 0 {
        bit -= 1;
        result = result * result % m;
        if (exp >> bit) & 1 == 1 {
            result = result * b % m;
        }
    }
    result as u64
}

/// BREADTH (this is the tier CBMC cannot reach): the Bézout identity over a
/// dense 0..512 grid, every consecutive-Fibonacci worst-case pair up to
/// `u64::MAX` (the inputs that maximize Euclid's iteration count), and 200k
/// full-range pseudo-random `u64` pairs. CBMC bit-blasts the `i128`
/// recurrence to a 51M-clause instance and OOMs past ~2^8 inputs, so the
/// wide domain is covered concretely here instead of symbolically.
#[test]
fn extended_gcd_bezout_wide() {
    let check = |a: u64, b: u64| {
        let (d, s, t) = extended_gcd(a, b);
        // Bézout (s*a + t*b == d) together with d | a and d | b pin d to
        // gcd(a,b) *independently* of the `gcd` function: gcd | d via Bézout,
        // d | gcd via the divisibilities, both non-negative, so d == gcd.
        assert_eq!(
            s * i128::from(a) + t * i128::from(b),
            i128::from(d),
            "bezout a={a} b={b}"
        );
        if d != 0 {
            assert_eq!(a % d, 0, "d does not divide a={a} (b={b})");
            assert_eq!(b % d, 0, "d does not divide b={b} (a={a})");
        } else {
            assert!(
                a == 0 && b == 0,
                "d == 0 only when a == b == 0 (a={a} b={b})"
            );
        }
    };
    for a in 0..512u64 {
        for b in 0..512u64 {
            check(a, b);
        }
    }
    // consecutive-Fibonacci pairs to u64::MAX — Euclid's worst case.
    let (mut x, mut y) = (0u64, 1u64);
    while let Some(z) = x.checked_add(y) {
        check(y, x);
        (x, y) = (y, z);
    }
    check(y, x);
    // 200k full-range pseudo-random pairs.
    let mut state = 0x2545_F491_4F6C_DD1D;
    for _ in 0..200_000 {
        check(lcg(&mut state), lcg(&mut state));
    }
}

/// PROOF (cofactor minimality — the bound the no-overflow safety rests on):
/// `extended_gcd` returns the *reduced* Bézout cofactors `|s| <= max(1, b/2)`
/// and `|t| <= max(1, a/2)`, and for `a, b >= 2` the products `s*a`, `t*b`
/// have opposite signs. The Bézout identity alone does NOT pin the cofactor
/// magnitude — any `(s + k*b/d, t - k*a/d)` also satisfies `s*a + t*b == d` —
/// so this is an *independent* property, and `|s*a|, |t*b| <= a*b/2 < 2^127`
/// is exactly what keeps the `i128` accumulators from overflowing. Same
/// domain as the Bézout breadth test: dense grid + Fibonacci worst cases +
/// 200k full-range pairs.
#[test]
fn extended_gcd_cofactor_bounds() {
    let check = |a: u64, b: u64| {
        let (_, s, t) = extended_gcd(a, b);
        assert!(
            s.unsigned_abs() <= u128::from((b / 2).max(1)),
            "|s| exceeds max(1, b/2) for a={a} b={b}: s={s}"
        );
        assert!(
            t.unsigned_abs() <= u128::from((a / 2).max(1)),
            "|t| exceeds max(1, a/2) for a={a} b={b}: t={t}"
        );
        // The products fit i128 (the safety claim) and oppose in sign once
        // both operands exceed 1 (the gcd is small, so the terms cancel).
        let (sa, tb) = (s * i128::from(a), t * i128::from(b));
        if a >= 2 && b >= 2 && sa != 0 && tb != 0 {
            assert!(
                (sa > 0) != (tb > 0),
                "s*a, t*b share a sign for a={a} b={b}"
            );
        }
    };
    for a in 0..512u64 {
        for b in 0..512u64 {
            check(a, b);
        }
    }
    let (mut x, mut y) = (0u64, 1u64);
    while let Some(z) = x.checked_add(y) {
        check(y, x);
        (x, y) = (y, z);
    }
    check(y, x);
    let mut state = 0x9E37_79B9_7F4A_7C15;
    for _ in 0..200_000 {
        check(lcg(&mut state), lcg(&mut state));
    }
}

/// BREADTH for modpow: cross-check right-to-left `modpow` against the
/// independent left-to-right oracle over a dense grid and 200k full-range
/// pseudo-random `(base, exp, modulus)` triples — the symbolic slice CBMC
/// cannot discharge (two 128-bit `% m` per exponent bit), covered here.
#[test]
fn modpow_differential_wide() {
    for base in 0..64u64 {
        for exp in 0..64u64 {
            for modulus in 1..64u64 {
                assert_eq!(
                    modpow(base, exp, modulus),
                    modpow_l2r(base, exp, modulus),
                    "grid base={base} exp={exp} mod={modulus}",
                );
            }
        }
    }
    let mut state = 0x9E37_79B9_7F4A_7C15;
    for _ in 0..200_000 {
        let (base, exp, modulus) = (lcg(&mut state), lcg(&mut state), lcg(&mut state));
        assert_eq!(
            modpow(base, exp, modulus),
            modpow_l2r(base, exp, modulus),
            "random base={base} exp={exp} mod={modulus}",
        );
    }
}

/// Brute-force modular inverse: scan `[0, modulus)` for `x` with
/// `a*x ≡ 1 (mod modulus)` — an oracle independent of extended Euclid.
fn brute_modinv(a: u64, modulus: u64) -> Option<u64> {
    if modulus <= 1 {
        return None;
    }
    let am = u128::from(a % modulus);
    let m = u128::from(modulus);
    (0..modulus).find(|&x| am * u128::from(x) % m == 1)
}

/// PROOF: `modinv` equals the brute-force inverse over 0..200 × 1..200, plus
/// the degenerate-modulus and non-invertible (`gcd != 1`) edges.
#[test]
fn modinv_matches_brute() {
    assert_eq!(modinv(0, 7), None, "0 has no inverse");
    assert_eq!(modinv(5, 1), None, "degenerate modulus 1");
    assert_eq!(modinv(5, 0), None, "degenerate modulus 0");
    assert_eq!(modinv(2, 4), None, "gcd(2,4) != 1");
    assert_eq!(modinv(3, 7), Some(5), "3*5 = 15 ≡ 1 (mod 7)");
    for modulus in 1..200u64 {
        for a in 0..200u64 {
            assert_eq!(
                modinv(a, modulus),
                brute_modinv(a, modulus),
                "a={a} mod={modulus}",
            );
        }
    }
}

/// PROOF that the result truly inverts: cross-checked against Fermat
/// (`a^{-1} ≡ a^(p-2) mod p` for prime `p`, via the independent `modpow`),
/// and verified as a genuine inverse over 100k full-range pseudo-random
/// `(a, modulus)` pairs.
#[test]
fn modinv_inverts_wide() {
    for &p in &[2u64, 3, 5, 7, 13, 101, 7919, 1_000_000_007] {
        for a in 1..50u64 {
            if a % p != 0 {
                assert_eq!(
                    modinv(a, p),
                    Some(modpow(a, p - 2, p)),
                    "fermat a={a} p={p}"
                );
            }
        }
    }
    let mut state = 0x1234_5678_9abc_def1;
    for _ in 0..100_000 {
        let modulus = lcg(&mut state).max(2);
        let a = lcg(&mut state);
        if let Some(inv) = modinv(a, modulus) {
            assert!(inv < modulus, "inv {inv} out of range mod {modulus}");
            assert_eq!(
                u128::from(a % modulus) * u128::from(inv) % u128::from(modulus),
                1,
                "a={a} mod={modulus} inv={inv}",
            );
        }
    }
}

/// Sieve of Eratosthenes — an independent primality oracle (no Miller–Rabin).
fn sieve(limit: usize) -> Vec<bool> {
    let mut is_p = vec![true; limit];
    for slot in is_p.iter_mut().take(2.min(limit)) {
        *slot = false; // 0 and 1 are not prime
    }
    let mut i = 2;
    while i * i < limit {
        if is_p[i] {
            let mut j = i * i;
            while j < limit {
                is_p[j] = false;
                j += i;
            }
        }
        i += 1;
    }
    is_p
}

/// PROOF: `is_prime` matches the sieve exhaustively over `0..1_000_000`.
#[test]
fn is_prime_matches_sieve() {
    let s = sieve(1_000_000);
    for (n, &prime) in s.iter().enumerate() {
        assert_eq!(is_prime(n as u64), prime, "n={n}");
    }
}

/// PROOF (known answers): large `u64` primes pass; Carmichael numbers (which
/// fool Fermat) and the `ψ_t` strong pseudoprimes (which fool Miller–Rabin
/// run on too few bases) are all rejected — each is composite and < `ψ_12`.
#[test]
fn is_prime_known_answers() {
    for &p in &[
        2u64,
        3,
        5,
        37,
        41,
        7919,
        4_294_967_291, // largest prime < 2^32
        2_147_483_647, // 2^31 - 1 (Mersenne M31)
        1_000_000_007,
        2_305_843_009_213_693_951, // 2^61 - 1 (Mersenne M61)
    ] {
        assert!(is_prime(p), "{p} should be prime");
    }
    for &c in &[
        0u64,
        1,
        4,
        9,
        561,
        1105,
        1729,
        41041,
        825_265,             // Carmichael numbers — fool Fermat, not MR
        2047,                // ψ_1: strong pseudoprime to base 2
        1_373_653,           // ψ_2
        25_326_001,          // ψ_3
        3_215_031_751,       // ψ_4
        2_152_302_898_747,   // ψ_5
        3_474_749_660_383,   // ψ_6
        341_550_071_728_321, // ψ_7
    ] {
        assert!(!is_prime(c), "{c} should be composite");
    }
}

/// PROOF (no false positive on composites): every product of two integers
/// `>= 2` is composite, so `is_prime` must reject all 100k full-range
/// products. Also exercises the single-base test directly.
#[test]
fn is_prime_rejects_random_composites() {
    let mut state = 0xC0FF_EE12_3456_789A;
    for _ in 0..100_000 {
        let a = (lcg(&mut state) >> 32).max(2); // [2, 2^32)
        let b = (lcg(&mut state) >> 32).max(2); // [2, 2^32)
        let n = a * b; // < 2^64, composite
        assert!(!is_prime(n), "false positive on composite {n} = {a}*{b}");
    }
    // a base-2 strong witness must reject a base-2 strong pseudoprime's
    // neighbours; 2047 is a strong *liar* for base 2 but composite (caught by
    // is_prime via other bases) — here just confirm the single-round contract.
    assert!(
        is_strong_probable_prime(2047, 2),
        "2047 is a base-2 strong liar"
    );
    assert!(
        !is_strong_probable_prime(2047, 3),
        "3 is a strong witness for 2047"
    );
}

/// PROOF (single-round Miller–Rabin semantics): the `ψ_t` strong pseudoprimes
/// pin down exactly what one round can and cannot do. `ψ_t` is the *smallest*
/// composite that is a strong pseudoprime to all of the first `t` prime
/// bases, so for each `ψ_t` the first `t` primes are strong *liars* and the
/// `(t+1)`-th prime is a strong *witness*. These values are from the
/// literature (independent of our code) and justify why `is_prime` runs
/// twelve bases.
#[test]
fn is_strong_probable_prime_liars_and_witnesses() {
    const PRIMES: [u64; 6] = [2, 3, 5, 7, 11, 13];
    // (ψ_t, t): the first t primes are liars; prime #(t+1) is a witness.
    let psi: [(u64, usize); 4] = [
        (2047, 1),          // 23 · 89
        (1_373_653, 2),     // 829 · 1657
        (25_326_001, 3),    // 2251 · 11251
        (3_215_031_751, 4), // 151 · 751 · 28351
    ];
    for (n, t) in psi {
        assert!(!is_prime(n), "ψ value {n} must be composite");
        for &a in &PRIMES[..t] {
            assert!(
                is_strong_probable_prime(n, a),
                "base {a} should be a strong liar for ψ={n}",
            );
        }
        let witness = PRIMES[t];
        assert!(
            !is_strong_probable_prime(n, witness),
            "prime {witness} should witness ψ={n} composite",
        );
    }
}

/// PROOF (a prime is a strong probable prime to every base; trivial bases
/// witness nothing): for each small prime `p`, every base `a in 0..2p` is a
/// liar — exhaustive over the base. The bases `a ≡ 0, 1, -1 (mod n)` pass
/// for any `n`, and `n < 2`, even `n > 2` are rejected (with `n = 2, 3` the
/// special cases).
#[test]
fn is_strong_probable_prime_edges_and_primes() {
    for &p in &[3u64, 5, 7, 11, 13, 97, 7919] {
        for a in 0..2 * p {
            assert!(
                is_strong_probable_prime(p, a),
                "prime {p} failed for base {a}"
            );
        }
    }
    for n in [9u64, 15, 21, 91, 561] {
        assert!(is_strong_probable_prime(n, 0), "base 0 must pass n={n}");
        assert!(is_strong_probable_prime(n, 1), "base 1 must pass n={n}");
        assert!(
            is_strong_probable_prime(n, n - 1),
            "base n-1 must pass n={n}"
        );
    }
    assert!(is_strong_probable_prime(2, 5), "n=2 is prime");
    assert!(is_strong_probable_prime(3, 7), "n=3 is prime");
    assert!(!is_strong_probable_prime(1, 5), "n<2 is not prime");
    assert!(!is_strong_probable_prime(0, 5), "n<2 is not prime");
    assert!(!is_strong_probable_prime(100, 7), "even n>2 is composite");
    // Even n > 2 is composite for EVERY base — crucially the trivial bases
    // a ≡ 0, 1, n-1 (mod n), which the MR machinery passes via the
    // `a <= 1 || a == n-1` short-circuit. This pins the `n & 1 == 0`
    // even-check: with it disabled (e.g. mutated to `n | 1`), those bases
    // would wrongly return `true` for an even composite.
    for n in [4u64, 6, 8, 10, 22, 100, 1000] {
        for a in [0u64, 1, 2, n - 1, n + 1, 3 * n + 1] {
            assert!(
                !is_strong_probable_prime(n, a),
                "even {n} must be composite for base {a}"
            );
        }
    }
}

/// PROOF (the round links to the deterministic test): over `3..20_000`, a
/// composite is below `ψ_6` ≈ `3.47e12`, so it cannot be a strong liar to all
/// six of `{2, 3, 5, 7, 11, 13}`. Hence "passes all six bases" ⟺ prime — an
/// independent re-derivation of `is_prime` from the single-round test.
#[test]
fn is_strong_probable_prime_consistent_with_is_prime() {
    const BASES: [u64; 6] = [2, 3, 5, 7, 11, 13];
    for n in 3..20_000u64 {
        let liar_to_all = BASES.iter().all(|&a| is_strong_probable_prime(n, a));
        assert_eq!(liar_to_all, is_prime(n), "n={n}");
    }
}

/// PROOF (BPSW correctness over `u64`): Baillie–PSW agrees with the independent
/// 12-base `is_prime` oracle on *every* input. Dense over `[0, 200_000)` — which
/// covers the first strong Lucas pseudoprimes (5459, 5777, 10877, …) and the
/// base-2 strong pseudoprimes (2047, 3277, 4033, …), the two error sets BPSW
/// composes away — plus 100k large pseudo-random `u64`. (No composite below
/// `2⁶⁴` passes both legs — Baillie–Fiori–Wagstaff 2021 — so this agreement is
/// total, not sampled, over the bounded range and holds beyond it.)
#[test]
fn bpsw_matches_is_prime() {
    for n in 0..200_000u64 {
        assert_eq!(bpsw(n), is_prime(n), "bpsw disagrees at n={n}");
    }
    let mut state = 0xB95E_1234_5678_9ABCu64;
    for _ in 0..100_000 {
        let n = lcg(&mut state);
        assert_eq!(bpsw(n), is_prime(n), "bpsw disagrees at n={n}");
    }
}

/// PROOF (Pocklington certificates are correct and re-checkable): the
/// constructive test agrees with the independent `is_prime` oracle over a dense
/// range — and, decisively, every certificate it *constructs* for a prime is
/// *re-verified* from scratch by `verify_certificate` (modular exponentiation
/// only, no trust in the constructor). An *acceptance* is therefore a genuine
/// proof, independent of `DET12` / Feitsma.
#[test]
fn is_prime_certified_matches_and_round_trips() {
    for n in 0..20_000u64 {
        assert_eq!(
            is_prime_certified(n),
            is_prime(n),
            "certified disagrees at n={n}"
        );
        if let Some(cert) = prime_certificate(n) {
            assert!(
                is_prime(n),
                "constructed a certificate for the composite {n}"
            );
            assert!(
                verify_certificate(n, &cert),
                "certificate for {n} fails re-verification"
            );
        } else {
            assert!(!is_prime(n) || n < 2, "no certificate for the prime {n}");
        }
    }
    // Larger primes: construct and independently re-verify the proof. The last,
    // 2·3·5·…·53 + 1 (the product of fifteen distinct primes, plus one — so n − 1
    // has the maximal `ω = 15` that fills a node's `children: [u8; 15]`),
    // exercises the child-loop bound at full width: an off-by-one there indexes
    // past the array.
    for &p in &[
        1_000_003u64,
        1_000_033,
        32_416_190_071,
        2_305_843_009_213_693_951,
        693_386_350_578_511_591,
    ] {
        let cert = prime_certificate(p).unwrap_or_else(|| panic!("no certificate for prime {p}"));
        assert!(verify_certificate(p, &cert), "round-trip failed for {p}");
    }
    // The tree is deduplicated post-order: each distinct prime appears once, the
    // root last. For 7 (7 − 1 = 2·3) that is exactly [2, 3, 7] — a duplicated
    // atom would betray a broken dedup lookup.
    let cert7 = prime_certificate(7).expect("7 is prime");
    assert_eq!(
        cert7.primes().collect::<Vec<_>>(),
        [2, 3, 7],
        "certificate tree is not deduplicated post-order"
    );
    // A certificate is bound to its n: it must not verify a different modulus.
    let cert = prime_certificate(1_000_003).unwrap_or_else(|| panic!("1_000_003 is prime"));
    assert!(
        !verify_certificate(1_000_005, &cert),
        "certificate verified the wrong n"
    );
}

/// Completeness of the certificate *builder* over a sweep. The `WITNESS_CAP` /
/// `MAX_NODES` headroom is an **empirical** claim (no theorem covers it — the
/// unconditional least-primitive-root bound is far too weak, and enumerating all
/// `u64` primes is infeasible; see PAPERS.md), and this test is the in-repo
/// artifact pinning it: every prime in a dense range gets a certificate, and so
/// does a deterministic sample of large (up to 64-bit) primes. A failure here
/// would be a completeness *miss* (`None` for a true prime — still sound, never
/// a wrong accept) and grounds to revisit the caps.
#[test]
fn certificate_completeness_sweep() {
    // Dense range, extending the round-trip test's [0, 20k).
    for n in 20_000..100_000u64 {
        if is_prime(n) {
            assert!(
                prime_certificate(n).is_some(),
                "completeness miss at prime {n}"
            );
        }
    }
    // Deterministic large sample: the next prime after each LCG draw, spanning
    // the full u64 range (so n − 1 exercises factorize's Pollard tier too).
    let mut x = 0x9E37_79B9_7F4A_7C15u64;
    for _ in 0..200 {
        x = x
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        if let Some(p) = next_prime(x) {
            let cert =
                prime_certificate(p).unwrap_or_else(|| panic!("completeness miss at prime {p}"));
            assert!(verify_certificate(p, &cert), "round-trip failed for {p}");
        }
    }
}

/// The strong Lucas leg on its own: every prime passes, and it independently
/// rejects the base-2 strong pseudoprimes (2047 = 23·89, 3277 = 29·113,
/// 4033 = 37·109) that fool a single base-2 Miller–Rabin round — the disjointness
/// that makes the BPSW composition strong. (A Lucas *pseudoprime* like 5459 the
/// Lucas leg alone accepts; the base-2 round rejects it, which `bpsw_matches`
/// covers.)
#[test]
fn is_lucas_probable_prime_basics() {
    // 2 and 3 pin the lower boundary: they pass through the `n < 2` guard into the
    // small-prime trial division and return true at `n == p`, so a guard widened to
    // `n == 2` or `n <= 2` would wrongly reject them.
    for p in [2u64, 3, 97, 101, 7919, 1_000_003, 2_305_843_009_213_693_951] {
        assert!(is_lucas_probable_prime(p), "prime {p} fails the Lucas test");
    }
    for c in [91u64, 561, 1105, 2047, 3277, 4033] {
        assert!(
            !is_lucas_probable_prime(c),
            "composite {c} passes the Lucas test"
        );
    }
    // is_lucas agrees with is_prime except on Lucas pseudoprimes; over a dense
    // range those are sparse, so any disagreement must be a composite (never a
    // missed prime — the Lucas test has no anti-pseudoprimes).
    for n in 3..20_000u64 {
        if is_lucas_probable_prime(n) != is_prime(n) {
            assert!(!is_prime(n), "Lucas test rejected the prime {n}");
        }
    }
}

/// PROOF (single-round MR over *arbitrary* `(n, a)` — the quadrant the
/// 12-base `is_prime` path never reaches): with an independent sieve as the
/// primality oracle, every prime is a strong probable prime to EVERY base
/// (completeness — the one-directional MR guarantee); soundness, "a witness
/// proves compositeness", is its contrapositive. The bases exercised include
/// the trivial cases `0, 1, n-1`, bases `>= n` (testing the `a mod n`
/// reduction), bases outside the fixed witness set `{2..37}`, and
/// pseudo-random `u64` — none of which `is_prime` ever feeds to a round.
/// Large known primes carry it past the sieve.
#[test]
fn spp_matches_independent_oracle() {
    let bound = 60_000usize;
    let prime = sieve(bound);
    let mut state = 0xA5A5_1234_5678_9ABC;
    for (n, &is_p) in prime.iter().enumerate().skip(3) {
        let nu = n as u64; // usize -> u64 is widening; never truncates
        let bases = [
            0u64,
            1,
            2,
            3,
            37,
            101,
            nu - 1,
            nu,
            nu + 1,
            nu / 2 + 1,
            lcg(&mut state),
        ];
        for &a in &bases {
            let spp = is_strong_probable_prime(nu, a); // exercise the path for all n
            if is_p {
                assert!(
                    spp,
                    "prime {nu} must be a strong probable prime to base {a}"
                );
            }
        }
    }
    // Past the sieve: a large prime is a strong probable prime to every base,
    // including bases far above the fixed witness set.
    for &p in &[
        2_147_483_647u64,
        4_294_967_291,
        1_000_000_007,
        2_305_843_009_213_693_951,
    ] {
        let mut st = p ^ 0xDEAD_BEEF;
        for _ in 0..2_000 {
            let a = lcg(&mut st);
            assert!(
                is_strong_probable_prime(p, a),
                "prime {p} failed SPRP to base {a}"
            );
        }
    }
}

/// Brute-force CRT: scan `[0, m1*m2)` for the solution — an oracle
/// independent of Garner's method.
#[allow(clippy::cast_possible_truncation)] // result < m1*m2 (small here) fits u64
fn brute_crt2(r1: u64, m1: u64, r2: u64, m2: u64) -> Option<u64> {
    let (m1u, m2u) = (u128::from(m1), u128::from(m2));
    let (a1, a2) = (u128::from(r1 % m1), u128::from(r2 % m2));
    (0..m1u * m2u)
        .find(|&x| x % m1u == a1 && x % m2u == a2)
        .map(|x| x as u64)
}

/// PROOF: `crt2` matches the brute-force solution over all coprime modulus
/// pairs `m1 in 1..40`, `m2 in 2..40`, and every residue pair.
#[test]
fn crt2_matches_brute() {
    for m1 in 1..40u64 {
        for m2 in 2..40u64 {
            if gcd(m1, m2) != 1 {
                continue;
            }
            for r1 in 0..m1 {
                for r2 in 0..m2 {
                    assert_eq!(
                        crt2(r1, m1, r2, m2),
                        brute_crt2(r1, m1, r2, m2),
                        "r1={r1} m1={m1} r2={r2} m2={m2}",
                    );
                }
            }
        }
    }
}

/// PROOF: `crt2` returns `None` for non-coprime moduli, a degenerate modulus,
/// and a product that overflows `u64`.
#[test]
fn crt2_none_cases() {
    for m1 in 2..40u64 {
        for m2 in 2..40u64 {
            if gcd(m1, m2) > 1 {
                assert_eq!(crt2(1, m1, 1, m2), None, "non-coprime m1={m1} m2={m2}");
            }
        }
    }
    assert_eq!(crt2(0, 5, 0, 1), None, "m2 = 1 is degenerate");
    assert_eq!(crt2(0, 1 << 33, 0, 1 << 33), None, "product overflows u64");
}

/// PROOF (defining property): whenever `crt2` returns `Some(x)`, `x` is in
/// range and satisfies both congruences — over 100k coprime near-`2^31`
/// modulus pairs (product fits `u64`), independent of Garner.
#[test]
fn crt2_satisfies_congruences_wide() {
    let mut state = 0x5151_5151_A5A5_A5A5;
    for _ in 0..100_000 {
        let m1 = (lcg(&mut state) >> 33).max(2); // [2, 2^31)
        let m2 = (lcg(&mut state) >> 33).max(2);
        if gcd(m1, m2) != 1 {
            continue;
        }
        let r1 = lcg(&mut state) % m1;
        let r2 = lcg(&mut state) % m2;
        let x = crt2(r1, m1, r2, m2).expect("coprime moduli, product < 2^62");
        assert_eq!(x % m1, r1, "x not ≡ r1 (m1={m1} m2={m2})");
        assert_eq!(x % m2, r2, "x not ≡ r2 (m1={m1} m2={m2})");
        assert!(
            u128::from(x) < u128::from(m1) * u128::from(m2),
            "x out of range"
        );
    }
}

/// PROOF (known answers + the k-modulus fold): textbook CRT instances, and
/// `crt` agreeing with iterated `crt2`.
#[test]
fn crt_known_answers() {
    assert_eq!(crt2(2, 3, 3, 5), Some(8), "8 ≡ 2 (3), 8 ≡ 3 (5)");
    assert_eq!(crt2(1, 2, 2, 3), Some(5), "5 ≡ 1 (2), 5 ≡ 2 (3)");
    assert_eq!(crt(&[(2, 3), (3, 5), (2, 7)]), Some(23), "mod 105");
    assert_eq!(crt(&[]), Some(0), "no constraints");
    assert_eq!(crt(&[(4, 9)]), Some(4), "single congruence");
    assert_eq!(crt(&[(2, 4), (3, 6)]), None, "moduli 4, 6 not coprime");
}

/// Deterministic LCG (Knuth MMIX constants) as a test [`Rng`].
struct LcgRng(u64);
impl Rng for LcgRng {
    fn next_u64(&mut self) -> u64 {
        self.0 = self
            .0
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        self.0
    }
}

/// PROOF: `next_prime(n)` is the least prime `>= n`, matched against the
/// sieve over `0..100_000`.
#[test]
fn next_prime_matches_sieve() {
    let is_p = sieve(100_200); // headroom past 100_000 for the next prime
    for n in 0..100_000usize {
        let want = (n..is_p.len()).find(|&m| is_p[m]).map(|m| m as u64);
        assert_eq!(next_prime(n as u64), want, "n={n}");
    }
}

/// PROOF (known answers, including the very top of the `u64` range).
#[test]
fn next_prime_known_answers() {
    assert_eq!(next_prime(0), Some(2));
    assert_eq!(next_prime(2), Some(2));
    assert_eq!(next_prime(3), Some(3));
    assert_eq!(next_prime(8), Some(11));
    assert_eq!(next_prime(14), Some(17));
    // 18_446_744_073_709_551_557 is the largest u64 prime.
    assert_eq!(
        next_prime(18_446_744_073_709_551_557),
        Some(18_446_744_073_709_551_557)
    );
    assert_eq!(next_prime(18_446_744_073_709_551_558), None);
}

/// PROOF: `random_prime` yields a certified prime of exactly the requested
/// bit length, for every `bits` in `2..=64`; out-of-range `bits` give `None`.
#[test]
fn random_prime_produces_primes() {
    assert_eq!(random_prime(&mut LcgRng(1), 1), None);
    assert_eq!(random_prime(&mut LcgRng(1), 65), None);
    let mut rng = LcgRng(0xDEAD_BEEF_CAFE_F00D);
    for bits in 2..=64u32 {
        for _ in 0..40 {
            let p = random_prime(&mut rng, bits).expect("2..=64 always yields a prime");
            assert!(is_prime(p), "not prime: {p} (bits={bits})");
            assert_eq!(
                64 - p.leading_zeros(),
                bits,
                "wrong bit length: {p} (bits={bits})"
            );
            assert_eq!(p & 1, 1, "even prime returned: {p}");
        }
    }
}

/// Brute-force least common multiple — the smallest positive multiple of both,
/// by scan. An oracle independent of the gcd-based formula.
fn brute_lcm(a: u64, b: u64) -> Option<u64> {
    if a == 0 || b == 0 {
        return Some(0);
    }
    let mut m = a;
    loop {
        if m.is_multiple_of(b) {
            return Some(m);
        }
        m = m.checked_add(a)?;
    }
}

/// PROOF: `lcm` equals the brute-force least common multiple over 0..200², is
/// divisible by both inputs, and satisfies the defining identity
/// `gcd(a,b) * lcm(a,b) == a * b`, plus the `u64`-overflow edge.
#[test]
fn lcm_matches_brute_and_identity() {
    assert_eq!(lcm(0, 0), Some(0));
    for a in 0..200u64 {
        for b in 0..200u64 {
            assert_eq!(lcm(a, b), brute_lcm(a, b), "a={a} b={b}");
            if a != 0 && b != 0 {
                let l = lcm(a, b).expect("0..200 product fits u64");
                assert_eq!(l % a, 0, "lcm {l} not divisible by a={a}");
                assert_eq!(l % b, 0, "lcm {l} not divisible by b={b}");
                assert_eq!(
                    u128::from(gcd(a, b)) * u128::from(l),
                    u128::from(a) * u128::from(b),
                    "gcd*lcm == a*b for a={a} b={b}",
                );
            }
        }
    }
    assert_eq!(lcm(1u64 << 62, (1u64 << 62) + 2), None, "overflow");
    assert_eq!(lcm(u64::MAX, u64::MAX), Some(u64::MAX), "lcm(n,n) = n");
}

/// Legendre symbol via Euler's criterion `a^((p-1)/2) mod p` — an oracle for
/// PRIME `p`, independent of the reciprocity algorithm (it uses `modpow`).
fn legendre_euler(a: u64, p: u64) -> i8 {
    let a = a % p;
    if a == 0 {
        return 0;
    }
    match modpow(a, (p - 1) / 2, p) {
        1 => 1,
        x if x == p - 1 => -1,
        _ => 0, // only reachable for composite p (Euler's criterion assumes prime)
    }
}

/// PROOF: for every odd prime `p`, `jacobi(a, p)` equals the Legendre symbol
/// by Euler's criterion over all `a` in `0..3p` — the Legendre case against an
/// independent `modpow`-based oracle — plus the even-`n` / zero edges.
#[test]
fn jacobi_matches_euler_for_primes() {
    for &p in &[3u64, 5, 7, 11, 13, 17, 19, 23, 97, 101] {
        for a in 0..3 * p {
            assert_eq!(jacobi(a, p), Some(legendre_euler(a, p)), "a={a} p={p}");
        }
    }
    // Large numerators near u64::MAX: the symbol is periodic in `a` (period n),
    // so `jacobi` must REDUCE `a` mod n and agree with Euler's criterion on the
    // reduced residue. A reduction that *added* n instead of reducing would
    // overflow here — this pins the `a %= n` step against that mutation.
    for &p in &[3u64, 5, 7, 97, 101, 7919] {
        for a in [
            u64::MAX,
            u64::MAX - 1,
            u64::MAX - 2,
            1u64 << 63,
            (1u64 << 63) + 1,
        ] {
            assert_eq!(
                jacobi(a, p),
                Some(legendre_euler(a, p)),
                "large a={a} p={p}"
            );
        }
    }
    assert_eq!(jacobi(3, 4), None, "even n undefined");
    assert_eq!(jacobi(3, 0), None, "zero n undefined");
    assert_eq!(jacobi(0, 9), Some(0), "gcd(0,9) != 1");
}

/// PROOF (composite `n`): `jacobi(a, n)` equals the product of the Legendre
/// symbols over `n`'s prime factors (the defining multiplicativity), and is
/// multiplicative in `a`: `jacobi(a*b mod n, n) == jacobi(a,n)*jacobi(b,n)`.
#[test]
fn jacobi_multiplicative_composite() {
    // jacobi(a,n) == ∏ legendre(a, p) over prime factors p (with multiplicity)
    // — an independent factor-based oracle for odd n.
    let factor_legendre = |a: u64, mut n: u64| -> i8 {
        let mut prod = 1i8;
        let mut d = 3u64;
        while d * d <= n {
            while n.is_multiple_of(d) {
                prod *= legendre_euler(a, d); // d is prime: smaller factors removed
                n /= d;
            }
            d += 2;
        }
        if n > 1 {
            prod *= legendre_euler(a, n);
        }
        prod
    };
    for n in (3..400u64).step_by(2) {
        for a in 0..n {
            assert_eq!(jacobi(a, n), Some(factor_legendre(a, n)), "a={a} n={n}");
        }
    }
    for n in [15u64, 21, 35, 45, 91] {
        for a in 0..n {
            for b in 0..n {
                assert_eq!(
                    jacobi(a * b % n, n).unwrap(),
                    jacobi(a, n).unwrap() * jacobi(b, n).unwrap(),
                    "multiplicative a={a} b={b} n={n}",
                );
            }
        }
    }
}

/// PROOF: for every prime `p` (covering `p ≡ 3` and `p ≡ 1 (mod 4)`, i.e. the
/// closed form and the Tonelli–Shanks loop) and every `a` in `[0, p)`,
/// `sqrtmod` returns `Some(x)` with `x*x ≡ a` and `x <= p/2` exactly when `a`
/// is a quadratic residue, and `None` exactly otherwise — adjudicated by an
/// independent residue-set scan.
#[test]
fn sqrtmod_self_verifies_over_primes() {
    for &p in &[2u64, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 97, 101] {
        let mut qr = std::collections::HashSet::new();
        for y in 0..p {
            qr.insert(y * y % p);
        }
        for a in 0..p {
            match sqrtmod(a, p) {
                Some(x) => {
                    assert_eq!(x * x % p, a, "sqrtmod({a},{p})={x}: x² ≢ a");
                    assert!(x <= p / 2, "non-canonical root {x} (p={p})");
                    assert!(qr.contains(&a), "root for non-residue a={a} p={p}");
                }
                None => assert!(!qr.contains(&a), "None for residue a={a} p={p}"),
            }
        }
    }
}

/// KAT: `sqrtmod` on large primes of both residue classes — `10^9+7`,
/// `2^31-1`, `2^61-1` (`≡ 3 mod 4`, closed form) and `998244353` (`≡ 1 mod 4`,
/// Tonelli–Shanks) — self-verified `x*x ≡ a (mod p)`.
#[test]
fn sqrtmod_large_primes() {
    for &p in &[
        1_000_000_007u64,
        2_147_483_647,
        998_244_353,
        2_305_843_009_213_693_951,
    ] {
        for a in [2u64, 3, 5, 10, 12_345, 999_999] {
            let a = a % p;
            if let Some(x) = sqrtmod(a, p) {
                assert_eq!(
                    u128::from(x) * u128::from(x) % u128::from(p),
                    u128::from(a),
                    "p={p} a={a} x={x}",
                );
            }
        }
    }
}

/// PROOF (misuse-safety): `sqrtmod` is total and bounded — every non-prime
/// modulus returns `None` rather than entering the `O(p)` non-residue search.
/// Covers `0`, `1`, small composites, a large odd strong-pseudoprime
/// (`3_215_031_751`, composite but a base-2 MR liar), and large even moduli
/// that would hang an unguarded search.
#[test]
fn sqrtmod_rejects_nonprime() {
    for p in [
        0u64,
        1,
        4,
        6,
        9,
        15,
        21,
        25,
        100,
        1000,
        3_215_031_751, // = 151·751·28351, composite
        1u64 << 40,    // large even
        (1u64 << 62) + 4,
    ] {
        for a in [0u64, 1, 2, 3, 7] {
            assert_eq!(sqrtmod(a, p), None, "non-prime p={p} a={a} must be None");
        }
    }
}

/// PROOF (deep Tonelli–Shanks, bidirectional): over FFT/NTT primes
/// `p ≡ 1 (mod 4)` with large 2-adic valuation `s = v₂(p-1)` (up to 30) — the
/// primes that drive the TS descent through its nested `t^(2^i)` loop, the
/// `i + 1 < m` bound, and the `1 << (m-i-1)` shift at large `m`, the path the
/// small-prime tests never reach — `sqrtmod` returns `Some(x)` with `x² ≡ a`
/// and `x ≤ p/2` exactly when `a` is a residue, `None` exactly otherwise. The
/// residue verdict is adjudicated independently by Euler's criterion
/// (`legendre_euler`, the `a^((p-1)/2)` route), a *different* computation from
/// the quadratic reciprocity `sqrtmod` uses internally.
#[test]
#[allow(clippy::cast_possible_truncation)] // b² mod p < p < 2^32 fits u64
fn sqrtmod_deep_tonelli_shanks() {
    for &p in &[
        998_244_353u64, // 119·2^23 + 1
        167_772_161,    // 5·2^25 + 1
        469_762_049,    // 7·2^26 + 1
        2_013_265_921,  // 15·2^27 + 1
        3_221_225_473,  // 3·2^30 + 1
    ] {
        assert!(is_prime(p), "{p} must be prime");
        assert!(
            (p - 1).trailing_zeros() >= 23,
            "{p} should have a deep 2-adic valuation"
        );
        let mut state = p ^ 0x5DEE_CE66_1357_9BDF;
        for _ in 0..3_000 {
            let a = lcg(&mut state) % p;
            let is_residue = a == 0 || legendre_euler(a, p) == 1;
            match sqrtmod(a, p) {
                Some(x) => {
                    assert!(is_residue, "sqrtmod Some for non-residue a={a} p={p}");
                    assert!(x <= p / 2, "non-canonical root {x} (p={p})");
                    assert_eq!(
                        u128::from(x) * u128::from(x) % u128::from(p),
                        u128::from(a),
                        "x² ≢ a for a={a} p={p} x={x}",
                    );
                }
                None => assert!(!is_residue, "sqrtmod None for residue a={a} p={p}"),
            }
        }
        // Guaranteed residues a = b² mod p heavily exercise the Some-arm at
        // large s without needing an O(p) scan to certify residue-ness.
        let mut bstate = p ^ 0xF00D_CAFE_1234_5678;
        for _ in 0..2_000 {
            let b = lcg(&mut bstate) % p;
            let a = (u128::from(b) * u128::from(b) % u128::from(p)) as u64;
            let x = sqrtmod(a, p).expect("b² is always a residue");
            assert_eq!(
                u128::from(x) * u128::from(x) % u128::from(p),
                u128::from(a),
                "perfect-square root wrong: p={p} b={b} a={a} x={x}",
            );
            assert!(x <= p / 2, "non-canonical root {x} (p={p})");
        }
    }
}

/// PROOF (k-ary CRT, wide + order-independent): over 50k random subsets of
/// distinct primes (pairwise coprime by construction, product kept `< 2^63`),
/// the `crt` fold returns a value satisfying *every* congruence and lying in
/// `[0, ∏ mᵢ)` — the defining property, independent of Garner — and the same
/// value regardless of fold order (a shuffled list yields the identical unique
/// solution). Closes the order-independence property the matrix only implied.
#[test]
#[allow(clippy::cast_possible_truncation)] // LCG-mod shuffle index fits usize
fn crt_satisfies_congruences_wide() {
    const PRIMES: [u64; 9] = [3, 5, 7, 11, 13, 17, 19, 23, 29];
    let mut state = 0x0C27_1234_ABCD_5151;
    for _ in 0..50_000 {
        let mask = lcg(&mut state);
        let mut pairs: Vec<(u64, u64)> = Vec::new();
        let mut prod: u128 = 1;
        for (i, &m) in PRIMES.iter().enumerate() {
            if (mask >> i) & 1 == 1 && prod * u128::from(m) < (1u128 << 63) {
                pairs.push((lcg(&mut state) % m, m));
                prod *= u128::from(m);
            }
        }
        let x = crt(&pairs).expect("distinct primes are coprime, product < 2^63");
        for &(r, m) in &pairs {
            assert_eq!(x % m, r, "x={x} not ≡ {r} (mod {m})");
        }
        assert!(u128::from(x) < prod, "x={x} out of range (∏ = {prod})");
        // Order-independence: a Fisher–Yates shuffle yields the same solution.
        let mut shuffled = pairs.clone();
        for i in (1..shuffled.len()).rev() {
            let j = (lcg(&mut state) % (i as u64 + 1)) as usize;
            shuffled.swap(i, j);
        }
        assert_eq!(
            crt(&shuffled),
            Some(x),
            "crt not order-independent: {pairs:?}"
        );
    }
}

/// PROOF (overflow boundary, independent oracle): the true lcm computed in
/// `u128` (`a/gcd · b`, exact since `gcd | a`) decides Some vs None — `lcm`
/// returns `Some(value)` exactly when it fits `u64`, `None` exactly when it
/// exceeds `u64::MAX`. Swept over 200k full-range pairs (most land in the
/// overflow regime the brute scan and the two KAT points cannot reach), the
/// `u64::MAX` corners, and a band straddling `u64::MAX / a`.
#[test]
fn lcm_overflow_boundary() {
    // Stein's binary GCD (shift/subtract, no division) — independent of the
    // production Euclid `gcd`, and O(log) so it runs at the 2^64 boundary
    // where the brute divisor scan cannot. Keeps the overflow oracle fully
    // independent (anti-tautology): a `gcd` bug cannot corrupt both the oracle
    // and `lcm` into agreement.
    fn binary_gcd(mut a: u64, mut b: u64) -> u64 {
        if a == 0 {
            return b;
        }
        if b == 0 {
            return a;
        }
        let shift = (a | b).trailing_zeros();
        a >>= a.trailing_zeros();
        loop {
            b >>= b.trailing_zeros();
            if a > b {
                std::mem::swap(&mut a, &mut b);
            }
            b -= a;
            if b == 0 {
                break;
            }
        }
        a << shift
    }
    let lcm_true = |a: u64, b: u64| -> u128 {
        if a == 0 || b == 0 {
            return 0;
        }
        u128::from(a / binary_gcd(a, b)) * u128::from(b) // exact: gcd | a
    };
    let check = |a: u64, b: u64| {
        let t = lcm_true(a, b);
        match lcm(a, b) {
            Some(l) => {
                assert!(
                    t <= u128::from(u64::MAX),
                    "lcm Some({l}) but true {t} > u64::MAX (a={a} b={b})"
                );
                assert_eq!(u128::from(l), t, "lcm {l} != true {t} (a={a} b={b})");
            }
            None => assert!(
                t > u128::from(u64::MAX),
                "lcm None but true {t} fits u64 (a={a} b={b})"
            ),
        }
    };
    let mut state = 0x1CEB_00DA_5A5F_00D5;
    for _ in 0..200_000 {
        check(lcg(&mut state), lcg(&mut state));
    }
    check(u64::MAX, u64::MAX);
    check(1u64 << 63, 1u64 << 63);
    check((1u64 << 63) + 1, 2);
    check(u64::MAX, 2);
    // A band straddling the boundary b ≈ u64::MAX / a.
    for k in 0..64u64 {
        let a = (1u64 << 32) + k;
        let b = u64::MAX / a;
        check(a, b);
        check(a, b + 1);
    }
}

/// Reference Miller–Rabin round as the full squaring sequence with no early
/// return — a separately-structured decision procedure used to cross-check
/// production `is_strong_probable_prime` on the SOUNDNESS quadrant (a strong
/// witness for a composite returns `false`) over arbitrary `(n, a)`, which the
/// `is_prime` path never drives (it only ever uses bases `{2..37}`).
#[allow(clippy::many_single_char_names)] // n, a, s, d, x — canonical MR names
fn reference_sprp(n: u64, a: u64) -> bool {
    if n == 2 || n == 3 {
        return true;
    }
    if n < 2 || n & 1 == 0 {
        return false;
    }
    let a = a % n;
    if a <= 1 || a == n - 1 {
        return true;
    }
    let s = (n - 1).trailing_zeros();
    let d = (n - 1) >> s;
    let mut x = modpow(a, d, n);
    if x == 1 {
        return true; // a^d ≡ 1
    }
    for _ in 0..s {
        if x == n - 1 {
            return true; // a^(2^j·d) ≡ -1 for some 0 <= j < s
        }
        x = modpow(x, 2, n);
    }
    false
}

/// PROOF (single-round MR soundness over arbitrary (n, a)): production
/// `is_strong_probable_prime` agrees with the reference MR for every `n` in
/// `3..30_000` (composite and prime) across bases inside and outside the
/// witness set `{2..37}` plus pseudo-random `u64`, and over large known
/// composites (incl. the ψ strong pseudoprimes) — pinning the witness ⇒
/// `false` direction the `is_prime` path never exercises.
#[test]
fn spp_soundness_vs_reference_mr() {
    let mut state = 0x50DD_1234_9999_ABCD;
    for n in 3..30_000u64 {
        let bases = [
            2u64,
            3,
            5,
            37,
            38,
            40,
            100,
            n - 1,
            n + 1,
            lcg(&mut state),
            lcg(&mut state),
        ];
        for &a in &bases {
            assert_eq!(
                is_strong_probable_prime(n, a),
                reference_sprp(n, a),
                "spp disagreement at n={n} a={a}",
            );
        }
    }
    for &n in &[
        3_215_031_751u64,      // ψ_4 = 151·751·28351
        2_152_302_898_747,     // ψ_5
        1_000_003 * 1_000_033, // product of two primes ≈ 10^12
        4_294_967_291 * 3,     // 3 · (largest prime < 2^32)
    ] {
        assert!(!is_prime(n), "{n} must be composite");
        let mut st = n ^ 0x1357_2468_ACE0_BDF9;
        for _ in 0..3_000 {
            let a = lcg(&mut st);
            assert_eq!(
                is_strong_probable_prime(n, a),
                reference_sprp(n, a),
                "spp disagreement at composite n={n} a={a}",
            );
        }
    }
}

/// PROOF (candidate-construction reachability — coprime's own `(rng & mask) |
/// top | 1` logic, not the RNG's job): for every band small enough to
/// enumerate, each odd prime in `[2^(bits-1), 2^bits)` is reachable — feeding
/// the word `w = q` reconstructs the candidate `q`, so `random_prime` returns
/// it — and the construction lands in-band, odd, and of the exact bit length
/// for arbitrary words. (2 is intentionally unreachable: `| 1` forces odd.)
#[test]
#[allow(clippy::cast_possible_truncation)] // band primes index the sieve losslessly on 64-bit
fn random_prime_reaches_every_band_prime() {
    struct Replay {
        words: Vec<u64>,
        i: usize,
        tail: u64,
    }
    impl Rng for Replay {
        fn next_u64(&mut self) -> u64 {
            let v = self.words.get(self.i).copied().unwrap_or(self.tail);
            self.i += 1;
            v
        }
    }
    let top_bits = 18u32;
    let is_p = sieve(1usize << top_bits);
    for bits in 2..=top_bits {
        let lo = 1u64 << (bits - 1);
        let hi = 1u64 << bits;
        for q in (lo..hi).filter(|&n| (n & 1) == 1 && is_p[n as usize]) {
            let mut rng = Replay {
                words: vec![q],
                i: 0,
                tail: q,
            };
            assert_eq!(
                random_prime(&mut rng, bits),
                Some(q),
                "prime {q} not reachable at bits={bits}",
            );
        }
    }
    // The construction maps arbitrary words into the band, odd, exact length.
    let mut state = 0xBEEF_0000_1111_2222;
    for bits in 2..=20u32 {
        let top = 1u64 << (bits - 1);
        let mask = if bits == 64 {
            u64::MAX
        } else {
            (1u64 << bits) - 1
        };
        for _ in 0..1_000 {
            let cand = (lcg(&mut state) & mask) | top | 1;
            assert!(
                cand >= top && cand < (top << 1),
                "candidate {cand} out of band (bits={bits})"
            );
            assert_eq!(cand & 1, 1, "candidate {cand} not odd");
            assert_eq!(
                64 - cand.leading_zeros(),
                bits,
                "candidate {cand} wrong bit length"
            );
        }
    }
}

/// Integer square root by a computationally distinct route: a hardware `f64`
/// sqrt seed corrected by a bounded integer walk. The comparisons use
/// division (`r > n/r ⟺ r*r > n`; `r+1 <= n/(r+1) ⟺ (r+1)² <= n`) so no
/// product ever overflows `u64` — independent of the bit-recurrence.
#[allow(
    clippy::cast_precision_loss,
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss
)]
fn isqrt_oracle(n: u64) -> u64 {
    if n == 0 {
        return 0;
    }
    let mut r = (n as f64).sqrt() as u64;
    while r > 0 && r > n / r {
        r -= 1; // f64 overshoot: r*r > n
    }
    while (r + 1) <= n / (r + 1) {
        r += 1; // (r+1)² <= n
    }
    r
}

/// PROOF: `isqrt` satisfies `r*r <= n < (r+1)²` (in the `u64`-safe form
/// `n - r*r <= 2*r`, since `(r+1)²` overflows at `n = u64::MAX`) and matches
/// the independent `f64`-seed oracle. Exhaustive over the boundary lattice
/// (every `k²-1`, `k²`, `k²+2k` for `k < 2^20`, where the answer changes),
/// the `2^32` overflow corner, KATs, and 200k full-range pairs.
#[test]
fn isqrt_matches_postcondition_and_oracle() {
    let check = |n: u64| {
        let r = isqrt(n);
        assert!(r <= 0xFFFF_FFFF, "isqrt({n})={r} exceeds 2^32-1");
        assert!(r * r <= n, "isqrt({n})={r}: r² > n"); // r*r <= n, so safe
        assert!(n - r * r <= 2 * r, "isqrt({n})={r}: n ≥ (r+1)²"); // u64-safe upper bound
    };
    for (n, want) in [
        (0u64, 0u64),
        (1, 1),
        (2, 1),
        (3, 1),
        (4, 2),
        (8, 2),
        (9, 3),
        (15, 3),
        (16, 4),
        (u64::MAX, 4_294_967_295),
    ] {
        assert_eq!(isqrt(n), want, "isqrt({n})");
        check(n);
    }
    // Exhaustive over every place the floor changes, k in 0..2^20.
    for k in 0..(1u64 << 20) {
        let sq = k * k;
        assert_eq!(isqrt(sq), k, "perfect square k={k}");
        if k > 0 {
            assert_eq!(isqrt(sq - 1), k - 1, "k²-1 k={k}");
        }
        assert_eq!(isqrt(sq + 2 * k), k, "k²+2k (last before (k+1)²) k={k}");
    }
    // The 2^32 corner where (r+1)² overflows u64.
    for k in 0xFFFF_FFF0u64..=0xFFFF_FFFF {
        check(k * k);
        check(k * k + 2 * k);
    }
    // Wide random vs the independent f64-seed oracle.
    let mut state = 0x15A1_2345_6789_ABCD;
    for _ in 0..200_000 {
        let n = lcg(&mut state);
        let r = isqrt(n);
        assert_eq!(r, isqrt_oracle(n), "oracle disagrees at n={n}");
        check(n);
    }
}

#[test]
fn isqrt_matches_std() {
    // u64::isqrt (std) is a second, fully independent ⌊√n⌋ — a different algorithm,
    // exact, and ceiling-free. coprime's binary digit-recurrence must agree with it
    // everywhere; the point is to witness coprime's own implementation, not adopt
    // std's. Exhaustive over a low range, then the perfect-square boundaries (where
    // the floor steps) across the full 2^32 root range, then full-u64 random.
    for n in 0..3_000_000u64 {
        assert_eq!(isqrt(n), n.isqrt(), "isqrt({n})");
    }
    let mut state = 0x5151_5151_2727_2727u64;
    for _ in 0..300_000 {
        let r = lcg(&mut state) >> 32; // candidate root in [0, 2^32)
        let rr = r * r; // < 2^64
        assert_eq!(isqrt(rr), rr.isqrt(), "isqrt(r²) r={r}");
        assert_eq!(
            isqrt(rr + 2 * r),
            (rr + 2 * r).isqrt(),
            "isqrt((r+1)²−1) r={r}"
        );
        if r > 0 {
            assert_eq!(isqrt(rr - 1), (rr - 1).isqrt(), "isqrt(r²−1) r={r}");
        }
    }
    for _ in 0..2_000_000 {
        let n = lcg(&mut state);
        assert_eq!(isqrt(n), n.isqrt(), "isqrt({n})");
    }
    assert_eq!(isqrt(u64::MAX), u64::MAX.isqrt());
}

/// Integer k-th root by a distinct route: a hardware `f64` `powf` seed
/// corrected by a bounded `checked_pow` walk — independent of the binary
/// search, and overflow-safe (an overflowing `r^k` counts as `> n`).
#[allow(
    clippy::cast_precision_loss,
    clippy::cast_possible_truncation,
    clippy::cast_sign_loss
)]
fn iroot_oracle(n: u64, k: u32) -> u64 {
    if k <= 1 || n <= 1 {
        return n;
    }
    let mut r = (n as f64).powf(1.0 / f64::from(k)) as u64;
    while r > 0 && r.checked_pow(k).is_none_or(|p| p > n) {
        r -= 1; // f64 overshoot (or overflow): r^k > n
    }
    while (r + 1).checked_pow(k).is_some_and(|p| p <= n) {
        r += 1; // (r+1)^k <= n
    }
    r
}

/// PROOF: `iroot` satisfies `r^k <= n < (r+1)^k` (checked via `checked_pow`,
/// so the `(r+1)^k` upper bound never overflows), agrees with `isqrt` at
/// `k = 2`, and matches the independent `f64`-`powf` oracle — over a `k × n`
/// grid, every perfect-power boundary `r^k ± 1`, KATs, and 100k random pairs.
#[test]
fn iroot_matches_isqrt_postcondition_and_oracle() {
    // k = 2 must agree with the dedicated isqrt.
    for n in 0..50_000u64 {
        assert_eq!(iroot(n, 2), isqrt(n), "iroot(.,2) vs isqrt at n={n}");
    }
    for (n, k, want) in [
        (0u64, 3u32, 0u64),
        (1, 3, 1),
        (8, 3, 2),
        (26, 3, 2),
        (27, 3, 3),
        (63, 3, 3),
        (64, 3, 4),
        (1_000_000, 3, 100),
        (u64::MAX, 2, 4_294_967_295),
        (255, 8, 1),
        (256, 8, 2),
        (123, 1, 123), // 1st root
        (123, 0, 123), // degenerate k = 0
    ] {
        assert_eq!(iroot(n, k), want, "iroot({n}, {k})");
    }
    let check = |n: u64, k: u32| {
        let r = iroot(n, k);
        assert_eq!(r, iroot_oracle(n, k), "oracle: n={n} k={k}");
        assert!(
            r.checked_pow(k).is_some_and(|p| p <= n),
            "r^k > n: n={n} k={k} r={r}"
        );
        assert!(
            (r + 1).checked_pow(k).is_none_or(|p| p > n),
            "(r+1)^k <= n: n={n} k={k} r={r}"
        );
    };
    for k in 2..=10u32 {
        for n in 0..2_000u64 {
            check(n, k);
        }
        // every perfect-power boundary r^k and r^k - 1, where the floor changes
        for r in 0..1_500u64 {
            if let Some(p) = r.checked_pow(k) {
                check(p, k);
                if p > 0 {
                    check(p - 1, k);
                }
            }
        }
    }
    let mut state = 0x1234_ABCD_5678_EF90;
    for _ in 0..100_000 {
        let n = lcg(&mut state);
        let k = 2 + (lcg(&mut state) % 31) as u32; // k in 2..=32
        check(n, k);
    }
}

/// Kronecker symbol by an independent factor-based route: pull the factors of
/// 2 (each contributing the `(a|2)` supplement) then multiply the Legendre
/// symbols of the odd prime factors via Euler's criterion (`legendre_euler`,
/// a `modpow` route) — never the quadratic reciprocity `kronecker` uses.
#[allow(clippy::many_single_char_names)]
fn kronecker_oracle(a: u64, n: u64) -> i8 {
    if n == 0 {
        return i8::from(a == 1);
    }
    let mut m = n;
    let mut prod = 1i8;
    while m.is_multiple_of(2) {
        m /= 2;
        prod *= if a.is_multiple_of(2) {
            0
        } else if a % 8 == 1 || a % 8 == 7 {
            1
        } else {
            -1
        };
    }
    // m is now odd; factor it and multiply Legendre symbols (m's prime factors).
    let mut d = 3u64;
    while d * d <= m {
        while m.is_multiple_of(d) {
            prod *= legendre_euler(a, d);
            m /= d;
        }
        d += 2;
    }
    if m > 1 {
        prod *= legendre_euler(a, m);
    }
    prod
}

/// PROOF: `kronecker` equals the independent factor-based oracle over all
/// `n` in `0..400` (including even `n`, `n = 0, 1, 2`) and `a` in `0..3n+3`,
/// is a faithful superset of `jacobi` (equal for every odd `n`), and matches
/// a KAT table of boundary values.
#[test]
fn kronecker_matches_factor_oracle_and_jacobi() {
    for (a, n, want) in [
        (3u64, 2u64, -1i8),
        (7, 2, 1),
        (1, 0, 1),
        (2, 0, 0),
        (0, 0, 0),
        (5, 12, -1),
        (15, 16, 1),
        (3, 4, 1),
        (0, 1, 1),
    ] {
        assert_eq!(kronecker(a, n), want, "KAT ({a}|{n})");
    }
    for n in 0..400u64 {
        for a in 0..3 * n + 3 {
            assert_eq!(kronecker(a, n), kronecker_oracle(a, n), "({a}|{n})");
        }
    }
    // Faithful extension of jacobi: equal for every odd n.
    for n in (1..400u64).step_by(2) {
        for a in 0..n {
            assert_eq!(kronecker(a, n), jacobi(a, n).unwrap(), "odd ({a}|{n})");
        }
    }
}

/// PROOF (reconstruction oracle, distinct from the rho search): `factorize`
/// returns ascending `(prime, exponent)` pairs whose product is `n` and whose
/// every base is prime — checked over `0..200_000` and hard KATs (the
/// `u64::MAX` factorization, a Mersenne prime, large semiprimes, and prime
/// powers). Reconstruction multiplies + primality-tests; it never factors.
#[test]
fn factorize_reconstructs_and_each_factor_prime() {
    let check = |n: u64| {
        let f = factorize(n);
        let s = f.as_slice();
        for w in s.windows(2) {
            assert!(w[0].0 < w[1].0, "factors not ascending/distinct for n={n}");
        }
        let mut prod = 1u128;
        for &(p, e) in s {
            assert!(is_prime(p), "non-prime factor {p} of n={n}");
            assert!(e >= 1, "zero exponent for {p} (n={n})");
            for _ in 0..e {
                prod *= u128::from(p);
            }
        }
        if n < 2 {
            assert!(s.is_empty(), "n={n} (<2) must have no prime factors");
        } else {
            assert_eq!(prod, u128::from(n), "reconstruction mismatch for n={n}");
        }
    };
    for n in 0..200_000u64 {
        check(n);
    }
    // Hard KATs: the cases trial division alone cannot resolve.
    check(u64::MAX); // 3·5·17·257·641·65537·6700417
    check(2_305_843_009_213_693_951); // 2^61 - 1, prime
    check(1_000_003 * 1_000_033); // semiprime ≈ 10^12
    check(4_294_967_291 * 4_294_967_279); // two large primes near 2^32
    check(1u64 << 63); // 2^63
    check(3u64.pow(40)); // 3^40
    check((1u64 << 31) - 1); // 2^31 - 1 (M31), prime
    assert_eq!(
        factorize(u64::MAX).as_slice(),
        &[
            (3, 1),
            (5, 1),
            (17, 1),
            (257, 1),
            (641, 1),
            (65_537, 1),
            (6_700_417, 1)
        ],
        "u64::MAX = 2^64 - 1 factorization",
    );
}

/// Brute totient: the count of `a` in `[1, n]` with `gcd(a, n) = 1` — the
/// *definition*, independent of the prime-factorization formula.
#[allow(clippy::cast_possible_truncation)] // count <= n fits u64
fn brute_totient(n: u64) -> u64 {
    if n == 0 {
        return 0;
    }
    (1..=n).filter(|&a| gcd(a, n) == 1).count() as u64
}

/// PROOF: `totient` equals the coprime-count definition over `1..=20_000`,
/// plus KATs and the largest `u64` prime (`φ(p) = p − 1`).
#[test]
fn totient_matches_coprime_count() {
    assert_eq!(totient(0), 0);
    for n in 1..=20_000u64 {
        assert_eq!(totient(n), brute_totient(n), "φ({n})");
    }
    for (n, want) in [(1u64, 1u64), (36, 12), (561, 320), (1_000_003, 1_000_002)] {
        assert_eq!(totient(n), want, "φ({n})");
    }
    let p = 18_446_744_073_709_551_557u64; // largest u64 prime
    assert_eq!(totient(p), p - 1);
}

/// Brute Carmichael: the `lcm` of the multiplicative orders of all units mod
/// `n`, each order found by repeated multiplication — independent of the
/// prime-power formula.
#[allow(clippy::cast_possible_truncation)] // x*a % n < n fits u64
fn brute_carmichael(n: u64) -> Option<u64> {
    if n == 0 {
        return None;
    }
    let mut lam = 1u64;
    for a in 1..n {
        if gcd(a, n) != 1 {
            continue;
        }
        let mut x = a;
        let mut k = 1u64;
        while x != 1 {
            x = (u128::from(x) * u128::from(a) % u128::from(n)) as u64;
            k += 1;
        }
        lam = lcm(lam, k)?;
    }
    Some(lam)
}

/// PROOF: `carmichael` equals the universal-exponent definition over
/// `1..=500`, plus the `λ(2^k)` family and Carmichael-number KATs.
#[test]
fn carmichael_matches_universal_exponent() {
    assert_eq!(carmichael(0), None);
    for n in 1..=500u64 {
        assert_eq!(carmichael(n), brute_carmichael(n), "λ({n})");
    }
    for (n, want) in [
        (1u64, 1u64),
        (2, 1),
        (4, 2),
        (8, 2),
        (16, 4),
        (15, 4),
        (561, 80),
        (1729, 36),
    ] {
        assert_eq!(carmichael(n), Some(want), "λ({n})");
    }
}

/// PROOF (closes the Euler-form implied-coverage notes for `modinv`/`modpow`):
/// over every modulus `m` in `2..=3_000` and every unit `a`, Euler's theorem
/// `a^φ(m) ≡ 1`, Carmichael's `a^λ(m) ≡ 1` with `λ(m) | φ(m)`, and the
/// inverse-as-power identity `modinv(a, m) == a^(φ(m)-1)`, asserted directly
/// via `totient`/`carmichael`. A large prime modulus carries the laws past the
/// brute range.
#[test]
fn euler_carmichael_laws() {
    for m in 2..=3_000u64 {
        let phi = totient(m);
        let lam = carmichael(m).expect("m >= 1 has a unit-group exponent");
        assert_eq!(phi % lam, 0, "λ({m}) does not divide φ({m})");
        for a in 1..m {
            if gcd(a, m) != 1 {
                continue;
            }
            assert_eq!(modpow(a, phi, m), 1, "Euler a^φ for a={a} m={m}");
            assert_eq!(modpow(a, lam, m), 1, "Carmichael a^λ for a={a} m={m}");
            assert_eq!(
                modinv(a, m),
                Some(modpow(a, phi - 1, m)),
                "inverse-as-power a={a} m={m}",
            );
        }
    }
    let m = 1_000_000_007u64; // prime ⇒ φ = m − 1
    assert_eq!(totient(m), m - 1);
    for a in [2u64, 3, 5, 123_456, m - 1] {
        assert_eq!(modpow(a, m - 1, m), 1, "Euler large m a={a}");
        assert_eq!(modinv(a, m), Some(modpow(a, m - 2, m)), "inverse a={a}");
    }
}

/// Brute multiplicative order: repeated multiplication until `a^k ≡ 1` — a
/// route independent of the factor-and-divide algorithm.
#[allow(clippy::cast_possible_truncation)] // x*a % m < m fits u64
fn brute_order(a: u64, m: u64) -> Option<u64> {
    if m <= 1 {
        return None;
    }
    let ar = a % m;
    if gcd(ar, m) != 1 {
        return None;
    }
    let mut x = ar;
    let mut k = 1u64;
    while x != 1 {
        x = (u128::from(x) * u128::from(ar) % u128::from(m)) as u64;
        k += 1;
    }
    Some(k)
}

/// PROOF: `order` matches the repeated-multiplication definition over all
/// `(a, m)` with `m in 1..=500`, plus KATs and the non-unit `None` cases.
#[test]
fn order_matches_brute() {
    for m in 1..=500u64 {
        for a in 0..m {
            assert_eq!(order(a, m), brute_order(a, m), "ord({a}, {m})");
        }
    }
    assert_eq!(order(2, 7), Some(3));
    assert_eq!(order(3, 7), Some(6));
    assert_eq!(order(1, 5), Some(1));
    assert_eq!(order(0, 5), None);
    assert_eq!(order(5, 1), None);
    // Large `a` (a >= m): `order` reduces via `a % m`, so the answer depends only
    // on the residue class. The overflow case a = u64::MAX pins that reduction —
    // an `a % m -> a + m` mutation wraps/overflows here and is caught.
    assert_eq!(order(u64::MAX, 3), None); // MAX % 3 == 0 — not a unit; a + m overflows
    assert_eq!(order(u64::MAX, 2), Some(1)); // MAX % 2 == 1 — the unit 1 mod 2
}

/// Brute primitive root: the least unit whose order equals `φ(m)` (counted,
/// not formula-derived), or `None` if no unit generates the whole group.
fn brute_primitive_root(m: u64) -> Option<u64> {
    if m == 0 {
        return None;
    }
    if m == 1 {
        return Some(0);
    }
    let phi = brute_totient(m);
    (1..m).find(|&g| gcd(g, m) == 1 && brute_order(g, m) == Some(phi))
}

/// PROOF: `primitive_root` matches the brute generator search over `0..=500`
/// (returning `None` for exactly the non-cyclic moduli), matches reference KATs, and
/// every generator it returns has order `φ(m)`.
#[test]
fn primitive_root_matches_brute_and_cyclicity() {
    for m in 0..=500u64 {
        assert_eq!(primitive_root(m), brute_primitive_root(m), "primroot({m})");
    }
    for (m, want) in [(2u64, 1u64), (4, 3), (7, 3), (18, 5), (25, 2)] {
        assert_eq!(primitive_root(m), Some(want), "primroot({m})");
    }
    for non_cyclic in [8u64, 12, 15, 16, 21, 24, 30] {
        assert_eq!(
            primitive_root(non_cyclic),
            None,
            "{non_cyclic} is not cyclic"
        );
    }
    // Every returned generator genuinely has order φ(m), including large m.
    for &m in &[3u64, 5, 9, 23, 25, 18, 22, 1_000_000_007] {
        if let Some(g) = primitive_root(m) {
            assert_eq!(order(g, m), Some(totient(m)), "generator order m={m}");
        }
    }
}

// ── Continued fractions ───────────────────────────────────────────────────
//
// `continued_fraction` is the Euclidean quotient sequence; its convergents are the
// best rational approximations of a/b. Two properties pin it, NEITHER a re-run of
// the production Euclid:
//   • reconstruction — the last convergent is a/b in lowest terms, i.e.
//       h_n·gcd(a,b) == a and k_n·gcd(a,b) == b;
//   • the determinant identity h_k·k_{k-1} − h_{k-1}·k_k == ±1 for every
//       consecutive pair — the cofactor identity the production code never forms,
//       checked in u128 since the cross-products reach 2^128.
// plus the best-approximation bound |a·k_k − b·h_k|·k_{k+1} <= b, the very fact
// Wiener's small-d attack exploits on e/N.

// Independent oracle: a fresh Euclid for the quotients, and the convergent
// recurrence carried in u128 — a different width and a separate code path from the
// u64 production iterator, so agreement is a genuine cross-check, not a tautology.
#[allow(clippy::many_single_char_names)] // x, y, r, q, h, k — Euclid / convergent names
fn cf_oracle(a: u64, b: u64) -> (Vec<u64>, Vec<(u128, u128)>) {
    let mut qs = Vec::new();
    let (mut x, mut y) = (a, b);
    if y == 0 {
        qs.push(x);
    } else {
        loop {
            qs.push(x / y);
            let r = x % y;
            if r == 0 {
                break;
            }
            (x, y) = (y, r);
        }
    }
    let mut convs = Vec::new();
    let (mut hm2, mut hm1) = (0u128, 1u128);
    let (mut km2, mut km1) = (1u128, 0u128);
    for &q in &qs {
        let q = u128::from(q);
        let h = q * hm1 + hm2;
        let k = q * km1 + km2;
        convs.push((h, k));
        (hm2, hm1) = (hm1, h);
        (km2, km1) = (km1, k);
    }
    (qs, convs)
}

#[test]
fn continued_fraction_known_answers() {
    // Classic expansions (Knuth TAOCP §4.5.3; the 649/200 Pell example).
    assert_eq!(continued_fraction(415, 93).as_slice(), &[4, 2, 6, 7]);
    assert_eq!(continued_fraction(649, 200).as_slice(), &[3, 4, 12, 4]);
    assert_eq!(continued_fraction(5, 3).as_slice(), &[1, 1, 2]);
    assert_eq!(continued_fraction(1, 2).as_slice(), &[0, 2]); // a < b opens with 0
    assert_eq!(continued_fraction(13, 8).as_slice(), &[1, 1, 1, 1, 2]); // Fibonacci ratio
    // Integers and degenerate edges.
    assert_eq!(continued_fraction(7, 1).as_slice(), &[7]);
    assert_eq!(continued_fraction(1, 1).as_slice(), &[1]);
    assert_eq!(continued_fraction(0, 5).as_slice(), &[0]);
    assert_eq!(continued_fraction(7, 0).as_slice(), &[7]); // a/0 ↦ [a] by convention
    assert_eq!(continued_fraction(0, 0).as_slice(), &[0]);
}

// A `u64/u64` expansion has at most 92 convergents (Lamé). Tests cap iteration of
// the production `Convergents` here, so a bug or mutation that breaks `next()`'s
// termination is CAUGHT fast by the assertions below (count/value mismatch) rather
// than collected unbounded into OOM or spun to a timeout. Correct code yields
// `<= 92` then `None`, so the cap is transparent.
const CONV_GUARD: usize = 128;

#[test]
fn continued_fraction_convergents_known() {
    let convs: Vec<(u64, u64)> = continued_fraction(415, 93)
        .convergents()
        .take(CONV_GUARD)
        .collect();
    assert_eq!(convs, vec![(4, 1), (9, 2), (58, 13), (415, 93)]);
    // a < b opens at 0/1 and still terminates at a/b in lowest terms.
    let convs: Vec<(u64, u64)> = continued_fraction(2, 5)
        .convergents()
        .take(CONV_GUARD)
        .collect();
    assert_eq!(convs.first(), Some(&(0, 1)));
    assert_eq!(convs.last(), Some(&(2, 5)));
}

#[test]
#[allow(clippy::many_single_char_names)] // a, b, g, h, k, q, w — CF / convergent names
fn continued_fraction_reconstructs_and_determinant_wide() {
    let check = |a: u64, b: u64| {
        let cf = continued_fraction(a, b);
        let (oracle_qs, oracle_convs) = cf_oracle(a, b);
        // (1) quotients agree with the independent Euclid oracle.
        assert_eq!(cf.as_slice(), oracle_qs.as_slice(), "quotients a={a} b={b}");
        let convs: Vec<(u64, u64)> = cf.convergents().take(CONV_GUARD).collect();
        assert_eq!(convs.len(), oracle_convs.len(), "len a={a} b={b}");
        // (2) production convergents equal the u128-oracle convergents (no overflow).
        for (i, (&(h, k), &(oh, ok))) in convs.iter().zip(&oracle_convs).enumerate() {
            assert_eq!(
                (u128::from(h), u128::from(k)),
                (oh, ok),
                "conv {i} a={a} b={b}"
            );
        }
        let last = *convs.last().unwrap();
        if b >= 1 {
            // (3) reconstruction: the last convergent is a/b in lowest terms.
            let g = u128::from(gcd(a, b));
            assert_eq!(u128::from(last.0) * g, u128::from(a), "recon h a={a} b={b}");
            assert_eq!(u128::from(last.1) * g, u128::from(b), "recon k a={a} b={b}");
        }
        // (4) determinant identity h_k·k_{k-1} − h_{k-1}·k_k = ±1, in u128.
        for w in convs.windows(2) {
            let (h0, k0) = (u128::from(w[0].0), u128::from(w[0].1));
            let (h1, k1) = (u128::from(w[1].0), u128::from(w[1].1));
            assert_eq!((h1 * k0).abs_diff(h0 * k1), 1, "det identity a={a} b={b}");
        }
        // (5) best-approximation bound |a·k_k − b·h_k|·k_{k+1} <= b (for k < n).
        if b >= 1 {
            for w in convs.windows(2) {
                let (h, k) = (u128::from(w[0].0), u128::from(w[0].1));
                let knext = u128::from(w[1].1);
                let err = (u128::from(a) * k).abs_diff(u128::from(b) * h);
                assert!(err * knext <= u128::from(b), "approx bound a={a} b={b}");
            }
        }
        // (6) canonical structure: a_i >= 1 for i >= 1; denominators non-decreasing.
        for &q in &cf.as_slice()[1..] {
            assert!(q >= 1, "partial quotient >=1 a={a} b={b}");
        }
        for w in convs.windows(2) {
            assert!(w[1].1 >= w[0].1, "denominators non-decreasing a={a} b={b}");
        }
    };

    // Small exhaustive grid.
    for a in 0..120u64 {
        for b in 0..120u64 {
            check(a, b);
        }
    }
    // Consecutive Fibonacci pairs — the worst case for length (the largest
    // representable pair is checked before checked_add overflows).
    let (mut p, mut q) = (1u64, 1u64);
    loop {
        check(p, q);
        check(q, p);
        match q.checked_add(p) {
            Some(n) => (p, q) = (q, n),
            None => break,
        }
    }
    // 100k full-range pairs.
    let mut state = 0x1234_5678_9abc_def1u64;
    for _ in 0..100_000 {
        check(lcg(&mut state), lcg(&mut state));
    }
}

#[test]
fn continued_fraction_length_bound() {
    // No u64/u64 expansion exceeds 92 partial quotients (Lamé), and the worst case
    // (consecutive Fibonacci) reaches it exactly — so the cap of 93 never panics.
    let mut max_len = 0usize;
    let mut state = 0xfeed_face_dead_beefu64;
    for _ in 0..500_000 {
        let n = continued_fraction(lcg(&mut state), lcg(&mut state))
            .as_slice()
            .len();
        max_len = max_len.max(n);
    }
    let (mut p, mut q) = (1u64, 1u64);
    loop {
        // p <= q, so cf(p, q) is the a < b orientation: a leading 0 then the
        // all-ones tail — the 92-term worst case at the largest Fibonacci pair.
        max_len = max_len.max(continued_fraction(p, q).as_slice().len());
        match q.checked_add(p) {
            Some(n) => (p, q) = (q, n),
            None => break,
        }
    }
    assert_eq!(
        max_len, 92,
        "max CF length over u64 is 92 (Fibonacci worst case)"
    );
}

// ── Adversarial KATs above the 0..10^6 sieve / n<=500 brute ceilings ──────────
//
// is_prime_matches_sieve is exhaustive only over 0..10^6 and the all-u64 Feitsma
// sweep is #[ignore]d; carmichael's brute oracle reaches only n<=500. Above those
// ceilings coverage otherwise rests on theorem-entailment, so these known-answer
// tests pin the sharpest adversarial points directly. Constants independently
// verified.

#[test]
fn is_prime_rejects_psi8_9_champion_and_needs_base_37() {
    // 3825123056546413051 = 149491 · 747451 · 34233211 (< 2^64) is a strong
    // pseudoprime to the first NINE primes — in fact a strong *liar* to the first
    // ELEVEN witness bases {2..31}; only base 37 witnesses it. The single sharpest
    // point of coprime's twelve-base theorem, far above the 10^6 sieve.
    const PSI: u64 = 3_825_123_056_546_413_051;
    assert_eq!(
        149_491u64 * 747_451 * 34_233_211,
        PSI,
        "factorization of the champion"
    );
    assert!(!is_prime(PSI), "the psi8/9 champion is composite");

    // Base-37 necessity: a strong liar to every base below 37, a witness at 37.
    for a in [2u64, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] {
        assert!(
            is_strong_probable_prime(PSI, a),
            "base {a} must be a strong LIAR for the champion"
        );
    }
    assert!(
        !is_strong_probable_prime(PSI, 37),
        "base 37 must WITNESS the champion — the twelfth base is load-bearing"
    );
}

#[test]
fn carmichael_two_adic_branch_full_range() {
    // λ(2)=1, λ(4)=2, and λ(2^e)=2^(e-2) for e >= 3 — the `1 << (e-2)` branch,
    // which the n<=500 brute oracle only ever reaches at e<=8. Sweep e up to 63.
    assert_eq!(carmichael(2), Some(1));
    assert_eq!(carmichael(4), Some(2));
    for e in 3u32..=63 {
        assert_eq!(
            carmichael(1u64 << e),
            Some(1u64 << (e - 2)),
            "λ(2^{e}) = 2^(e-2)"
        );
    }
}

#[test]
fn is_prime_count_checkpoint_above_sieve_ceiling() {
    // π(2·10^6) − π(10^6) = 70435 (OEIS A006880). The segment [10^6, 2·10^6) is
    // disjoint from is_prime_matches_sieve's 0..10^6 domain — an independent
    // known-answer count over a region the everyday oracle never touches.
    let count = (1_000_000u64..2_000_000).filter(|&n| is_prime(n)).count();
    assert_eq!(count, 70_435, "primes in [10^6, 2·10^6)");
}

// ── Adversarial inputs (assume a bug; construct the input that finds it) ───────
//
// Each targets a NAMED weakness a subtly-wrong implementation could pass
// everything else and fail here. Every constant is derived independently of
// coprime (an independent reference / by construction).

#[test]
fn is_prime_adversarial_near_and_above_2_64() {
    // Kills: an MR squaring-loop bug on a Fermat-fooler above the n<2^32 composite
    // fuzz. Chernick Carmichael 1171·2341·3511 fools Fermat for every coprime base,
    // but base 2 strong-WITNESSES it (v2(n-1)=3 drives the squaring loop).
    assert_eq!(1_171u64 * 2_341 * 3_511, 9_624_742_921);
    assert!(!is_prime(9_624_742_921));
    assert!(!is_strong_probable_prime(9_624_742_921, 2));
    // Kills: a truncated MR squaring loop (false-negative on a prime needing a deep
    // descent). Prime in (2^63, 2^64) with v2(p-1)=28 — a deeper squaring descent
    // than any other prime KAT (which have v2 ≤ 2).
    assert_eq!((16_188_946_917_747_589_121u64 - 1) >> 28, 60_308_526_895);
    assert!(is_prime(16_188_946_917_747_589_121));
    assert!(is_strong_probable_prime(16_188_946_917_747_589_121, 2));
}

#[test]
fn factorize_primorial15_fills_capacity() {
    // Kills: an off-by-one in Factorization::push at the [(u64,u32);15] capacity.
    // The primorial of the first 15 primes is the unique smallest u64 with 15
    // distinct factors — it fills the array exactly (the 16th prime overflows 2^64).
    assert_eq!(
        factorize(614_889_782_588_491_410).as_slice(),
        &[
            (2, 1),
            (3, 1),
            (5, 1),
            (7, 1),
            (11, 1),
            (13, 1),
            (17, 1),
            (19, 1),
            (23, 1),
            (29, 1),
            (31, 1),
            (37, 1),
            (41, 1),
            (43, 1),
            (47, 1)
        ]
    );
}

#[test]
fn factorize_rho_prime_power_and_batch_collision() {
    // Kills: rho/coalesce bugs on factors above the d<1024 trial-division ceiling.
    // 1031^6 reaches Brent-rho on a prime POWER, six occurrences coalescing to e=6
    // (3^40 is peeled entirely by trial division, never reaching rho).
    assert_eq!(
        factorize(1_201_024_845_477_409_681).as_slice(),
        &[(1031, 6)]
    );
    // Kills: removal of brent()'s `g == n` batch-collision recovery branch — the
    // path that almost never runs; 1031·1039 drives it on the first rho constant.
    assert_eq!(factorize(1_071_209).as_slice(), &[(1031, 1), (1039, 1)]);
}

#[test]
fn modpow_u128_product_at_ceiling() {
    // Kills: a u64 (instead of u128) product in modpow — a bug the all-u128
    // differential oracles would SHARE and miss. base = 2^64-2, modulus = the
    // largest u64 prime: every square-and-multiply step nears 2^128. Constant from
    // an independent big-integer reference.
    assert_eq!(
        modpow(
            18_446_744_073_709_551_614,
            u64::MAX,
            18_446_744_073_709_551_557
        ),
        2_012_073_826_774_673_798
    );
}

#[test]
fn extended_gcd_tightest_i128_cofactor() {
    // Kills: i64 (instead of i128) cofactor accumulators. The tightest corner over
    // all u64 is NOT a Fibonacci pair (cofactors ~2^65) but the gcd-1 pair
    // (2^64-1, 2^64-3), whose |t·b| is 127-bit — the documented "tight at the
    // ceiling" bound, pinned deterministically.
    let (d, s, t) = extended_gcd(u64::MAX, u64::MAX - 2);
    assert_eq!(
        (d, s, t),
        (1, -9_223_372_036_854_775_806, 9_223_372_036_854_775_807)
    );
    assert_eq!(s * i128::from(u64::MAX) + t * i128::from(u64::MAX - 2), 1);
    assert_eq!(
        (t * i128::from(u64::MAX - 2)).unsigned_abs(),
        170_141_183_460_469_231_685_570_443_531_610_226_691
    );
}

#[test]
fn iroot_overflow_and_large_k_edges() {
    // Kills: a wrapping (not checked) (r+1)^k. k=3 at the largest cube in u64;
    // 2642246^3 overflows and must count as > n.
    assert_eq!(iroot(18_446_724_184_312_856_125, 3), 2_642_245);
    assert_eq!(iroot(18_446_724_184_312_856_124, 3), 2_642_244);
    assert_eq!(iroot(u64::MAX, 3), 2_642_245);
    // Kills: a large-k bracket/overflow bug (k > the k≤32 sampling never reaches).
    assert_eq!(iroot(1u64 << 63, 63), 2);
    assert_eq!(iroot((1u64 << 63) - 1, 63), 1);
    assert_eq!(iroot(u64::MAX, 64), 1);
    assert_eq!(iroot(u64::MAX, 63), 2);
}

#[test]
fn sqrtmod_deep_tonelli_near_2_64() {
    // Kills: a Tonelli–Shanks shift/depth bug at large 2-adic valuation. p near
    // 2^64, p ≡ 1 (mod 4), v2(p-1) = 40 — a deeper descent than the p < 2^32, s ≤ 30
    // regime any other deep-TS KAT reaches. Self-verifying (x² ≡ 5) plus the canonical root.
    let p = 18_446_742_974_197_923_841u64;
    assert_eq!(sqrtmod(5, p), Some(6_350_389_231_041_696_576));
    assert_eq!(
        6_350_389_231_041_696_576u128 * 6_350_389_231_041_696_576 % u128::from(p),
        5
    );
    assert_eq!(sqrtmod(19, p), None); // a quadratic non-residue mod p
}

#[test]
fn sqrtmod_rejects_composite_with_perfect_square() {
    // Kills: removal/weakening of the prime-modulus guard. n is composite and a is a
    // perfect square, so without the guard sqrtmod would spuriously find a root — a
    // non-square a fails closed anyway, so only a perfect-square a exercises the guard.
    let n = 18_446_743_979_220_271_189u64; // 4294967291 · 4294967279
    assert_eq!(sqrtmod(4, n), None);
    assert_eq!(sqrtmod(25, n), None);
}

#[test]
fn carmichael_lambda_is_attained() {
    // Kills: carmichael returning a MULTIPLE of λ. The universal-exponent check
    // (a^λ ≡ 1) passes for any multiple too; tightness needs an element of order
    // EXACTLY λ. n = 5²·7·11·13 has λ=60 ≪ φ=14400, above the n≤500 brute oracle.
    assert_eq!(carmichael(25_025), Some(60));
    assert_eq!(totient(25_025), 14_400);
    assert_eq!(order(10_727, 25_025), Some(60)); // λ attained by a concrete unit
    assert_ne!(modpow(10_727, 30, 25_025), 1); // order is not 60/2
    assert_ne!(modpow(10_727, 20, 25_025), 1); // not 60/3
    assert_ne!(modpow(10_727, 12, 25_025), 1); // not 60/5
}

#[test]
fn order_descent_over_rich_totient() {
    // Kills: a per-prime mistake in order's factor-and-divide descent over a highly
    // composite φ. p = M31, φ = p−1 = 2·3²·7·11·31·151·331; the order collapses via
    // order(gʲ) = order(g)/gcd(order(g), j) to independently-known small values.
    let p = 2_147_483_647u64;
    assert_eq!(order(7, p), Some(2_147_483_646)); // 7 is a primitive root (order = φ)
    assert_eq!(order(modpow(7, (p - 1) / 7, p), p), Some(7));
    assert_eq!(order(modpow(7, (p - 1) / 3, p), p), Some(3));
}
