In a previous post, I showed how to compute $\zeta(3)$ using polynomial ODEs — but the construction required transcendental initial values like $e$ and $1/(e-1)$. Here’s a different approach: start from a generating function whose coefficients encode $\zeta(3)$, derive the ODE it satisfies, and see what goes right and what goes wrong.

The Generating Function

The coefficients $a_n = \frac{(-1)^{n-1}}{n^3\binom{2n}{n}}$ decay like $\frac{\sqrt{\pi n}}{n^3 \cdot 4^n}$ (since $\binom{2n}{n} \sim 4^n/\sqrt{\pi n}$), so the radius of convergence is $R = 4$.

This is a classical identity — it can be derived from the beta function integral $\frac{1}{\binom{2n}{n}} = \frac{n}{2} B(n, n) = \frac{n}{2}\int_0^1 t^{n-1}(1-t)^{n-1}\,dt$, then summing over $n$ and evaluating the resulting integrals. The key point for us is: $\zeta(3) = \frac{5}{2}F(1)$, and $x = 1$ is well inside the convergence disk.

So if $F(x)$ satisfies a nice ODE, we can compute $\zeta(3)$ by solving it from $x = 0$ to $x = 1$.

Step 1: The Recurrence

Since this ratio is a rational function of $n$, the generating function $F(x)$ is holonomic (D-finite): it satisfies a linear ODE with polynomial coefficients. The degree of the denominator (3) predicts a third-order equation.

Step 2: From Recurrence to ODE via the Euler Operator

This is the heart of the derivation. We convert the recurrence into a differential equation by multiplying by $x^n$ and summing.

Powers of $\theta$ multiply by powers of $n$: $\theta^k\left[\sum a_n x^n\right] = \sum n^k a_n x^n$. Explicitly:

(You can verify these by induction: $\theta^{k+1}(f) = x\frac{d}{dx}[\theta^k(f)]$.)

The other ingredient: multiplication by $x$ corresponds to the backward shift on coefficients ($a_n \mapsto a_{n-1}$), and division by $x$ corresponds to the forward shift ($a_n \mapsto a_{n+1}$, valid when $a_0 = 0$).

Summing the recurrence

Second sum. This is just $\theta^3 F$ (since $a_0 = 0$, the $n=0$ term contributes nothing).

Assembling the equation

Expanding in terms of derivatives

Substitute the expressions for $\theta^k$:

Adding:

The equation $(4\theta^3 - 2\theta^2 + x\theta^3)F = x$ becomes:

Dividing both sides by $x$:

A third-order linear ODE with integer polynomial coefficients and constant right-hand side.

Step 3: Initial Conditions

Wait — let me be careful. $a_2 = \frac{(-1)^1}{2^3\binom{4}{2}} = \frac{-1}{8\cdot 6} = \frac{-1}{48}$. So $F''(0) = 2a_2 = -\frac{1}{24}$.

All rational. Compare this to the integral-based construction from the previous post, where the initial values were $e$, $1/(e-1)$, $e/(e-1)$ — transcendental numbers that themselves require computation. Here, we start from $0$, $1/2$, $-1/24$. That’s a dramatic improvement.

Step 4: Converting to a Dynamical System

To compute $\zeta(3) = \frac{5}{2}F(1)$, we need to “walk” from $x = 0$ to $x = 1$ along the solution. The standard trick: substitute $x = 1 - e^{-t}$ so that $x$ increases from $0$ to $1$ as $t$ goes from $0$ to $\infty$. This gives $\dot{x} = 1 - x$.

Define state variables $y_0 = F$, $y_1 = F'$, $y_2 = F''$. The chain rule converts $x$-derivatives to $t$-derivatives: $\dot{y}_k = y_{k+1}\cdot\dot{x} = y_{k+1}(1-x)$. From the ODE, $F''' = \frac{1 - x(10+3x)y_2 - (2+x)y_1}{x^2(4+x)}$, so:

Numerical integration confirms: all trajectories are bounded, $y_0(t) \to F(1) = \frac{2}{5}\zeta(3)$, and the convergence rate is $e^{-t}$ — first-floor, the fastest possible.

So this system has:

1. Rational initial conditions.
2. Integer coefficients in the ODE.
3. First-floor convergence.
4. Bounded trajectories.

Everything you could want — except one thing.

The Obstacle: A Regular Singular Point

The equation for $\dot{y}_2$ contains $\frac{1}{x^2(4+x)}$, which diverges at $x = 0$ — exactly where we start. The right-hand side is not polynomial.

This is a regular singular point (Fuchsian singularity) of the ODE. The solution $F(x)$ is perfectly analytic at $x = 0$ — the singularity is an artifact of the coefficient structure, not a true blow-up. Concretely, at $x = 0$:

• Numerator of $F'''$: $1 - 0\cdot(\ldots)\cdot y_2 - 2\cdot y_1 = 1 - 2\cdot\frac{1}{2} = 0$.
• Denominator: $x^2(4+x) = 0$.

The $0/0$ limit exists: $F'''(0) = 6a_3 = 6\cdot\frac{1}{540} = \frac{1}{90}$. Numerically, starting from $x = \varepsilon$ works perfectly. But formally, the system is not a polynomial ODE.

Can We Remove the Singularity?

Operator factorization

This has only a simple zero at $x = 0$ (instead of the double zero $x^2$). Progress — but the relationship $H = F'$ and $G = xH' + H$ means the full system for $(F, H, G)$ still contains $\dot{H} = \frac{(G-H)(1-x)}{x}$, a division by $x$.

Time reparametrization

Replace $t$ with a new time variable satisfying $d\tau/dt = 1/x$. This absorbs the $1/x$ and makes the right-hand side polynomial:

(Here $w = 1/(4+x)$ is an auxiliary variable making the rational function polynomial.)

This system is polynomial, has rational initial conditions, and is bounded. But $x = 0$ is now a fixed point of $\dot{x} = x(1-x)$. The trajectory never leaves the origin. The system is frozen.

The indicial equation

The obstruction runs deeper than any particular trick. The indicial equation of the original ODE at $x = 0$ is $I(\rho) = 2\rho^2(2\rho - 1)$, with roots $\rho = 0$ (multiplicity 2) and $\rho = 1/2$.

The double root at $\rho = 0$ is the culprit. It means the ODE has a logarithmic solution $F_{\log}(x) = F_1(x)\log x + \tilde{F}(x)$ alongside the analytic solution. Attempts to resolve the singularity by blow-up (defining new variables to absorb the $1/x$ factors) lead to an infinite regression: each new variable introduces a new $0/0$ form at $x = 0$, requiring yet another blow-up. The chain never terminates.

If the indicial roots were all simple, a single blow-up would suffice. The multiplicity is the structural obstruction.

Composing Two Systems: The PIVP Compiler

The time reparametrization above gives a polynomial system with rational initial conditions — but frozen at $x = 0$. Meanwhile, the DNA32 integral approach gives a polynomial system that converges to $\zeta(3)$ — but requires transcendental initial conditions $e$ and $1/(e-1)$.

If we could run this “Block A” to completion first, then feed $u = e$ as the initial condition into “Block B” (the DNA32 system), we’d have a polynomial PIVP for $\zeta(3)$ with rational initial conditions.

But in a single ODE system, both blocks run concurrently from $\tau = 0$. Block B starts before Block A has converged. The PIVP compiler problem is: can we compile a sequential composition of PIVPs into a single concurrent system with the correct limit?

Gate chains

One approach: delay Block B using a chain of polynomial gates. Each gate variable $g_k' = g_{k-1}^2(1 - g_k)$ waits for the previous one to saturate before it begins to rise. Block B’s evolution rate is multiplied by $g_k^2$, so it stays nearly frozen until the full chain opens — by which time Block A has approximately converged.

