infsup

by Zinan Huang 🌸

Computing π Without an Inverter: Three Direct Designs

2026-04-25


The Chudnovsky design in the previous post computes $1/\pi$ first and then inverts. The inverter is a low-pass filter: its decay rate is fixed by the upstream modular constant and cannot be enlarged by spending more states. This post takes a different route. We build three polynomial PIVPs whose output state $P(\tau)$ targets $\pi$ directly, without any inversion step. Two of them (Leibniz, Machin) compute π exactly; the third (BBP) hits a structural obstruction at $z = 0$ and ends up as an obstruction example, not a valid construction. The trade for the two valid designs is convergence rate for honesty: each is slower per unit time than Chudnovsky-with-inverter, but the fixed-point readout is exactly $\pi$, not its reciprocal.

Why direct readout matters

A PIVP is “polynomial” because each right-hand side is a polynomial in the state variables. To compile it onto a chemical reaction network we need every coefficient and initial value to be rational. The output is read off as the limit value of one designated state.

When the natural series for a constant $L$ has shape

$$ \frac{1}{L} \;=\; \sum_{k\ge 0} a_k \, c^k, $$

a CRN can sum the $a_k c^k$ terms straightforwardly into a state $S$ that converges to $1/L$. To extract $L$ itself we add an inverter $\dot P = c_0 - L_0 P$ whose fixed point is $c_0/L_0 = L$. The inverter’s decay rate at the fixed point is exactly $-L_0$ — the value of the upstream modular constant. Spending extra states inside the inverter cannot speed it up; the previous post discusses a rate-$k$ variant that mitigates the cap but does not remove it.

The clean fix is to replace the series for $1/L$ with one whose partial sums converge to $L$ directly. Then no inversion is needed and the fixed-point readout is the constant itself.

Three direct-π series, three PIVPs

We pick three classical π series whose partial sums approach $\pi$ — not $1/\pi$ — and write each one as a polynomial PIVP with rational coefficients.

Design 1: Leibniz / arctan(1)

Leibniz’s $\pi = 4(1 - 1/3 + 1/5 - \cdots) = 4\arctan(1)$ is the slowest of the three, but the cleanest in continuous time. The generating function $4\arctan(z) = \sum_{k\ge 0} (-1)^k 4 z^{2k+1}/(2k+1)$ has derivative $4/(1+z^2)$, so we encode states $z$, $w = 1/(1+z^2)$, $u = 4\arctan(z)$, drive $z$ from 0 to 1 along $\dot z = 1 - z$, and compute $\dot w$ by the chain rule.

$$ \begin{cases} \dot z \;=\; 1 - z, \\ \dot w \;=\; -2\,z\,(1-z)\,w^{2}, \\ \dot u \;=\; 4\,(1-z)\,w. \end{cases} $$

Initial values

$$ \begin{cases} z(0) \;=\; 0, \\ w(0) \;=\; 1, \\ u(0) \;=\; 0. \end{cases} $$

State count: 3. Output: $P = u$. Convergence rate near the fixed point: linearizing gives $\dot z \approx -(z - 1)$, so $|1 - z(t)|$ decays like $e^{-t}$ and the readout error inherits this rate.

The integrand $1/(1+z^2)$ stays bounded and analytic on $[0,1]$, so nothing dramatic happens at the boundary; the $\arctan$ series' notoriously slow $\sim 1/n$ per-term convergence is replaced, in the continuous-time form, by clean exponential decay set by the drive.

Design 2: Machin / Stormer

Machin’s identity

$$ \pi \;=\; 16\arctan\!\bigl(\tfrac{1}{5}\bigr) - 4\arctan\!\bigl(\tfrac{1}{239}\bigr) $$

evaluates two arctangents at rational interior points and combines them. The PIVP runs two parallel arctan channels. For $j = 1, 2$ with target $z_j^* \in \{1/5,\ 1/239\}$:

$$ \begin{cases} \dot z_j \;=\; z_j^{*} - z_j, \\ \dot w_j \;=\; -2\,z_j\,(z_j^{*} - z_j)\,w_j^{2}, \\ \dot u_j \;=\; (z_j^{*} - z_j)\,w_j. \end{cases} $$

Initial values

$$ \begin{cases} z_j(0) \;=\; 0, \\ w_j(0) \;=\; 1, \\ u_j(0) \;=\; 0, \end{cases} $$

with output

$$ P \;=\; 16\,u_1 - 4\,u_2. $$

State count: 6. Each channel converges with rate 1 in $\tau$. The per-term series error in the discrete sum $\arctan(1/m) = \sum (-1)^k m^{-(2k+1)}/(2k+1)$ falls geometrically at ratio $m^{-2}$; that geometric structure is what makes Machin quantitatively serious historically (Shanks computed 707 digits with it). In the continuous PIVP form, this fast per-term decay disappears behind the rate-1 drive, just as in Leibniz; the relevant question is state count, not series rate.

