Implement OrthogonalPolynomialMomentFunctionalSequence and its fractional-Riccati subclass

[crowlogic]

Issue: Implement OrthogonalPolynomialMomentFunctionalSequence and its Fractional-Riccati Subclass

What this object is (the math, not the application)

Given a linear functional 𝓛 : R[x] → R on polynomials (where R is any
coefficient ring the existing infrastructure supports — Real,
RealPolynomial, Complex, ComplexPolynomial, etc.), presented as its
moment sequence

m(k) := 𝓛[x^k],     k = 0, 1, 2, …

the bilinear form ⟨p, q⟩ := 𝓛[p·q] endows R[x] with a notion of pairing.
When the form is nondegenerate on every R[x]_{≤n} — equivalently, when every
leading Hankel minor det[m(i+j)]_{i,j≤n} is a unit in R — Chihara calls
the functional quasi-definite (or regular), and R[x] admits a
unique distinguished basis: the monic polynomials P(n, x) obtained by
Gram–Schmidt of 1, x, x², … against ⟨·,·⟩.

That basis is the orthogonal polynomial sequence of 𝓛. It is the basis
of R[x] in which the moment functional's bilinear form is diagonal. From
it the following are reconstructable:

• The Jacobi matrix J(n) — the tridiagonal matrix of multiplication-
  by-x in the OPS basis. The three-term recurrence (α(n), β(n)) are
  its entries.
• The Stieltjes transform C(z) = Σ_k m(k) / z^{k+1} — its
  J-fraction expansion has (α(n), β(n)) as partial quotients.
• The representing measure dμ (positive-definite case only) — the
  spectral measure of J at e_0; supp(dμ) = spec(J). Recovered by
  Stieltjes inversion.
• The diagonal Padé approximants of the moment generating series —
  denominators are the OPS polynomials reciprocal-flipped.

Structural facts. Once the OPS exists:

• ⟨P(n,·), P(m,·)⟩ = 0 for n ≠ m — the basis diagonalizes the pairing.
• ⟨P(n,·), P(n,·)⟩ = h(n) ∈ R records the functional's scale at degree n.
• Multiplication by x is self-adjoint with respect to ⟨·,·⟩. Its matrix
  in the P(n, ·) basis is therefore tridiagonal — the Jacobi matrix
  J. The three-term recurrence
  x · P(n, x) = P(n+1, x) + α(n) · P(n, x) + β(n) · P(n-1, x)
  is just J read off entry-wise.
• Favard's theorem is the converse: any sequence of monic polynomials
  satisfying a three-term recurrence with β(n) units arises from a
  unique quasi-definite moment functional.

So moment functionals and three-term recurrences carry the same content,
and this class is the constructor of one from the other.

What to call this class

The class is the OPS of a moment functional. The name is

OrthogonalPolynomialMomentFunctionalSequence

Concrete applications subclass it by supplying a momentSequence().

Existing arb4j infrastructure to build on

This issue does not introduce the Müntz–Padé / fractional-Riccati
machinery from scratch — most of the parameter-handling and recurrence
plumbing already lives in arb.functions.complex:

• MuntzPadeFunctional — holds α ∈ (0,1) and a
  ComplexPolynomialSequence a (the Müntz coefficients k ↦ v ↦ a(k, v)),
  returns a MuntzPadeApproximant on evaluation.
• RiccatiMuntzPadeFunctional extends MuntzPadeFunctional — takes
  ComplexPolynomialNullaryFunction P, Q, R (caller-supplied parameter
  factories), registers ComplexPolynomial p, q, r as variables in the
  Context, compiles the Müntz coefficient recurrence
  a:k➔v➔when(k=1, p(v)/Γ(μ+1), else, (Γ((k-1)μ+1)/Γ(kμ+1)) · (q(v)·a(k-1)(v) + r(v)·∑ a(j)(v)·a(k-1-j)(v){j=1..k-2}))
  as a ComplexPolynomialSequence, exposes the discriminant
  q(v)² − 4·p(v)·r(v), supports invalidateCache() on parameter changes,
  and computes the Jacobian ∂y/∂v symbolically via
  ComplexPolynomial.derivative().

