\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{lmodern}
\usepackage{amsmath,amssymb,amsthm,mathtools,bm}
\usepackage{enumitem,microtype}
\usepackage[hidelinks]{hyperref}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\DeclareMathOperator{\dist}{dist}
\DeclareMathOperator{\Var}{Var}
\DeclareMathOperator{\TV}{TV}
\let\epsilon\varepsilon
\newcommand{\E}{\mathbb E}
\renewcommand{\P}{\mathbb P}
\allowdisplaybreaks[2]
\title{An FPRAS for Cell-Bounded Contingency Tables}
\author{OpenAI}
\date{September 24, 2026}

\begin{document}
\maketitle

\begin{abstract}
We give a fully polynomial randomized approximation scheme (FPRAS)
for counting nonnegative integer matrices with prescribed row sums,
column sums, and individual entry bounds.  Both dimensions vary, all
numerical data are encoded in binary, and zero bounds (forbidden
entries) are allowed.  The algorithm uses only independent unbiased
random bits, and a single fixed polynomial in the input length, in
$1/\epsilon$, and in $\log(1/\delta)$ bounds its total bit operations
\emph{on every execution}, including executions whose numerical output
falls outside the requested relative interval.  The same argument
yields an FPRAS and an almost-uniform sampler for bounded integral
single-commodity flows in arbitrary explicitly listed networks.
\end{abstract}

\section*{How to read this paper}
\addcontentsline{toc}{section}{How to read this paper}

The paper has one theorem and one corollary, and one long proof of
each.  A reader meeting the problem for the first time should read
this guide, then the introduction (Section~\ref{sec:intro}), then the
statement of the main estimate in Section~\ref{sec:scales} together with
the contamination bound at its end, then the annealing algorithm in
Section~\ref{sec:algorithm} and the proof of the main theorem there.
Those four places contain the whole idea.  The remaining sections are
the technical engines that make the idea rigorous:
Section~\ref{sec:scales} also builds the small/large decomposition;
Section~\ref{sec:integer-signatures} proves the summation lemma and the
near-boundary estimate; Section~\ref{sec:signatures} constructs the two
factorial lifts and their quadratic signatures;
Section~\ref{sec:transport} proves the transport, variance, and trace
bounds; and Section~\ref{sec:evaluation} implements the approximate
evaluation of the completion weights.

Three conventions recur.  Every count is an integer; every capacity is
written in binary; every randomized procedure draws only independent
unbiased bits (``fair coins'') and runs for a deterministic number of
steps chosen in advance.  A quantity of the form ``$1-\epsilon$'' is a
\emph{relative} error, not an absolute one.

\section{Introduction}\label{sec:intro}

\paragraph{The problem.}
A \emph{cell-bounded contingency table} is a nonnegative integer matrix
with prescribed row sums, prescribed column sums, and prescribed
individual entry bounds.  Formally, for positive integers $m,n$, margins
$r\in\mathbb Z_{\ge0}^m$ and $c\in\mathbb Z_{\ge0}^n$ of equal total, and
bounds $b\in\mathbb Z_{\ge0}^{m\times n}$, let
\begin{equation}\label{eq:feasible-set}
\begin{aligned}
 \Omega(r,c,b)&=\Bigl\{X\in\mathbb Z_{\ge0}^{m\times n}:\;
   \sum_{j=1}^n X_{ij}=r_i,\quad
   \sum_{i=1}^m X_{ij}=c_j,\quad
   X_{ij}\le b_{ij}\ \text{for all } i,j\Bigr\},\\
 Z(r,c,b)&=|\Omega(r,c,b)|.
\end{aligned}
\end{equation}
Each table is counted once.  A bound $b_{ij}=0$ forbids the entry; no
positivity is imposed on the support or the margins.  In the binary
encoding model a capacity may be exponentially larger than its
description, and both $m$ and $n$ are part of the input.  The problem is
to compute $Z(r,c,b)$, or to approximate it well in this model.

