Add trace_moment: faster tr(H^k) via identity-multiset enumeration (#80)

[thedaemon-wizard]

Add trace_moment: faster tr(H^k) via identity-multiset enumeration

Closes #80.

Summary

Adds trace_moment, which computes moments tr(o^k) (and tr(A^k B^l)) by enumerating only
the Pauli-string multisets whose product is the identity and summing each ordering phase
analytically — without ever constructing o^(k/2), exactly as requested in #80.

This is a complete implementation covering:

• trace_moment(o::Operator, k; scale=0) — the issue's example (plain Operator).
• trace_moment(A::Operator, k, B::Operator, l; scale=0) — the tr(A^k B^l) generalisation.
• trace_moment(o::Operator{<:PauliStringTS}, k; scale=0) — translation-symmetric operators,
  exploiting translation invariance by anchoring the first factor to a stored representative.
• Multithreading for all of the above (independent search seeds, deterministic reduction).

trace_product is left unchanged; trace_moment is added as a new, non-breaking function.

How it works

From tr(H^k) = Σ c_{l₁}…c_{lₖ} tr(P_{l₁}…P_{lₖ}), a term contributes only when the XOR of the
chosen Pauli strings is the identity — a property of the multiset, not the order.

1. Enumerate identity multisets. A pruned DFS over the terms tracks the running XOR and keeps
   only branches that can still reach the target (support-reachability via a precomputed suffix
   union, plus a ceil(active_sites / max_support) lower bound; terms sorted by lowest active
   site). This removes the k! factor of the naive ordered sum.
2. Analytic per-multiset phase. The ordering sum factors as Σ_orderings σ = K_ref(M)·D̃(M),
   where K_ref is the sign of one canonical ordering (O(k)) and D̃ depends only on the
   anticommutation graph: terms commuting with all others peel off as binomial coefficients,
   leaving a small "anticommuting core" handled by a mixed-radix multiplicity DP.
3. tr(A^k B^l). The cheaper side is tabulated by XOR target R (summing its multiset
   contributions); the other side is enumerated for each R, and the two blocks combine through
   a (-1)^{ycount(R)} cross phase: tr(A^k B^l)=scale·Σ_R fA(R)·fB(R)·(-1)^{ycount(R)}.
4. Translation symmetry. For OperatorTS, anchoring the first factor to a representative and
   multiplying by the number of translations folds the search by ~1/|G| (the stored-vs-physical
   coefficient relation makes the stabilizer factors cancel exactly).
5. Multithreading. The search is split into independent seeds (each in its own workspace) and
   the partial sums reduced in a fixed order, so the result is independent of the thread count.

The i^ycount Clifford phases are carried inside the stored coefficients, so the bookkeeping is
real ±1, exactly as in trace_product.

Performance

Issue example — plain Operator, TFIM, N=20, k=1:14, scale=1:

	trace_product	trace_moment (1 thr)	trace_moment (8 thr)
total k=1:14	76.2 s	40.5 s (1.9×)	20.4 s (3.7×)
k=11	3.26 s	0.065 s (50×)	0.029 s (111×)
k=13	36.8 s	0.39 s (95×)	0.23 s (158×)

tr(H^k · O), H=TFIM(N=20), O=X₁, k=1:14 (the regime the moment method wins most):

	trace_product	trace_moment (8 thr)
total k=1:14	119.8 s	5.6 s (21×)
k=14	71.0 s	0.39 s (182×)

Odd / high moments are dramatically faster (the existing path builds H^{(k+1)/2}; the multiset
method gets odd moments almost for free and never materialises a high power). Even moments are
matched/beaten once multithreading is enabled. The translation-symmetric method gives the same
wins for OperatorTS inputs.

Results agree with trace_product and with the dense operator power (trace(o^k),
trace(A^k*B^l)) to machine precision.

What's included

• src/moments.jl: trace_moment (plain Operator, 4-arg, and OperatorTS) + self-contained
  helpers (pruned DFS, peel + mixed-radix-DP phase, XOR-table for the 4-arg, threaded seeds).
• src/PauliStrings.jl: export trace_moment.
• docs/src/docstrings.md: added under "Power and moments".
• test/moments.jl and test/moments_ts.jl: test sets wired into runtests.jl.
• benchmark/benchmarks.jl: trace_moment vs trace_product entries (plain, 4-arg, 1D/2D TS).

Test plan

• Full suite passes on Julia 1.10 and 1.6 (the CI matrix).
• trace_moment vs trace_product and vs the dense operator power (trace(o^k),
  trace(A^k*B^l)), to machine precision — plain & TS, 1D/2D/mixed periodicity, k = 0…8,
  complex / non-Hermitian coefficients, the scale keyword, k=0/l=0, error handling.
• serial == multithreaded for every method.

Relation to #91

PR #91 also implements the moment method for Operator. This PR was developed independently and
additionally provides the translation-symmetric OperatorTS method and multithreading.
Happy to consolidate with #91 in whatever way is most useful to the maintainers.

Limitations

• Integer phase bookkeeping is exact up to k ≈ 20, beyond which the moment value itself exceeds
  the Float64 range.

Theory & references that informed the implementation

• The moment method itself — the PauliStrings.jl paper, N. Loizeau et al.,
  Quantum many-body simulations with PauliStrings.jl (arXiv:2410.09654),
  states Tr H^k = Σ h_{i₁}…h_{iₖ} Tr(τ_{i₁}…τ_{iₖ}) and that this is tractable because the operators
  are sparse in the Pauli basis and the expansion contains large commuting sets — the basis for the
  whole trace_moment approach.
• Binary symplectic Pauli algebra — Dehaene & De Moor, Clifford group, stabilizer states, and
  linear and quadratic operations over GF(2), Phys. Rev. A 68, 042318 (2003) — the (v,w)
  representation (already used by this package) and its quadratic-form commutation phase. The whole
  phase computation reduces to XOR (p₁ ⊻ p₂) and the parity count_ones(p₁.v & p₂.w) (mod 2),
  so all per-multiset phases are exact integer ±1.
• Signed ordering sums = q-binomials at q = −1 — the anticommuting-core ordering sum D̃(M)
  is a (-1)-weighted (Mahonian / Gaussian q-binomial at q=-1) enumeration; terms commuting
  with everything factor out as ordinary binomials, the rest are summed by a small mixed-radix DP.
  (Standard enumerative combinatorics — e.g. Stanley, Enumerative Combinatorics Vol. 1, §1.7,
  Gaussian binomial coefficients.)
• Why high moments matter (motivation) — Parker, Cao, Avdoshkin, Scaffidi & Altman,
  A Universal Operator Growth Hypothesis, Phys. Rev. X 9, 041017 (2019)
  (arXiv:1812.08657); and the numerical moment/Lanczos study
  arXiv:2203.00533: Lanczos coefficients (hence Krylov/operator
  complexity) are obtained directly from the Hamiltonian moments.
• Moments → Lanczos → Green's functions (application) — Greene-Diniz et al., Quantum Computed
  Green's Functions using a Cumulant Expansion of the Lanczos Method, Quantum 8, 1383 (2024)
  (arXiv:2309.09685).
• High-performance Pauli-string processing (design corroboration) — Krötz, PauLIB: A
  High-Performance Library for Processing Pauli Strings (2026,
  arXiv:2605.25974), which independently argues for a bit-packed
  symplectic representation, a sorted-array layout instead of hash maps, and multi-threaded
  merging — matching the choices here (alphabet sorted by lowest active site, bitwise XOR algebra,
  thread-local search seeds with a fixed-order reduction).

AI assistance disclosure

Per the unitaryHACK Ethical AI policy: I used Claude (Anthropic) as a coding co-pilot — for
drafting, the algebra of the phase factorisation / translation folding, and test/benchmark
scaffolding. I reviewed and understand every part and can explain it; correctness was verified by
me locally against two independent oracles (trace_product and the dense operator power) with the
full suite passing on Julia 1.6 and 1.10. No unverified output was submitted.