Relevant typed sequence/functional API already in place:

• Functional<DOMAIN, CODOMAIN, FUNC> and its specializations.
• Sequence<R> and the full family: RealSequence, ComplexSequence,
  RealPolynomialSequence, ComplexPolynomialSequence,
  RealSequenceSequence, ComplexSequenceSequence, etc.
• RecurrentlyGeneratedOrthogonalPolynomialSequence<R, V, E> and its
  Real… / Complex… specializations.
• All Sequence instances are auto-cached.

Codomain note: complex throughout. The characteristic-function /
Heston application is intrinsically complex (u ∈ ℝ, but P(u), Q(u), R(u) ∈ ℂ[u]). This issue's classes are therefore in the
arb.functions.polynomials.orthogonal.complex package, with codomain
ComplexPolynomial for the OPS itself, and a recurrence-coefficient
ring R = ComplexPolynomial (rather than Complex) — α(n) and
β(n) are polynomials in the Fourier parameter u, not scalars.

Relationship to this issue: the new
OrthogonalPolynomialMomentFunctionalSequence is a consumer of the
Müntz coefficient sequence produced by MuntzPadeFunctional.muntzBasis(),
not a reimplementation. Specifically,
OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence should
accept either (μ, P, Q, R) directly (constructing its own
RiccatiMuntzPadeFunctional internally) or an existing
RiccatiMuntzPadeFunctional instance whose muntzBasis() is fed in as
the moment sequence. The latter form lets a caller share one Riccati
functional between its Müntz–Padé evaluation path and its OPS / Jacobi-
matrix / Padé-denominator path without recomputing the Müntz coefficients.

The polynomial Chebyshev/Wheeler S(n, z) machinery, the (A, B, C)
extraction, and the reciprocal-polynomial flip are the genuinely new
pieces; everything upstream of momentSequence() is already written.

Class hierarchy

RecurrentlyGeneratedOrthogonalPolynomialSequence<P, V, R>                          (existing)
  └── ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence (NEW, R = ComplexPolynomial)
        └── OrthogonalPolynomialMomentFunctionalSequence                            (NEW, abstract)
              └── OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence     (NEW, concrete)

The new parent class
ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence
exists because the existing Complex…RecurrentlyGenerated… peer uses
Complex (scalars) as the recurrence-coefficient ring, while the
characteristic-function setting forces the recurrence coefficients to be
elements of ℂ[u].

Every classical OPS already in the codebase (LaguerrePolynomials,
HermitePolynomials, LegendrePolynomials, Type1ChebyshevPolynomials,
JacobiPolynomialSequence, etc.) is implicitly a moment-functional OPS
that happens to skip the abstract base because its (α(n), β(n)) are
tabulated as elementary functions of n and the weight parameters. The
new abstract class is the genus; those are explicit-form specimens; the
Riccati subclass is the recursive-form specimen.

Construction — polynomial Chebyshev/Wheeler recursion (revised)

This is the core new contribution and replaces the bivariate-scalar
S(n, k) formulation in the prior draft of this issue. Full statement
and proof: docs/SolutionMethodologyForConstantCoefficientFractionalRiccatiEquations.md,
Theorems 4.1–5.1 and Remark 5.2.

For each n ≥ −1, define the generating polynomial of the OPS

S(n, z) := Σ_{k ≥ 0} 𝓛[x^k · P(n, x)] · z^k

— a univariate polynomial in z whose coefficients are in ℂ[u]. By
orthogonality, 𝓛[x^k · P(n, x)] = 0 for k < n, so S(n, z) is
divisible by z^n. Base cases:

S(−1, z) = 1,    S(0, z) = Σ_{k ≥ 0} m(k, u) · z^k.