\paragraph{A tiny example.}
Take $m=n=2$, $r=(2,1)$, $c=(1,2)$, and every $b_{ij}=2$.  Writing
$X=\begin{psmallmatrix}a&b\\c&d\end{psmallmatrix}$, the constraints read
$a+b=2$, $c+d=1$, $a+c=1$, $b+d=2$.  Solving gives $a=d$, $c=1-a$,
$b=2-a$ with $a\in\{0,1\}$, so
\[
 \begin{pmatrix}0&2\\1&0\end{pmatrix},
 \qquad
 \begin{pmatrix}1&1\\0&1\end{pmatrix},
\]
and $Z=2$.  Real inputs have many rows and columns and astronomically
large capacities.

\paragraph{Why it is hard.}
Deciding whether a feasible table exists is easy: build a directed
network with a source, the rows, the columns, and a sink; give the
source--row, row--column, and column--sink arcs capacities $r_i$,
$b_{ij}$, and $c_j$; then a feasible table is exactly an integral flow
of value $N=\sum_i r_i=\sum_j c_j$, and Edmonds--Karp decides
feasibility in time polynomial in the number of \emph{bits} of the
input (see \cite{EK1972}).  Counting is different.  Dyer, Kannan, and
Mount \cite[Theorem~1]{DKM1997} proved that exact counting is
$\#\mathrm P$-complete already for two rows, so no exact
polynomial-time algorithm is expected.  An approximation is still
possible; that is the object of an FPRAS.

\paragraph{What an FPRAS is.}
An FPRAS for $Z$ receives the data together with rational
$\varepsilon,\delta\in(0,1)$ and returns a nonnegative rational
$\widehat Z$ with
\[
 \P\bigl((1-\varepsilon)Z\le\widehat Z\le(1+\varepsilon)Z\bigr)\ge1-\delta.
\]
``Fully polynomial'' means the running time is bounded by a fixed
polynomial in the number of bits $L$ of the input (including the
encodings of $\varepsilon$ and $\delta$), in $1/\varepsilon$, and in
$\log(1/\delta)$ -- never in $1/\delta$ or in the numeric capacities
themselves.  Bit complexity is the strict form of this demand: an input
that names a cap of $2^{100}$ must cost polynomial time in the roughly
$100$ digits, not $2^{100}$ steps.  Most naive approaches fail exactly
here.

\paragraph{The main result.}
The paper proves the following.

\begin{theorem}\label{thm:main}
There is one randomized algorithm with the following guarantee.  Given
$(r,c,b)$ as above and rational $\varepsilon,\delta\in(0,1)$, it outputs
a nonnegative rational $\widehat Z$ such that
\[
 \P\bigl((1-\varepsilon)Z(r,c,b)\le\widehat Z
                 \le(1+\varepsilon)Z(r,c,b)\bigr)\ge1-\delta.
\]
If $\Omega(r,c,b)$ is empty, the output is $0$ on every execution.  The
algorithm uses independent unbiased bits, and one fixed polynomial in
$L$, $\varepsilon^{-1}$, and $\log(\delta^{-1})$ bounds its total bit
operations on every execution, where $L$ is the total binary input
length, including the accuracy parameters.  No counting or sampling
oracle is required.
\end{theorem}

The theorem allows arbitrary cell bounds and margins: no fixed number of
rows or columns, no lower bound on positive capacities, no density or
balance assumption.  Its running-time guarantee covers preprocessing,
random-bit generation, rational arithmetic, and output, even on
executions in the failure event.  The polynomial exponents are
deliberately generous; the content is a uniform complexity guarantee,
not a practical running-time estimate.

\paragraph{A second result: flows.}
Section~\ref{sec:flows} records a consequence for integer flows.  For an
explicitly given directed multigraph with finite integer arc bounds and
prescribed vertex balances, a single bijection to a cell-bounded table
transfers Theorem~\ref{thm:main}, yielding both an FPRAS for the number
of feasible integral flow vectors and a sampler that is within total
variation $\tau$ of uniform at cost polynomial in the input length and
$\tau^{-1}$.

