Two families, one recursion. Forward: {𝔼[Wrpκ(W)]}. Inverse: {𝔼[W−rpκ(W−1)]}. The same Stein/Haff identity generates both, in all three division algebras at once — real (β = 1), complex (β = 2), quaternion (β = 4) — deformed through the single constant cβ = 2/β − 1 ∈ {1, 0, −½}. But it generates them in opposite directions, and that one asymmetry is what separates the two halves of this page: forward, each level is a linear rewrite of the level below, with coefficients polynomial in n; inverse, each level is a linear system to be solved, with coefficients rational in ñ = n − m + 1 − α whose poles are exactly the existence thresholds. Use the switch below to move between them.
Every state Φr,κ = 𝔼[Wrpκ(W)] of the central family, arranged by level N = r + |κ|. Solid arrows are the three terms of the central master recursion — each state is a linear rewrite of states one level below, identically in β (only two weights change). Dashed arrows are the scalar trace reduction, which stays within a level. Press Play to watch the recursion fill in level by level, or click any state to see its instantiated recursion.
Part 1 shows the structure. This section computes the actual central tables — the full family {𝔼[Wr pκ(W)] : r + |κ| ≤ N} — live, in your browser, in exact arithmetic, for all three cases: real W = Z′Z (β = 1) and complex W = Z*Z (β = 2) — the two cases treated by Hillier & Kan (2024) — and quaternion W = ZDZ (β = 4, twin constant c₄ = −½; coefficients are signed dyadic rationals, exact in float64). Every table is a combination of words in Σ alone with polynomial-in-n coefficients; the stored-term counts are identical across the three cases (no exact cancellations occur centrally) — only the coefficients differ. The computation runs in a background Web Worker, is capped at N ≤ 14 (the float64-exact domain), and is gated: the Self-test button recomputes all three cases to N = 6 and checks them, term by term via SHA-256 of a canonical serialization, against digests generated from an independent exact-rational implementation of the same recursion.
Plug a numeric Σ (m×m) and n into the exact table above (float64; the identity is polynomial in Σ and n, so any symmetric — Hermitian in the complex case, quaternion-Hermitian with entries a+bi+cj+dk in the quaternion case — matrix is valid; trace factors are evaluated as Re tr in the quaternion case, as in the stored tables). Enter the upper triangle; the lower mirrors automatically.
Timing note: central tables are small — level 14 stores 157,354 terms across all states, and a typical laptop computes the whole family in well under a second. The largest coefficient at level N occurs in the real case and is exactly 2N−1(N−1)! — 5.1×1013 at N = 14, against the float64-exact ceiling 253 = 9.0×1015. The cap therefore keeps a full level of margin; beyond it, exact bignum engines take over.
The inverse family is Φr,κ = 𝔼[W−rpκ(W−1)], arranged by the same level N = r + |κ|. The Stein/Haff identity applies to it verbatim, but it comes out pointing the other way: what it gives is Σ−1Φr,κ at level N−1 as a combination of states at level N. So the arrows below run from an equation to the unknowns it contains, and the unknowns of a level are coupled: each column is a square linear system — the matrix D−(N−1) — that has to be solved. Everything that distinguishes the inverse side follows from that one fact.
The same three cases, and the same exactness, but a different computation. Because each level is a solve rather than a rewrite, the coefficients are rational functions of ñ = n − m + 1 − α instead of polynomials in n, and every level shares one denominator. The arithmetic below is exact rational arithmetic on BigInts. Each level needs D−k−1, and getting it by elimination is ruinous; instead the page uses the fact that Mk = Dk − nI is diagonalisable with spectrum among the α-contents of the addable cells, so p(x) = ∏(x − c) annihilates it and D−k−1B = (Σj cj(ñ)MkjB) / p(ñ) — products of a scalar matrix against a polynomial matrix, with no polynomial × polynomial multiplication anywhere. (This is what dkimapr.m does in the Matlab package.) The level denominator is then reduced to denpoly's by an exact division, which doubles as a check.
Even so the cap here is N ≤ 11, against 14 on the forward side, and the reason is structural rather than a matter of tuning: the forward array is sparse — level 12 stores 43,954 nonzero terms across every state — while the inverse array is completely dense. Every one of the 195 × 195 entries of C̃12 is a nonzero rational function, 774,704 stored coefficients in all. Inverting a sparse matrix does not give a sparse matrix, and no amount of cleverness in the solve changes that.
Everything on this page depends on Σ only through Σ−1 and its power sums, and any Hermitian Σ is unitarily conjugate to a real diagonal one, so a real symmetric Σ covers every case here — in the complex and quaternion cases the answer for a general Hermitian Σ is the conjugate of this one. Enter the upper triangle; the lower mirrors automatically. Evaluation is float64; the table above is exact.
Where the poles come from: the denominator of level k is ∏(ñ − c) over the α-contents c = αj − i of the cells (i,j) of the partitions of k, each with the largest multiplicity it attains in a single partition. Its largest root is (k−1)α — which is exactly the existence threshold (n−m+1)/α > k for the inverse moments of degree k. The run above prints the denominator and its roots at every level, so you can watch the threshold move.