Theorem 5.1 (Polynomial Chebyshev–Wheeler). For n ≥ 1, write
S(n, z) = z^n · S̃(n, z) with S̃(n, 0) = h(n) ∈ ℂ[u]ˣ (a unit by
quasi-definiteness). The sequence {S(n, z)} satisfies the closed
polynomial recurrence

┌──────────────────────────────────────────────────────┐
│                                                      │
│   S(n+1, z) = S(n, z) / z                            │
│                − α(n) · S(n, z)                      │
│                − β(n) · S(n−1, z)                    │
│                                                      │
└──────────────────────────────────────────────────────┘

in which α(n), β(n) ∈ ℂ[u] are the unique elements making the right-
hand side divisible by z^{n+1}. Equivalently, matching the two relevant
orders of z in this divisibility constraint yields, as identities in
ℂ[u][[z]],

β(n) = S(n, z) / ( z · S(n−1, z) ),
α(n) = − S(n+1, z) / S(n, z).

In implementation, neither identity is evaluated as a power-series
division — β(n) is determined by matching the z^{n−1} coefficient
of the recurrence, and α(n) by matching the z^n coefficient; both
produce elements of ℂ[u] directly.

The three-term recurrence coefficients in the standard form
P(n+1, x) = (A(n) · x + B(n)) · P(n, x) − C(n) · P(n−1, x) are then

A(n) = 1,
B(n) = −α(n) = S(n+1, z) / S(n, z),
C(n) =  β(n) = S(n, z) / ( z · S(n−1, z) ).

Seeds: P(0, x) = 1, P(1, x) = x − m(0, u) / 1 (degenerate n = 0
case handled by the divisibility constraint cancelling the single pole
of S(0, z) / z).

Why this form, not the bivariate-scalar form. The prior draft of
this issue stated the construction as a bivariate-scalar array
S(n, k) := 𝓛[x^k · P(n, x)], with α(n), β(n) recovered as ratios of
specific S(n, k) entries. That formulation requires scalar coefficient
extractions from ℂ[u] at every step, which is forbidden — α(n) and
β(n) must remain opaque elements of ℂ[u] throughout, accessible only
via the ring operations of ℂ[u]. The polynomial form above achieves
this: every step is one polynomial division, one polynomial linear
combination, and exactly two z-coefficient reads. No scalar
coefficient extractions in the parameter ring ℂ[u] are performed —
the two coefficient reads are at the outer z level, which is the
truncation level of the moment series, not the inner u level.

Cost: Θ(M²) polynomial ops in ℂ[u] (one acb_poly_div_series and
one acb_poly_add / _sub / scalar-mul per step; one
acb_poly_get_coeff for each of α(n), β(n)).

The motivating application: fractional Riccati

For μ ∈ (0,1], P(u), Q(u), R(u) ∈ ℂ[u] — functions of the Fourier
parameter u, constant in t — solve

D^μ y(t,u) = P(u) + Q(u)·y(t,u) + R(u)·y(t,u)²,
I^{1−μ} y(0,u) = 0.

The Müntz–Tau ansatz y(t,u) = Σ_{k≥1} a(k, u) · t^{kμ} gives the
recurrence

a(1, u) = P(u) / Γ(μ + 1)
a(k+1, u) = [Γ(kμ + 1) / Γ((k+1)μ + 1)]
            · ( Q(u) · a(k, u)
                + R(u) · Σ_{j+ℓ=k, j,ℓ≥1} a(j, u) · a(ℓ, u) )

The substitution z = t^μ gives g(z, u) = Σ_{k≥1} a(k, u) · z^k,
meromorphic on ℂ for each u (Theorem 9.18 of §9 of the companion PDF).
Setting m(k, u) := a(k+1, u) makes this the moment sequence of a
quasi-definite (but signed, never positive-definite) functional over
the coefficient ring ℂ[u]; the OPS produced is, after reciprocal-
polynomial flip, the diagonal [M/M] Padé denominator of g(z, u).

