Add an x86_64 Pasta field assembly backend

[ValarDragon]

What changed

This adds an opt-in pasta_curves/x86_64-asm feature for Pasta field
multiplication and squaring on x86-64.

• The normal Rust implementation remains the portable fallback.
• Generic x86-64 builds perform one cached CPUID check for ADX and BMI2, then
  use checked-in Montgomery multiplication and squaring assembly when both are
  available.
• Builds compiled explicitly with both target features call the assembly path
  directly.
• Other architectures and unsupported targets continue using portable Rust.

The feature is non-default and adds no C or C++ code.

Why

Field multiplication is a large part of Orchard proof generation. The ADX/BMI2
instructions provide two carry chains and a non-destructive multiply, which the
portable Rust arithmetic does not reliably expose to LLVM from its current
serial u128 dataflow.

On the dedicated Linux/x86-64 benchmark host, against exact main at
bcf22f33, a genuine one-Action Orchard proof at k = 11 improved from a
2.155071707 s geometric center to 2.068251146 s:

• 4.028662% lower proving latency
• 1.041978x throughput
• paired reductions of 3.891614% and 4.165515%
• control drift +0.028306%; candidate drift -0.256767%

The balanced A/B/B/A run used Rust 1.88, locked dependencies, byte-identical
harnesses, RAYON_NUM_THREADS=1, and CPU pinning. Every leg recorded zero CPU
steal, zero swap activity, and no competing benchmark/compiler process. Both
binaries generated and verified the proof before timing.

Implementation and provenance

The four object-format-specific assembly files are generated from Semolina
v0.1.4 (supranational/semolina@13ffc780…) under Apache-2.0. They retain the
upstream copyright and SPDX header and document the local symbol-prefixing and
routine-pruning adaptations.

The Rust boundary is private. It caches the public CPU-capability result in an
atomic byte and calls only two prefixed native entry points. The ELF symbols are
hidden, Mach-O symbols are private externs, and helpers are local/private where
the object format permits. The ELF object also marks a non-executable stack.

Validation

• Rust 1.88 Fp and Fq differentials on Linux for both runtime CPUID dispatch
  and compile-time +adx,+bmi2: 36 boundary pairs plus 1,024 deterministic
  random pairs per field, covering multiplication and squaring.
• Warning-denied Linux checks for both dispatch modes.
• Apple AArch64 portable-fallback checks with x86_64-asm enabled, including
  coexistence with aarch64-asm.
• Generated-object audit: expected hidden/local symbols, no undefined symbols,
  non-executable stack, and exact mulx/adcx/adox instruction census.
• cargo package --list, formatting, diff hygiene, and focused Orchard proof
  generation/verification.

Commit e386dce adds CI coverage only; the benchmarked production sources are
the byte-identical parent commit c91e923. CI explicitly covers generic and
static-feature Linux builds, Windows MSVC and GNU objects, Intel macOS, and the
portable fallback on Apple AArch64.

API surface

The only downstream-visible change is the new non-default Cargo feature
x86_64-asm. There are no new or changed Rust pub or pub(crate) items,
function signatures, trait methods, enum variants, type aliases, constants, or
constructors. The implementation adds two private prefixed native link symbols
and confines unsafe code to the private CPUID/assembly binding module.

Benchmark artifacts:
pasta-field-x86-asm-main-c91e923-benchmark-20260815.tar.gz, SHA-256
d2379f5bf01a142efc30da6e2ffbd5c2add4fa76835a4f8189747ea4a16adc30.

Current final-MSM stack benchmark (2026-08-19)

The x86 backend was replayed without source changes (apart from additive
changelog composition) over the accepted GLV/XYZZ/multicore/c=10 stack at PR
#93. The exact incremental production boundary was a6b10fa -> 12e8971.

Standard balanced A/B/B/A geometric centers of means on the same Xeon 8358
host were:

Workload	Portable	x86 ASM	Reduction
prover, one worker	1365.259 ms	1298.929 ms	4.858404%
verifier B1, one worker	25.202 ms	24.209 ms	3.940504%
verifier B2, one worker	26.992 ms	25.978 ms	3.755019%
verifier B16, one worker	49.401 ms	47.875 ms	3.088564%
verifier B64, one worker	126.243 ms	124.624 ms	1.281914%
verifier B64, eight workers	41.848 ms	40.555 ms	3.090010%

The serial B64 mean had 1.20% candidate drift; its median center was a cleaner
1.830772% reduction, with both paired directions positive. The eight-worker
B64 paired mean reductions were 3.584677% and 2.592805%.

This run used the authenticated 64-proof corpus, native Fp/Fq differentials,
genuine proof preflights, Rust 1.88, locked dependencies, and the standard
telemetry gates. The combined x86-backend plus #56 specialized-square ratios
are 5.956128% for proving, 2.239169% for serial B64 verification, and
4.165871% for eight-worker B64 verification over PR #93.

Current-stack archive SHA-256:
f35c6c98d9e1aa88f2280c4310ebef573b79ae6f621f59be0f47bd306013521f.

[v12-auditor]

Note

Complete: Audit complete. V12 found two issues worth reviewing.

Open the full results here.

Finding	Severity	Details
F-228964	🟡 Medium	MIT option excludes vendored assembly All four newly vendored Semolina-derived assembly files carry Supranational LLC copyright and an Apache-2.0-only SPDX identifier. The package manifest nevertheless advertises the complete crate as MIT OR Apache-2.0, and both packaged licensing documents tell recipients they may choose either license. The bundled MIT license is granted by The Electric Coin Company and does not grant MIT terms over Supranational's files. Because the package has no include or exclude restriction, Cargo source packages contain these files independently of whether the non-default feature is enabled. Per-file headers disclose the narrower terms to a reader but do not correct the package-level expression consumed by registries, SBOM generators, and license scanners.
F-228965	🔵 Low	MSVC cross-builds require unavailable MASM The build selector maps every x86-64 Windows/MSVC target to the MASM-syntax pasta_mulx-x86_64-win64.asm solely from target properties, then passes that file to cc. The resolved cc dependency treats an .asm input for an MSVC target as Microsoft macro-assembler input and invokes ml64.exe; configuring a standard Linux/macOS MSVC cross-toolchain around clang-cl does not replace that command. The selected source confirms its MASM-only dialect through OPTION, SEGMENT, PROC, and ENDP directives. Standard cargo-xwin/Clang cross-builders do not provide ml64.exe, so feature-enabled MSVC cross-builds terminate in the build script even though the repository already includes a GAS-syntax COFF implementation with the same Win64 ABI adaptation.