Higher-order Stormer formulas (e.g. the four-term $48\arctan(1/49) + 128\arctan(1/57) - 20\arctan(1/239) + 48\arctan(1/110443)$) give 12-state PIVPs by the same construction. They do not improve the PIVP rate; they are useful only when each $\arctan$ is computed discretely from its series.

Design 3: BBP via the $G_j$ generating function (obstruction example)

Caveat up front. Unlike Leibniz and Machin, the BBP construction below does not compute π exactly. The cubic drive forces a non-zero starting point $z(0) = \varepsilon$ and a truncated rational IV, so the continuous-time limit is $\pi + \delta$ with $\delta \neq 0$. We include the design because it is a clean illustration of the recurring $z = 0$-fixed-point obstruction in BBP-style holonomic encodings, not as a valid CRN that computes π.

The Bailey–Borwein–Plouffe formula

$$ \pi \;=\; \sum_{k\ge 0} \frac{1}{16^{k}}\!\left(\frac{4}{8k+1} - \frac{2}{8k+4} - \frac{1}{8k+5} - \frac{1}{8k+6}\right) $$

discovered in 1995 is famous for letting one compute the $n$-th hexadecimal digit of $\pi$ without computing the previous digits. It also gives a direct-π construction in PIVP form, but with a small obstruction we have to walk around.

The natural integral form uses the substitution $\xi^8 = u$ with $u^{*} = 1/16$, giving evaluation point $\xi^{*} = 1/\sqrt 2$. That $\sqrt 2$ is irrational and has no place in a $\mathbb Q$-coefficient drive equation. The bypass is to keep the parameter as $z = u$ rather than $\xi$, evaluate at $z^* = 1/16$ (rational), and use the first-order ODE for the four series

$$ G_j(z) \;=\; \sum_{k\ge 0} \frac{z^{k}}{8k+j}, \qquad j \in \{1, 4, 5, 6\}. $$

Multiplying the recurrence $(8k+j)\,a_k = 1$ by $z^k$ and summing gives the rational-coefficient ODE

$$ 8\,z\,(1-z)\,G_j'(z) \;+\; j\,(1-z)\,G_j(z) \;=\; 1. $$

To make this a polynomial PIVP we choose a drive that kills the $1/(z(1-z))$ factor. The cubic drive $\dot z = z(1-z)(1/16 - z)$ does so cleanly, with all three roots $z = 0,\ 1/16,\ 1$ as fixed points. Integration starts at $z = \varepsilon$ small and rational and flows toward $z = 1/16$.

$$ \begin{cases} \dot z \;=\; z\,(1-z)\,(1/16 - z), \\ \dot G_{1} \;=\; \tfrac{1}{8}\,(1/16 - z)\,\bigl(1 - (1-z)\,G_{1}\bigr), \\ \dot G_{4} \;=\; \tfrac{1}{8}\,(1/16 - z)\,\bigl(1 - 4\,(1-z)\,G_{4}\bigr), \\ \dot G_{5} \;=\; \tfrac{1}{8}\,(1/16 - z)\,\bigl(1 - 5\,(1-z)\,G_{5}\bigr), \\ \dot G_{6} \;=\; \tfrac{1}{8}\,(1/16 - z)\,\bigl(1 - 6\,(1-z)\,G_{6}\bigr). \end{cases} $$

Initial values are taken at a small rational $\varepsilon$ with each $G_j(0)$ given by the finite, rational truncated power series

$$ G_{j}(0) \;=\; \sum_{k=0}^{K-1} \frac{\varepsilon^{k}}{8k+j}, \qquad j \in \{1, 4, 5, 6\}. $$

Concretely we use $\varepsilon = 1/10000$ and $K = 12$, so each $G_j(0)$ is a sum of 12 rational terms. The first few read

$$ \begin{cases} G_{1}(0) \;=\; 1 \;+\; \frac{1}{9 \cdot 10^{4}} \;+\; \frac{1}{17 \cdot 10^{8}} \;+\; \frac{1}{25 \cdot 10^{12}} \;+\; \cdots \;+\; \frac{1}{89 \cdot 10^{44}}, \\[4pt] G_{4}(0) \;=\; \frac{1}{4} \;+\; \frac{1}{12 \cdot 10^{4}} \;+\; \frac{1}{20 \cdot 10^{8}} \;+\; \cdots \;+\; \frac{1}{92 \cdot 10^{44}}, \\[4pt] G_{5}(0) \;=\; \frac{1}{5} \;+\; \frac{1}{13 \cdot 10^{4}} \;+\; \frac{1}{21 \cdot 10^{8}} \;+\; \cdots \;+\; \frac{1}{93 \cdot 10^{44}}, \\[4pt] G_{6}(0) \;=\; \frac{1}{6} \;+\; \frac{1}{14 \cdot 10^{4}} \;+\; \frac{1}{22 \cdot 10^{8}} \;+\; \cdots \;+\; \frac{1}{94 \cdot 10^{44}}. \end{cases} $$