\paragraph{History and significance.}
The enumeration of fixed-margin arrays is classical in statistics and
combinatorics; Diaconis and Gangolli \cite{DG1995} survey its methods.
Ordinary contingency tables are the special case of \eqref{eq:feasible-set}
in which the margins force $X_{ij}\le\min\{r_i,c_j\}$; even these are
$\#\mathrm P$-complete to count exactly for two rows \cite{DKM1997}.
Approximate counting is nevertheless efficient under broader conditions
than exact counting.  For ordinary tables, Dyer, Kannan, and Mount
\cite{DKM1997} gave geometric sampling and counting for sufficiently
large margins, and Morris \cite{Morris2002} weakened the required margin
bounds.  Dyer and Greenhill \cite{DG2000} handled arbitrary two-row
margins.  Cryan and Dyer \cite{CD2003}, and later Dyer \cite{Dyer2003},
obtained an FPRAS for a fixed number of rows, and Gopalan, Klivans,
Meka, \v{S}tefankovi\v{c}, Vempala, and Vigoda \cite{GKM2011} gave a
deterministic FPTAS for a fixed number of rows.  These results already
accept binary-encoded margins; their restriction is on the number of
rows or on the sizes of the margins, not on the encoding.

Individual caps connect the table problem to several other models.  When
every cap is $0$ or $1$ the table is a subgraph of the allowed bipartite
graph with prescribed degrees; Jerrum, Sinclair, and Vigoda
\cite{JSV2004} give an FPRAS for this whole case, and all-one margins
give perfect matchings.  Bez\'akov\'a, Bhatnagar, and Vigoda
\cite{BBV2007} gave a direct annealing algorithm for binary tables.
More generally, capped tables are integral flows in a bipartite network:
Cryan, Dyer, and Randall \cite{CDR2010} gave an FPRAS for integral flows
with sufficiently large tight capacities, and for cell-bounded tables
with a fixed number of rows and arbitrary positive margins and caps.
They identified unrestricted approximation as an open problem.  Their
tight-capacity reduction and maximum-capacity spanning-tree
representation are the direct antecedents of the large-cell construction
here; a forest and soft constraints let us accommodate arbitrary
mixtures of small and large bounds, with both dimensions variable and
zero bounds permitted.

Convex optimization gives complementary estimates.  Barvinok
\cite[Theorem~1.3]{Barvinok2009} bounds ordinary weighted table counts,
including forbidden cells, within a factor $N^{O(m+n)}$ where $N$ is the
total margin; Barvinok, Luria, Samorodnitsky, and Yong \cite{BLSY2010}
obtain relative approximations in complexity $N^{O(\log N)}$ for
specified smooth-margin classes.  This scale is quasipolynomial in the
numeric total, not in its binary length.  Br\"and\'en, Leake, and Pak
\cite[Theorem~2.1]{BLP2021} give capacity bounds for arbitrary
individual caps with explicit multiplicative factors in the margins and
dimensions.  Theorem~\ref{thm:main}, by contrast, delivers any requested
relative tolerance in polynomial binary time with no smoothness
assumption.  The geometric object here is the set of all integer points
of a transportation polytope, not only its vertices (for instance the
$2\times2$ tables with every margin equal to two form a segment with
three integer points and two vertices); Guo and Jerrum
\cite[Proposition~4.2]{GJ2023} treat vertices of promised $0/1$
polytopes from network matrices, a different scope.  Counting also
yields sampling by self-reduction \cite[Theorem~6.3]{JVV1986}, which
Section~\ref{sec:flows} carries out by bisecting cell intervals with
exact feasibility tests at every branch.