And two more auto-invalidated findings.

Analyzed nine files, diff bcf22f3...e386dce.

[TalDerei]

left some observations, cc @ebfull for visibility.

[TalDerei]

ACK runtime CPUID dispatch which is what i'd expect for ADX/BMI2. i'm on aarch64-apple-darwin so couldn't actually execute assembly and compare against rust without a docker env, so dispatched agents to do a static analysis instead.

[TalDerei]

actually also did differential executions between assembly / rust on my linux box with ryzen cpu that has both ADX/BMI2. i know x86_64_asm_matches_portable_arithmetic already runs in CI, but wanted to scrutinize it more.

[TalDerei]

few CI commands still use pasta_curves, but current main renamed it to zakura-pasta-curves.

[TalDerei]

oh supranational's code and the assembly was written before AI :)

[TalDerei]

something i noticed; CI does't deterministically execute ADX/BMI2. CI compiles the assembly path, but runtime tests can detect that the runner lacks those CPU instructions and use portable rust instead, so like green CI isn't actually a guarantee the assembly was executed. this is a test coverage gap imo.

[TalDerei]

Tested on my linux box:

The variants match — everywhere I could push them

I checked out PR #53 into a worktree (/home/margulus/zcash/wt-pr53, branch pr-53), confirmed this CPU has ADX+BMI2 so the assembly path genuinely executes, and ran the assembly and portable Rust against each other far past the PR's own tests. Zero divergence: every result was bit-for-bit identical, across runtime CPUID dispatch and compile-time +adx,+bmi2 builds, debug and release, on rustc 1.90 (the default 1.86 toolchain is below the crate's 1.88 MSRV). I also verified the linked objects: symbols are local/hidden as claimed and the stack is non-executable.

Scrutiny of x86_64_asm_matches_portable_arithmetic

The test is honestly constructed (it correctly uses the inherent Fp::mul/Fp::square as the portable reference, and raw-limb equality makes the comparison bit-exact), but it has real gaps:

1. The big one — untested non-canonical inputs that production actually uses. from_u512 (behind Field::random and hash_to_curve) computes d0 * R2 + d1 * R3 through the dispatched operator, where d0/d1 are raw 512-bit-hash halves that can be ≥ p, up to 2²⁵⁶−1. The asm binding's doc says it takes "canonical Montgomery residues", and every input in the PR's test goes through from_raw (which reduces), so this reachable domain was never covered. It happens to be mathematically fine — with one factor < p the Montgomery intermediate stays below 2p, and the assembly keeps a fifth accumulator limb with a full 5-limb conditional subtract — and my new tests now prove it empirically, but the doc comment on x86_64_asm::mul understates the real precondition and should be fixed.
2. A boundary value is a duplicate. Fp::from_raw([1, 0, 0, 0]) converts to Montgomery form, so it is Fp::one() — the boundary set has 5 distinct values, not 6, and the raw limb pattern [1,0,0,0] (field element R⁻¹) was never tested. Likely the author meant Fp([1, 0, 0, 0]).
3. Randoms are uniform-only. 1,024 from_raw samples don't target carry-chain stressors: all-ones limbs, single bits at word boundaries, p±1, INV, 32-bit half-words, R2±1, (p±1)/2, etc.
4. Differential-only. With no implementation-independent oracle, a bug shared by both variants is invisible in principle.
5. No asm self-consistency checks — commutativity, and sqrx vs mulx(a,a) (which enter the shared reduction through different stubs) — and no iterated feedback of asm outputs into asm inputs.

New adversarial coverage (all passing)

I added a x86_64_asm_edge_cases module to both fp.rs and fq.rs (~440 lines each, identical apart from the field alias and fuzz seed) — 5 new tests per field, roughly 79,000 assembly invocations per field:

• Edge-pattern grid: 215 raw limb patterns (edge words, single bits ±1 at every word boundary, offsets of p, R, R², R³, (p±1)/2, 2p±1…) → 209 distinct canonical values, all-pairs multiplied (21,945 pairs, both argument orders, plus the dispatched operator) and squared.
• Non-canonical domain: the 97 patterns at/above p multiplied against R², R³ and other canonical values in both orders — exactly the from_u512 envelope — every result oracle-verified; plus an end-to-end from_u512 check against a portable-only reconstruction.
• Independent oracle: a bit-serial shift-and-subtract checker asserting canonicity and the defining property out·2²⁵⁶ ≡ a·b (mod p), sharing no algorithm with either implementation.
• Structured fuzzing: 2,048 canonical pairs + 512 extended pairs with limbs biased toward edge words and sparse/dense bit soups.
• Iterated feedback: an 8,192-step walk feeding outputs back as inputs so any divergence compounds.

The full 44-test suite passes in every mode, cargo fmt --check is clean, clippy adds zero new warnings (36 pre-existing upstream lints, unchanged), and warning-denied builds succeed.

