An arcsine ghost in the binomial channel: numerically chasing a √(n·log log n) lower bound

Paper. An Improved Lower Bound on Support Size of Capacity-Achieving Inputs for the Binomial Channel: Extended version, M. Baniasadi, L. Barletta, A. Dytso, arXiv:2605.12472 (cs.IT).

The binomial channel is the toy that refuses to be tame: input a probability $x\in[0,1]$, get back $Y \sim \mathrm{Bin}(n,x)$. The capacity-achieving input is a finite mixture of point masses, but how many of them is a question that has been re-opened every few years for two decades. Baniasadi, Barletta and Dytso just nudged the lower bound from $\sqrt n$ to $\sqrt{n\log\log n}$. This post implements the three numerical pillars that hold their proof up — Blahut–Arimoto for the capacity asymptote, the Beta(1/2,1/2) reference output, and a $\chi^2$-rigidity sweep — in roughly three hundred lines of Python that you can paste into a notebook.

Why Python

I usually reach for Julia or Rust on numerics-heavy posts, but this paper lives almost entirely inside scipy.special. Every quantity I need — gammaln, betaln, the Beta-binomial PMF, the digamma function — is in there, log-stable, and SIMD-vectorised already. Blahut–Arimoto is one matrix multiply per iteration; the bottleneck is iteration count, not language. Going to Julia would buy a 2× constant and cost me the day; going to Rust would mean wiring the special functions by hand. The interesting code is what I evaluate, not how fast. Python it is.

Restating the result in my own words

The binomial channel takes a real number $x \in [0,1]$ and returns the count of heads in $n$ independent flips of a coin biased to land heads with probability $x$. The output alphabet is $\{0,\dots,n\}$, finite. The input alphabet is uncountable, but a classical result of Witsenhausen (1980) and Burshtein–Della Pietra–Kanevsky (1992) says the capacity-achieving input distribution $P_{X^*}$ is discrete and supported on finitely many points. Call that count $K(n)$. Two questions have been open ever since:

1. How fast does $K(n)$ grow with $n$?
2. What does the capacity look like as $n \to \infty$?

For (1) the standard upper bound is $K(n) \le n/2 + 1$ (Witsenhausen). Lower bounds were stuck at $K(n) \gtrsim \sqrt n$ for a long time. The new paper improves this to

with an explicit constant $c$. That extra $\sqrt{\log\log n}$ factor sounds tiny — and analytically it is — but it required a genuinely new tool: a quantitative $\chi^2$-rigidity bound for how well a $K$-component binomial mixture can imitate a particular continuous mixture.

For (2) the paper sharpens existing capacity expansions to

(in nats; bits is the same formula divided by $\ln 2$). The $o(1)$ here is what makes the proof go: the gap in capacity, not the value, controls how close a candidate input distribution can be to optimal.

The proof has three legs. None of them stands alone:

• Leg A. Tighten the capacity to $C(n) = \tfrac12\log(n\pi/(2e)) + o(1)$.
• Leg B. Show that the input $X_r \sim \mathrm{Beta}(1/2, 1/2)$ — the arcsine prior, also known as Jeffreys' prior for the Bernoulli mean — induces a Beta-binomial output distribution that is asymptotically optimal: it converges to the capacity-achieving output $P_{Y^*}$ in relative entropy, and (after a comparison step) in $\chi^2$ divergence.
• Leg C. Prove a $\chi^2$ lower bound: any $K$-atom input induces an output whose $\chi^2$ distance from the Beta-binomial reference cannot be smaller than something growing in $n$, unless $K \gtrsim \sqrt{n\log\log n}$.

Combine: the optimal input must reach a Beta-binomial-like output, but $\chi^2$-rigidity forbids a small number of atoms from doing so.

The implementation in this post hits all three legs numerically. I cannot reproduce the analytic $o(1)$ constants from a laptop, but I can show:

• the capacity gap shrinks at the rate Leg A predicts,
• the Beta-binomial KL/χ² distance to the BA-optimal output shrinks (Leg B), and
• the smallest $K$ that can imitate the Beta-binomial in $\chi^2$ tracks $\sqrt{n\log\log n}$ almost perfectly (Leg C).

The channel and its log-PMF

A row of the channel transition matrix is a binomial:

For numerical work the only sane thing is to live in log-space: a tail $W(y\mid x)$ for large $n$ and small $x$ vanishes faster than float64 can hold. Here's the entire channel, in five lines after the bookkeeping:

from scipy.special import gammaln, betaln, xlogy
import numpy as np
LOG2 = np.log(2.0)

def log_binom_pmf(n: int, x: np.ndarray) -> np.ndarray:
    """log P(Y=k | X=x) on the full output alphabet, for a vector of x."""
    x = np.clip(np.asarray(x, dtype=np.float64), 1e-12, 1 - 1e-12)
    k = np.arange(n + 1)
    log_binom = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1)
    return log_binom[None, :] \
         + k[None, :] * np.log(x)[:, None] \
         + (n - k)[None, :] * np.log1p(-x)[:, None]

gammaln keeps the binomial coefficient stable up to $n \sim 10^6$. The clip at 1e-12 is what saved this script: at machine epsilon, 1.0 - 1e-300 == 1.0 exactly, and log1p(-1.0) is $-\infty$. I lost an iteration of the script on that one and put a comment in the source to remind me.

The Beta-binomial PMF — the marginal of $Y$ when $X \sim \mathrm{Beta}(\alpha,\beta)$ — is just as compact:

def betabin_log_pmf(n: int, a: float, b: float) -> np.ndarray:
    k = np.arange(n + 1)
    log_binom = gammaln(n + 1) - gammaln(k + 1) - gammaln(n - k + 1)
    return log_binom + betaln(k + a, n - k + b) - betaln(a, b)

For $\alpha=\beta=1/2$ this is the arcsine Beta-binomial — a discrete U-shape with most of the mass at $y=0$ and $y=n$. It is the paper's reference output throughout.

Blahut–Arimoto: a finite-grid stand-in for the optimum

The capacity is

Blahut–Arimoto (1972) is the textbook fixed-point iteration for the discrete-input case. To use it for a continuous input, we discretise the input alphabet: pick a fine uniform grid on $[0,1]$, run BA, and let the iteration concentrate the mass on whatever subset of the grid actually wants to carry probability. The capacity-achieving distribution being supported on finitely many points means BA on a fine grid will almost sparsify itself; near each true atom it will smear mass across a few adjacent cells. We post-process by clustering.

def blahut_arimoto(n: int, grid_size: int = 1001,
                   max_iter: int = 4000, tol: float = 1e-12):
    x = np.linspace(0.0, 1.0, grid_size)
    log_W = log_binom_pmf(n, x)         # (M, n+1)
    W = np.exp(log_W)                   # rows sum to 1
    M = grid_size
    p = np.full(M, 1.0 / M)             # initial input: uniform on grid

    I_prev = -np.inf
    for it in range(1, max_iter + 1):
        Q = p @ W                       # induced output distribution
        Q = np.maximum(Q, 1e-300)
        log_Q = np.log(Q)
        # d(x) = sum_y W(y|x) [log W(y|x) - log Q(y)]
        d = (W * (log_W - log_Q[None, :])).sum(axis=1)
        # multiplicative update on the simplex of input PMFs
        log_p = np.log(np.maximum(p, 1e-300)) + d
        log_p -= log_p.max()
        p = np.exp(log_p)
        p /= p.sum()
        I = (p * d).sum()               # mutual information in nats
        if abs(I - I_prev) < tol and it > 50:
            break
        I_prev = I

    Q = p @ W
    return {"x_grid": x, "p_x": p, "p_y": Q / Q.sum(),
            "capacity": I / LOG2, "iters": it}

The update is the classical one:

d is the relative entropy from the row $W(\cdot\mid x)$ to the current induced output $P_Y^{(t)}$. Inputs whose conditional sticks out from the current output get amplified; inputs whose conditional already looks like the average get damped. Convergence to the true capacity is monotone and linear, which is great theoretically and slow practically — for $n=256$ on a 1201-point grid, BA needs a few thousand iterations to settle.

A clustering helper for support extraction:

def extract_support(x_grid, p_x, rel_thresh=1e-3, gap=2):
    thresh = rel_thresh * p_x.max()
    active = np.where(p_x > thresh)[0]
    if active.size == 0: return []
    clusters, cur = [], [active[0]]
    for j in active[1:]:
        if j - cur[-1] <= gap:
            cur.append(j)
        else:
            clusters.append(cur); cur = [j]
    clusters.append(cur)
    out = []
    for cl in clusters:
        idx = np.array(cl); w = p_x[idx]; m = w.sum()
        out.append((float((w * x_grid[idx]).sum() / m), float(m)))
    return out

What BA reports