Why this functional is never positive-definite

For a positive-definite moment functional over ℝ, the OPS has all real
simple zeros — proof: P(n, ·) is the characteristic polynomial of the
leading principal submatrix J(n) of the Jacobi matrix, which is real
symmetric tridiagonal with positive off-diagonals, hence has real simple
spectrum. Real zeros of P(n, ·) mean real poles of g(z, u) after the
flip.

But for the fractional Riccati, g(z, u) has no real poles for any
u in the §9.5 contraction regime. The Volterra contraction estimate
(§9.5 of the companion PDF) bounds g(·, u) on every real interval and
forces all poles off the real axis. Equivalently, the Stahl compact
Δ_g(u) on which the poles condense is a complex set disjoint from ℝ
for every such u.

So the Riccati moment functional is forced to live in the
quasi-definite, signed regime — the squared norms h(n, u) ∈ ℂ[u]
are sign-indefinite as functions of u, the Hankel determinants vanish
on real codimension-1 strata of u-space, and the OPS zeros are
genuinely complex. The class should not expose an isPositiveDefinite()
predicate; it should expose quasi-definiteness diagnostics
(firstHankelVanishingLocus() over u-space) instead.

Class skeletons

Parent class — polynomial-coefficient recurrence

package arb.functions.polynomials.orthogonal.complex;

import arb.Complex;
import arb.ComplexPolynomial;
import arb.functions.polynomials.orthogonal
       .RecurrentlyGeneratedOrthogonalPolynomialSequence;

/**
 * Common parent for any complex-codomain OPS whose three-term recurrence
 * coefficients live in ℂ[u] rather than ℂ. Differs from the existing
 * ComplexRecurrentlyGeneratedOrthogonalPolynomialSequence only in the
 * recurrence-coefficient ring R.
 */
public abstract class
ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence
    extends RecurrentlyGeneratedOrthogonalPolynomialSequence
            <ComplexPolynomial, Complex, ComplexPolynomial>
{
    protected ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence(
        int bits, String A, String B, String C)
    {
        super(bits, A, B, C);
    }
}

Abstract base — polynomial Chebyshev/Wheeler in ℂ[u]

package arb.functions.polynomials.orthogonal.complex;

import arb.ComplexPolynomial;
import arb.functions.integer.ComplexPolynomialSequence;

/**
 * The orthogonal polynomial sequence (Chihara, 1978, Ch. I) of a
 * quasi-definite moment functional over ℂ[u], constructed from its
 * moment sequence via the polynomial Chebyshev/Wheeler recurrence
 * (Theorem 5.1 of docs/SolutionMethodologyForConstantCoefficientFractionalRiccatiEquations.md).
 *
 * Subclasses supply momentSequence(): ℤ_{≥0} → ℂ[u].
 *
 * Makes no assumption of positive-definiteness — the squared norms h(n)
 * may be sign-indefinite as elements of ℂ[u]. For classical OPS families
 * the (A, B, C) are usually elementary functions of n and bypass this
 * class entirely. This class is for the case where the moments are
 * known (closed-form or recursive) but the OPS recurrence is not
 * elementary in n.
 *
 * Implementation: maintains the generating polynomial sequence S(n, z)
 * with S(n, z) ∈ ℂ[u][z]. At step n: S(n, z) and S(n-1, z) are known;
 * β(n) is the z^{n-1}-coefficient of S(n, z) / (z · S(n-1, z)) read off
 * after one acb_poly_div_series in ℂ[u][[z]] truncated at z^{n}; α(n) is
 * the z^n-coefficient of the recurrence; S(n+1, z) is then the
 * polynomial linear combination of S(n, z) / z, S(n, z), S(n-1, z) with
 * coefficients (1, -α(n), -β(n)).
 *
 * No scalar coefficient extractions in ℂ[u] are performed — α(n) and
 * β(n) remain opaque elements of ℂ[u] throughout. The only coefficient
 * reads are on the outer (z) level.
 */
