diff --git a/research/pr4554-roots-recovery-check/.gitattributes b/research/pr4554-roots-recovery-check/.gitattributes
new file mode 100644
index 0000000000..cf808dabba
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/.gitattributes
@@ -0,0 +1 @@
+PR4554-ROOTS-RECOVERY.md -text
diff --git a/research/pr4554-roots-recovery-check/PR4554-ROOTS-RECOVERY.md b/research/pr4554-roots-recovery-check/PR4554-ROOTS-RECOVERY.md
new file mode 100644
index 0000000000..5791526c61
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/PR4554-ROOTS-RECOVERY.md
@@ -0,0 +1,223 @@
+# Native root seeds and exact product recovery
+
+Prepared by Dugongue with substantial OpenAI Codex assistance. Apache-2.0.
+Every algorithm here acts directly on bounded words, arrays and arithmetic.
+There is no Turing-machine primitive, simulator, conversion or multiplication
+black box on an entire variable-length input. Ordinary MUL/DIV on word-sized
+operands are allowed, with explicit sufficient width below.
+
+## 1. A generated primitive-root seed with certified precision
+
+For every integer D>=1 and requested p>=1, root_seed.py constructs Gaussian
+dyadic coordinates x/2^p,y/2^p satisfying
+
+    max(|x/2^p - cos(2pi/D)|, |y/2^p - sin(2pi/D)|) < 2^-p.
+
+Take P=8(p+16), S=2^P, K=128(p+1). It suffices that every word has at least
+2P+20 bits and that D and all array/control addresses fit in one word. Signed
+integers are a sign bit and an unsigned magnitude; signed arithmetic expands
+to bounded comparisons, ordinary ADD/SUB and word MUL/DIV. This is a uniform
+finite loop algorithm, not a stored transcendental constant or root oracle.
+For p=Theta(log N) and D<=N^C with fixed C, the sufficient word width remains
+Theta(log N). Constants are deliberately conservative.
+
+### Construct pi rather than assuming it
+
+Use Machin's identity pi=16 arctan(1/5)-4 arctan(1/239). Its identity can be
+verified without a numerical pi supplier: tan(2 arctan(1/5))=5/12 and
+ tan(4 arctan(1/5))=120/119. The tangent subtraction formula then gives
+ tan(4 arctan(1/5)-arctan(1/239))=1; the angle is in (0,pi/2), so it is pi/4.
+
+For q=5 and 239, sum the alternating terms
+
+    (-1)^j floor(S / ((2j+1) q^(2j+1)))
+
+while q^(2j+1)<=S. Generate successive powers by multiplication by the fixed
+integer q^2. There are at most P+1 terms. Each floored term differs from its
+exact scaled term by less than 1, and the alternating remainder after stopping
+is at most 1/S. Consequently the constructed integer A obeys
+
+    |A/S - pi| <= 20(P+2)/S.
+
+Every denominator is at most S(2P+3); the next q-power is at most 239^2 S.
+These are bounded O(P)-bit integers. Compute the integer angle t=floor(2A/D).
+Then 0<=t<=8S and
+
+    |t/S - 2pi/D| <= [40(P+2)+1]/S.
+
+No search for a prime or primitive root is involved.
+
+### Fixed-point complex Taylor evaluation
+
+Start the Taylor term and sum at (S,0). For k=1,...,K, replace the current
+term (a,b) by
+
+    (trunc(-b*t/(S*k)), trunc(a*t/(S*k))),
+
+and add it to the sum. Truncation is toward zero, implemented on sign/magnitude
+words. Let theta=t/S. Each coordinate's rounding error in term k satisfies
+E_k <= (8/k) E_(k-1)+1/S. Expanding this recurrence and using
+j!/k!<=1/(k-j)! yields E_k<exp(8)/S<4096/S. Thus the summed rounding error in
+each coordinate is at most 4096K/S.
+
+For completeness, exp(8)<4096 follows from the elementary estimate
+exp(1)<11/4 and (11/4)^8<4096. The first estimate follows by summing through
+1/3! and bounding the tail from 1/4! by (1/24)sum_j 5^-j.
+
+The omitted Taylor tail is at most 2*8^(K+1)/(K+1)!, since the subsequent ratio
+is below 1/2. The last K/2+1 factorial factors are each at least 128, so this
+is at most 2^(-K/2-3), which is <=32/S for all p>=1. The cosine and sine
+functions are 1-Lipschitz on the real line. Combining the angle, rounding and
+tail errors gives the explicit coordinate bound
+
+    E_internal <= [40(P+2)+1+4096K+32]/S.
+
+The displayed numerator times 2^(p+1) is strictly less than S. Finally round
+the two fixed-point sums to nearest at p fractional bits. This adds at most
+2^(-p-1), proving the stated strict 2^-p error.
+
+Term and partial-sum magnitudes are below 2^14 S: the exact absolute Taylor
+mass is below exp(8), and the accumulated rounding envelope is smaller than S.
+Products with t are below 2^(2P+17), denominators and all Machin intermediates
+are smaller still. Thus 2P+20 bits is a sufficient native arithmetic bound.
+The implementation charges a conservative
+
+    100(P+1)+40(K+1)
+
+word-operation envelope. Both P and K are O(p), so setup is O(p) word
+operations and O(1) live arithmetic words, besides caller storage. It is not
+necessary to form factorials or unbounded exact rational denominators.
+
+## 2. Rounding and CRT recover exact convolution coefficients
+
+Let two nonnegative base-beta digit arrays have lengths a,b, beta=2^h>=2.
+Every digit is below beta. Their ordinary linear convolution has length
+a+b-1 and coefficients
+
+    c_j = sum_(i+k=j) u_i v_k,
+    0 <= c_j <= C=min(a,b)(beta-1)^2.
+
+**Rounded numerical route.** If a computed dyadic real value Y_j/2^P satisfies
+|Y_j/2^P-c_j|<1/2 for every j, nearest-integer rounding recovers every c_j
+exactly. Range-checking 0<=c_j<=C is useful, but does not itself prove the
+error premise. The numerical transform must establish that premise from its
+actual precision and prefix-norm bounds; it cannot be inferred from good
+sample outputs. The recovery loop costs O(a+b) word operations. Its input
+integers fit when 2^P(C+1), plus fixed rounding slack, fits in a word.
+
+**Exact modular route.** For distinct certified primes p_1,...,p_s, let
+Q=product p_i. If Q>C, the canonical residue in [0,Q) is c_j itself. For
+signed coefficients |c_j|<=C, use Q>2C and the representative in (-Q/2,Q/2).
+
+The supplied CRT constructor checks each prime by literal trial division and
+constructs inverses of preceding products by modular repeated squaring with
+exponent p_i-2. Primality and inverse work is charged setup. Trial division
+costs O(sum sqrt(p_i)) word operations; modular inversion costs
+O(sum log p_i). This is a constructive positive algorithm, but these setup
+costs are not asserted to fit a fast FFT bound. A fast application must provide
+its own fully costed modulus construction or choose the numerical route above.
+
+For each coefficient maintain its reconstructed x in [0,M), initially x=0,
+M=1. The precomputed next modulus data give
+
+    t = ((r_i - (x mod p_i)) mod p_i) * inv(M mod p_i) mod p_i,
+    x <- x + M*t,   M <- M*p_i.
+
+Compute the subtraction as r_i+p_i-(x mod p_i), then reduce; it is nonnegative.
+Induction proves all processed congruences and 0<=x<M. At the end uniqueness
+proves exact recovery. The cost is O(s) native word operations per coefficient,
+with no per-coefficient root/inverse search. It suffices that Q, 2Q, every
+p_i^2 and all actual buffer indices, with constant arithmetic slack, fit in
+a word. If Q is polynomial in input length and s is fixed, this is a genuine
+logarithmic-word O(a+b) recovery stage. No such conclusion is made if Q grows
+beyond the declared word width.
+
+## 3. Linear native carry propagation and exact product
+
+Append c_(a+b-1)=0. Set carry=0. For j=0,...,a+b-1 perform
+
+    t = c_j + carry;
+    digit_j = t AND (beta-1);
+    carry = t >> h.
+
+If carry<=min(a,b)(beta-1), then
+
+    t <= min(a,b) beta(beta-1),
+
+and the next carry obeys the same bound. Thus this quantity, together with
+buffer addresses, is a sufficient word bound. No large input-sized integer
+is ever multiplied by the algorithm.
+
+The induction invariant is
+
+    sum_(j<i) digit_j beta^j + beta^i carry_i
+      = sum_(j<i) c_j beta^j.
+
+At completion the right side is the exact product of the two original digit
+arrays. That product is below beta^(a+b), so after a+b emitted digits the carry
+is zero. The output digits therefore represent the exact integer product.
+Each iteration takes a constant number of native operations, for O(a+b) time
+and O(1) arithmetic scratch. Separate source/result arrays give an immediate
+frame/preservation property; ordinary streaming array stores are sufficient.
+
+If the original inputs had exactly n bits and their limb arrays are padded to
+ceil(n/h) words, the output may contain high zero limbs. The bound product<2^(2n)
+justifies removing them and masking any partial top output limb in O(a+b)
+time; adding zeros to an MSB-first representation would be a different action.
+
+## 4. Direct native multiplication assembly
+
+Take h=Theta(log n), so each input has L=Theta(n/log n) digits. A directly
+implemented approximate DFT/NTT convolution on N=Theta(L) entries, with its
+own paid roots, normalization, bounded-word arithmetic and error<1/2, combines
+with the recovery above in O(L) additional word operations. The numerical
+root seed costs O(log n) when its requested precision is Theta(log n).
+
+Thus a direct O(N log N) word-operation convolution gives a direct O(n)
+Word RAM integer multiplication algorithm. If a fully implemented stronger
+native convolution has cost O(N log^(1-kappa) N), then the same explicit
+assembly gives O(n/log^kappa n) word operations. The latter sentence is a
+sufficient composition theorem, not an assertion that this packet supplies
+that stronger convolution. It does not use a TM black box anywhere.
+
+## 5. Evidence, library reuse and next questions
+
+`root_seed.py` checks 80 seeds, p=1,...,128 and orders through 2^40+87, against
+a separate high-precision Decimal implementation. Its transcendental reference
+is validation-only and is not called by the native constructor. An undersized
+word rejects; the exact order-four result is checked.
+
+`recover.py` checks 960 products, 10,752 CRT coefficients and 10,752 rounded
+coefficients, including extremal digit arrays, unequal lengths, signed CRT,
+composite/duplicate moduli, an invalid carry input and an actual half-unit
+rounding corruption. These are finite tests, not universal proof certificates.
+
+The existing Problem 1.8 sources provide directly relevant precedents:
+`MpLimbs20261007.lv_dig`, `dig_add_high`, `addCode` and `add_run` establish
+literal limb representation and carries in bounded IMP; `Frame` describes
+preserved arrays and IO. `Problem18IMPPrimeInverse20261001.inverseRoutine_run`
+constructs literal modular inversion from a prime/nonzero residue and charges
+84*size(p)+34 bounded instructions. Their source was inspected, not freshly
+kernel-audited here. This packet uses these as precedents, not as claims that
+the newly supplied Python/Word RAM algorithms have been formalized.
+
+The federated DAG search for Machin also found the existing source-only
+`FamousTheorems.machin_formula` public-catalog record and a Gerver-sofa
+certificate declaration; neither hit was re-elaborated or imported. Actual
+inspected project precedents include Dugongue's PR #3877,
+`research/gaussian-expansion-headroom/PROOF.md`, whose interval backend is
+credited to Mason Hale and whose retained lineage includes Jose Manuel's
+physical chirp/Gaussian work. It already uses Machin arctangent remainders and
+outward Taylor enclosures. The registered OpenAI preprint
+*The maximum number of mutually unbiased bases in dimension six*
+(2026-09-24), `build/sections/rounding.tex`, already writes the tangent proof
+and explicit alternating arctangent tail. These are method precedents, not
+freshly checked proofs or new claims of invention. Our new component is the
+explicit bounded-word seed algorithm, global precision envelope and its
+integration contract. The sources and exact inspected hashes are recorded.
+
+Concrete next tasks are to bind the recovery threshold to the actual native
+transform's row-norm/error proof, emit the signed fixed-point/root routines in
+the actual RamStd/IMP language, and audit their chosen word/arena width. A fast
+prime-based variant additionally needs a fast fully paid prime/root supplier.
+None of those requirements is hidden inside the executable root constructor.
diff --git a/research/pr4554-roots-recovery-check/README.md b/research/pr4554-roots-recovery-check/README.md
new file mode 100644
index 0000000000..d9abbba784
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/README.md
@@ -0,0 +1,60 @@
+# Independent check of #4554: root seeds and exact recovery (`ROOTS-RECOVERY.md`)
+
+**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.**
+
+This completes my coverage of the native multiplier in Dugongue's #4554 (head `4f46c6a9`). #4557 checked the FFT error chain, assuming a root seed with coordinate error `< 2^−p`. This package checks that seed construction, and the CRT and carry recovery, against `ROOTS-RECOVERY.md`. The text is vendored as `PR4554-ROOTS-RECOVERY.md` (blob `5791526c`, SHA256-pinned). Nothing is imported from #4554's `root_seed.py` or `recover.py`.
+
+## Results
+
+All checked claims are confirmed.
+
+**1. Root seed.** I implemented the construction from the spec text:
+- `P = 8(p+16)`, `S = 2^P`, `K = 128(p+1)`;
+- Machin π from floored alternating terms while `q^(2j+1) ≤ S` (at most `P+1` terms);
+- `t = ⌊2A/D⌋`;
+- `K` fixed-point complex Taylor steps with truncation toward zero;
+- nearest rounding to `p` bits.
+
+The checks:
+- **Per-term rounding (12 runs, `p ≤ 4`).** Every term is compared with the exact rational Taylor term at `θ = t/S`. Each term's error stays within the recurrence `E_k ≤ (8/k)E_(k−1) + 1/S` and below `e^8/S`; the worst observed error is `6.9/S`. The summed error is at most `4096K/S`, with the worst run at 7.9·10⁻⁵ of that bound.
+- **Machin and angle.** Against a separate AGM π, `|A/S − π| ≤ 20(P+2)/S` and `|t/S − 2π/D| ≤ [40(P+2)+1]/S` hold for `p ≤ 64` and `D` up to `2^40+87`.
+- **Final error `< 2^−p`.**
+  - 112 checks against exact algebraic values for `D ∈ {1, 2, 3, 4, 6, 8, 12}` and `p ≤ 128`. Values involving `√2/2` and `√3/2` are decided by integer square comparisons, so no transcendental code is involved.
+  - 36 checks against a separate decimal reference for other `D` up to `2^40+87`.
+- **Word budget.** Every intermediate stays below `2^(2P+20)`. The observed peak is about 283 bits for `p = 1`, inside the stated `2^(2P+17)` product bound of 289 bits.
+
+**2. Lemmas, exactly.**
+- The Machin tangent identity: `tan 2a = 5/12`, `tan 4a = 120/119` and `tan(4a − b) = 1`.
+- `e < 11/4`, from the partial sum through `1/3!` and the `(1/24)Σ5^−j` tail, and `(11/4)^8 < 4096`.
+- For every `p ≤ 4096`:
+  - the tail ratio is `< 1/2`;
+  - the last `K/2+1` factorial factors are all `≥ 128`;
+  - `2^(−K/2−3) ≤ 32/S`;
+  - `[40(P+2)+1+4096K+32]·2^(p+1) < S`.
+- The factorial tail `2·8^(K+1)/(K+1)! ≤ 2^(−K/2−3)` itself is checked exactly at `K = 256`, 384 and 512.
+
+**3. Recovery.**
+- **CRT.** Garner CRT with Fermat inverses recovers all 23,800 coefficients exactly, both unsigned and signed. It runs on 300 random and extremal digit-array pairs with `β ≤ 2^16`, and the invariant `0 ≤ x < M` holds throughout.
+- **Carries.** `t ≤ min(a,b)β(β−1)` and `carry ≤ min(a,b)(β−1)` hold. The carry invariant holds at every step, the final carry is zero, and every product is exact.
+
+**4. Binding.** 7 statements are bound verbatim to the vendored text.
+
+## Verify
+
+```text
+python3 -B research/pr4554-roots-recovery-check/verify.py
+```
+
+This runs in about 26 s with the standard library only.
+
+## Limitations
+
+- **Not checked:** the operation envelope `100(P+1)+40(K+1)` and the CRT setup costs. Only the structural facts behind them are checked: at most `P+1` Machin terms and `K` Taylor steps.
+- **Finite checks.** The per-term error check uses exact rationals only for `p ≤ 4`; larger `p` are checked on the final error. The universal statement rests on the written recurrence, which these checks corroborate.
+- **Out of scope:** the Gaussian adapters other than the banded inverse (#4556 covers that), the phase and exp tables, and the outer assembly.
+
+## Open next questions
+
+- An independent check of `PHASE-TABLES.md` and `NATIVE-EXP-TABLE.md`, the paid native phase and exp tables used by the Gaussian adapters. First step: replay their stated table precision against exact rational or algebraic references.
+
+Credit: Dugongue (#4554). Machin (1706) and the standard Taylor and CRT methods; the method precedents the PR cites (#3877 and Mason Hale's interval backend). Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance.
diff --git a/research/pr4554-roots-recovery-check/SOURCES.json b/research/pr4554-roots-recovery-check/SOURCES.json
new file mode 100644
index 0000000000..a2a5ea8bbd
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/SOURCES.json
@@ -0,0 +1,18 @@
+{
+  "schema_version": 1,
+  "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1",
+  "files": [
+    {
+      "path": "certificates/selected-result.json",
+      "usage": "provenance",
+      "git_blob": "ed33caef5c5827693e9572d252a66f99cbb73d40",
+      "sha256": "5732707cc1e460b090b53f2463162749ad860de878a18f7408938671d295a9d9"
+    },
+    {
+      "path": "docs/research/agent-task-catalog.json",
+      "usage": "provenance",
+      "git_blob": "64bfc5af6709ed9811171e0a2108f05e1219cb65",
+      "sha256": "352a6184713b8843d96fd2c2071eb7aeaa655906b7c51d45ad444912666c8a96"
+    }
+  ]
+}
diff --git a/research/pr4554-roots-recovery-check/TASK.json b/research/pr4554-roots-recovery-check/TASK.json
new file mode 100644
index 0000000000..03b4efc954
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/TASK.json
@@ -0,0 +1,100 @@
+{
+  "schema_version": 1,
+  "task_id": "independent-moment",
+  "priority_revision": 7,
+  "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1",
+  "original_card": {
+    "id": "independent-moment",
+    "question": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?",
+    "starting_points": [
+      "docs/research/research-family-status.md",
+      "certificates/selected-result.json",
+      "research/staged544-mixed-fourth-entry/README.md",
+      "research/staged-fourth-mixture-ceiling/README.md",
+      "research/r437-two-window-carrier/README.md"
+    ],
+    "methods": [
+      "Independently reconstruct one changed paid ledger or pooled dual with exact source and model binding.",
+      "Check a genuinely uncovered moment implication with a different enclosure method and semantic controls."
+    ],
+    "acceptance_gate": "Demonstrate independence, source binding, positive cases and semantic failure controls. Name the existing checks and the missing coverage. More digits or another equivalent evaluation of an unchanged root is not by itself new verification coverage.",
+    "first_step": {
+      "command": [
+        "python3",
+        "-B",
+        "scripts/research_baseline.py"
+      ],
+      "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY",
+      "inspect": "Choose one uncovered implication in #4088 literal mixed fill, #4091 pooled dual/model binding or #4067 changed frame/host bill. #4091 pins an earlier #4088 head: validate word, capacity and stock equivalence before applying that ceiling to the later fill.",
+      "scope": "Source-binding orientation only; no scientific replay, theorem acceptance or supplier admission. Budget any heavy replay separately."
+    },
+    "first_experiment": "Independently reconstruct one changed ledger or pooled certificate with a positive case and a semantic failure control. Compare exact sources and name what existing checks miss. Keep this separate from supplier-boundary-audit (ordinary/outer dispatch) and uniform-transfer-audit (all-size proof). Stop after the chosen implication; more digits of an unchanged root are not additional mechanism coverage.",
+    "bottleneck_link": "Independent pricing prevents unsupported frontier claims and can reveal genuine accounting slack.",
+    "allowed_variations": [
+      "Reconstruct either a source ledger or a numeric enclosure with a different method."
+    ],
+    "scope_change_triggers": [
+      "Importing the primary function and calling its output independent verification.",
+      "Calling an aggregate histogram or source archive an independently accepted physical supplier."
+    ],
+    "overlap_terms": [
+      "independent moment",
+      "rank histogram",
+      "exact enclosure",
+      "2979",
+      "rank7 joint endframes",
+      "paired extraction",
+      "joint lowering",
+      "4088",
+      "4091",
+      "4067",
+      "staged544-mixed-fourth-entry",
+      "staged-fourth-mixture-ceiling",
+      "r437-two-window-carrier"
+    ],
+    "existing_checks": [
+      {
+        "scope": "PR #4088: mixed local/owner/phase ledger and moment reported.",
+        "evidence": "research/staged544-mixed-fourth-entry/README.md",
+        "remaining": "Independent binding of changed construction to paid histogram, distinct from full admission.",
+        "status": "author-checked",
+        "source_commit": "72d7efc181edb21806a878f174847c96dd410ab0",
+        "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129"
+      },
+      {
+        "scope": "PR #4091: sharp pooled ceiling proposed for fixed THREE/FOUR words and inventory.",
+        "evidence": "research/staged-fourth-mixture-ceiling/README.md",
+        "remaining": "Validate dual, exhaustive pattern scope and exact source equivalence before applying to newer fills.",
+        "status": "author-checked",
+        "source_commit": "a29bdd8f56cd18b805917775cdec09ff935eddc4",
+        "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129"
+      },
+      {
+        "scope": "PR #4067: paid (2,2) to (1,3) exchange and fixed OQ4 host reported.",
+        "evidence": "research/r437-two-window-carrier/README.md",
+        "remaining": "Independent full cost delta; remaining caller/precision/stock obligations are separate.",
+        "status": "author-checked",
+        "source_commit": "a95f11c52e3ca1e4c5171b3c20d8f4b0815aee2e",
+        "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129"
+      }
+    ]
+  },
+  "scope_status": "unchanged",
+  "actual_scope": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?",
+  "scope_change_reason": null,
+  "bottleneck_link": "Completes independent coverage of #4554's native multiplier: the root-seed premise assumed by its FFT error bound, and exact recovery.",
+  "budget_minutes": 45,
+  "stop_rule": "Stop at ROOTS-RECOVERY.md; operation envelopes not re-counted.",
+  "outcome": "completed",
+  "evidence": "Own root-seed implementation: per-term rounding vs exact rational Taylor terms (12 runs), Machin/angle bounds vs AGM pi, final < 2^-p on 112 algebraic and 36 reference checks, word budget; exact lemmas for p<=4096; 23,800 CRT coefficients and carry invariants exact; 7 statements bound.",
+  "author": "Ryanson Jonathan, with Claude Code (Claude Opus 5.5) assistance",
+  "ai_assistance": "Substantial Claude Code (Claude Opus 5.5) assistance; no independent human review.",
+  "open_questions": [
+    {
+      "question": "Independent check of #4554 PHASE-TABLES.md and NATIVE-EXP-TABLE.md.",
+      "why_it_matters": "Paid native tables used by the Gaussian adapters.",
+      "first_step": "Replay their stated table precision against exact references.",
+      "blocked_by": "None."
+    }
+  ]
+}
diff --git a/research/pr4554-roots-recovery-check/VERIFY.json b/research/pr4554-roots-recovery-check/VERIFY.json
new file mode 100644
index 0000000000..8a77ada4cb
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/VERIFY.json
@@ -0,0 +1,18 @@
+{
+  "schema_version": 1,
+  "status": "ready",
+  "python": "3.12",
+  "requirements": [],
+  "inputs": [],
+  "commands": [
+    {
+      "argv": [
+        "{python}",
+        "-B",
+        "{package}/verify.py"
+      ],
+      "timeout_seconds": 300,
+      "checks": "Own root-seed implementation: per-term rounding recurrence vs exact terms, Machin/angle bounds, final < 2^-p vs algebraic values and decimal reference, word budget; exact lemmas p<=4096; Garner CRT and carry recovery; binding to vendored spec (sha256)."
+    }
+  ]
+}
diff --git a/research/pr4554-roots-recovery-check/verify.py b/research/pr4554-roots-recovery-check/verify.py
new file mode 100644
index 0000000000..6f14c384f8
--- /dev/null
+++ b/research/pr4554-roots-recovery-check/verify.py
@@ -0,0 +1,290 @@
+#!/usr/bin/env python3
+"""Independent check of PR #4554 (Dugongue) ROOTS-RECOVERY.md, written from the vendored text only.
+
+Nothing is imported from #4554's root_seed.py or recover.py.
+  1. Root seed. My own implementation of the stated construction (P = 8(p+16), S = 2^P, K = 128(p+1); Machin pi by
+     floored alternating terms while q^(2j+1) <= S; t = floor(2A/D); K-term fixed-point complex Taylor with truncation
+     toward zero; nearest rounding to p bits). For small p every term is compared with the exact rational Taylor term
+     at theta = t/S: per-term error <= the recurrence E_k <= (8/k)E_(k-1) + 1/S < e^8/S, summed error <= 4096K/S;
+     |A/S - pi| <= 20(P+2)/S and |t/S - 2pi/D| <= [40(P+2)+1]/S against a separate AGM pi; the final coordinate error
+     < 2^-p against exact algebraic values (D = 1, 2, 3, 4, 6, 8, 12, by integer square comparisons) and against a
+     separate high-precision decimal reference for other D up to 2^40+87; observed magnitudes stay below 2^(2P+20).
+  2. Lemmas, exactly: Machin's tangent identity (tan 2a = 5/12, tan 4a = 120/119, tan(4a - b) = 1); e < 11/4 from the
+     stated partial sum and geometric tail; (11/4)^8 < 4096; the Taylor tail 2*8^(K+1)/(K+1)! <= 2^(-K/2-3) <= 32/S
+     and [40(P+2)+1+4096K+32]*2^(p+1) < S for every p <= 4096.
+  3. Recovery: Garner CRT with Fermat inverses (unsigned and signed representatives) and the linear carry invariant,
+     including the bound carry <= min(a,b)(beta-1), on random and extremal digit arrays.
+Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0.
+"""
+import decimal
+from fractions import Fraction as Fr
+import hashlib
+import json
+import math
+from pathlib import Path
+import random
+import sys
+
+HERE = Path(__file__).resolve().parent
+SPEC = HERE / "PR4554-ROOTS-RECOVERY.md"  # #4554 head 4f46c6a9, research/native-word-fft/ROOTS-RECOVERY.md, blob 5791526c
+SPEC_SHA256 = "5e115c3ded1734dbcb119c889755a13ce9e9357598ec20a3eb1d3cdd727638f4"
+
+
+def need(condition, message):
+    if not condition:
+        raise AssertionError(message)
+
+
+def trunc_div(a, b):
+    """Integer division truncating toward zero (b > 0)."""
+    return a // b if a >= 0 else -((-a) // b)
+
+
+class Track:
+    def __init__(self, P):
+        self.limit = 1 << (2 * P + 20)
+        self.peak = 0
+
+    def __call__(self, v):
+        self.peak = max(self.peak, abs(v))
+        need(abs(v) < self.limit, "magnitude below 2^(2P+20)")
+        return v
+
+
+def machin(P, T):
+    S = 1 << P
+
+    def arctan_inv(q):
+        total, j, power = 0, 0, q
+        while power <= S:
+            term = T(S // T((2 * j + 1) * power))
+            total += -term if j % 2 else term
+            power = T(power * q * q)
+            j += 1
+        need(j <= P + 1, "at most P+1 terms")
+        return total
+
+    return 16 * arctan_inv(5) - 4 * arctan_inv(239)
+
+
+def seed(D, p, exact_terms=False):
+    P = 8 * (p + 16)
+    S = 1 << P
+    K = 128 * (p + 1)
+    T = Track(P)
+    A = machin(P, T)
+    t = T(2 * A // D)
+    need(0 <= t <= 8 * S, "0 <= t <= 8S")
+    a, b = S, 0
+    sx, sy = S, 0
+    errs = []
+    if exact_terms:
+        theta = Fr(t, S)
+        ex = (Fr(1), Fr(0))
+    for k in range(1, K + 1):
+        a, b = T(trunc_div(T(-b * t), S * k)), T(trunc_div(T(a * t), S * k))
+        sx, sy = T(sx + a), T(sy + b)
+        if exact_terms:
+            ex = (-ex[1] * theta / k, ex[0] * theta / k)
+            errs.append(max(abs(Fr(a, S) - ex[0]), abs(Fr(b, S) - ex[1])))
+    shift = P - p
+    x = (sx + (1 << (shift - 1))) >> shift  # nearest at p bits
+    y = (sy + (1 << (shift - 1))) >> shift
+    return x, y, A, t, S, K, errs, T.peak
+
+
+def agm_pi(prec):
+    decimal.getcontext().prec = prec + 10
+    D = decimal.Decimal
+    a, b, t, pw = D(1), D(1) / D(2).sqrt(), D(1) / D(4), D(1)
+    for _ in range(int(math.log2(prec)) + 4):
+        an = (a + b) / 2
+        b = (a * b).sqrt()
+        t -= pw * (a - an) ** 2
+        a = an
+        pw *= 2
+    return (a + b) ** 2 / (4 * t)
+
+
+def ref_cos_sin(D, digits):
+    decimal.getcontext().prec = digits + 20
+    pi = agm_pi(digits + 20)
+    x = 2 * pi / D
+    while x > pi:
+        x -= 2 * pi
+    c, s, term, k = decimal.Decimal(0), decimal.Decimal(0), decimal.Decimal(1), 0
+    eps = decimal.Decimal(10) ** (-(digits + 15))
+    while abs(term) > eps or k < 4:
+        if k % 4 == 0:
+            c += term
+        elif k % 4 == 1:
+            s += term
+        elif k % 4 == 2:
+            c -= term
+        else:
+            s -= term
+        k += 1
+        term = term * x / k
+    return c, s, pi
+
+
+def within_sqrt(x, p, num, den):
+    """|x/2^p - sqrt(num/den)| < 2^-p for x >= 1, via integer comparisons: (x-1)^2 den < num 4^p < (x+1)^2 den."""
+    return (x - 1) ** 2 * den < num * 4 ** p < (x + 1) ** 2 * den
+
+
+def check_seeds():
+    out = {"exact_term_runs": 0, "algebraic_checks": 0, "reference_checks": 0, "worst_term_err_times_S": 0.0,
+           "worst_sum_err_over_4096K": 0.0, "peak_bits_over_2P20": 0.0}
+    for p in (1, 2, 4):
+        for D in (1, 3, 7, 16):
+            x, y, A, t, S, K, errs, peak = seed(D, p, exact_terms=True)
+            P = 8 * (p + 16)
+            E = Fr(0)
+            for k, e in enumerate(errs, start=1):
+                E = Fr(8, k) * E + Fr(1, S)
+                need(e <= E and e < Fr(4096, S), "per-term rounding recurrence")
+            out["worst_term_err_times_S"] = max(out["worst_term_err_times_S"], float(max(errs) * S))
+            need(sum(errs) <= Fr(4096 * K, S), "summed rounding error")
+            out["worst_sum_err_over_4096K"] = max(out["worst_sum_err_over_4096K"], float(sum(errs) * S / (4096 * K)))
+            out["peak_bits_over_2P20"] = max(out["peak_bits_over_2P20"], peak.bit_length() / (2 * P + 20))
+            out["exact_term_runs"] += 1
+    # Machin and angle bounds against a separate AGM pi.
+    for p in (1, 8, 32, 64):
+        P = 8 * (p + 16)
+        S = 1 << P
+        pi = agm_pi(P // 3 + 60)
+        A = machin(P, Track(P))
+        need(abs(decimal.Decimal(A) / decimal.Decimal(S) - pi) <= decimal.Decimal(20 * (P + 2)) / S, "Machin bound")
+        for D in (1, 3, 1000, (1 << 40) + 87):
+            t = 2 * A // D
+            need(abs(decimal.Decimal(t) / S - 2 * pi / D) <= decimal.Decimal(40 * (P + 2) + 1) / S, "angle bound")
+    # Final error: exact algebraic values.
+    for p in (1, 2, 3, 8, 16, 33, 64, 128):
+        P2 = 1 << p
+        for D, (cx, cy) in {1: ((1, 1), (0, 1)), 2: ((-1, 1), (0, 1)), 4: ((0, 1), (1, 1)), 8: ("r2", "r2"),
+                            3: ((-1, 2), "r3"), 6: ((1, 2), "r3"), 12: ("r3", (1, 2))}.items():
+            x, y = seed(D, p)[:2]
+            for v, c in ((x, cx), (y, cy)):
+                if c == "r2":
+                    need(within_sqrt(v, p, 1, 2), "sqrt(2)/2 within 2^-p")
+                elif c == "r3":
+                    need(within_sqrt(v, p, 3, 4), "sqrt(3)/2 within 2^-p")
+                else:
+                    need(abs(Fr(v, P2) - Fr(*c)) < Fr(1, P2), "rational value within 2^-p")
+                out["algebraic_checks"] += 1
+    # Final error: decimal reference for general orders.
+    for p in (1, 5, 17, 40, 64, 128):
+        for D in (5, 7, 10, 97, 1 << 20, (1 << 40) + 87):
+            x, y = seed(D, p)[:2]
+            c, s, _ = ref_cos_sin(D, int(p * 0.302) + 30)
+            P2 = decimal.Decimal(1 << p)
+            need(abs(decimal.Decimal(x) / P2 - c) < 1 / P2 and abs(decimal.Decimal(y) / P2 - s) < 1 / P2, "final < 2^-p")
+            out["reference_checks"] += 1
+    return out
+
+
+def check_lemmas():
+    def tan_add(u, v):
+        return (u + v) / (1 - u * v)
+
+    a = Fr(1, 5)
+    t2 = tan_add(a, a)
+    t4 = tan_add(t2, t2)
+    need(t2 == Fr(5, 12) and t4 == Fr(120, 119), "tan 2a, tan 4a")
+    need((t4 - Fr(1, 239)) / (1 + t4 * Fr(1, 239)) == 1, "tan(4a - b) = 1")
+    partial = sum(Fr(1, math.factorial(k)) for k in range(4))  # through 1/3!
+    tail = Fr(1, 24) * Fr(1, 1 - Fr(1, 5))  # (1/24) sum_j 5^-j bounds sum_{k>=4} 1/k!
+    need(sum(Fr(1, math.factorial(k)) for k in range(4, 40)) < tail, "tail bound for e")
+    need(partial + tail < Fr(11, 4) and Fr(11, 4) ** 8 < 4096, "e < 11/4, (11/4)^8 < 4096")
+    checked = 0
+    for p in range(1, 4097):
+        P, K = 8 * (p + 16), 128 * (p + 1)
+        need(8 * 2 < K + 2, "tail ratio 8/(K+2) < 1/2")
+        need(K // 2 + 1 <= K + 1 - 127, "last K/2+1 factors are at least 128")
+        need(K // 2 + 3 >= P - 5, "2^(-K/2-3) <= 32/S")
+        need((40 * (P + 2) + 1 + 4096 * K + 32) << (p + 1) < 1 << P, "numerator * 2^(p+1) < S")
+        checked += 1
+    for K in (256, 384, 512):  # the factorial tail itself, exactly
+        need(Fr(2 * 8 ** (K + 1), math.factorial(K + 1)) <= Fr(1, 2 ** (K // 2 + 3)), "factorial tail")
+    return checked
+
+
+def check_recovery(rng):
+    primes = [1000003, 1000033, 1000037, 1000039]
+    for q in primes:
+        need(all(q % d for d in range(2, int(q ** 0.5) + 1)), "trial-division primality")
+    coeffs = 0
+    for trial in range(300):
+        h = rng.choice((1, 2, 4, 8, 16))
+        beta = 1 << h
+        a, b = rng.randint(1, 40), rng.randint(1, 40)
+        extremal = trial % 5 == 0
+        u = [beta - 1 if extremal else rng.randrange(beta) for _ in range(a)]
+        v = [beta - 1 if extremal else rng.randrange(beta) for _ in range(b)]
+        c = [0] * (a + b - 1)
+        for i, x in enumerate(u):
+            for k, y in enumerate(v):
+                c[i + k] += x * y
+        C = min(a, b) * (beta - 1) ** 2
+        need(all(0 <= x <= C for x in c), "coefficient range")
+        mods = []
+        Q = 1
+        for q in primes:
+            mods.append(q)
+            Q *= q
+            if Q > 2 * C:
+                break
+        need(Q > 2 * C, "enough moduli")
+        for j, cj in enumerate(c):
+            for signed in (False, True):
+                val = -cj if signed and j % 2 else cj
+                x, M = 0, 1
+                for q in mods:
+                    r = val % q
+                    inv = pow(M % q, q - 2, q)
+                    tt = ((r + q - (x % q)) % q) * inv % q
+                    x, M = x + M * tt, M * q
+                    need(0 <= x < M, "Garner invariant")
+                rec = x - M if signed and x > M // 2 else x
+                need(rec == val, "CRT recovery")
+                coeffs += 1
+        # Carry propagation.
+        cc = c + [0]
+        carry, digits = 0, []
+        for j in range(len(cc)):
+            tt = cc[j] + carry
+            need(tt <= min(a, b) * beta * (beta - 1), "t bound")
+            digits.append(tt & (beta - 1))
+            carry = tt >> h
+            need(carry <= min(a, b) * (beta - 1), "carry bound")
+            need(sum(d << (h * i) for i, d in enumerate(digits)) + (carry << (h * (j + 1)))
+                 == sum(x << (h * i) for i, x in enumerate(cc[:j + 1])), "carry invariant")
+        need(carry == 0, "final carry")
+        U = sum(x << (h * i) for i, x in enumerate(u))
+        V = sum(x << (h * i) for i, x in enumerate(v))
+        need(sum(d << (h * i) for i, d in enumerate(digits)) == U * V, "exact product")
+    return coeffs
+
+
+def check_binding():
+    raw = SPEC.read_bytes()
+    need(hashlib.sha256(raw).hexdigest() == SPEC_SHA256, "vendored ROOTS-RECOVERY.md bytes")
+    text = raw.decode("utf-8")
+    for s in ("Take P=8(p+16), S=2^P, K=128(p+1).", "|A/S - pi| <= 20(P+2)/S.", "|t/S - 2pi/D| <= [40(P+2)+1]/S.",
+              "E_k <= (8/k) E_(k-1)+1/S", "E_internal <= [40(P+2)+1+4096K+32]/S.", "x <- x + M*t,   M <- M*p_i.",
+              "t <= min(a,b) beta(beta-1),"):
+        need(s in text, f"spec states: {s}")
+    return 7
+
+
+def main():
+    rng = random.Random(45542)
+    report = {"status": "PASS_PR4554_ROOTS_RECOVERY_CHECK", "seeds": check_seeds(), "lemma_p_values": check_lemmas(),
+              "crt_coefficients_recovered": check_recovery(rng), "spec_statements_bound": check_binding(),
+              "new_kappa_claimed": False}
+    print(json.dumps(report, indent=1))
+
+
+if __name__ == "__main__":
+    sys.exit(main())