Every term is a rational number with denominator dividing $94 \cdot 10^{44}$, so each $G_j(0) \in \mathbb Q$ and so is $z(0) = 1/10000$. The truncation tail $\bigl| G_j(\varepsilon) - G_j(0) \bigr| \le \varepsilon^{K}/(1-\varepsilon) < 2 \cdot 10^{-48}$ is below the working precision floor and contributes no observable error.

Output: $P = 4\,G_{1} - 2\,G_{4} - G_{5} - G_{6}$. State count: 5.

Linearizing the cubic drive at $z = 1/16$ gives rate

$$ \partial_z\bigl[z(1-z)(\tfrac{1}{16}-z)\bigr]\Big|_{z=1/16} \;=\; -\tfrac{1}{16}\!\cdot\!\tfrac{15}{16} \;=\; -\tfrac{15}{256} \;\approx\; -0.0586. $$

So $|1/16 - z(\tau)|$ decays at rate $\sim 0.06$, ten times slower than the rate-1 drives of Leibniz and Machin. The cost of avoiding $\sqrt 2$ is this slow rate — a polynomial drive whose roots include both endpoints is necessarily slower at the interior root than a linear drive at the same target. A linear drive $\dot z = (1/16 - z)$ would restore rate 1 but reintroduces the $1/(z(1-z))$ factor in $\dot G_j$, requiring an auxiliary state $W = 1/(8z(1-z))$ and sensitive to numerical drift in $W$. We stay with the cubic drive for the cleaner construction.

Numerical comparison

A 4th-order Runge–Kutta integration from each design’s initial condition gives Figure 1.

Direct-π PIVPs: Leibniz, Machin, BBP

Panel (a) shows $P(\tau)$ rising to $\pi$ for each design and stopping at $\pi$ on the dashed line. Panel (b) shows $\log_{10}|P(\tau) - \pi|$ versus $\tau$ on a common axis. Leibniz and Machin run at rate 1 in $\tau$ and bottom out near machine precision around $\tau \approx 35$. BBP runs at rate $0.06$ and reaches the same precision around $\tau \approx 600$, ten times longer in continuous time. All three hit floating-point limits with no plateau, no drift, no inverter.

The asymptotic slopes match the linearization predictions: Leibniz and Machin slope $\approx -0.434$ (which is $-1/\ln 10$), BBP slope $\approx -0.025$ (which is $-15/(256\ln 10)$). The drive’s eigenvalue at the fixed point is the only thing that sets the rate.

A theoretical distinction we have to be honest about

The three designs above are not on equal footing once we separate the GPAC’s idealized continuous-time limit from any numerical approximation. The cleanest way to see this is to ask: does the continuous-time evolution of the GPAC, with rational coefficients and rational initial values, give $P(\infty) = \pi$ exactly?

For Leibniz the drive is $\dot z = 1 - z$. At $z = 0$, $\dot z = 1 \neq 0$, so $z = 0$ is not a fixed point of the drive. We start at $z(0) = 0$ exactly, with $w(0) = 1, u(0) = 0$ — every initial value is rational and exact. The trajectory follows $z(\tau) = 1 - e^{-\tau}$ exactly, $u(\tau) = 4\arctan(z(\tau))$ exactly, and $u(\infty) = 4\arctan(1) = \pi$ exactly. Same story for Machin: $\dot z_j = z_j^* - z_j$ at $z_j = 0$ gives $\dot z_j = z_j^* \neq 0$, so we start exactly at $z_j(0) = 0$ with exact rational IV, and the limit is exactly $16\arctan(1/5) - 4\arctan(1/239) = \pi$.

The BBP cubic-drive design is different. The drive $\dot z = z(1-z)(1/16 - z)$ has $z = 0$ as a fixed point: $\dot z|_{z=0} = 0$. If we set $z(0) = 0$, the trajectory stays at $z = 0$ forever and we never reach the readout point $z^* = 1/16$. So we are forced to start at $z(0) = \varepsilon > 0$. But that means we have to give an initial value for $G_j$ at $z = \varepsilon$, and the exact value $G_j(\varepsilon) = \sum_{k \ge 0} \varepsilon^k/(8k+j)$ is an infinite sum — irrational for almost any rational $\varepsilon$. The rational IV has to be a truncation of this sum, and the truncation introduces a finite error in the limit.