public abstract class OrthogonalPolynomialMomentFunctionalSequence extends
    ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence
{
    public OrthogonalPolynomialMomentFunctionalSequence(int bits)
    {
        super(bits,
              "1",                                // A(n)
              "S(n+1, z)/S(n, z)",                // B(n) = -α(n) = z^n-coeff of recurrence
              "S(n, z)/(z·S(n-1, z))");           // C(n) =  β(n) = z^{n-1}-coeff of recurrence

        context.register("m", momentSequence());

        // S(n, z) as ComplexPolynomialSequence, indexed by n.
        // Each S(n, z) is a polynomial in z whose coefficients are in ℂ[u],
        // truncated at degree 2M+2 (the working order).
        context.register("S", ComplexPolynomialSequence.express(
            "n ➔ when(n = -1, 1,"
          + "        n =  0, sum(j ➔ m(j)·z^j {j=0..2*M+1}),"
          + "        else,   S(n-1)/z - α(n-1)·S(n-1) - β(n-1)·S(n-2))",
            context));

        p0.set("1");
        p1.set("x - m(0)");
    }

    protected abstract ComplexPolynomialSequence momentSequence();
}

Fractional Riccati subclass

package arb.functions.polynomials.orthogonal.complex;

import arb.ComplexPolynomial;
import arb.Real;
import arb.functions.complex.RiccatiMuntzPadeFunctional;
import arb.functions.complex.ComplexPolynomialNullaryFunction;
import arb.functions.integer.ComplexPolynomialSequence;

/**
 * OPS of the moment functional
 *   𝓛_Riccati[x^k; u] := a(k+1, u)
 * where a(k, u) are the Müntz–Tau coefficients of the constant-coefficient
 * fractional Riccati equation
 *   D^μ y(t,u) = P(u) + Q(u)·y(t,u) + R(u)·y(t,u)²,    I^{1-μ} y(0,u) = 0,
 * with P(u), Q(u), R(u) ∈ ℂ[u] and μ ∈ (0,1].
 *
 * The functional is quasi-definite and signed over ℂ[u]; never positive-
 * definite at any u ∈ ℝ for which the §9.5 contraction hypotheses hold
 * (Volterra contraction → no real poles of g(·, u)). OPS zeros are complex;
 * their reciprocals are the poles of the meromorphic generating function
 * g(z, u), condensing on the Stahl compact Δ_g(u) ⊂ ℂ \ ℝ.
 */
public class OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence
    extends OrthogonalPolynomialMomentFunctionalSequence
{
    private final RiccatiMuntzPadeFunctional muntz;

    public OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence(
        int bits, Real μ,
        ComplexPolynomialNullaryFunction P,
        ComplexPolynomialNullaryFunction Q,
        ComplexPolynomialNullaryFunction R)
    {
        this(bits, new RiccatiMuntzPadeFunctional(μ, P, Q, R));
    }

    public OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence(
        int bits, RiccatiMuntzPadeFunctional muntz)
    {
        super(bits);
        this.muntz = muntz;
    }

    @Override
    protected ComplexPolynomialSequence momentSequence() {
        // m(k, u) = a(k+1, u), shift the Müntz coefficient sequence by one.
        return muntz.muntzBasis().shift(1);
    }
}

Padé numerator sibling

The Padé numerators satisfy the same three-term recurrence with
different seeds:

P^(1)(-1, x) = -1,    P^(1)(0, x) = 0.

Add a sibling subclass reusing the same α(n), β(n) machinery with
these shifted seeds. With the reciprocal-polynomial flip
Q̂(M, z) = z^M · P(M, 1/z) and
P̂(M, z) = z^M · P^(1)(M, 1/z), this assembles the full diagonal
[M/M] Padé approximant

R̂(M, z, u) = z · P̂(M, z) / Q̂(M, z).

Required native binding