The resulting system is fully polynomial with rational initial conditions. But the limit carries a residual error from Block A’s transient phase. Numerically, the error decays as roughly $O(3^{-k})$ for $k$ gate levels: measurably better with more gates, but never zero for any fixed depth.

Integrals vs fixed points

The gate chain works perfectly when Block B computes via a fixed point — convergence to an attractor that “forgets” its initial conditions. But the $\zeta(3)$ computation is of integral type: Block B accumulates

along a trajectory, and this integral depends on the entire history of $u(\tau)$, not just its limit. During the early transient, $u(\tau) \neq e$, and Block B has already begun sweeping through part of its integration range. The accumulated error

is bounded but nonzero for any fixed gate depth.

This integral-vs-fixed-point distinction is the central structural issue: fixed-point computations are robust to transient parameter errors, while integral computations are not.

Finite-time blow-ups

An alternative compilation strategy uses explosions ($v' = v^2$, blow-up at $t = 1$) to create “infinity at a point.” After DNA32-style resolution, Block A has infinite effective time by $t = 1$, so $u$ reaches $e$ exactly at that instant. But $u(t) = e^t$ passes through $e$ at $t = 1$ without freezing — and freezing it would require a polynomial vanishing at $e$, impossible since $e$ is transcendental. Meanwhile, Block B’s integration variable sweeps its range during the explosion, reading $u(t)$ at varying times rather than the snapshot $u(1) = e$.

Both obstructions — unfreezing and the snapshot problem — stem from the same root cause: in a concurrent ODE, one block cannot take a “snapshot” of another block’s output at a specific instant.

What This Means

Here’s the updated scorecard:

Three polynomial systems, each failing in a different way: transcendental initial conditions, a fixed point at the origin, or residual transient error. The approaches are complementary but not yet combinable.

The PIVP compiler problem now reduces to a concrete question: can a Lyapunov function $V$ be constructed for the coupled estimator-plus-main system such that $V' \le 0$ is achieved by polynomial dynamics? If so, the compiled system would converge to the correct value regardless of transient estimator errors, unifying fixed-point and integral-type computations. Alternatively, reformulating $\zeta(3)$ as an attracting fixed point (rather than an accumulated integral) would sidestep the issue entirely.

The Derivation in Retrospect

Stepping back, the method itself is general and worth knowing:

1. Start with a series $F(x) = \sum a_n x^n$ whose value at some point gives the constant you want.
2. Find a recurrence for the coefficients $a_n$ — usually by computing the ratio $a_{n+1}/a_n$ and checking that it’s rational in $n$.
3. Multiply by $x^n$ and sum, using the Euler operator $\theta = xD$ to translate $n^k a_n$ into differential operators.
4. Simplify to get a linear ODE with polynomial coefficients.

This is a special case of the theory of holonomic functions: a power series is D-finite (satisfies a linear ODE with polynomial coefficients) if and only if its coefficient sequence is P-recursive (satisfies a linear recurrence with polynomial coefficients). The Euler operator $\theta$ is the bridge between the two worlds.

For our series, the ratio $a_{n+1}/a_n$ was rational of degree 3 in the denominator, which predicted (and delivered) a third-order ODE. For other number-theoretic series — Catalan numbers, Delannoy numbers, Apéry numbers for $\zeta(2)$ — the same machine produces their differential equations, each with its own singularity structure and its own story about whether the singularity can be resolved.

This is part of an ongoing project on bounded analog complexity — specifically, the question of whether $\zeta(3) \in R_\mathrm{RTCRN}$. The generating function approach reveals the exact obstruction (a double indicial root at a regular singular point), while the PIVP compiler attempts show that the deeper issue is the integral-vs-fixed-point distinction in how $\zeta(3)$ is computed. The open directions are: Lyapunov-based compilation, a fixed-point characterization of $\zeta(3)$, or a fundamentally different ODE.