I ran the sweep $n \in \{4,8,16,32,64,128,256\}$ with grid sizes scaling from 401 to 1201 and iteration budgets from 1500 to 5000. Total wall time: about ten seconds. The headline numbers:

The gap halves roughly every doubling of $n$ — i.e. it goes like $n^{-1/2}$ in this regime, consistent with an $o(1)$ correction that includes the next-order Stirling-type term. By $n=256$ the BA capacity is within 0.11 bits of the asymptotic prediction, and by $n=1024$ a back-of-envelope extrapolation puts the gap below 0.05 bits.

The recovered support locations for small $n$ are the part I find the most pleasing:

n = 4,  K = 3
   x = 0.000000   p = 0.369969
   x = 0.500000   p = 0.256408
   x = 1.000000   p = 0.369969

n = 8,  K = 4
   x = 0.000000   p = 0.299576
   x = 0.281666   p = 0.197885
   x = 0.718334   p = 0.197885
   x = 1.000000   p = 0.299576

n = 16, K = 6
   x = 0.000000   p = 0.234222
   x = 0.144268   p = 0.146068
   x = 0.379891   p = 0.113481
   x = 0.620109   p = 0.113481
   x = 0.855732   p = 0.146068
   x = 1.000000   p = 0.234222

n = 32, K = 7
   x = 0.000000   p = 0.178910
   x = 0.074544   p = 0.112548
   x = 0.213261   p = 0.094258
   x = 0.500000   p = 0.208233
   x = 0.786739   p = 0.094258
   x = 0.925456   p = 0.112548
   x = 1.000000   p = 0.178910

Every one of these is symmetric about $x=1/2$, every one puts heavy mass at the endpoints, and the spacing of the interior atoms looks suspiciously like Chebyshev nodes — the same nodes one would get from the arcsine measure. That visual coincidence is the entire intuition behind the paper's choice of $\mathrm{Beta}(1/2,1/2)$ as the reference input. The BA optimum is trying to be arcsine-distributed; the discrete atoms are samples of that continuous reference, density-equalised.

A caveat I will not paper over

For $n \ge 128$, BA on a finite grid stops being trustworthy as a support counter. With grid spacing $\Delta x = 1/(M-1)$, the binomial conditionals $W(\cdot\mid x)$ and $W(\cdot\mid x+\Delta x)$ become numerically near-collinear when $x(1-x)$ is small — the channel manifold gets ill-conditioned. BA then leaves a thin smear of mass across the entire grid instead of fully concentrating to the true atoms, and any clustering threshold becomes arbitrary. Inspecting the top 25 grid masses for $n=128$ confirms this:

n=128: top 25 grid masses cover 37.69% of total
   x=0.00000  p=9.88e-02
   x=0.01800  p=1.81e-02
   x=0.01900  p=2.64e-02     <-- one true atom near 0.019
   x=0.02000  p=1.20e-02
   x=0.05100  p=3.33e-03
    .   .   .                 <-- a smeared cluster around 0.054
   x=0.05800  p=3.18e-03
   x=0.94200  p=3.18e-03      <-- mirror cluster around 0.946
    .   .   .
   x=0.98100  p=2.64e-02
   x=1.00000  p=9.88e-02

So at $n=128$ I can confidently say there are atoms at 0, ≈0.019, ≈0.054, and (by symmetry) their mirror images, but mass somewhere between 0.06 and 0.94 is small and noisy. Counting it accurately is a separate research project — the ground-truth $K(n)$ for moderate $n$ is computed in the literature with cutting-plane methods that I am not implementing here. What I can trust from BA at any $n$ is the capacity value, since that is a smooth functional of the input distribution and not sensitive to where the small-mass atoms sit.

Leg B: the Beta-binomial as a reference output

The paper picks $X_r \sim \mathrm{Beta}(1/2,1/2)$ and looks at the induced output $P_{Y_r}$, which is the Beta-binomial:

Its claim — Leg B in our breakdown — is that this $P_{Y_r}$ is asymptotically optimal: the relative entropy $D(P_{Y_r}\,\|\,P_{Y^*})$ vanishes, and after a comparison argument so does the $\chi^2$ divergence. We can check both directly:

Both columns shrink monotonically, with the KL halving roughly every two doublings of $n$ (so $D \sim n^{-1/4}$ in the visible range, slower than the $\chi^2$). This is the empirical face of Leg B. The Beta-binomial is not the capacity-achieving output for any finite $n$ — it cannot be, since $P_{Y^*}$ is induced by a discrete input — but it is a smooth continuous-input approximant that lives close enough to converge.