The only native-side gap initially identified was
arb_poly_compose / acb_poly_compose, needed by an earlier
formulation that manipulated m(k, u) as a polynomial in u
inside the S(n, k) scalar array. The revised polynomial
formulation above does not hit this gap — every step is a single
acb_poly_* call at the z-level and the u-level is carried opaquely
through ComplexPolynomial arithmetic. If the binding is added later
for unrelated reasons, the OPS construction here will not need to
change.

ComplexPolynomial.reciprocalFlip() (coefficient reversal of degree M)
is still needed to assemble the Padé denominator in one line.

Validation oracles

1. Polynomial vs scalar Chebyshev/Wheeler. Pick a fixed u_0 ∈ ℂ,
   substitute into the moment sequence to get a scalar m(k, u_0) ∈ ℂ,
   run a scalar Chebyshev/Wheeler recursion to produce reference
   (α^scalar(n), β^scalar(n)) ∈ ℂ, and compare with α(n)(u_0),
   β(n)(u_0) from the polynomial recurrence above. Verified
   numerically: the two paths agree to ~10^{-47} at 128-bit
   precision.

2. Padé identity in the orthogonality block. Compute
   r(z) := g(z) · Q̂(M, z) − z · P̂(M, z) at a numerical u_0. The
   coefficients [z^{M+1}], …, [z^{2M+1}] of r(z) must vanish to
   working precision. Verified numerically: residue ~10^{-53}.
   (Low-order coefficients [z^1], …, [z^M] are the expected free part
   of the Padé identity matched by the numerator construction and are
   not part of the orthogonality block.)

3. R(u) ≡ 0 Mittag-Leffler oracle (Remark 9.14 of the PDF).
   Moments collapse to
   m(k, u) = P(u) · Q(u)^k / Γ((k+1) μ + 1); the [M/M] Padé must
   reproduce y(t, u) = (P(u) / Q(u)) · (E_{μ, 1}(Q(u) · t^μ) − 1) to
   working precision.

4. μ ≡ 1 classical Riccati oracle (Remark 9.15 of the PDF).
   Reduces to the classical Riccati Taylor expansion; the [M/M] Padé
   is exact for M ≥ deg of the rational closed-form solution.

5. Volterra cross-check. Pick fixed u_0; compute y(t, u_0) by
   direct high-precision Volterra discretization, compare with
   R̂(M, t^μ, u_0) for M = 8, 16, 32. Expect geometric convergence.

6. Pole locus on Stahl compact. Eigenvalues of the non-symmetric
   tridiagonal Jacobi matrix J(M, u_0) (subdiagonal β(n),
   superdiagonal 1, diagonal α(n); see Theorem 7.2 of the
   methodology doc) should condense on Δ_g(u_0) ⊂ ℂ \ ℝ as M → ∞.
   The symmetrization β(n) ↦ √β(n) of the positive-definite case is
   unavailable because β(n) is complex-valued.

7. Quasi-definiteness audit. Compute h(n, u) = ⟨P(n, ·), P(n, ·)⟩ ∈ ℂ[u] for n ≤ M. Locate (symbolically or numerically) the real
   codimension-1 strata in u-space where h(n, u) = 0 and verify the
   working u_0 stays off them.

Future work: O(n) specialization for the Riccati case

The polynomial Chebyshev/Wheeler is Θ(M²) ring ops in ℂ[u] and is
universal — works for any moment sequence. For the Riccati moment
functional specifically, additional structure is available that should
permit an O(n) direct recurrence on (α(n), β(n)):

The R(u) ≡ 0 reduction is sharp. Hankel factorization
H(n, u) = P(u) · D(u) · G(u) · D(u) with
D(u) = diag(1, Q(u), Q(u)², …, Q(u)^n) and
G(u) = [γ(i+j, u)] (the pure-Gamma Hankel matrix with
γ(k, u) = 1 / Γ((k+1) μ + 1)) gives

α(n, u) = Q(u)  · α^Γ(n, μ)
β(n, u) = Q(u)² · β^Γ(n, μ)