\paragraph{Organization.}
Section~\ref{sec:scales} performs the deterministic reductions, splits
cells into small and large, builds the maximum-capacity forest and the
soft penalty $D$, and defines the softened counts $h_j$ and $C(j)$, with
the contamination bound as its goal.  Section~\ref{sec:integer-signatures}
proves the integer summation lemma and the log-concavity of tree
marginals, and thereby completes the proof of the contamination bound.
Section~\ref{sec:signatures} introduces the paired binary encoding, the
two factorial lifts $f_j$ and $\widetilde f_j$, their quadratic
signature properties, and the defect-mass inequalities.
Section~\ref{sec:transport} builds the enlarged Metropolis chain and
proves the one-hole transport, transversal-variance, defect-mean,
observable-variance, time-average, and trace-mixing bounds.
Section~\ref{sec:evaluation} gives the bounded weight-evaluation
procedure via bins, the Pr\'ekopa--Leindler interpolation inequality, the
bin conductance bound, and the finite-bit uniform-draw lemma.
Section~\ref{sec:algorithm} assembles the annealing algorithm, proves
the reference-experiment accuracy, couples the bounded implementation to
it, and proves Theorem~\ref{thm:main} by amplification.
Section~\ref{sec:flows} proves the flow corollary.

\paragraph{Conventions.}
All logarithms with subscript $2$ are binary; $\log$ without a subscript
in exponent estimates is natural.  A symmetric matrix has
\emph{at most one positive eigenvalue} when its positive inertia index
is at most one, counting multiplicity; this property of quadratic forms
is preserved by restriction and by limits.  A sum over edges in a
Dirichlet form is over unordered edges unless stated otherwise.  Every
$\lceil\cdot\rceil$ and $\lfloor\cdot\rfloor$ of a binary logarithm used
for parameters or bit counts is computable by integer shifts and
comparisons.  The large matrices, conditioning trees, and transport
flows that appear in proofs are analytic devices: the algorithm never
forms them.

\section{Tight capacities and two scales}\label{sec:scales}

This section performs the deterministic reductions and sets every
polynomial parameter that the rest of the proof uses.  The plan is
fixed once and for all: first make each cell bound \emph{tight} (its
two endpoints attained), then split the remaining cells by size into a
polynomially bounded \emph{small} part and a \emph{large} part handled
by a spanning forest, then attach a soft penalty that measures how far
a candidate completion is from being a genuine table.  The use of tight
capacities and a maximum-capacity spanning forest follows the
representation of integral flows in \cite{CDR2010}, but the estimates
below allow both dimensions and all capacities to vary.

\subsection{Feasibility, tightening, and parameters}

\paragraph{Feasibility is a flow question.}
Form a directed network with a source $\sigma$, the row vertices, the
column vertices, and a sink $\tau$.  Give the source--row arcs
capacities $r_i$, the row--column arcs capacities $b_{ij}$, and the
column--sink arcs capacities $c_j$.  An integral flow of value
$N=\sum_i r_i=\sum_j c_j$ is exactly a feasible table: at that value
every arc incident to the source or the sink is saturated, and the
row--column values are the table entries.  The shortest-augmenting-path
algorithm of Edmonds and Karp \cite[Theorem~1]{EK1972} uses $O(VE)$
augmentations, where $V$ and $E$ are the numbers of vertices and arcs;
the present network is simple, so $E=O(V^2)$ and the augmentation count
is $O(V^3)$, independent of the numerical capacities.  Each augmenting
path is found by breadth-first search in $O(V+E)$ operations, and
starting from the zero flow, augmentation by the minimum residual
capacity preserves integrality.  All flow and residual-capacity values
have polynomial binary length, so this is a deterministic
polynomial-bit feasibility test for the present network.  If the input
is infeasible the final algorithm returns zero; we now suppose it is
feasible.

\paragraph{Tightening.}
For a feasible input let $\ell_s$ and $u_s$ be the minimum and the
maximum feasible integer values of cell $s=(i,j)$.  Each extremum is
computable by binary search using flow feasibility with one additional
upper or lower bound: a lower bound is subtracted from its cell and
from the two incident margins before the flow test, and a negative
# ^C interrupted by user!