Committed as 570f351 on local branch pr-53 — not pushed. Two caveats: I set a repo-local git identity Tal Bowe <tal@bowe.tech> (amend if you'd like different details), and the commit is unsigned because your ~/.ssh/tachyon signing key needs a passphrase I can't enter — git commit --amend -S will fix that. If you want, I can push it to the PR branch, and/or prepare the small follow-up fixing the x86_64_asm::mul doc contract and the duplicate boundary value in the original test.

git am-able patch (2 files, +880 lines)

From 570f351f91beab2799c3f57e6cff3f77ed79bfbf Mon Sep 17 00:00:00 2001
From: Tal Bowe <tal@bowe.tech>
Date: Sun, 16 Aug 2026 11:41:34 -0700
Subject: [PATCH] Add adversarial differential tests for the x86_64 field
 assembly

Extend the x86_64-asm coverage in Fp and Fq with a shared edge-case
module that checks the assembly, the portable arithmetic, and an
independent bit-serial oracle against each other:

- an all-pairs grid over ~209 canonicalized limb patterns built from
  carry-stressing words, word-boundary bits, and offsets of p, R, R2,
  and R3, in both argument orders plus the dispatched operator;
- direct coverage of the non-canonical operands (values at and above p)
  that from_u512 feeds the assembly through d0 * R2 + d1 * R3, plus an
  end-to-end from_u512 reconstruction check;
- structured-random fuzzing biased toward edge words and sparse or
  dense bit soups, and an 8192-step iterated-feedback walk;
- a shift-and-subtract oracle asserting canonicity and the defining
  Montgomery property out * 2^256 = a * b (mod p) without sharing any
  algorithm with either implementation.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 pasta_curves/src/fields/fp.rs | 440 ++++++++++++++++++++++++++++++++++
 pasta_curves/src/fields/fq.rs | 440 ++++++++++++++++++++++++++++++++++
 2 files changed, 880 insertions(+)

diff --git a/pasta_curves/src/fields/fp.rs b/pasta_curves/src/fields/fp.rs
index 869ae22..2825f31 100644
--- a/pasta_curves/src/fields/fp.rs
+++ b/pasta_curves/src/fields/fp.rs
@@ -1060,6 +1060,446 @@ fn x86_64_asm_matches_portable_arithmetic() {
     }
 }
 