where α^Γ(n, ·), β^Γ(n, ·) are the recurrence coefficients of the
Mittag-Leffler moment functional with moments
γ(k, μ) = 1 / Γ((k+1) μ + 1) — a one-parameter family in μ known in
the literature on Padé approximation of Mittag-Leffler functions
(Stahl 2003; Wimp 1981). P(u) drops out of (α(n), β(n)) entirely
(global scale on the functional); the u-dependence enters only
through Q(u).

General R(u) ≢ 0. Magnus's framework (Magnus 1995; Van Assche
2018) for semi-classical orthogonal polynomials predicts that
(α(n, u), β(n, u)) satisfy a finite-order discrete nonlinear
recurrence in n — a Laguerre–Freud / string equation — whose form
is dictated by the (fractional) Riccati ODE on the moments. Cost
reduction: Θ(M²) → O(M) ring operations in ℂ[u]. Tracked
separately as a follow-up.

The polynomial Chebyshev/Wheeler in this issue stays as:

• the reference oracle that the fast path will be validated against;
• the universal fallback for any other moment sequence the codebase
  ever wants to support.

Open questions / decision points

1. Division semantics in ℂ[u]. B(n), C(n) are ratios of
   polynomials in u. Default: RationalFunction<ComplexPolynomial>,
   clear denominators at evaluation. Switch to Geronimus division-free
   form if denominator growth bites.

2. Precision propagation. Chebyshev cancellation in S(n, z) is
   catastrophic at high n if bits is too low; bits should grow at
   least linearly with M. Provide recommendedBits(int M).

3. Truncation level of S(n, z). S(n, z) only needs to be known
   modulo z^{2M+2} for an [M/M] Padé. Provide truncationLevel(M)
   and store S(n, z) truncated at that level.

Deliverables checklist

• Add
  ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence
  parent class.
• Add OrthogonalPolynomialMomentFunctionalSequence abstract class
  with the polynomial Chebyshev/Wheeler recurrence on S(n, z).
• Add OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence
  concrete class, accepting either (μ, P, Q, R) directly or an
  existing RiccatiMuntzPadeFunctional. In the former case the
  constructor builds a RiccatiMuntzPadeFunctional internally and
  reuses its muntzBasis(); in the latter the Müntz coefficient
  sequence is shared with the caller.
• Add Padé-numerator sibling class (same α, β; different seeds).
• Add ComplexPolynomial.reciprocalFlip().
• Unit test: polynomial Chebyshev/Wheeler agrees with a scalar
  Chebyshev/Wheeler reference at a fixed u_0 to working precision.
• Unit test: Padé identity g · Q̂ − z · P̂ = O(z^{2M+2}) holds in
  the orthogonality block z^{M+1}, …, z^{2M+1}.
• Unit test: R(u) ≡ 0 Mittag-Leffler oracle.
• Unit test: μ ≡ 1 classical Riccati oracle.
• Unit test: Volterra cross-check on a Heston-style
  (P(u), Q(u), R(u), μ).
• Unit test: J(M, u_0) has no real eigenvalues at working precision.
• Unit test: J(M, u_0) eigenvalues condense on Δ_g(u_0) as
  M → ∞.
• Decide division semantics in ℂ[u].
• Document precision-vs-M scaling in Javadoc.
• Follow-up issue: derive the Laguerre–Freud / string equation on
  (α(n, u), β(n, u)) for the Riccati moment functional; implement
  FastOrthogonalPolynomialFractionalRiccatiMomentFunctionalSequence
  at O(M) cost.

§9.5 Lemma A corrigendum

High-precision maximum of the scalar function
h(μ) = Γ(μ + 1) / Γ(2μ + 1) on the parameter range μ ∈ (0, 1] is

μ* = 0.14503474316667…,    h(μ*) = 1.03967510771…

(not μ* ≈ 0.1339, h(μ*) ≈ 1.0428 as written in the companion PDF).
The 21/20 uniform bound and the role of C_μ in the majorant ODE are
unchanged.