$$ G_j'(z) + \frac{j}{8z}\,G_j(z) \;=\; \frac{1}{8z(1-z)}, $$$$ \bigl|G_j^{\mathrm{trunc}}(1/16) - G_j(1/16)\bigr| \;\le\; |C| \cdot 16^{j/8} \;=\; \varepsilon^{K} \cdot (16\varepsilon)^{j/8}. $$

For $\varepsilon = 10^{-4}, K = 12$ this is at most $2 \times 10^{-48}$ in the BBP readout $P = 4G_1 - 2G_4 - G_5 - G_6$. The factor $(16\varepsilon)^{j/8} < 1$ comes from $z^{-j/8}$ being a decreasing function of $z$ on $(0, 1]$ for $j > 0$ — perturbations introduced at $z = \varepsilon$ are damped (not amplified) as $z$ flows toward $z^* = 1/16$. But they are not driven to zero: the limit is $\pi + \delta$ with $\delta \neq 0$ for any finite $\varepsilon, K$.

The criterion is binary: limit equals π exactly, or the construction does not compute π. A family parametrized by $(\varepsilon, K)$ whose limit is $\pi + \delta_{\varepsilon, K}$ is not a GPAC that computes π, even if $\delta \to 0$ as $\varepsilon \to 0, K \to \infty$. By that criterion:

The decisive feature is whether the drive has a fixed point at $z = 0$. If not, $z(0) = 0$ with exact rational IV gives an exact-π limit. If so, we are forced into truncation territory and the construction no longer hits π exactly. Most BBP-like generating functions inherit a regular singular point of their generating-function ODE at $z = 0$, which forces the drive to vanish there — so this obstruction is not specific to our particular cubic-drive choice; it is a structural feature of BBP-type holonomic-series encodings.

Open question. Is there a polynomial PIVP that hits π exactly at BBP-style per-term rate? Concretely: a drive without a fixed point at the IV starting point, with a generating-function ODE regular at the starting point, rational indicial value, and continuous-time limit exactly π at the readout point. Leibniz and Machin escape the obstruction with their linear drives but at the cost of slow per-term rate; whether the BBP-rate analogue exists is open.

A caveat on convergence rate is also in order. The “rate 1 in $\tau$” figures for Leibniz and Machin look misleadingly competitive because the continuous-time parametrization $z(\tau) = z^* (1 - e^{-\tau})$ erases the per-term rate of the underlying series. What survives in physical complexity is the per-term rate, because that controls how many states (or how much working precision) one needs for $n$-digit accuracy:

So Leibniz, Machin, and BBP are all in the same “$\sim 1$ digit per term” regime — the geometric prefactors differ but the order of magnitude is the same. Chudnovsky’s modular acceleration is in a different league, and it is the only one of the four that lives there. Removing the inverter is one axis (the topic of this post); reaching modular-arithmetic per-term rates is a separate axis (the topic of the open question below).

What this rules out, and what it doesn’t

Two clean direct-π PIVPs now exist:

  1. Direct π readout in 3 states (Leibniz). Slow but irreducibly small. Limit is exactly π.
  2. Direct π readout at rational interior points (Machin), 6 states. Limit is exactly π.

The BBP design we wrote down sits in a different category: it is an obstruction example, not a valid construction. Its cubic drive vanishes at $z = 0$, forcing a non-zero starting point and a truncated rational IV; the continuous-time limit is therefore $\pi + \delta$ with $\delta \neq 0$, not π. The structural barrier — drives that vanish at $z = 0$ — recurs across BBP-type holonomic encodings, and finding a BBP-rate exact-π PIVP is one of the open questions below.

What is not yet possible is direct π readout at Chudnovsky rate. Chudnovsky’s series outputs $1/\pi$, by construction; there is no known holonomic decomposition

$$ A_n - \pi B_n \;=\; \int\!\!\int\!\!\int (\cdots)^{n} $$

where $A_n, B_n \in \mathbb Q$ are individually holonomic and the right-hand side gives a Beukers-style integral driving the residual to zero at Chudnovsky rate. For Apéry’s $\zeta(3)$ such a triple integral exists (Beukers 1979) and turns the $\zeta(3)$ readout into a no-inverter two-channel design. For $\pi$ the analogous identity is open — see Adamczewski–Rivoal and the Zudilin bibliography for partial progress.

This is the cleanest stated open problem at the intersection of analog complexity, irrationality theory, and modular forms: find rationals $A_n, B_n$ such that $A_n - \pi B_n$ has an explicit ODE-encodable integral form converging at modular-arithmetic rate. A positive answer would let us build a Chudnovsky-rate inverter-free π-PIVP. A negative answer would explain why π is harder than $\zeta(3)$ at the level of continuous-time computation.


Code: code/direct_pi_candidates.py in the HolonomicCRN paper repository generates the figure.


问渠那得清如许?为有源头活水来。

— 朱熹《观书有感》

How does the canal stay so clear? Because fresh water keeps arriving from the source.