Spigot algorithms: streaming π one digit at a time

In Python, you can write an infinite, lazy sequence of π’s digits in a handful of lines and consume it like any other iterator:

g = pi_digits()
next(g)                    # 3
list(islice(g, 20))        # 20 more digits, computed on demand

That’s a strange claim once you stop to look at it. π is irrational; every digit depends, in principle, on infinite series that don’t terminate. How can a program emit digit 50 without first deciding how many digits in total it’s going to compute?

The answer is a class of algorithms called spigot algorithms — they let π “drip” out one digit at a time, with the digit emitted as soon as it’s provably fixed, no matter how many more terms of the series remain unread.

What “spigot” means

A normal high-precision π routine picks a target precision N, computes π to N digits in fixed-point or fixed-precision floats using something like Machin’s formula or Chudnovsky, then prints the result. To get one more digit, you have to start over with a bigger N.

A spigot algorithm flips the dependency: instead of “given N, produce N digits,” it maintains some bounded state and asks at every step “is the next digit pinned down?” If yes, it emits the digit and updates state. If not, it absorbs another term of an underlying series and tries again.

There are two flavours:

• Bounded spigot (Rabinowitz & Wagon, 1995): you tell it up front how many digits you want; it allocates an array proportional to that count and grinds through it. Each output digit is produced once and never revised. Famous as a short, surprising program — it fits in a tweet — but it’s not truly streaming: the N is baked in.
• Unbounded spigot (Gibbons, 2006): a true lazy stream. The state size grows with the digits produced, but the algorithm never needs to know N in advance. You can pipe it into head, islice, or anything else and stop whenever you like.

Both rely on the same trick at heart: keeping the uncertainty interval for π narrow enough that the next digit’s identity is provable from the current state.

Gibbons’ unbounded spigot

The Gibbons algorithm threads a single iteration loop that’s either in produce mode (emit a digit) or consume mode (pull in another term of the underlying series). The state is six integers (q, r, t, k, n, l). The full Python is:

def pi_digits():
    q, r, t, k, n, l = 1, 0, 1, 1, 3, 3
    while True:
        if 4*q + r - t < n*t:
            yield n
            q, r, t, k, n, l = (
                10*q,
                10*(r - n*t),
                t,
                k,
                (10*(3*q + r)) // t - 10*n,
                l,
            )
        else:
            q, r, t, k, n, l = (
                q*k,
                (2*q + r)*l,
                t*l,
                k + 1,
                (q*(7*k + 2) + r*l) // (t*l),
                l + 2,
            )

The condition 4*q + r - t < n*t is the producibility check: it tests whether the current state has pinned the next digit to a single value n. If yes, emit n and shift the state by a factor of 10 (so the next digit becomes the next most-significant). If no, absorb another series term — q, r, t all grow, the bounds tighten, and we try again on the next iteration.

Try it

Press step to advance one iteration at a time. Each step is either a produce (the < branch fires, a digit is appended on the right) or a consume (the ≥ branch fires, no digit emitted but q, r, t get larger). Press play to let it run; the first ten digits come out within ~30 steps, but the gap between digits widens as the integers get larger.

Watch what happens to the bit-lengths of q, r, t. After 100 digits they’re a few hundred bits each; by 1000 digits they’re thousands. That’s the cost of being lazy — you’re paying in space what a bounded algorithm pays in commitment.

What the state actually represents

The six numbers aren’t arbitrary. The pair (q, r, t, k, n, l) encodes a Möbius transformation x ↦ (q·x + r) / (t·x + ?) composed with the partial series for π expanded so far. Each k indexes one term of a series for π:

π = Σ ( (2k+1) · k!² · 2^(k+1) ) / (2k+1)!     [for k ≥ 0]

(equivalent to Euler’s transformed series; the exact form isn’t crucial here). The consume step absorbs the k-th term into the accumulator. The produce step extracts one decimal digit from the front of the accumulator without disturbing what’s behind it.

The clever part is the producibility check: by computing the Möbius image of the unit interval [3, 4] (the loosest possible bound for the remainder of the series), the algorithm can ask “does the entire remaining series map to numbers that all agree on the next digit?” If yes — produce; the digit is locked in. If no — consume more.

This is why you never see a digit get revised. The algorithm only ever emits a digit once it’s provably the right one given the state, regardless of what the unread series terms turn out to be.

A few practical notes

• Integer arithmetic only. Everything is int — no floats, no rounding, no precision loss. Python’s arbitrary-precision integers are doing the heavy lifting; the same algorithm in C needs a bignum library (GMP).
• Cost per digit grows. Roughly: digit n requires state with O(n) digits and a multiplication on each step. So generating N digits is O(N²) bit-operations. Fine for the first thousand; painful past a million.
• Faster algorithms exist, but not as a stream. Chudnovsky (used by mpmath, y-cruncher, and every record-setting π computation) is much faster — but it’s a bounded algorithm: you pick N, allocate the precision, compute. There’s no Chudnovsky generator.
• The BBP formula is a different shape entirely. Bailey–Borwein–Plouffe lets you compute the n-th hexadecimal digit of π directly, without computing any earlier digit, in O(n log² n) time. That’s not a spigot — it’s random-access. Different superpower; same family of “weird things you can do with π’s series representations.”

Why this is worth knowing

The Gibbons spigot is a small but striking demonstration that laziness composes with mathematics. The Python generator protocol expects a function that can be paused and resumed, advancing a finite step at a time. π is an infinite object. The spigot is the bridge: it turns an infinite series into a finite state machine with a clear “I can emit one more digit now” criterion.

You’ll see the same pattern in interval arithmetic, continued fraction expansions, exact real arithmetic (e.g. Edalat–Potts), and Boehm’s “constructive reals” — all of them are different ways of arranging an infinite computation so that a finite consumer can ask for “the next bit” and get a guaranteed-correct answer.

Source

Jeremy Gibbons, Unbounded Spigot Algorithms for the Digits of Pi (American Mathematical Monthly, 2006) — the original paper, eight pages, very readable. It also covers a streaming spigot for e and explains the Möbius-transformation framework in full.

For the bounded version, Stanley Rabinowitz & Stan Wagon, A Spigot Algorithm for the Digits of Pi (AMM, 1995) — older, shorter, and the one most people mean when they say “the spigot algorithm” without qualifying it.