References

• Chihara, T. S. (1978). An Introduction to Orthogonal Polynomials.
  Ch. I: the OPS of a quasi-definite moment functional; Favard's theorem.
• Akhiezer, N. I. (1965). The Classical Moment Problem.
• Gautschi, W. (1970, 1982). On generating orthogonal polynomials.
  Chebyshev / Wheeler recursion analysis.
• Gragg, W. B. (1974). Matrix interpretations and applications of
  the continued fraction algorithm. Rocky Mountain J. Math. 4,
  213–225. OPS ↔ Padé denominator equivalence.
• Brent, R. P., Gustavson, F. G., Yun, D. Y. Y. (1980). Fast
  solution of Toeplitz systems of equations and computation of Padé
  approximants. J. Algorithms 1, 259–295.
• Deift, P. (1999). Orthogonal Polynomials and Random Matrices: A
  Riemann–Hilbert Approach.
• Magnus, A. P. (1995). Painlevé-type differential equations for
  the recurrence coefficients of semi-classical orthogonal polynomials.
  J. Comput. Appl. Math. 57, 215–237.
• Van Assche, W. (2018). Orthogonal Polynomials and Painlevé
  Equations. Cambridge University Press.
• Stahl, H. (2003). Spurious poles in Padé approximation of
  algebraic and Mittag-Leffler functions.
• Wimp, J. (1981). Computation with Recurrence Relations.
• §9 of docs/ThePadeMuntzSpectralTauMethodForConstantCoeffecientFractionalRiccations.pdf
  — Müntz–Tau Puiseux series, meromorphy of g, Riemann–Hilbert
  representation, Stahl compact Δ_g ⊂ ℂ \ ℝ.
• docs/SolutionMethodologyForConstantCoefficientFractionalRiccatiEquations.md
  — full statement and proof of the polynomial Chebyshev/Wheeler
  recurrence (Theorem 5.1) used in this issue.

[crowlogic]

Scope reduction (no multivariate dependency)

#1022 (acb_mpoly bindings) is not needed for this work. The evaluator boils any rough Heston parameterisation down to:

• P(v) — polynomial in v, degree ≤ 2
• Q(v) — polynomial in v, degree ≤ 1
• R(v) — polynomial in v, degree 0 (scalar)

All three are univariate ComplexPolynomial instances in the single variable v (the Fourier / Riccati parameter). The Wheeler recurrence is therefore expressed entirely over ComplexPolynomial (no multivariate ring), with sequences S, β, T, α : ℤ → ℂ[v]. The OPS implementation in OrthogonalPolynomialMomentFunctionalSequence already lives in this univariate-coefficient world via ComplexPolynomialSequence and ComplexPolynomialCoefficientRecurrentlyGeneratedOrthogonalPolynomialSequence.

This is also more efficient: no multivariate polynomial arithmetic overhead, smaller bytecode, no acb_mpoly allocation churn.

Current state (wip/1021-polynomial-recurrence)

Full test suite passes 805/805 after the four commits:

• 056e4c4 — thread shared Context through OPS so m's expression survives
• 4ffea52 — filter propagateContextVariablesByReference by target's declared fields
• a0bf6ad — per-emission guard in linkSharedVariableToReferencedFunction
• d5d7478 — implement ComplexPolynomial.set(Fraction) (was an empty stub)

The OPS smoke test (OrthogonalPolynomialFractionalRiccatiMomentFunctionalSequenceTest.testConstructionAndFirstFewTerms) now constructs the full Wheeler cluster (m, β, T, α, S) over the shared muntz context and evaluates m(0..2), S(0..2) without faulting.

Remaining work per the original TODO:

• Validation oracles for the Wheeler β(n), α(n) coefficients against a hand-computed small-n table
• Smoke test that exercises the parent class's P(n) three-term recurrence (B(n) = −α(n), C(n) = β(n)) end-to-end