Trace reconstruction from deletion traces
A complete and readable account of the quasipolynomial sample bound
Companion to trace_reconstruction_anchor_induction_original.txt
Version 2. Plain text. ASCII only.
TABLE OF CONTENTS
0. What this document is
1. The problem
2. The answer
3. How this compares with earlier work
4. Notation and the tools we borrow
5. The shape of the proof in eight steps
6. Doubling, anchors, probes, and the difference polynomial
7. The induction goal
8. The base case: a period-avoiding word
9. The step: a Gaussian window and an arc test
10. Shuffle: merging the copies
11. Center Fourier expansion: an admissible mask
12. Radial compactification
13. Stopping: from masks to phase parameters
14. Interpolation and the trace statistic
15. Counting traces
16. General alphabets
17. What is and is not proved
18. Glossary
19. A reading guide
20. References
========================================================================
0. WHAT THIS DOCUMENT IS
========================================================================
This is a readable companion to a research paper about trace
reconstruction. The paper proves an upper bound on how many independent
deletion traces are needed to determine an unknown string of length n.
The original is written for specialists. This version keeps every step
and every constant, but explains each one in plain words. A reader who
knows what a sequence is, what a probability is, and what a polynomial
is should be able to follow the whole argument.
Two warnings.
First, this is a proof of an information-theoretic statement. It says:
this many traces contain enough information to identify the string. It
does not give a fast algorithm for extracting the string from the
traces. The paper is explicit about this and so is this document.
Second, the constants are enormous (10^100 appears). They are
bookkeeping, chosen so that many inequalities are obviously true, not
because anybody believes they are tight. Read them as "some fixed
number large enough". The interesting content is in the exponents.
========================================================================
1. THE PROBLEM
========================================================================
Imagine a long string of symbols: a piece of DNA, or a line of text
copied by a careless scribe.
Nature makes a copy by walking along the string and, independently at
each position, keeping the symbol with probability p and deleting it
with probability q = 1 - p. What is left is called a trace. A trace is
a subsequence of the original: the surviving symbols in their original
order, with no record of where the deletions happened, and with no
labels on the survivors saying which original positions they came from.
You are given many independent traces. You are also told p and the
original length n. You must recover the original string.
Question. How many traces do you need?
Two words to fix.
Worst case means the bound must hold for every possible original
string, not merely for typical ones. An adversary picks the hardest
string.
Sample complexity means the number of traces. It does not mean the
running time of whatever program processes them.
A tiny example. Take n = 4 and x = 0110, with p = 3/4. The retention
pattern (keep, keep, delete, keep) has probability (3/4)^3 (1/4) and
produces the trace 010. The pattern (delete, keep, keep, delete)
produces 11. The pattern (keep, keep, keep, keep) produces 0110. The
pattern (delete, delete, delete, delete) produces the empty trace.
If p is small, most traces are short, and each trace tells you very
little. The whole difficulty of the subject is that p may be as small
as 1/n, in which case a typical trace has length about 1.
Two features of the problem make it hard.
The traces do not tell you where the deletions were. If they did, one
trace would trivially suffice.
Two different strings can produce traces that are easy to confuse. A
string with many repeated symbols is harder than a random string.
========================================================================
2. THE ANSWER
========================================================================
Throughout, log means natural logarithm unless another base is stated,
and exp(t) means e to the power t.
THEOREM 1 (main bound). Let n >= 2 be an integer and let 0 < p <= 1.
For every binary string x of length n, there is a reconstruction rule
that recovers x from independent deletion traces with probability at
least 2/3, provided the number of traces is at least
exp( O( p^{-1} (log n)^3 (1 + log log(2n))^6 ) ).
The O(.) hides an absolute constant that does not depend on n, on p, or
on x. The same number of traces works for all x simultaneously. The
rule may depend on n and on p (both of which are known) but not on x.
THEOREM 2 (finer form). Put q = 1 - p and
ell_* = 1 + log(1/q)
X = 1 + log(n) / ell_*.
Then the logarithm of the number of traces needed is at most
O( p^{-1} (log n) X^2 (1 + log X)^6 ).
THEOREM 3 (low deletion). Suppose q <= n^{-epsilon} for some fixed
epsilon > 0. Then either q = 0, in which case one trace suffices, or
ell_* >= 1 + epsilon log n, so X <= 1 + 1/epsilon, while p >= 1/2 for
large n. Hence polynomially many traces suffice; the degree of the
polynomial may depend on epsilon.
THEOREM 4 (general alphabets). The same asymptotic bounds hold for
strings over an arbitrary alphabet, assuming the surviving symbols are
observed exactly.
Sanity check on sizes. Fix p = 1/2 and n = 10^9. Then log n is about
20.7, (log n)^3 is about 8900, log log(2n) is about 3.05,
(1 + log log(2n))^6 is about 1250, and the product is about 1.1 * 10^7.
So the theorem promises fewer than exp(2.2 * 10^7) traces. That is a
finite number, and the shape of the exponent is the point. A trivial
bound is 2^n, which here is 10^{3 * 10^8}: vastly larger.
Vocabulary. A function of n is polynomial if it is at most n^C for some
fixed C. It is quasipolynomial if it is at most exp((log n)^C) for some
fixed C. Quasipolynomial growth sits strictly between polynomial and
exponential: for every fixed C and every fixed c > 0, exp((log n)^C)
grows faster than n^c and slower than exp(n^c).
Three things the theorem does not say.
It is not a matching lower bound. The true number of traces needed
might be much smaller.
It is not an efficient algorithm. See section 17.
For fixed p in (0,1) the bound stays quasipolynomial. It does not
become polynomial in n. Getting polynomial for any fixed p is open.
========================================================================
3. HOW THIS COMPARES WITH EARLIER WORK
========================================================================
The previous quantitative quasipolynomial estimate is due to
Burudgunte, Valiant, and Wang, cited as BVW. Their Theorem 43 states a
quasipolynomial bound with an unspecified log-power c and a p^{-7/3} in
the exponent.
If you unwind their proof (their Section 8) you find concrete
prescriptions:
b_old = ceil( log_{1.1} log_2 n )
p_base = 3^{-b_old} min(p, 1/4)
pair-test bound = 4 tau_base^{-2 C^{b_old}} for some C >= 1
log(1/tau_base) = constant * d_base^{1/3} p_base^{-1/3}
d_base = max( p_base^{-6}, O(1) )
Even if you ignore the C^{b_old} multiplier entirely, the inverse
retention term alone forces a log-power of at least
(7/3) log_{1.1}(3) > 26
for fixed p. The present result improves on that proved estimate. It is
not a claim that no other optimization of BVW's approach could do
better, and it does not go below quasipolynomial order for fixed
0 < p < 1.
What is genuinely new.
The core idea is kept from BVW: test relations between phases using
products of Fourier coefficients.
The induction is different. It runs directly on generating functions
of the source string, instead of on the pair-test structure.
The phase argument is weighted and uses up to four factors at once.
Masks are made spectrally compact: extra low-dimensional factors are
propagated with a bounded frequency footprint. That lets the proof
relax the angle budget during the induction and spend it at the very
last step.
The base case keeps the period-avoiding-word idea of Chase and
Rubinstein.
For polynomial-time bounds in the low deletion regime, see Chen, De,
Lee, Servedio, and Sinha. See also the references at the end.
========================================================================
4. NOTATION AND THE TOOLS WE BORROW
========================================================================
4.1 Notation.
log natural logarithm
exp(t) e to the power t
a^b a to the power b
|z| modulus of the complex number z
Re z, Im z real and imaginary parts
sum, prod sum and product
inf infinity
O(g), o(g) see 4.3
|x|_2 Euclidean length of a vector
|x|_inf largest coordinate in absolute value
||mu|| total variation of a complex measure
sup, inf supremum, infimum
1_E indicator of the event E
gcd(a,b) greatest common divisor
A polynomial here is always an ordinary polynomial in one complex
variable z. A quantity is analytic in a region if it is given there by
a convergent power series. A measure is a rule assigning a complex
number to each nice set, additively; its total variation is the largest
possible value of the sum of absolute values of the pieces, over all
finite partitions.
4.2 Seven standard facts.
FACT 1 (geometric series). If |z| < 1 then sum_{m>=0} z^m = 1/(1-z),
and sum_{m>=0} |z|^m = 1/(1-|z|).
FACT 2 (Cauchy-Schwarz). |sum_i x_i y_i|^2 <= (sum_i |x_i|^2)
(sum_i |y_i|^2).
FACT 3 (union bound). For any events A and B, P(A or B) <= P(A) +
P(B).
FACT 4 (Hoeffding). Let S be the average of N independent random
variables, each in [0,1]. Then for every t > 0,
P( |S - E S| >= t ) <= 2 exp(-2 N t^2).
So to resolve a gap of size t you need about log(1/delta)/t^2
samples, where delta is the error you are willing to tolerate.
FACT 5 (Poisson). If f is analytic and nonzero on a disk, then log|f|
is harmonic there, so its value at the centre is a weighted average
of its boundary values, with the Poisson weight
(1 - r^2) / (2 pi (1 + r^2 - 2 r cos(angle)))
at radius r. We only use one consequence: a large value at one
interior point forces a comparable value somewhere on any arc that
carries enough Poisson weight. Boundary zeros are harmless; apply the
inequality to log(|f| + delta) and let delta tend to 0.
FACT 6 (Hadamard three-lines). If G is analytic and |G| <= 1 in the
vertical strip a <= Re s <= b, then log sup_{Re s = t} |G| is a
convex function of t in [a,b]. Consequence: a lower bound at one
interior line plus an upper bound at both edges gives a lower bound
at the other edge, and the price is the ratio of the distances.
FACT 7 (Fine and Wilf). If a word w has periods p and q and
|w| >= p + q - gcd(p,q), then w has period gcd(p,q).
A period of a word of length L is a positive integer delta such that
every pair of positions at distance delta carries the same symbol,
for as long as both positions lie inside the word. Every word
trivially has every period delta >= L.
4.3 Big-O and little-o.
"f(n) = O(g(n))" means there is a constant C, depending only on
quantities that are explicitly fixed and not on n, such that
|f(n)| <= C g(n) for all sufficiently large n.
"f(n) = o(g(n))" means f(n)/g(n) tends to 0 as n tends to infinity.
In this paper, every O(.) and o(.) statement is uniform over the
remaining parameters: the channel probabilities in the stated range,
the induction level, the choice of pair, and so on.
"For all sufficiently large n" means there is an absolute threshold
n_0 such that the conclusion holds for n >= n_0. These thresholds are
never computed, and they are always far smaller than the actual size
of n in any application.
========================================================================
5. THE SHAPE OF THE PROOF IN EIGHT STEPS
========================================================================
Here is the whole argument at a glance. Each step is expanded in a
later section.
STEP 1. Double the input. Prepend a long block of zeros that
everybody knows. This makes room to work and puts the first mismatch
far from the left end. Section 6.
STEP 2. Look at several positions at once. A probe picks indices
i_1 < ... < i_k, tests a predicate at each, and may force some
adjacent pairs to be consecutive. Weighting by z to the power of the
last index minus the mismatch position gives a polynomial H(z).
Section 6.
STEP 3. Set up an induction. At stage j the polynomial is large at a
point near z = 1, at radius exp(-a_j). Each stage moves a much closer
to 0, i.e. closer to the unit circle. Section 7.
STEP 4. Base case. Build a probe out of a period-avoiding word whose
occurrences are far apart, so H has the shape 1 - eta z^u plus a
small tail; the Poisson inequality then supplies a pivot. Section 8.
STEP 5. Step. Multiply by a Gaussian window, then run a phase test on
the arc of radius exp(-a_j): either a product of two, three, or four
Fourier coefficients is already large, or the deficit budget would be
exceeded. The budget is a Cauchy-kernel vertex-cover statement.
Section 9.
STEP 6. Shuffle. Expand the products. Each copy contributes its own
tuple; merge all copies into one big tuple with ties allowed. The
phases collapse into one outer phase and r-1 inner phases in the
differences between latest offsets. Section 10.
STEP 7. Make the mask admissible. A Fourier expansion in the outer
offset converts the Gaussian window into a compactly supported
frequency measure, giving a legitimate mask of the new level.
Sections 11 and 12.
STEP 8. Stop and count. Expand all masks, pick one frequency
combination, interpolate the loose-edge parameters with Hadamard
three-lines, build a bounded trace statistic whose expectation
differs between the two candidates, and apply Hoeffding plus a union
bound over all 2^n rivals. Sections 13, 14, 15.
========================================================================
6. DOUBLING, ANCHORS, PROBES, AND THE DIFFERENCE POLYNOMIAL
========================================================================
6.1 Small cases.
n = 1. About 1/p traces suffice, because each trace either shows the
single symbol or not.
2 <= n <= some fixed constant. Collect ceil(p^{-n} log 3) traces.
With probability at least 2/3 at least one trace has full length, and
then you can simply read off x.
p < 1/n with n >= 2. The same full-length strategy costs at most
exp(O(p^{-1} log n)). Indeed, with s = n p <= 1,
n log(1/p) = p^{-1} s (log n + log(1/s)) = O(p^{-1} log n).
q < 1/n with n >= 2. Then p^n >= (1 - 1/n)^n >= 1/4, so a constant
number of traces suffices by the same strategy.
So assume from now on that p >= 1/n and q >= 1/n.
6.2 Doubling.
Prepend a known block of zeros of length n to each candidate. A trace
of the extended input is simulated from a trace of the original input
by prepending an independent deletion trace of the zero block. From
here on the letter n denotes the doubled length. At the very end switch
back to the original length; the exponent changes only by an absolute
factor.
Fix two distinct extended binary candidates x and y. They share the
known zero prefix. Let d >= n/2 be their first mismatch, using 0-based
indices.
6.3 The constants.
All of these are absolute constants or simple expressions in n, p, q.
They look arbitrary. They are chosen once and never changed.
D = 10^14 a safety constant
epsilon = 10^{-5} the mask decay rate
ell_* = 1 + log(1/q)
b = ceil( log_2( 1 + log n / ell_* ) )
L_h = 2^{b-h} for 0 <= h <= b
M_h = 4^{b-h} for 0 <= h <= b
a_0 = exp(-10^16) / (D L_0)
t_h = a_{h+1} = D q L_h a_h^2
P_0 = ceil(1 / a_0)
R_h = D q (h+1)^2 M_h a_h^2
B_j = 10^100 ell_* 8^j (b+1)^3
f(u) = exp(u^2 / 20)
What each one is for.
ell_* measures how much a single trace carries. It is large exactly
when p is small. It appears everywhere as the cost of seeing one
symbol.
b is the number of induction stages. For fixed p and large n it
grows like log log n; more precisely, 2^b is about log(n)/ell_*.
a_h is the radius parameter at stage h. It shrinks super-
exponentially: setting v_h = D q L_h a_h gives v_0 = q exp(-10^16)
and v_{h+1} = v_h^2 / 2, so v_h = 2 (v_0/2)^{2^h}, and a_h decays
even faster. In particular a_b n < .0001.
t_h = a_{h+1} is the decay rate that level-h masks may use. Since a
is decreasing in the index, t_h >= a_j for every h < j.
P_0 is the length parameter of the period-avoiding word. It is
O(2^b), so log(2 P_0) = O(b).
R_h is the decay rate of the compact representation of a level-h
mask, and R_h / t_h = (h+1)^2 M_h / L_h >= 1.
B_j is the log-scale accuracy demanded at stage j. It multiplies by 8
at each stage, but 8^b is only about X^3.
Easy consequences, used repeatedly:
b >= 1;
ell_* (2^b - 1) >= log n, hence ell_* (b+1) -> infinity, hence
B_0 -> infinity uniformly;
a_{h+1}/a_h = v_h <= v_0 = q exp(-10^16), which is far below .001;
log(1/a_h) <= 10^18 (b + 2^h) ell_*;
P_0 = O(2^b), so log(2 P_0) = O(b), independent of ell_*.
6.4 Probes.
A tuple probe of order k >= 1 consists of
indices 0 <= i_1 < ... < i_k < n,
predicates chi_e : {0,1} -> {0,1}, for e = 1,...,k,
edges the pairs (i_e, i_{e+1}) for e = 1,...,k-1.
The value of the probe at that tuple is the product of the numbers
chi_e(x_{i_e}) over e.
Some edges are designated exact. An exact edge carries the extra
factor 1_{i_{e+1} - i_e = 1}; it forces the two indices to be
consecutive and contributes nothing otherwise. All other edges are
called loose.
The block count A is one plus the number of loose edges. Equivalently,
A is the number of maximal runs of consecutive forced indices. Note
that different blocks need not be separated by gaps of size greater
than one; they can touch.
6.5 Masks.
A probe at stage j also carries finitely many masks of levels h < j.
A mask of level h uses r-1 differences, with 2 <= r <= 4:
g_s = i_{b_s} - i_{c_s} for s = 1,...,r-1,
for some chosen indices 1 <= c_s <= b_s <= k. All of them are
nonnegative by the ordering, and one may be zero if the two endpoints
coincide.
A mask is a factor
K(g) * prod_{s=1}^{r-1} exp(i beta_s g_s),
where the beta_s are real with |beta_s| <= 10^4 a_h, and K may be
complex. K must satisfy, for all g_s >= 0,
(M1) |K(g)| <= 2 exp( - epsilon t_h (g_1 + ... + g_{r-1}) ),
and must also have a compact representation
(M2) K(g) = exp( - R_h (g_1 + ... + g_{r-1}) )
* int_{R^{r-1}} exp(i xi . g) mu(d xi),
with mu a complex measure supported on the cube |xi|_inf <= C R_h and
total variation
||mu|| <= exp( C (R_h / t_h) B_h ).
Here C is a large absolute constant, the same at every level. Its
existence follows from the construction, and one value serves for all
stages.
How to read (M1) and (M2).
(M1) says that a mask is negligible unless all of its differences are
small, of order 1/t_h. This is what keeps the probe's sums bounded.
(M2) says that a mask is an average of pure phases with bounded
frequencies. In other words, a mask behaves like a finite combination
of honest phase factors, so that later steps can treat it as a tuple
with extra phases. The total variation is the price of that
approximation, and it is exponential in (R_h/t_h) B_h.
A difference g_s covers exactly the edges (i_e, i_{e+1}) inside its
interval, counting each edge once; a zero difference covers nothing.
The masks' difference intervals are how the induction remembers
things. At each stage, the loose edges of the new probe must be covered
by masks, either old ones or the one new mask created at that stage.
6.6 Latest-anchor arrays and the difference polynomial.
Let W_x(i_1,...,i_k) be the product of all bit, exact-gap, and mask
factors of a probe. Define
s_x(m) = sum over 0 <= i_1 < ... < i_k < n with i_k = d + m
of W_x(i_1,...,i_k),
and define s_y(m) the same way with y. Invalid latest positions
contribute zero. Now define the difference polynomial
H(z) = sum_m ( s_x(m) - s_y(m) ) z^m.
Key observation. For m < 0 the condition i_k = d + m forces every index
to be less than d. On that region x and y agree, so each term cancels
against its counterpart. Hence H has no negative powers. It is an
ordinary polynomial.
This is the whole point of the latest-anchor convention. The latest
index is either at the mismatch or to the right of it; anything to the
left has already been decided and cancels.
========================================================================
7. THE INDUCTION GOAL
========================================================================
At stage j, with 0 <= j <= b, we want a probe with the following
properties.
(G1) Order and blocks.
1 <= k <= 2 P_0 4^j and 1 <= A <= 4^j.
(G2) Masks. For every level h < j there are at most 4^{j-h-1} masks of
that level, and every loose edge is covered by at least one of
their differences.
(G3) A pivot. There is a real theta with |theta| <= 5 a_j such that
|H( exp(-a_j + i theta) )| >= exp( - B_j / f(theta/a_j) ).
Stop as soon as a_j n <= .001. This happens by j = b, because
a_b n < .0001. At the moment of stopping, the radius parameter a_j is at
most .001/n, i.e. very far inside the unit disk.
Interpretation of (G3). The point exp(-a_j + i theta) sits just inside
the unit circle, at distance about a_j from it. The factor
1/f(theta/a_j) is between exp(-25/20) and 1 when |theta| <= 5 a_j, so it
is a mild penalty for choosing a direction theta away from the positive
real axis.
Size of the sums. Every probe satisfying (G1) and (G2) has
sup_m |s_x(m)| <= Z_j, sup_m |s_y(m)| <= Z_j,
where
Z_j = 2^{4^j} * ( 2 / (epsilon a_j) )^{4^j}.
Why. Each loose edge contributes a decay factor at least
exp(-epsilon a_j Delta), where Delta is its gap, because every level
h < j has t_h >= a_j and the coverage condition forces at least one
applicable mask. Summing over Delta >= 1 gives at most
1/(1 - exp(-epsilon a_j)) <= 2/(epsilon a_j). There are at most
A <= 4^j loose edges. Bit predicates are at most 1. Each mask factor has
modulus at most 2 by (M1), and there are at most 4^j of them, giving
2^{4^j}. Exact edges are fixed at distance 1 and contribute 1.
Consequence, used many times below: the auxiliary magnitudes are tiny
compared to the accuracy budget. For 0 <= j < b,
log Z_j + log Z_{j+1} + log(1/a_j) + log(1/t_j)
<= 10^21 ell_* 4^j (b + 2^j)
<= 10^{-60} B_j.
Proof: each term is at most about 4^j log(1/a_j), and
log(1/a_h) <= 10^18 (b + 2^h) ell_*. So the left side is at most
10^18 * 4^j * (b + 2^j) * ell_* times a small constant. And
4^j (b + 2^j) <= 8^j (b + 1) <= 8^j (b + 1)^3.
Since B_j contains 10^100 * 8^j * (b + 1)^3 * ell_*, there is enormous
room.
========================================================================
8. THE BASE CASE: A PERIOD-AVOIDING WORD
========================================================================
We must produce a probe at stage 0. It will have all edges exact, no
masks, and order k = 2 P_0.
8.1 Finding the word.
By the padding, d >= 2 P_0 for all sufficiently large n. Let
w' = x_{d-2P_0+1} ... x_{d-1} = y_{d-2P_0+1} ... y_{d-1}
be the common last 2 P_0 - 1 symbols before the mismatch.
CLAIM. At least one of the two words w'0 and w'1 has no period strictly
smaller than P_0.
Proof. Suppose w'0 has period p_1 < P_0 and w'1 has period p_2 < P_0.
Then w', being a prefix of each, has both periods. Since
|w'| = 2 P_0 - 1 >= p_1 + p_2 - gcd(p_1, p_2)
(because p_1, p_2 <= P_0 - 1), Fine and Wilf gives that w' has period
g = gcd(p_1, p_2).
Look at positions 2 P_0 - p_1 and 2 P_0 - p_2 inside w'. They differ by
p_2 - p_1, which is a multiple of g, so they carry the same symbol.
But period p_1 of w'0 forces (w'0)_{2P_0 - p_1} to equal (w'0)_{2P_0},
which is 0. And period p_2 of w'1 forces (w'1)_{2P_0 - p_2} to equal
(w'1)_{2P_0}, which is 1. So the common symbol would be both 0 and 1.
Contradiction.
Choose w in {w'0, w'1} with no period below P_0.
8.2 Building the base probe.
Let k = 2 P_0. Take all edges exact. Take predicates chi_e that test
the e-th symbol of w. There are no masks, so A = 1.
The number of occurrences of w ending exactly at d differs between x
and y: by construction one of them has w there and the other does not.
Hence s_x(0) - s_y(0) is +1 or -1, so the two strings already disagree
at offset 0.
8.3 Occurrences are far apart.
Suppose w occurs in x ending at d + m_1 and again ending at d + m_2,
with 0 < m_2 - m_1 < P_0. Call delta = m_2 - m_1.
The two occurrences have length 2 P_0 and are shifted by delta. At every
position covered by both, the symbol is determined twice: once from the
first occurrence and once from the second. Reading off the two
determinations, we get w_j = w_{j+delta} whenever both indices lie in
[1, 2 P_0]. That is exactly the statement that w has period delta.
Since delta < P_0, this contradicts the choice of w. The same argument
applies in y. So occurrence endings in each individual string are
P_0-separated.
8.4 The shape of H.
Putting the two previous paragraphs together, up to an overall sign,
H(z) = 1 - eta z^u + sum_{m >= P_0} c_m z^m
with
1 <= u < P_0, eta in {0,1}, |c_m| <= 1,
and if eta = 0 then u is arbitrary, since the term -eta z^u is then
absent for every choice of u.
Explanation. s_x(0) - s_y(0) = +-1 gives the constant term. The other
string has its nearest occurrence ending at d + u for some 1 <= u < P_0,
or none at all, in which case eta = 0; this gives the -eta z^u term.
Every further occurrence, in either string, ends at least P_0 away from
d, giving the tail.
8.5 The pivot by the Poisson inequality.
Let v = 10 log(P_0) / P_0. At the positive real point z = exp(-v):
The first two terms have modulus at least
1 - exp(-u v) >= 1 - exp(-v) >= v/2,
since v is tiny.
The tail is at most
sum_{m >= P_0} exp(-v m) = exp(-v P_0)/(1 - exp(-v))
<= 2 P_0^{-10} / v <= v/6,
because exp(-v P_0) = P_0^{-10} and 1 - exp(-v) >= v/2.
Therefore |H(exp(-v))| >= v/2 - v/6 = v/3.
On the circle |z| = exp(-a_0) the crude bound |H(z)| <= 2/a_0 holds.
There are at most about 1/a_0 contributing offsets, each coefficient at
most 2; since 1/a_0 is astronomically large, this is trivial.
The Poisson inequality on the disk of radius exp(-a_0) centred at 0,
evaluated at the interior point exp(-v) which is much closer to 1, then
gives
max_{|theta| <= a_0} |H(exp(-a_0 + i theta))|
>= exp( - O( (v/a_0) log(1/a_0) ) ).
Indeed, the harmonic weight of that arc is at least a constant times
a_0/v (because v >= 2 a_0 and v < .01), and log|H| is bounded above by
log(2/a_0) elsewhere. Boundary zeros are removed by the usual
regularization log(|H| + delta), delta -> 0.
Size of the exponent. We have
log(1/a_0) = 10^16 + log(D L_0) <= 10^16 + 35 + b,
v/a_0 <= 10 log(2 P_0) <= 10 (10^16 + 36 + b),
so the exponent is at most about 10 (10^16 + 36 + b)^2, which is under
1.1 * 10^33 for any conceivable n. On the other hand B_0 is at least
10^100. So the exponent is at most 10^{-66} B_0, comfortably below the
requirement B_0 / f(theta/a_0) >= 0.94 B_0.
This gives (G1), (G2), and (G3) at j = 0: order 2 P_0 <= 2 P_0 4^0,
block count 1 <= 4^0, no masks, and |theta| <= a_0 <= 5 a_0.
========================================================================
9. THE STEP: A GAUSSIAN WINDOW AND AN ARC TEST
========================================================================
Suppose (G1), (G2), (G3) hold at stage j < b. Write
a = a_j, t = a_{j+1}, B = B_j, E = B / f(theta/a_j).
Recall t/a <= v_0, far below .01; a <= a_0, far below .001; and, during
a genuine step (one not at the stopping point), a > .001/n.
9.1 Step one: a window.
Define the window
sigma^2 = 2.4 B / t^2
ell = floor( 10^{-10} E )
v(m) = exp(-y^2/2) * sum_{k=0}^{ell} (y^2/2)^k / k!,
with y = m/sigma throughout. The function v is a polynomial in m^2 of
degree 2 ell; it satisfies 0 <= v <= 1; and it is a smooth window, equal
to 1 inside a window and 0 far outside.
CLAIM 9.1. For 0 <= m <= 3 E/a we have 1 - v(m) <= exp(-2E).
Proof. On this range,
y^2/2 <= 2 E (t/a)^2 <= 2 E exp(-2 * 10^16),
because t/a <= exp(-10^16). So w = y^2/2 is astronomically small, and
the Taylor remainder of the exponential series after ell terms,
multiplied by exp(-w), is at most (e w/(ell+1))^{ell+1}, which is far
below exp(-2E) because ell is a positive fraction of E.
Interpretation: the window radius 3E/a is astronomically smaller than
the Gaussian scale sigma, so within the window the correction is
negligible.
Project to real or imaginary part. Choose c(s) = Re(s) or c(s) = Im(s),
whichever keeps at least half of the pivot. Since H = Re H + i Im H,
one of the two has modulus at least half of |H|.
Define, for 0 <= m with d + m < n, and zero elsewhere,
f_x(m) = ( 2 + c( s_x(m)/Z_j ) ) v(m) / Y
f_y(m) = ( 2 + c( s_y(m)/Z_j ) ) v(m) / Y
Y = 3 sum_m v(m).
Two immediate facts. Since |s_x|, |s_y| <= Z_j, we have
0 <= 2 + c(s/Z_j) <= 3, so 0 <= f_x, f_y and sum_m f_x(m), sum_m f_y(m)
<= 1. Also Y >= 3, because the term m = 0 has v(0) = 1; so there is no
division by zero.
The normalizing factor is harmless: log Y <= 10^{-8} E. To see this,
split at H_* = ceil(1000 B / t). The sum of v within that range is at
most about 2 H_*, and outside it v(m) is at most
2^ell exp(-m^2/(4 sigma^2)), which decays like a Gaussian. Using
log(1/t) <= 10^{-60} B, B <= 4E, and log B = o(E), the claim follows.
9.2 Step two: a polynomial on the unit circle.
Define
J(z) = sum_m ( f_x(m) - f_y(m) ) z^m.
Two bounds.
On the whole unit circle, sup_{|z| = 1} |J(z)| <= 2, because
sum_m |f_x - f_y| <= sum_m f_x + sum_m f_y <= 2.
At the pivot z_0 = exp(-a + i theta),
|J(z_0)| >= exp( -1.0005 E ).
Reason: masking costs at most exp(-2E) on the window (Claim 9.1) and
is exponentially small beyond it; 1/Y costs at most exp(-10^{-8} E);
log Z_j costs at most 10^{-60} B. Adding these up gives the stated
loss.
The truncation of J to |m| <= H_* is denoted J_diamond. Because all
negative powers of z vanish in the difference, and the omitted tail has
l_1 norm below exp(-200E), we also get
|J_diamond(z_0)| >= exp( -1.0006 E ).
9.3 Step three: a discrete grid on the arc.
Put
N = ceil( (1000 B / t)^6 ),
h = 2 pi / (N a),
v_0^* = theta / a,
gamma = floor( 10^4 / h ),
F(l) = sum_m f_x(m) exp( 2 pi i l m / N ),
G(l) = sum_m f_y(m) exp( 2 pi i l m / N ).
The grid points are z_l = exp(-a) exp(2 pi i l / N) = exp(-a + i a l h),
for |l| <= gamma. Note that |F|, |G| <= 1.
CLAIM 9.2 (phase test). Put Lambda = exp(-6B). Then there exist
r in {2,3,4} and integers l_1,...,l_r with
|l_i| <= gamma and |l_1 + ... + l_r| <= 1,
such that
| prod_i F(l_i) - prod_i G(l_i) | > Lambda.
Sketch of proof. Introduce the deficits
m_l = max( |F(l)|, |G(l)| ) (so 0 <= m_l <= 1)
g_l = min( 50, E^{-1} log(1/m_l) )
q_l = min( 50, E^{-1} log(2 m_l / |F(l) - G(l)|) ),
with the conventions g_l = 50 when m_l = 0, and q_l = 50 when the
denominator vanishes. Both are nonnegative and symmetric, because F and
G are real on the grid by construction.
The grid weights are
k_l = h / ( pi ( 1 + (l h - v_0^*)^2 ) ) for |l| <= gamma,
and we define
G_cal = sum_l k_l g_l, Q = sum_l k_l q_l.
The Poisson inequality, applied to the pivot bound on J_diamond, implies
the deficit budget
G_cal + Q <= 1.003.
The proof is a Riemann-sum approximation to the Poisson integral.
Concretely, for each l,
g_l + q_l <= min( 100, E^{-1} log( 2 / |J(exp(-a + i a l h))| ) )
up to the truncation correction, and the Poisson integral of that capped
quantity is at most 1.0025. The rest is bookkeeping. The Poisson density
in the rescaled variable x = l h, for the disk of radius exp(-a), is at
least (1 - 2a)/(pi (1 + (x - v_0^*)^2)), and the grid sum differs from
the integral by o(1) uniformly, by layer cake and bounded variation.
Now the combinatorial core. For every integer j'shift with
1 <= j'shift <= S/h, every fractional vertex cover of the graph on the
integer points in [-gamma, gamma] with edges (l, l + j'shift) has total
k_l weight at least M. Here S and M are the numbers in the table below.
Proof. Restrict to one residue class mod j'shift: it is a chain. Extend
the weights to the full Cauchy kernel on the whole line; the tail costs
less than .0001. A fractional cover assigns a number in [0,1] to each
vertex, with endpoint sums at least 1. It can be rounded to a random
binary cover without increasing expected weight: pick one uniform
threshold U_0 in [0,1], use "vertex number >= U_0" on even vertices and
"vertex number >= 1 - U_0" on odd vertices. So it is enough to bound
binary covers.
In each chain, take the central vertex l nearest to v_0^*/h. The
weights decrease outward along both arms. A binary cover either omits
the centre, in which case both first neighbours are forced and the
disjoint edges (2,3), (4,5), ... force the odd weights on both arms; or
it includes the centre, in which case the disjoint edges (1,2), (3,4),
... force the even weights. So the minimum is the lesser of the two
parity masses.
The parity masses have closed forms. With d' = j'shift * h,
x = l h - v_0^*, and K_Cauchy(x) = 1/(pi (1 + x^2)),
Xi(x) = sum_{nu in Z} K_Cauchy(x + 2 nu d')
= (1/(2d')) * (1 + 2 sum_{m>=1} exp(-m pi/d') cos(m pi x/d')),
using the Fourier transform identity
int K_Cauchy(x) exp(i xi x) dx = exp(-|xi|) and summing over periods.
The even mass is h Xi(x). The odd mass is h Xi(x + d'), which is never
larger: write Xi in closed form with r_0 = exp(-pi/d'),
Xi(x) = (1 - r_0^2) / (2 d' (1 + r_0^2 - 2 r_0 cos(pi x/d'))),
and cos is maximal at x = 0.
Summing the odd masses over the j'shift central representatives, up to
O(h) error, gives the integral
int_{-d'/2}^{d'/2} Xi(x + d') dx
= 1/2 - (2/pi) arctan( exp(-pi/d') ).
This is smallest at the largest possible d', which is S. So the cover
weight is at least M.
Second, directly by Riemann sum,
sum_{|l| <= S/h} k_l > J_*.
Table. Round |v_0^*| up to w in {1,2,3,4,5}. Set
w 1 2 3 4 5
S 1.5 2.1 2.8 3.6 4.4
M .419 .357 .297 .246 .208
J_* .524 .453 .380 .335 .292
These are conservative, in the sense of the four inequalities
M < 1/2 - (2/pi) arctan( exp(-pi/S) ) - .0002
J_* < (1/pi)( arctan(S+w) + arctan(S-w) ) - .0002
T_lim(v) = 6 f(v) - .1
T_lim(w-1) > max{ 2.006/M,
1.005 / (J_* - (1 + 2 J_*/M)^2 M/16) }.
All of them can be checked with rational arithmetic, without any
sensitive numerical integration:
use 3.14159 < pi < 3.14160;
bound exp(-pi/S) from above by a truncated exponential series;
write arctan(S+w) = pi/2 - arctan(1/(S+w));
the remaining arctangent arguments have modulus at most .6, so
evaluating their alternating series to six terms costs at most
.6^{15}/15;
lower-bound exp(u) in f(w-1) by 1 + u + u^2/2.
The resulting worst-endpoint values are
w 1 2 3 4 5
max{...} 4.96 5.79 6.94 8.26 9.71
while
T_lim(w-1) = 6 exp((w-1)^2/20) - .1
equals about 5.90, 6.52, 7.41, 8.73, 10.51. Each exceeds the
corresponding maximum.
Now the contradiction. Put R = G_cal/M and T_0 = T_lim(v_0^*), which is
at least T_lim(w-1). Then 2R < T_0, and there exists an integer s with
|s| <= S/h, 2R + q_s < T_0, 2g_s + q_s < T_0.
Otherwise, pointwise on the whole interval,
q_s >= (T_0 - 2R) max( 0, 1 - 2 g_s/T_0 ),
so by the cover bound
Q >= (T_0 - 2 G_cal/M)(J_* - 2 G_cal/T_0),
and minimizing this quadratic in G_cal gives
G_cal + Q >= T_0 ( J_* - (1 + 2 J_*/M)^2 M/16 ) > 1.004,
contradicting G_cal + Q <= 1.003. Choose the witness with s >= 0 by
symmetry.
For each 0 <= j'shift <= s, choose a reference pair (e_j, e_j +
j'shift) whose end deficits sum to at most R. Such a pair exists: if
the minimum were R' > R, the clipped deficits min(1, g_l/R') would be a
fractional cover of cost G_cal/R' < M, contradicting the cover bound.
For j'shift = 0, minimize 2 g_l and use the total kernel mass > .999
together with M < .499.
Suppose the claim fails. At every pair (l, -l), the failure bounds the
modulus of F(l) - G(l) by Lambda/m_l. At all relevant indices
2 g_l < T_0 < 50, so the cap is inactive and m_l = exp(-g_l E). Write
U(l) = arg F(l) - arg G(l) modulo 2 pi.
By reality of F and G, U(-l) = -U(l). Set
d_j = U(e_j + j) - U(e_j), so d_0 = 0.
The failed quadruple product at frequencies
(e_j + j, -e_j, -e_{j-1} - j + 1, e_{j-1})
sums to 1, so it bounds d_j - d_{j-1} mod 2 pi by 100 Lambda exp(2 R E):
both products have modulus at least exp(-2RE)/16, and two complex
numbers differing by at most Lambda, both of modulus at least v, have
angular distance at most pi Lambda/v.
Summing over j gives
|d_s| <= 100 N Lambda exp(2 R E) modulo 2 pi.
Similarly, the failed triple product at (e_s + s, -e_s, -s) sums to 1
and gives
|U(s) - d_s| <= 100 Lambda exp( (R + g_s) E ).
But R + g_s + q_s <= max(2R + q_s, 2g_s + q_s) < T_0. Therefore
|F(s) - G(s)| / m_s
<= Lambda/m_s^2
+ 100 Lambda ( exp((R+g_s)E) + N exp(2RE) )
< exp(-q_s E),
using log N <= 10^{-8} E and R + g_s + q_s < T_0 < 50. That contradicts
the uncapped definition of q_s, which was supposed to stay below 50. The
claim holds.
Consequence. We now have r in {2,3,4} and angles theta_1,...,theta_r
with
|theta_s| <= 10^4 a and |theta_1 + ... + theta_r| <= .01 t,
and a product discrepancy larger than Lambda.
========================================================================
10. SHUFFLE: MERGING THE COPIES
========================================================================
The discrepancy is a difference of products. We expand it.
10.1 Components.
Each array copy contributing a factor (one F or one G) is a bounded
combination of three component types:
(C1) a constant "latest valid" probe with a single wildcard position
(value 1),
(C2) the stage-j probe,
(C3) the complex conjugate of the stage-j probe (which also
conjugates the masks).
Each copy also carries the Gaussian weight v(m_s), the phase
exp(i theta_s m_s), and a scalar coefficient of modulus at most 1 (from
the floor, from 1/Y, from 1/Z_j, and from the projection c).
10.2 Categories.
Group the expanded products into categories. A category records
for each copy, its type (C1), (C2), or (C3);
the total weak ordering of all selected source positions across all
copies, including all equalities, i.e. all ties.
The number P of nonempty categories satisfies log P = o(B), uniformly.
Why. Once the types are fixed and all exact adjacencies inside each copy
are honoured, each copy contributes at most 4^{j+1} contiguous blocks.
Specify the weak order of the block starts: at most m! 2^m possibilities
for m <= 4^{j+1} blocks. Then process blocks in start order. A new start
either lies strictly beyond the maximum end index seen so far, or its
offset from that maximum is an integer between 0 and k-1. This
determines the entire weak ordering. The total number of categories is
at most 81 (4k+1)^{4^{j+1}}, so
log P = O( 4^j log( 4k + 4^{j+1} + 2 ) ) = O(4^j b) = o(B),
including the case k = 1.
10.3 Merged tuples.
Fix a category. All of its exact adjacency constraints act across
adjacent distinct slots of the merged order: if an old exact edge
spanned an intervening distinct slot, the category would be empty.
Collect these as the merged exact edges, and declare the rest loose.
The merged order has k' <= 4k slots and at most as many blocks as the
sum of the old block counts (each old block is contained in one merged
block, and each merged block contains at least one old block). Tie
bit-predicates are multiplied. Retain every old mask unchanged; each
contributes at most 4 copies here. The wildcard component contributes
no mask.
Coverage. Every loose merged edge is covered either by an inherited
mask or by the new mask we are about to add. Consider a copy with a slot
at the left endpoint of a loose merged edge. If that slot is the copy's
latest, it is not the overall latest (since the copy is not the
component chosen last), so it contributes a new covering interval. If it
is not the copy's latest, then the copy has a consecutive pair
straddling the edge; that pair was loose in the copy, so an old mask
already covered it.
Latest component and differences. Choose the component whose latest
selected index is last in the category, with ties broken
deterministically, and relabel it as s = r. Let Y_0 be the overall
latest offset. Put
g_r = 0,
g_s = Y_0 - m_s >= 0 for s < r,
g_bar = (g_1 + ... + g_r)/r,
Q(g) = sum_s (g_s - g_bar)^2.
10.4 Phase bookkeeping.
The product of phases factorizes:
prod_{s=1}^r exp(i theta_s m_s)
= exp( i (theta_1 + ... + theta_r) Y_0 )
* prod_{s<r} exp( -i theta_s g_s ).
So the outer phase carries the tiny sum theta_1 + ... + theta_r, and the
inner phases live in the differences g_s. The inner phases become a mask
of level j. This is the only place where the relaxation of the angle
budget pays off.
========================================================================
11. CENTER FOURIER EXPANSION: AN ADMISSIBLE MASK
========================================================================
The Gaussian weight is not yet a mask in the sense of (M1) and (M2). We
convert it.
Momentarily vary Y_0 as a real variable and define
W_g(Y_0) = exp(t Y_0) * prod_{s=1}^r v(Y_0 - g_s).
CLAIM 11.1. There is a centred Gaussian probability nu(d omega) with
variance r / ((1 - zeta) sigma^2), where zeta = 10^{-5}, such that
W_g(Y_0) = exp(1.203 B) int_R exp(i omega Y_0) K(g, omega) nu(d omega),
where
K(g, omega) = exp(-1.203 B) * W_g_hat(omega) / (d nu / d omega),
and W_g_hat is the Fourier transform
W_g_hat(omega) = (2 pi)^{-1} int_R W_g(Y_0) exp(-i omega Y_0) dY_0.
Moreover, for all real omega and all g_s >= 0,
|K(g, omega)| <= exp( - epsilon t (g_1 + ... + g_r) ). (11.1)
Sketch of proof. Scale:
z = (Y_0 - g_bar) sqrt(r) / sigma
rho = t sigma / sqrt(r)
u' = omega sigma / sqrt(r)
d_s = (g_s - g_bar) / sigma
p_ell(x) = sum_{k=0}^{ell} (x^2/2)^k / k!.
Shifting the contour by -i u', which is legal because the integrand is a
polynomial times a Gaussian, gives
W_g_hat = (2 pi)^{-1/2}
* exp( (t - i omega) g_bar - Q(g)/(2 sigma^2)
+ rho^2/2 - (u')^2/2 - i rho u' )
* E_{X ~ N(rho,1)} prod_s p_ell( (X - i u')/sqrt(r) - d_s ).
The product of polynomials is at most zeta^{-r ell} times
exp( (zeta/2)( X^2 + (u')^2 + Q(g)/sigma^2 ) ). Integrating X gives
|W_g_hat(omega)| / (density of nu at omega)
<= zeta^{-r ell}/(1 - zeta)
* exp( t g_bar - (1 - zeta) Q(g)/(2 sigma^2)
+ 2.4 B / (2 r (1 - zeta)) ).
Now Q(g) >= r g_bar^2/(r-1), because g_r = 0. So the exponent above,
plus the epsilon t sum_s g_s term, which equals epsilon t r g_bar, is at
most
(1 + r epsilon)^2 (2.4 B) (r-1) / (2 (1 - zeta) r).
For r = 4 this is at most .9001 B. Add 2.4B/(2r(1-zeta)) <= .3001 B and
the negligible log term; the total is below 1.203 B. Hence (11.1).
Interpretation. K is a Gaussian bump in g. For now it decays like
exp(-epsilon t |g|_1), which is the (M1) behaviour, before we make it
compact.