Intrinsic primitive via discovered OPS — given only f, produce F via uniformly-convergent OPS expansion in the discovered measure dμ_f

[crowlogic]

Proposal

Compute the primitive F(x) = ∫ f(t) dt of an arbitrary integrand f (real, complex, or sign-changing) on its domain by intrinsically discovering — from f alone, with no externally-supplied measure, basis, or auxiliary data — the unique (up to the integration constant +C) measure dμ_f and associated orthogonal polynomial sequence (OPS) {p_n^{(f)}} such that the partial sums

F_N(x) = Σ_{n=0}^{N} c_n · p_n^{(f)}(x)

converge uniformly to F on the domain.

The construction is intrinsic to f: the only input is the integrand. The measure dμ_f is discovered, not selected. The OPS, the three-term recurrence coefficients (β_j, γ_j), the Jacobi matrix, the expansion coefficients c_n of F, and ultimately F itself, are all produced by a single Θ(M²) sweep driven by the existing MuntzPadeApproximant σ-sliding-window engine (issue #1016).

Mathematical setup

Given f on a compact (or compactified) domain D, the primitive F has its own moment sequence

m_k = ∫_D t^k · dν_F(t)        k = 0, 1, 2, …

where dν_F is the unique (up to +C) measure/distribution whose moments characterize F via the Hamburger–Nevanlinna correspondence. Equivalently, the Cauchy transform of dν_F,

ĝ_F(z) = ∫_D dν_F(t) / (z − t)

is determined by f alone (one integration of f against the geometric kernel Σ t^k / z^{k+1} produces the moments m_k, term by term).

On a compact domain, the moment problem for dν_F is determinate, so dν_F is unique. The OPS {p_n^{(f)}} orthogonal w.r.t. dν_F is then unique by Gram–Schmidt. The (β_j, γ_j) of that OPS is what the σ-window machinery already in MuntzPadeApproximant computes from the moment sequence.

The expansion coefficients of F in the OPS basis,

c_n = ⟨F, p_n^{(f)}⟩_{ν_F} / ⟨p_n^{(f)}, p_n^{(f)}⟩_{ν_F}

reduce to a finite expression in f and (β_j, γ_j) via integration by parts (F' = f is the integrand we started with — no extra data needed). The partial sums F_N converge uniformly to F on D by the standard OPS-expansion theorem upgraded by Jackson-type smoothness: F is one degree smoother than f, so its coefficients decay one power faster than the corresponding expansion of f — uniform convergence is automatic whenever f's expansion is L²-convergent, which on a compact domain is automatic for f ∈ L^1.

The signed / complex case (sign-changing f, complex-valued f) is handled identically — dν_F becomes signed or complex, the OPS exists formally via the same recurrence, and h_j ∈ ℂ instead of h_j > 0. The only new failure mode is h_j straddling zero, which is exactly the HankelDegeneracyException / PrecisionExhaustedException pair the σ engine already throws.

Architecture

One new class — call it PrimitiveByOPS (final name negotiable) — composed on top of the existing MuntzPadeApproximant σ engine. No new native code, no new SWIG, no new matrix machinery: every piece of computation is a ComplexPolynomial / Complex operation already in arb4j.

public final class PrimitiveByOPS implements ComplexFunction, AutoCloseable
{
  PrimitiveByOPS(ComplexFunction f, Complex domainStart, Complex domainEnd,
                 Complex basePoint, int workingBits) { … }

  @Override Complex evaluate(Complex x, int order, int bits, Complex result);
}

Pipeline (all driven from f alone)

1. Moment extraction. Compute m_k = ∫_{a..b} t^k · f(t) dt for k = 0..2M−1 to working precision. M is discovered adaptively, same ΦBound-style monotone-refinement loop as MuntzPadeApproximant.evaluate. The integration of t^k · f(t) is the only place f is touched as a callable; from this point forward everything is a finite computation on the moment sequence.

   • For analytic f given as an arb4j expression — moments are computable symbolically against monomials in many cases; fall back to Gauss-Legendre or Clenshaw–Curtis quadrature on the moment integrand otherwise.
   • For f given as a callable ComplexFunction — moments by adaptive quadrature at workingBits. The quadrature precision is the bottleneck; everything downstream is exact in the moment sequence.

2. σ recurrence. Hand the moment sequence to the existing σ engine. Out come (β_j, γ_j) for j = 1..M, the OPS {p_n^{(f)}} implicitly via the Jacobi matrix J = diag(β_j) + offdiag(√γ_j), and the monic OPS polynomials Q_M if needed. This is literally the code from MuntzPadeApproximant.growOne(). The plan is to factor that code out of MuntzPadeApproximant into a reusable ChebyshevRecurrenceState helper class that both MuntzPadeApproximant and PrimitiveByOPS instantiate. Internal-only refactor; the public API of MuntzPadeApproximant is unaffected.

3. Coefficient extraction c_n of F. Integration by parts in the OPS inner product, using F' = f:

   c_n · h_n = ⟨F, p_n⟩ = boundary terms involving F(a), F(b)
                            + ⟨f, P̃_n⟩
   

   where P̃_n is the antiderivative of p_n (a polynomial of degree n+1, expressible in the OPS basis via the (β, γ) recurrence). The boundary F(a) is the +C ambiguity — fixed by choosing the base point at which F is declared zero. F(b) is then Σ c_n · p_n(b), consistent with the expansion.

   Each c_n is a finite linear combination of the moments m_k and the recurrence coefficients — Θ(n) work per coefficient, Θ(M²) total.

4. Clenshaw evaluation. F(x) = Σ_{n=0}^{M} c_n · p_n^{(f)}(x) evaluated by Clenshaw recurrence using (β_j, γ_j) directly — no need to ever materialize the polynomials. Θ(M) per evaluate(x, …) call.

Adaptive convergence

Same shape as MuntzPadeApproximant.evaluate: monotone-refinement loop on M, with a ΦBound-style estimator measuring |F_M(x) − F_{M−1}(x)|² / (|F_{M−1}(x) − F_{M−2}(x)| − |F_M(x) − F_{M−1}(x)|) (Aitken acceleration of the convergence rate). Stop when below 2^{−bits}.

Spectral convergence (exponential in M for analytic f) means typical M ≈ 10–30 at 128 bits, same regime as the Müntz–Padé path.

Failure modes

Inherited from the σ engine:

• HankelDegeneracyException(j) — moment functional indefinite at order j for sign-changing f whose moments happen to make h_j = 0 exactly. Rare; signals the expansion has a singular pivot at that order. Caller can retry with a slightly perturbed base point or domain.
• PrecisionExhaustedException(j, bits) — h_j straddles zero at working precision. Caller retries with workingBits ← 2·workingBits.

The "+C ambiguity" is not a failure mode — it's a free parameter exposed as the basePoint constructor argument.

Net result

The user calls

try (PrimitiveByOPS F = new PrimitiveByOPS(f, a, b, basePoint, bits)) {
    Complex value = F.evaluate(x, 1, bits, null);   // F(x) − F(basePoint)
}

— no measure supplied, no basis chosen, no quadrature rule selected. The integrand f and the domain go in; the primitive comes out, evaluable at any point in the domain, with the OPS / measure / Jacobi matrix discovered intrinsically.

Out of scope (for this issue)

• Multi-dimensional integrands (would require multivariate OPS — separate construction).
• Singular integrands at the endpoints (would require a pre-extracted Jacobi-type weight; not "intrinsic to f" anymore — separate issue if it ever comes up).
• The Θ(M log² M) half-GCD recurrence path (out-of-scope for Replace Hankel bordered-inverse with Chebyshev orthogonal-polynomial recurrence (sliding-window σ-table) #1016, same reasoning here).

Implementation steps

1. Refactor MuntzPadeApproximant to extract the σ-window + Q/P state machine into a reusable ChebyshevRecurrenceState helper. No behavior change; MuntzPadeApproximant becomes a thin client of the helper. Full test suite re-run.
2. Add PrimitiveByOPS with the pipeline above. Quadrature for moment extraction lives behind a small interface so it can be swapped (initial impl: adaptive Gauss-Legendre at working precision).
3. Add the Clenshaw evaluate path driven by (β_j, γ_j).
4. Add an integration test: pick f(t) = sin(t) on [0, π], compare F(π) to 2 and F(π/2) to 1; pick f(t) = e^{−t²} on [−4, 4], compare F(∞) − F(−∞) to √π; pick a sign-changing analytic f and verify uniform error below 2^{−bits} at a grid of test points.
5. Document in docs/ how the construction discovers dμ_f from f alone.

References

• Stahl & Totik, General Orthogonal Polynomials (1992), Ch. 4 (Markov's theorem, uniform convergence of Padé convergents off the support).
• Akhiezer, The Classical Moment Problem (1965), Ch. 1–2 (Hamburger / Nevanlinna correspondence; moment-problem determinacy on compact support).
• Gautschi, Orthogonal Polynomials: Computation and Approximation (2004), Ch. 2 (modified Chebyshev algorithm = the σ-sliding-window engine).
• Crowley, qrh.tex §9.7 (the σ recurrence as implemented in MuntzPadeApproximant, issue Replace Hankel bordered-inverse with Chebyshev orthogonal-polynomial recurrence (sliding-window σ-table) #1016).

[crowlogic]

Intrinsic Primitive via Discovered OPS — Full Algorithm Derivation

Problem: Given only f, produce its antiderivative F via a uniformly-convergent OPS expansion in the discovered spectral measure dμ_f — without any prescribed weight or domain.

---

Part I — OPS and the Three-Term Recurrence

Let μ be a positive Borel measure on (a,b) ⊆ ℝ with finite moments of all orders. The associated OPS {Pₙ} satisfies

⟨Pₘ, Pₙ⟩_μ = ∫_a^b Pₘ(x) Pₙ(x) dμ(x) = hₙ δₘₙ, hₙ > 0

Every such family obeys Favard's three-term recurrence:

Pₙ₊₁(x) = (x − αₙ) Pₙ(x) − βₙ Pₙ₋₁(x),  n ≥ 1     (1)

with P₋₁ ≡ 0, P₀ ≡ 1, αₙ ∈ ℝ, βₙ > 0. The pair ({αₙ}, {βₙ}) uniquely determines μ when the moment problem is determinate.

The symmetric tridiagonal Jacobi matrix J has diagonal (α₀, α₁, α₂, …) and super/sub-diagonal (√β₁, √β₂, √β₃, …). Its spectral measure with respect to e₀ is precisely μ.

Theorem (Integration is banded): The indefinite integration operator 𝒥: f ↦ ∫_a^x f in OPS coefficient space is a banded matrix of bandwidth 2 — making integration O(N).

---

Part II — Measure Discovery: Stieltjes–Lanczos

The Stieltjes procedure recovers αₙ, βₙ from inner-product evaluations alone:

p₋₁ = 0,  p₀ = 1

for n = 0 to N−1:
    αₙ    = ⟨x·pₙ, pₙ⟩ / ⟨pₙ, pₙ⟩
    pₙ₊₁  = (x − αₙ)·pₙ − (‖pₙ‖²/‖pₙ₋₁‖²)·pₙ₋₁
    βₙ₊₁  = ‖pₙ₊₁‖² / ‖pₙ‖²

This is identical to Lanczos iteration on the multiplication operator Mₓ: g ↦ xg with initial vector f/‖f‖. The discovered measure is the spectral measure of the resulting Jacobi matrix J_N:

dμ_f ≈ Σₖ wₖ δ(λₖ)     (4)

where λₖ are eigenvalues of J_N and wₖ = (Sₖ₁)² are the squares of the first components of the normalized eigenvectors.

---

Part III — Expansion Coefficients

f = Σₙ cₙ Pₙ,   cₙ = ∫_a^b f(x) Pₙ(x) dμ_f(x)     (5)

Computed via Gauss quadrature adapted to dμ_f:

cₙ ≈ Σₖ f(λₖ) Pₙ(λₖ) wₖ     (6)

Exact for polynomials of degree ≤ 2N−1; exponentially convergent for analytic f.

---

Part IV — Term-by-Term Integration

Define 𝒬ₙ(x) = ∫_x₀^x Pₙ(t) dt. From the recurrence identity

x Pₙ = Pₙ₊₁ + αₙ Pₙ + βₙ Pₙ₋₁     (13)

integrating both sides from x₀ to x yields the integration recurrence:

𝒬ₙ₊₁ = (1/βₙ₊₁) · [ ∫_x₀^x t·Pₙ(t) dt  −  αₙ 𝒬ₙ  −  βₙ 𝒬ₙ₋₁ ]     (15)

The integral ∫ t·Pₙ dt on the right is resolved by integration by parts, involving only 𝒬ₙ, 𝒬ₙ₋₁, and boundary values Pₙ(x), Pₙ(x₀) — preserving the three-term structure. The antiderivative is:

F(x) = F(x₀) + Σₙ cₙ 𝒬ₙ(x)     (9)

Chebyshev special case (αₙ = 0, βₙ = 1/4, arcsine measure on [−1,1]):

∫ Tₙ dx = Tₙ₊₁ / (2(n+1))  −  Tₙ₋₁ / (2(n−1)),  n ≥ 2     (11)
bₙ = (cₙ₋₁ − cₙ₊₁) / (2n),  n ≥ 1     (12)

This is the Clenshaw–Curtis integration formula. Recurrence (15) reduces exactly to (12) when αₙ = 0, βₙ = 1/4.

---

Part V — Uniform Convergence

Key inequality: For partial sums S_N f = Σₙ₌₀ᴺ cₙ Pₙ and F_N = ∫_x₀^x S_N f,

‖F − F_N‖_∞ ≤ (b − a) · ‖f − S_N f‖_∞     (17)

Uniform convergence of F_N → F follows directly from uniform convergence of S_N f → f.

Sufficient conditions:

• f ∈ C([a,b])
• Lebesgue constants Λ_N = O(log N) — holds for Chebyshev; verified numerically for discovered measures
• For analytic f extending to a Bernstein ellipse of parameter ρ > 1: ‖f − S_N f‖_∞ = O(ρ⁻ᴺ), hence ‖F − F_N‖_∞ = O(ρ⁻ᴺ)

Theorem: Let f ∈ C([a,b]) and suppose S_N f → f uniformly. Then F_N → F uniformly.

Proof: |F(x) − F_N(x)| = |∫_x₀^x (f − S_N f) dt| ≤ (b−a) ‖f − S_N f‖_∞ → 0 uniformly in x. □

Certified error bound (arb4j): With cₙ computed as balls [cₙ_mid ± εₙ],

‖F − F_N‖_∞  ≤  (b−a) Σₙ₌ₙ₊₁^∞ |cₙ| ‖Pₙ‖_∞  +  Σₙ₌₀ᴺ εₙ ‖𝒬ₙ‖_∞     (20)

Both sums are computable: the first decays geometrically for analytic f; the second is the accumulated ball-arithmetic radius.

---

Part VI — Complete Algorithm

ALGORITHM: IntrinsicPrimitive(f, [a,b], x₀, N)

STEP 1 — MEASURE DISCOVERY (Stieltjes–Lanczos)
  p₋₁ = 0,  p₀ = 1
  inner products via high-order Gauss–Legendre in ball arithmetic
  for n = 0 to N−1:
    αₙ    = ⟨x·pₙ, pₙ⟩ / ⟨pₙ, pₙ⟩
    pₙ₊₁  = (x − αₙ)·pₙ − (‖pₙ‖²/‖pₙ₋₁‖²)·pₙ₋₁
    βₙ₊₁  = ‖pₙ₊₁‖² / ‖pₙ‖²
  Build Jacobi matrix J_N from (α₀…αₙ₋₁), (β₁…βₙ)

STEP 2 — EXPANSION COEFFICIENTS
  {λₖ, wₖ} = eigendecomposition of J_N
  cₙ = Σₖ f(λₖ) · Pₙ(λₖ) · wₖ    for n = 0..N

STEP 3 — ANTIDERIVATIVE COEFFICIENTS
  𝒬₋₁ = 0,  𝒬₀ = (x − x₀)
  for n = 0 to N−1:
    𝒬ₙ₊₁ = (1/βₙ₊₁) · [∫_x₀^x t·Pₙ dt  −  αₙ·𝒬ₙ  −  βₙ·𝒬ₙ₋₁]
  F_N(x) = F(x₀) + Σₙ₌₀ᴺ cₙ · 𝒬ₙ(x)

STEP 4 — CERTIFIED ERROR BOUND
  tail  = (b−a) · Σₙ₌ₙ₊₁^{2N} |cₙ|     (geometric decay estimate)
  arith = Σₙ₌₀ᴺ rad(cₙ) · ‖𝒬ₙ‖_∞
  err   = tail + arith

RETURN F_N, err

Complexity: O(N²) overall — Stieltjes dominates (N inner products each O(N)); Steps 2–3 are O(N).

---

Part VII — Edge Cases

• Determinacy: Carleman's condition Σₙ (β₁···βₙ)^(−1/2) = ∞ guarantees a unique dμ_f; holds automatically for compactly supported μ. For entire f with super-exponential growth, use the Friedrichs extension of J.
• Non-compact support: Three-term recurrence still holds on ℝ, but uniform convergence is only guaranteed on compact sub-intervals; bound (17) applies locally.
• Inner product choice: Natural choice is ⟨g, h⟩ = ∫ |f| g h dx; Lebesgue measure on [a,b] also works and relies on rapid coefficient decay from the structure of f.
• Stability: Ball arithmetic flags orthogonality loss via radius inflation. Periodic re-orthogonalization (cost: factor of 2 in radius) restores certified orthogonality.

---

Part VIII — arb4j Infrastructure Map

Component	arb4j class/method	Role
Ball arithmetic	RealPolynomial, arb_poly_*	All recurrences and inner products
Gauss quadrature	RealPolynomials.gaussianQuadrature	Inner product evaluation (Step 1)
Jacobi eigenvalues	Jacobi matrix eigensolver	Gauss nodes λₖ (Step 2)
Three-term recurrence	OrthogonalPolynomial.evaluate	Pₙ(λₖ) evaluation (Step 2)
Polynomial arithmetic	RealPolynomial.add/multiply	𝒬ₙ construction (Step 3)
Radius tracking	Ball midpoint/radius	Certified bound (Step 4)

Three new components needed:

1. Stieltjes driver — iterative inner-product loop producing J_N from f
2. Integration recurrence loop — Step 3, using αₙ, βₙ from Step 1
3. Tail bound estimator — fits geometric decay to {|cₙ|}, computes series remainder

---

Conclusion

The intrinsic primitive is a corollary of three classical results:

1. Favard's theorem — every OPS arises from a measure via its Jacobi matrix
2. Tridiagonality of the integration operator in OPS coefficient space
3. Term-by-term integration of uniformly convergent series

The discovered-measure strategy makes convergence optimal: f is by construction well-represented in L²(dμ_f), and for analytic f coefficients decay as O(ρ⁻ᴺ). All steps are compatible with arb4j's ball-arithmetic model, and the certified error bound (20) is computable from the expansion coefficients alone — making this a fully rigorous antiderivative oracle.