Leg C: $\chi^2$-rigidity and the heart of the new bound

Leg C is the genuine novelty. Suppose somebody hands you a $K$-atom input and asks: can the induced binomial mixture mimic the Beta-binomial $P_{Y_r}$? The paper answers no, quantitatively, and the rigidity strengthens with $n$.

Numerically I phrased it as an optimisation:

where $\chi^2(P\,\|\,Q) = \sum_y (P_y - Q_y)^2 / Q_y$. I solve it with mirror descent on the simplex of weights, with the locations initialised on Chebyshev / uniform nodes plus jitter, multiple restarts. Best-of-restarts is taken as a proxy for the global minimum.

def best_kpoint(n, K, n_starts=8, max_iter=500, lr=0.4, seed=0):
    rng = np.random.default_rng(seed + K * 17 + n)
    target = betabin_pmf(n, 0.5, 0.5)
    target_safe = np.maximum(target, 1e-300)
    best = np.inf
    cheb = 0.5 * (1 - np.cos(np.pi * (np.arange(K) + 0.5) / K))
    unif = (np.arange(K) + 0.5) / K
    for x0 in (cheb, unif):
        for trial in range(n_starts):
            x = np.clip(x0 + 1e-3 * rng.standard_normal(K), 1e-7, 1 - 1e-7)
            x.sort()
            log_W = log_binom_pmf(n, x); W = np.exp(log_W)
            w = np.full(K, 1.0 / K)
            for it in range(max_iter):
                Q = np.maximum(w @ W, 1e-300)
                grad = -2.0 * (W @ ((target - Q) / target_safe))
                w = w * np.exp(-lr * (grad - grad.min())); w /= w.sum()
            chi = chi2_div(target, w @ W)
            if chi < best: best = chi
    return best

Mirror descent on the simplex, also called the exponentiated-gradient update, is the right choice here because it preserves nonnegativity and normalisation for free. A subtraction of grad.min() before exponentiation keeps the iterate in float range. I keep locations fixed within each restart (so each inner solve is convex in the weights) and rely on multiple restarts to explore location-space.

The output is striking. Each row is one $n$; each cell shows the best $\chi^2$ achievable with $K$ atoms:

A few things jump out:

• The $\chi^2$ explodes when $K$ is too small. At $n=512$, two atoms get you a $\chi^2$ of $4 \cdot 10^{63}$. That is the rigidity: there is no way to mimic the Beta-binomial output with a small mixture, no matter where you put the atoms.
• The "elbow" — the smallest $K$ that pushes $\chi^2$ below a fixed target like $10^{-2}$ — moves with $n$ in a controlled way. Reading off:

The ratio in the rightmost column hovers between 1.0 and 1.3 across two orders of magnitude in $n$. Compare that to the ratio $K^*/\sqrt n$, which would grow without bound — by $n=512$, $K^*/\sqrt n \approx 1.4$ and rising. The data is genuinely consistent with the new $\sqrt{n\log\log n}$ rate and not with the previous $\sqrt n$ rate.

This is what the paper's third leg is, made finite. The $\chi^2$-rigidity argument they prove analytically says: for any $K$-atom input, the induced output is forced to differ from the Beta-binomial by a $\chi^2$ amount that scales like (modulo constants) $\exp(-cK^2/n)$ at best. For that to be small, you need $K \gg \sqrt n$; for it to be small enough to compete with the rate at which the Beta-binomial converges to $P_{Y^*}$, you need $K \gg \sqrt{n\log\log n}$. The double-log factor comes from chaining the $\chi^2$ rate to the $L^\infty$ approximation rate of the Beta-binomial output.

A worked example: $n=16$, the easy case

To show the pieces snapping together, let me chase $n=16$ end-to-end.

Capacity-achieving input from BA, six atoms:

x = 0.0000   p = 0.2342
x = 0.1443   p = 0.1461
x = 0.3799   p = 0.1135
x = 0.6201   p = 0.1135
x = 0.8557   p = 0.1461
x = 1.0000   p = 0.2342

The induced output $P_{Y^*}$ (length 17), computed as $\sum_i p_i\,\binom{16}{y}x_i^y(1-x_i)^{16-y}$:

y=0 :  P_Y* = 0.23448      P_BB = 0.20264
y=1 :  P_Y* = 0.04994      P_BB = 0.05776
y=2 :  P_Y* = 0.04144      P_BB = 0.04153
y=3 :  P_Y* = 0.03823      P_BB = 0.03551
y=4 :  P_Y* = 0.03457      P_BB = 0.03247
y=5 :  P_Y* = 0.03241      P_BB = 0.03077
y=6 :  P_Y* = 0.03113      P_BB = 0.02985
y=7 :  P_Y* = 0.03052      P_BB = 0.02948
y=8 :  P_Y* = 0.03037      P_BB = 0.02938
... (symmetric in y -> 16-y)
y=16:  P_Y* = 0.23448      P_BB = 0.20264

(I generated these with the script in binch.py above.)

The gap is largest at the endpoints, where the Beta-binomial under-estimates the true capacity-achieving mass — the discrete optimum piles a bit more probability on the endpoints than the smooth Beta-binomial does. Their KL distance is

a few percent of the capacity itself. By $n=256$ both have shrunk by roughly $4\times$, and the trend is clear.

Now ask: can a $K=4$ input mimic the Beta-binomial well at $n=16$? From the rigidity table, the best $\chi^2$ is $0.0365$ — about a sixth of the BA-induced gap. So four atoms is not yet enough to look indistinguishable from Beta-binomial; we need $K \ge 6$ to get below $10^{-3}$. At $n=16$ the predicted $\sqrt{n\log\log n}\approx 4$, so we need a constant-factor more atoms than the rate would suggest, which is consistent with the explicit constant in the bound being modestly larger than 1.

Putting the run on a laptop

A single bash line reproduces the headline tables:

PYTHONPATH=/tmp/pylib python3 run_capacity.py     # ~10 seconds
PYTHONPATH=/tmp/pylib python3 run_kpoint.py       # ~25 seconds

run_capacity.py is the BA sweep we just walked through:

import sys, time, math, json
import numpy as np
sys.path.insert(0, "/labs-output")
from binch import (blahut_arimoto, extract_support, betabin_pmf,
                   kl_bits, chi2_div, C_asymptotic_bits)

CONFIG = [(4, 401, 1500), (8, 401, 2000), (16, 401, 2500),
          (32, 601, 3000), (64, 801, 3500),
          (128, 1001, 4000), (256, 1201, 5000)]