+#[cfg(all(
+    test,
+    feature = "x86_64-asm",
+    target_arch = "x86_64",
+    any(unix, windows)
+))]
+mod x86_64_asm_edge_cases {
+    //! Adversarial differential coverage for the x86_64 assembly backend.
+    //!
+    //! Every check compares three implementations: the assembly, the portable
+    //! inherent arithmetic, and a bit-serial shift-and-subtract oracle that
+    //! shares no algorithm with either. The oracle verifies the defining
+    //! Montgomery property `out * 2^256 = a * b (mod p)` plus canonicity
+    //! `out < p`.
+
+    use super::Fp as F;
+    use super::{INV, MODULUS, R, R2, R3};
+    use crate::fields::x86_64_asm;
+
+    use ff::Field;
+    use rand::{RngCore, SeedableRng};
+    use rand_xorshift::XorShiftRng;
+    use std::vec::Vec;
+
+    type Limbs = [u64; 4];
+
+    const P: Limbs = MODULUS.0;
+
+    /// Returns whether `a < b` as 256-bit little-endian integers.
+    fn less_than(a: &Limbs, b: &Limbs) -> bool {
+        for (a, b) in a.iter().zip(b.iter()).rev() {
+            if a != b {
+                return a < b;
+            }
+        }
+        false
+    }
+
+    /// Subtracts the modulus from `a`, which must not underflow.
+    fn sub_p(a: &Limbs) -> Limbs {
+        let mut out = [0; 4];
+        let mut borrow = false;
+        for ((out, a), p) in out.iter_mut().zip(a).zip(&P) {
+            let (diff, underflow) = a.overflowing_sub(*p);
+            let (diff, chained) = diff.overflowing_sub(borrow as u64);
+            *out = diff;
+            borrow = underflow | chained;
+        }
+        assert!(!borrow, "oracle subtraction underflowed");
+        out
+    }
+
+    /// Reduces a little-endian 512-bit value modulo `p`, one bit at a time.
+    fn oracle_mod_p(wide: &[u64; 8]) -> Limbs {
+        let mut acc = [0u64; 4];
+        for bit in (0..512).rev() {
+            let mut carry = (wide[bit / 64] >> (bit % 64)) & 1;
+            for limb in acc.iter_mut() {
+                let shifted = (*limb << 1) | carry;
+                carry = *limb >> 63;
+                *limb = shifted;
+            }
+            assert_eq!(carry, 0, "oracle accumulator overflowed");
+            if !less_than(&acc, &P) {
+                acc = sub_p(&acc);
+            }
+        }
+        acc
+    }
+
+    /// Schoolbook 256x256 -> 512-bit multiplication.
+    fn oracle_wide_mul(a: &Limbs, b: &Limbs) -> [u64; 8] {
+        let mut wide = [0u64; 8];
+        for (i, a) in a.iter().enumerate() {
+            let mut carry = 0u128;
+            for (j, b) in b.iter().enumerate() {
+                let t = u128::from(wide[i + j]) + u128::from(*a) * u128::from(*b) + carry;
+                wide[i + j] = t as u64;
+                carry = t >> 64;
+            }
+            wide[i + 4] = carry as u64;
+        }
+        wide
+    }
+
+    /// Asserts that `out` is canonical and `out * 2^256 = a * b (mod p)`,
+    /// the defining property of a Montgomery product.
+    fn assert_montgomery_product(out: &Limbs, a: &Limbs, b: &Limbs, what: &str) {
+        assert!(less_than(out, &P), "{what}: {out:016x?} is not canonical");
+        let mut shifted = [0u64; 8];
+        shifted[4..].copy_from_slice(out);
+        assert_eq!(
+            oracle_mod_p(&shifted),
+            oracle_mod_p(&oracle_wide_mul(a, b)),
+            "{what}: {a:016x?} * {b:016x?} produced {out:016x?}"
+        );
+    }
+
+    fn add_one(mut v: Limbs) -> Limbs {
+        for limb in v.iter_mut() {
+            let (sum, carry) = limb.overflowing_add(1);
+            *limb = sum;
+            if !carry {
+                return v;
+            }
+        }
+        panic!("pattern overflowed 2^256");
+    }
+
+    fn sub_one(mut v: Limbs) -> Limbs {
+        for limb in v.iter_mut() {
+            let (diff, borrow) = limb.overflowing_sub(1);
+            *limb = diff;
+            if !borrow {
+                return v;
+            }
+        }
+        panic!("pattern underflowed zero");
+    }
+
+    fn halve(v: Limbs) -> Limbs {
+        [
+            (v[0] >> 1) | (v[1] << 63),
+            (v[1] >> 1) | (v[2] << 63),
+            (v[2] >> 1) | (v[3] << 63),
+            v[3] >> 1,
+        ]
+    }
+
+    fn double(v: Limbs) -> Limbs {
+        assert_eq!(v[3] >> 63, 0, "pattern overflowed 2^256");
+        [
+            v[0] << 1,
+            (v[1] << 1) | (v[0] >> 63),
+            (v[2] << 1) | (v[1] >> 63),
+            (v[3] << 1) | (v[2] >> 63),
+        ]
+    }
+
+    /// Maps any 256-bit value to its canonical residue; every 256-bit value
+    /// is below 4p because p exceeds 2^254.
+    fn canonicalize(mut v: Limbs) -> Limbs {
+        for _ in 0..3 {
+            if !less_than(&v, &P) {
+                v = sub_p(&v);
+            }
+        }
+        assert!(less_than(&v, &P));
+        v
+    }
+
+    /// Words chosen to stress the mulx partial products and both ADX carry
+    /// chains: all-zero and all-one words, bits at word boundaries,
+    /// alternating patterns, 32-bit halves, and the words of the modulus and
+    /// of the Montgomery constant.
+    fn edge_words() -> [u64; 20] {
+        [
+            0,
+            1,
+            2,
+            0x0000_0000_ffff_ffff,
+            0xffff_ffff_0000_0000,
+            0x5555_5555_5555_5555,
+            0xaaaa_aaaa_aaaa_aaaa,
+            (1 << 62) - 1,
+            1 << 62,
+            (1 << 63) - 1,
+            1 << 63,
+            u64::MAX - 1,
+            u64::MAX,
+            INV,
+            P[0] - 1,
+            P[0],
+            P[0] + 1,
+            P[1] - 1,
+            P[1],
+            P[1] + 1,
+        ]
+    }
+
+    /// Raw 256-bit input patterns, deliberately including values at and above
+    /// the modulus: edge-word grids, single bits and their neighbours at
+    /// every word boundary, and offsets of the modulus and the Montgomery
+    /// constants.
+    fn edge_patterns() -> Vec<Limbs> {
+        fn push(patterns: &mut Vec<Limbs>, value: Limbs) {
+            if !patterns.contains(&value) {
+                patterns.push(value);
+            }
+        }
+
+        let mut patterns = Vec::new();
+        for &word in edge_words().iter() {
+            push(&mut patterns, [word; 4]);
+            for position in 0..4 {
+                let mut sparse = [0; 4];
+                sparse[position] = word;
+                push(&mut patterns, sparse);
+                let mut dense = [u64::MAX; 4];
+                dense[position] = word;
+                push(&mut patterns, dense);
+            }
+        }
+        for bit in [
+            0usize, 1, 32, 63, 64, 65, 127, 128, 129, 191, 192, 193, 252, 253, 254, 255,
+        ] {
+            let mut power = [0u64; 4];
+            power[bit / 64] = 1 << (bit % 64);
+            push(&mut patterns, sub_one(power));
+            push(&mut patterns, power);
+            push(&mut patterns, add_one(power));
+        }
+        for constant in [P, R.0, R2.0, R3.0] {
+            push(&mut patterns, sub_one(constant));
+            push(&mut patterns, constant);
+            push(&mut patterns, add_one(constant));
+        }
+        push(&mut patterns, sub_one(sub_one(P)));
+        push(&mut patterns, halve(sub_one(P)));
+        push(&mut patterns, add_one(halve(sub_one(P))));
+        push(&mut patterns, sub_one(double(P)));
+        push(&mut patterns, double(P));
+        push(&mut patterns, add_one(double(P)));
+        patterns
+    }
+
+    /// Differentially checks one canonical multiplication pair through the
+    /// raw assembly bindings (in both argument orders) and the dispatched
+    /// operator, against the portable implementation.
+    fn check_mul_pair(a: Limbs, b: Limbs, oracle: bool) {
+        let portable = F::mul(&F(a), &F(b)).0;
+        assert_eq!(
+            (&F(a) * &F(b)).0,
+            portable,
+            "dispatched mul disagrees with portable: {a:016x?} * {b:016x?}"
+        );
+        if oracle {
+            assert_montgomery_product(&portable, &a, &b, "portable mul");
+        }
+        if x86_64_asm::is_available() {
+            let asm_ab = x86_64_asm::mul(&a, &b, &P, INV).unwrap();
+            let asm_ba = x86_64_asm::mul(&b, &a, &P, INV).unwrap();
+            assert_eq!(asm_ab, portable, "asm mul: {a:016x?} * {b:016x?}");
+            assert_eq!(asm_ba, portable, "asm mul swapped: {b:016x?} * {a:016x?}");
+        } else {
+            assert_eq!(x86_64_asm::mul(&a, &b, &P, INV), None);
+        }
+    }
+
+    /// Differentially checks one canonical squaring, including agreement
+    /// between the assembly square and the assembly self-multiplication,
+    /// which enter the shared reduction through different stubs.
+    fn check_square(a: Limbs, oracle: bool) {
+        let portable = F::square(&F(a)).0;
+        assert_eq!(
+            <F as Field>::square(&F(a)).0,
+            portable,
+            "dispatched square disagrees with portable: {a:016x?}"
+        );
+        if oracle {
+            assert_montgomery_product(&portable, &a, &a, "portable square");
+        }
+        if x86_64_asm::is_available() {
+            let asm_square = x86_64_asm::square(&a, &P, INV).unwrap();
+            let asm_self_mul = x86_64_asm::mul(&a, &a, &P, INV).unwrap();
+            assert_eq!(asm_square, portable, "asm square: {a:016x?}");
+            assert_eq!(
+                asm_self_mul, asm_square,
+                "asm mul(a, a) != asm square(a): {a:016x?}"
+            );
+        } else {
+            assert_eq!(x86_64_asm::square(&a, &P, INV), None);
+        }
+    }
+
+    /// Checks the extended domain that `from_u512` reaches in production:
+    /// one operand is an arbitrary 256-bit value, the other a canonical
+    /// residue. With one factor below p the Montgomery intermediate stays
+    /// below 2p, so the portable result is exact and the assembly must match
+    /// it bit-for-bit in both argument orders.
+    fn check_noncanonical_pair(raw: Limbs, canonical: Limbs) {
+        assert!(less_than(&canonical, &P));
+        let portable = F::mul(&F(raw), &F(canonical)).0;
+        assert_montgomery_product(&portable, &raw, &canonical, "portable mul, extended domain");
+        assert_eq!(
+            (&F(raw) * &F(canonical)).0,
+            portable,
+            "dispatched mul disagrees with portable: {raw:016x?} * {canonical:016x?}"
+        );
+        if x86_64_asm::is_available() {
+            let asm_ab = x86_64_asm::mul(&raw, &canonical, &P, INV).unwrap();
+            let asm_ba = x86_64_asm::mul(&canonical, &raw, &P, INV).unwrap();
+            assert_eq!(
+                asm_ab, portable,
+                "asm mul, non-canonical lhs: {raw:016x?} * {canonical:016x?}"
+            );
+            assert_eq!(
+                asm_ba, portable,
+                "asm mul, non-canonical rhs: {canonical:016x?} * {raw:016x?}"
+            );
+        }
+    }
+
+    /// All-pairs differential grid over canonicalized edge patterns. The
+    /// bit-serial oracle additionally spot-checks the grid and every square.
+    #[test]
+    fn x86_64_asm_matches_portable_on_edge_patterns() {
+        let mut values: Vec<Limbs> = Vec::new();
+        for pattern in edge_patterns() {
+            let value = canonicalize(pattern);
+            if !values.contains(&value) {
+                values.push(value);
+            }
+        }
+
+        for (i, &a) in values.iter().enumerate() {
+            check_square(a, true);
+            for (j, &b) in values.iter().enumerate().skip(i) {
+                check_mul_pair(a, b, (i + j) % 29 == 0);
+            }
+        }
+    }
+
+    /// `from_u512` (hash-to-field, `Field::random`) multiplies raw 256-bit
+    /// halves by R2 and R3 through the dispatched operator, so the assembly
+    /// receives operands at and above the modulus. Cover that domain
+    /// directly.
+    #[test]
+    fn x86_64_asm_matches_portable_on_noncanonical_inputs() {
+        let multipliers = [
+            R2.0,
+            R3.0,
+            [1, 0, 0, 0],
+            R.0,
+            sub_one(P),
+            canonicalize([u64::MAX; 4]),
+        ];
+        for pattern in edge_patterns() {
+            for &multiplier in &multipliers {
+                check_noncanonical_pair(pattern, multiplier);
+            }
+        }
+    }
+
+    /// Pins the whole `from_u512` conversion against a reconstruction that
+    /// uses only the portable inherent operations.
+    #[test]
+    fn x86_64_asm_from_u512_matches_portable_reconstruction() {
+        fn check(lo: Limbs, hi: Limbs) {
+            let dispatched = F::from_u512([lo[0], lo[1], lo[2], lo[3], hi[0], hi[1], hi[2], hi[3]]);
+            let portable = F::mul(&F(lo), &R2).add(&F::mul(&F(hi), &R3));
+            assert_eq!(
+                dispatched.0, portable.0,
+                "from_u512 diverged: lo={lo:016x?} hi={hi:016x?}"
+            );
+        }
+
+        for &pattern in edge_patterns().iter() {
+            check(pattern, [0; 4]);
+            check([0; 4], pattern);
+            check(pattern, pattern);
+            check(pattern, [u64::MAX; 4]);
+            check([u64::MAX; 4], pattern);
+        }
+    }
+
+    /// Biased fuzzing: words drawn from the edge set, sparse and dense random
+    /// bit soups, and uniform words, reduced when canonicity is required.
+    #[test]
+    fn x86_64_asm_matches_portable_on_structured_random_patterns() {
+        fn draw_word(words: &[u64], rng: &mut XorShiftRng) -> u64 {
+            match rng.next_u32() % 4 {
+                0 => words[rng.next_u32() as usize % words.len()],
+                1 => rng.next_u64() & rng.next_u64() & rng.next_u64(),
+                2 => rng.next_u64() | rng.next_u64() | rng.next_u64(),
+                _ => rng.next_u64(),
+            }
+        }
+
+        fn draw(words: &[u64], rng: &mut XorShiftRng) -> Limbs {
+            [
+                draw_word(words, rng),
+                draw_word(words, rng),
+                draw_word(words, rng),
+                draw_word(words, rng),
+            ]
+        }
+
+        let words = edge_words();
+        let mut rng = XorShiftRng::from_seed([0xa5; 16]);
+        for round in 0..2048 {
+            let a = canonicalize(draw(&words, &mut rng));
+            let b = canonicalize(draw(&words, &mut rng));
+            check_mul_pair(a, b, round % 16 == 0);
+            check_square(a, round % 16 == 0);
+        }
+        for _ in 0..512 {
+            let raw = draw(&words, &mut rng);
+            let canonical = canonicalize(draw(&words, &mut rng));
+            check_noncanonical_pair(raw, canonical);
+        }
+    }
+
+    /// Feeds assembly outputs back into assembly inputs for thousands of
+    /// steps so any divergence compounds, nudging the trajectory regularly
+    /// so it keeps crossing the final conditional-subtraction boundary.
+    #[test]
+    fn x86_64_asm_matches_portable_on_iterated_feedback() {
+        if !x86_64_asm::is_available() {
+            return;
+        }
+
+        let mut x = R2.0;
+        let mut y = sub_one(P);
+        for step in 0..8192u64 {
+            let product = F::mul(&F(x), &F(y)).0;
+            let square = F::square(&F(y)).0;
+            assert_eq!(
+                x86_64_asm::mul(&x, &y, &P, INV).unwrap(),
+                product,
+                "asm mul diverged at step {step}: x={x:016x?} y={y:016x?}"
+            );
+            assert_eq!(
+                x86_64_asm::square(&y, &P, INV).unwrap(),
+                square,
+                "asm square diverged at step {step}: y={y:016x?}"
+            );
+            if step % 64 == 0 {
+                assert_montgomery_product(&product, &x, &y, "walk mul");
+                assert_montgomery_product(&square, &y, &y, "walk square");
+            }
+            x = product;
+            y = square;
+            if step % 17 == 0 {
+                y = F(y).add(&F([step + 1, 0, 0, 0])).0;
+            }
+        }
+    }
+}
+
 #[test]
 fn test_inv() {
     // Compute -(r^{-1} mod 2^64) mod 2^64 by exponentiating
diff --git a/pasta_curves/src/fields/fq.rs b/pasta_curves/src/fields/fq.rs
index f6c8cf8..dba3ca9 100644
--- a/pasta_curves/src/fields/fq.rs
+++ b/pasta_curves/src/fields/fq.rs
@@ -1059,6 +1059,446 @@ fn x86_64_asm_matches_portable_arithmetic() {
     }
 }
 
+#[cfg(all(
+    test,
+    feature = "x86_64-asm",
+    target_arch = "x86_64",
+    any(unix, windows)
+))]
+mod x86_64_asm_edge_cases {
+    //! Adversarial differential coverage for the x86_64 assembly backend.
+    //!
+    //! Every check compares three implementations: the assembly, the portable
+    //! inherent arithmetic, and a bit-serial shift-and-subtract oracle that
+    //! shares no algorithm with either. The oracle verifies the defining
+    //! Montgomery property `out * 2^256 = a * b (mod p)` plus canonicity
+    //! `out < p`.
+
+    use super::Fq as F;
+    use super::{INV, MODULUS, R, R2, R3};
+    use crate::fields::x86_64_asm;
+
+    use ff::Field;
+    use rand::{RngCore, SeedableRng};
+    use rand_xorshift::XorShiftRng;
+    use std::vec::Vec;
+
+    type Limbs = [u64; 4];
+
+    const P: Limbs = MODULUS.0;
+
+    /// Returns whether `a < b` as 256-bit little-endian integers.
+    fn less_than(a: &Limbs, b: &Limbs) -> bool {
+        for (a, b) in a.iter().zip(b.iter()).rev() {
+            if a != b {
+                return a < b;
+            }
+        }
+        false
+    }
+
+    /// Subtracts the modulus from `a`, which must not underflow.
+    fn sub_p(a: &Limbs) -> Limbs {
+        let mut out = [0; 4];
+        let mut borrow = false;
+        for ((out, a), p) in out.iter_mut().zip(a).zip(&P) {
+            let (diff, underflow) = a.overflowing_sub(*p);
+            let (diff, chained) = diff.overflowing_sub(borrow as u64);
+            *out = diff;
+            borrow = underflow | chained;
+        }
+        assert!(!borrow, "oracle subtraction underflowed");
+        out
+    }
+
+    /// Reduces a little-endian 512-bit value modulo `p`, one bit at a time.
+    fn oracle_mod_p(wide: &[u64; 8]) -> Limbs {
+        let mut acc = [0u64; 4];
+        for bit in (0..512).rev() {
+            let mut carry = (wide[bit / 64] >> (bit % 64)) & 1;
+            for limb in acc.iter_mut() {
+                let shifted = (*limb << 1) | carry;
+                carry = *limb >> 63;
+                *limb = shifted;
+            }
+            assert_eq!(carry, 0, "oracle accumulator overflowed");
+            if !less_than(&acc, &P) {
+                acc = sub_p(&acc);
+            }
+        }
+        acc
+    }
+
+    /// Schoolbook 256x256 -> 512-bit multiplication.
+    fn oracle_wide_mul(a: &Limbs, b: &Limbs) -> [u64; 8] {
+        let mut wide = [0u64; 8];
+        for (i, a) in a.iter().enumerate() {
+            let mut carry = 0u128;
+            for (j, b) in b.iter().enumerate() {
+                let t = u128::from(wide[i + j]) + u128::from(*a) * u128::from(*b) + carry;
+                wide[i + j] = t as u64;
+                carry = t >> 64;
+            }
+            wide[i + 4] = carry as u64;
+        }
+        wide
+    }
+
+    /// Asserts that `out` is canonical and `out * 2^256 = a * b (mod p)`,
+    /// the defining property of a Montgomery product.
+    fn assert_montgomery_product(out: &Limbs, a: &Limbs, b: &Limbs, what: &str) {
+        assert!(less_than(out, &P), "{what}: {out:016x?} is not canonical");
+        let mut shifted = [0u64; 8];
+        shifted[4..].copy_from_slice(out);
+        assert_eq!(
+            oracle_mod_p(&shifted),
+            oracle_mod_p(&oracle_wide_mul(a, b)),
+            "{what}: {a:016x?} * {b:016x?} produced {out:016x?}"
+        );
+    }
+
+    fn add_one(mut v: Limbs) -> Limbs {
+        for limb in v.iter_mut() {
+            let (sum, carry) = limb.overflowing_add(1);
+            *limb = sum;
+            if !carry {
+                return v;
+            }
+        }
+        panic!("pattern overflowed 2^256");
+    }
+
+    fn sub_one(mut v: Limbs) -> Limbs {
+        for limb in v.iter_mut() {
+            let (diff, borrow) = limb.overflowing_sub(1);
+            *limb = diff;
+            if !borrow {
+                return v;
+            }
+        }
+        panic!("pattern underflowed zero");
+    }
+
+    fn halve(v: Limbs) -> Limbs {
+        [
+            (v[0] >> 1) | (v[1] << 63),
+            (v[1] >> 1) | (v[2] << 63),
+            (v[2] >> 1) | (v[3] << 63),
+            v[3] >> 1,
+        ]
+    }
+
+    fn double(v: Limbs) -> Limbs {
+        assert_eq!(v[3] >> 63, 0, "pattern overflowed 2^256");
+        [
+            v[0] << 1,
+            (v[1] << 1) | (v[0] >> 63),
+            (v[2] << 1) | (v[1] >> 63),
+            (v[3] << 1) | (v[2] >> 63),
+        ]
+    }
+
+    /// Maps any 256-bit value to its canonical residue; every 256-bit value
+    /// is below 4p because p exceeds 2^254.
+    fn canonicalize(mut v: Limbs) -> Limbs {
+        for _ in 0..3 {
+            if !less_than(&v, &P) {
+                v = sub_p(&v);
+            }
+        }
+        assert!(less_than(&v, &P));
+        v
+    }
+
+    /// Words chosen to stress the mulx partial products and both ADX carry
+    /// chains: all-zero and all-one words, bits at word boundaries,
+    /// alternating patterns, 32-bit halves, and the words of the modulus and
+    /// of the Montgomery constant.
+    fn edge_words() -> [u64; 20] {
+        [
+            0,
+            1,
+            2,
+            0x0000_0000_ffff_ffff,
+            0xffff_ffff_0000_0000,
+            0x5555_5555_5555_5555,
+            0xaaaa_aaaa_aaaa_aaaa,
+            (1 << 62) - 1,
+            1 << 62,
+            (1 << 63) - 1,
+            1 << 63,
+            u64::MAX - 1,
+            u64::MAX,
+            INV,
+            P[0] - 1,
+            P[0],
+            P[0] + 1,
+            P[1] - 1,
+            P[1],
+            P[1] + 1,
+        ]
+    }
+
+    /// Raw 256-bit input patterns, deliberately including values at and above
+    /// the modulus: edge-word grids, single bits and their neighbours at
+    /// every word boundary, and offsets of the modulus and the Montgomery
+    /// constants.
+    fn edge_patterns() -> Vec<Limbs> {
+        fn push(patterns: &mut Vec<Limbs>, value: Limbs) {
+            if !patterns.contains(&value) {
+                patterns.push(value);
+            }
+        }
+
+        let mut patterns = Vec::new();
+        for &word in edge_words().iter() {
+            push(&mut patterns, [word; 4]);
+            for position in 0..4 {
+                let mut sparse = [0; 4];
+                sparse[position] = word;
+                push(&mut patterns, sparse);
+                let mut dense = [u64::MAX; 4];
+                dense[position] = word;
+                push(&mut patterns, dense);
+            }
+        }
+        for bit in [
+            0usize, 1, 32, 63, 64, 65, 127, 128, 129, 191, 192, 193, 252, 253, 254, 255,
+        ] {
+            let mut power = [0u64; 4];
+            power[bit / 64] = 1 << (bit % 64);
+            push(&mut patterns, sub_one(power));
+            push(&mut patterns, power);
+            push(&mut patterns, add_one(power));
+        }
+        for constant in [P, R.0, R2.0, R3.0] {
+            push(&mut patterns, sub_one(constant));
+            push(&mut patterns, constant);
+            push(&mut patterns, add_one(constant));
+        }
+        push(&mut patterns, sub_one(sub_one(P)));
+        push(&mut patterns, halve(sub_one(P)));
+        push(&mut patterns, add_one(halve(sub_one(P))));
+        push(&mut patterns, sub_one(double(P)));
+        push(&mut patterns, double(P));
+        push(&mut patterns, add_one(double(P)));
+        patterns
+    }
+
+    /// Differentially checks one canonical multiplication pair through the
+    /// raw assembly bindings (in both argument orders) and the dispatched
+    /// operator, against the portable implementation.
+    fn check_mul_pair(a: Limbs, b: Limbs, oracle: bool) {
+        let portable = F::mul(&F(a), &F(b)).0;
+        assert_eq!(
+            (&F(a) * &F(b)).0,
+            portable,
+            "dispatched mul disagrees with portable: {a:016x?} * {b:016x?}"
+        );
+        if oracle {
+            assert_montgomery_product(&portable, &a, &b, "portable mul");
+        }
+        if x86_64_asm::is_available() {
+            let asm_ab = x86_64_asm::mul(&a, &b, &P, INV).unwrap();
+            let asm_ba = x86_64_asm::mul(&b, &a, &P, INV).unwrap();
+            assert_eq!(asm_ab, portable, "asm mul: {a:016x?} * {b:016x?}");
+            assert_eq!(asm_ba, portable, "asm mul swapped: {b:016x?} * {a:016x?}");
+        } else {
+            assert_eq!(x86_64_asm::mul(&a, &b, &P, INV), None);
+        }
+    }
+
+    /// Differentially checks one canonical squaring, including agreement
+    /// between the assembly square and the assembly self-multiplication,
+    /// which enter the shared reduction through different stubs.
+    fn check_square(a: Limbs, oracle: bool) {
+        let portable = F::square(&F(a)).0;
+        assert_eq!(
+            <F as Field>::square(&F(a)).0,
+            portable,
+            "dispatched square disagrees with portable: {a:016x?}"
+        );
+        if oracle {
+            assert_montgomery_product(&portable, &a, &a, "portable square");
+        }
+        if x86_64_asm::is_available() {
+            let asm_square = x86_64_asm::square(&a, &P, INV).unwrap();
+            let asm_self_mul = x86_64_asm::mul(&a, &a, &P, INV).unwrap();
+            assert_eq!(asm_square, portable, "asm square: {a:016x?}");
+            assert_eq!(
+                asm_self_mul, asm_square,
+                "asm mul(a, a) != asm square(a): {a:016x?}"
+            );
+        } else {
+            assert_eq!(x86_64_asm::square(&a, &P, INV), None);
+        }
+    }
+
+    /// Checks the extended domain that `from_u512` reaches in production:
+    /// one operand is an arbitrary 256-bit value, the other a canonical
+    /// residue. With one factor below p the Montgomery intermediate stays
+    /// below 2p, so the portable result is exact and the assembly must match
+    /// it bit-for-bit in both argument orders.
+    fn check_noncanonical_pair(raw: Limbs, canonical: Limbs) {
+        assert!(less_than(&canonical, &P));
+        let portable = F::mul(&F(raw), &F(canonical)).0;
+        assert_montgomery_product(&portable, &raw, &canonical, "portable mul, extended domain");
+        assert_eq!(
+            (&F(raw) * &F(canonical)).0,
+            portable,
+            "dispatched mul disagrees with portable: {raw:016x?} * {canonical:016x?}"
+        );
+        if x86_64_asm::is_available() {
+            let asm_ab = x86_64_asm::mul(&raw, &canonical, &P, INV).unwrap();
+            let asm_ba = x86_64_asm::mul(&canonical, &raw, &P, INV).unwrap();
+            assert_eq!(
+                asm_ab, portable,
+                "asm mul, non-canonical lhs: {raw:016x?} * {canonical:016x?}"
+            );
+            assert_eq!(
+                asm_ba, portable,
+                "asm mul, non-canonical rhs: {canonical:016x?} * {raw:016x?}"
+            );
+        }
+    }
+
+    /// All-pairs differential grid over canonicalized edge patterns. The
+    /// bit-serial oracle additionally spot-checks the grid and every square.
+    #[test]
+    fn x86_64_asm_matches_portable_on_edge_patterns() {
+        let mut values: Vec<Limbs> = Vec::new();
+        for pattern in edge_patterns() {
+            let value = canonicalize(pattern);
+            if !values.contains(&value) {
+                values.push(value);
+            }
+        }
+
+        for (i, &a) in values.iter().enumerate() {
+            check_square(a, true);
+            for (j, &b) in values.iter().enumerate().skip(i) {
+                check_mul_pair(a, b, (i + j) % 29 == 0);
+            }
+        }
+    }
+
+    /// `from_u512` (hash-to-field, `Field::random`) multiplies raw 256-bit
+    /// halves by R2 and R3 through the dispatched operator, so the assembly
+    /// receives operands at and above the modulus. Cover that domain
+    /// directly.
+    #[test]
+    fn x86_64_asm_matches_portable_on_noncanonical_inputs() {
+        let multipliers = [
+            R2.0,
+            R3.0,
+            [1, 0, 0, 0],
+            R.0,
+            sub_one(P),
+            canonicalize([u64::MAX; 4]),
+        ];
+        for pattern in edge_patterns() {
+            for &multiplier in &multipliers {
+                check_noncanonical_pair(pattern, multiplier);
+            }
+        }
+    }
+
+    /// Pins the whole `from_u512` conversion against a reconstruction that
+    /// uses only the portable inherent operations.
+    #[test]
+    fn x86_64_asm_from_u512_matches_portable_reconstruction() {
+        fn check(lo: Limbs, hi: Limbs) {
+            let dispatched = F::from_u512([lo[0], lo[1], lo[2], lo[3], hi[0], hi[1], hi[2], hi[3]]);
+            let portable = F::mul(&F(lo), &R2).add(&F::mul(&F(hi), &R3));
+            assert_eq!(
+                dispatched.0, portable.0,
+                "from_u512 diverged: lo={lo:016x?} hi={hi:016x?}"
+            );
+        }
+
+        for &pattern in edge_patterns().iter() {
+            check(pattern, [0; 4]);
+            check([0; 4], pattern);
+            check(pattern, pattern);
+            check(pattern, [u64::MAX; 4]);
+            check([u64::MAX; 4], pattern);
+        }
+    }
+
+    /// Biased fuzzing: words drawn from the edge set, sparse and dense random
+    /// bit soups, and uniform words, reduced when canonicity is required.
+    #[test]
+    fn x86_64_asm_matches_portable_on_structured_random_patterns() {
+        fn draw_word(words: &[u64], rng: &mut XorShiftRng) -> u64 {
+            match rng.next_u32() % 4 {
+                0 => words[rng.next_u32() as usize % words.len()],
+                1 => rng.next_u64() & rng.next_u64() & rng.next_u64(),
+                2 => rng.next_u64() | rng.next_u64() | rng.next_u64(),
+                _ => rng.next_u64(),
+            }
+        }
+
+        fn draw(words: &[u64], rng: &mut XorShiftRng) -> Limbs {
+            [
+                draw_word(words, rng),
+                draw_word(words, rng),
+                draw_word(words, rng),
+                draw_word(words, rng),
+            ]
+        }
+
+        let words = edge_words();
+        let mut rng = XorShiftRng::from_seed([0x5a; 16]);
+        for round in 0..2048 {
+            let a = canonicalize(draw(&words, &mut rng));
+            let b = canonicalize(draw(&words, &mut rng));
+            check_mul_pair(a, b, round % 16 == 0);
+            check_square(a, round % 16 == 0);
+        }
+        for _ in 0..512 {
+            let raw = draw(&words, &mut rng);
+            let canonical = canonicalize(draw(&words, &mut rng));
+            check_noncanonical_pair(raw, canonical);
+        }
+    }
+
+    /// Feeds assembly outputs back into assembly inputs for thousands of
+    /// steps so any divergence compounds, nudging the trajectory regularly
+    /// so it keeps crossing the final conditional-subtraction boundary.
+    #[test]
+    fn x86_64_asm_matches_portable_on_iterated_feedback() {
+        if !x86_64_asm::is_available() {
+            return;
+        }
+
+        let mut x = R2.0;
+        let mut y = sub_one(P);
+        for step in 0..8192u64 {
+            let product = F::mul(&F(x), &F(y)).0;
+            let square = F::square(&F(y)).0;
+            assert_eq!(
+                x86_64_asm::mul(&x, &y, &P, INV).unwrap(),
+                product,
+                "asm mul diverged at step {step}: x={x:016x?} y={y:016x?}"
+            );
+            assert_eq!(
+                x86_64_asm::square(&y, &P, INV).unwrap(),
+                square,
+                "asm square diverged at step {step}: y={y:016x?}"
+            );
+            if step % 64 == 0 {
+                assert_montgomery_product(&product, &x, &y, "walk mul");
+                assert_montgomery_product(&square, &y, &y, "walk square");
+            }
+            x = product;
+            y = square;
+            if step % 17 == 0 {
+                y = F(y).add(&F([step + 1, 0, 0, 0])).0;
+            }
+        }
+    }
+}
+
 #[test]
 fn test_inv() {
     // Compute -(r^{-1} mod 2^64) mod 2^64 by exponentiating
-- 
2.34.1

Authored with the help of Claude Fable.

[ValarDragon]

Need to redo with new data