for n, grid, mit in CONFIG:
    t0 = time.time()
    res = blahut_arimoto(n, grid_size=grid, max_iter=mit, tol=1e-12)
    bb = betabin_pmf(n, 0.5, 0.5)
    K = len(extract_support(res["x_grid"], res["p_x"],
                            rel_thresh=2e-3, gap=max(2, grid // 200)))
    print(f"n={n:>4d}  C_BA={res['capacity']:.6f}  "
          f"C_asym={C_asymptotic_bits(n):.6f}  "
          f"K={K:>3d}  KL_BB={kl_bits(bb, res['p_y']):.3e}  "
          f"chi2_BB={chi2_div(bb, res['p_y']):.3e}  "
          f"({time.time()-t0:.1f}s, {res['iters']} it)")

run_kpoint.py runs the rigidity sweep with the best_kpoint solver from earlier. Both scripts and binch.py are reproduced in this post in full; nothing of mine lives only in a tarball.

The KL/χ² helpers, for completeness:

def kl_bits(p, q):
    p = np.asarray(p); q = np.asarray(q); m = p > 0
    return float(xlogy(p[m], p[m] / q[m]).sum() / LOG2)

def chi2_div(p, q):
    p = np.asarray(p); q = np.asarray(q); m = q > 0
    return float(((p[m] - q[m]) ** 2 / q[m]).sum())

def C_asymptotic_bits(n):
    return 0.5 * np.log2(n * np.pi / (2 * np.e))

xlogy(p, p/q) is scipy.special.xlogy, which evaluates $p\log(p/q)$ with the right limit at $p=0$ — small thing, but it lets the KL be written as a one-liner.

What makes this paper different from prior work

The previous-best lower bound of $\sqrt n$ comes from a counting argument that bounds how many distinct PMFs over $\{0,\dots,n\}$ can be sufficiently far apart in some divergence — a packing-style argument with constants that are pleasant but inflexible. To improve on it, the new paper does two things in tandem:

1. It picks a specific reference output, the Beta-binomial induced by the arcsine prior, instead of arguing about all possible outputs. The arcsine is the unique input under which the Fisher information of the binomial channel is constant in $x$ — a fact that makes the chi-square rigidity argument tractable. (You can see this fingerprint in the Chebyshev-node spacing of the BA optima above.)
2. It plays the rigidity bound off against the capacity gap, not against a generic divergence. Because the asymptotic capacity is now pinned down to $\tfrac12\log(n\pi/(2e)) + o(1)$ (Leg A is also new in the paper, sharpening earlier constants), any candidate input must have $I(X;Y)$ exceed this benchmark up to vanishing error. That converts a "how close to optimal in capacity" question into a "how close to Beta-binomial in $\chi^2$" question, and the rigidity bound takes over from there.

The double-log factor falls out of the third step. To get the capacity gap closed, the candidate output must approximate the Beta-binomial in $L^\infty$ (or $\chi^2$) at a rate that scales with $\log\log n$. Combine that with the $K^2/n$-style $\chi^2$ rigidity, and the lower bound on $K$ picks up the $\sqrt{\log\log n}$ factor.

That's why the rate jumped: the proof gained access to a reference distribution whose specific structure could be exploited. In the rigidity table above, you can see the consequence directly — the elbow $K^*$ tracks $\sqrt{n\log\log n}$ much more cleanly than it tracks $\sqrt n$.

Limits of what I implemented

To be honest about scope:

• No proof of the constants. The paper gives explicit constants in front of $\sqrt{n\log\log n}$. I did not extract or numerically verify those — the empirical ratio $K^*/\sqrt{n\log\log n} \approx 1.0$–$1.3$ is suggestive but not a check on the analytic constant.
• BA loses support resolution at $n \ge 128$. Ground-truth $K(n)$ for the binomial channel above $n\approx 100$ is computed in the literature with cutting-plane / SDP-based methods that I am not implementing here. My BA support counts are honest to the grid; the capacity values they produce are honest.
• The K-point optimisation is local. Mirror descent on weights with multiple restarts is not a guarantee of global minimum. A determined adversary could possibly find a slightly smaller $\chi^2$ for a given $K$ than my best-of-restarts. The qualitative scaling — orders of magnitude of $\chi^2$ separating $K = c_1\sqrt n$ from $K = c_2\sqrt{n\log\log n}$ — is robust to that.
• Leg A (the capacity asymptote) is empirical only. I check $C_{\mathrm{BA}}(n) - C_{\mathrm{asym}}(n) \to 0$; I do not derive the next-order term.

If you want the analytic story end-to-end, read the paper. If you want a way to see that the story is plausible and to play with the channel yourself, this post is the playground.

Files used

uv pip install --target /tmp/pylib --quiet numpy scipy mpmath is the only setup. No GPU, no cluster, no proprietary data — the entire run finishes in well under a minute on one CPU core.

References

• M. Baniasadi, L. Barletta, A. Dytso. An Improved Lower Bound on Support Size of Capacity-Achieving Inputs for the Binomial Channel: Extended version. arXiv:2605.12472 [cs.IT], 2026. https://arxiv.org/abs/2605.12472
• H. S. Witsenhausen. Some aspects of convexity useful in information theory. IEEE Trans. Inf. Theory, 26(3):265–271, 1980. (The $K(n) \le n/2 + 1$ upper bound and discreteness.)
• J. G. Smith. The information capacity of amplitude- and variance-constrained scalar Gaussian channels. Inform. Control, 18:203–219, 1971. (The original use of the discreteness-of-capacity-achieving-input idea, applied to the AWGN channel under amplitude constraints.)
• R. E. Blahut. Computation of channel capacity and rate-distortion functions. IEEE Trans. Inf. Theory, 18(4):460–473, 1972. (The fixed-point iteration we use.)
• S. Arimoto. An algorithm for computing the capacity of arbitrary discrete memoryless channels. IEEE Trans. Inf. Theory, 18(1):14–20, 1972.
• scipy.special.gammaln, scipy.special.betaln, scipy.special.xlogy — the actual reason this post is in Python.

@article{BaniasadiBarlettaDytso2026,
  author  = {Mohammadamin Baniasadi and Luca Barletta and Alex Dytso},
  title   = {{An Improved Lower Bound on Support Size of Capacity-Achieving
             Inputs for the Binomial Channel: Extended version}},
  journal = {arXiv preprint arXiv:2605.12472},
  year    = {2026},
  eprint  = {2605.12472},
  archivePrefix = {arXiv},
  primaryClass  = {cs.IT},
  url     = {https://arxiv.org/abs/2605.12472}
}