Trace reconstruction from deletion traces
A complete and readable account of the quasipolynomial sample bound







========================================================================
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.06,
(1 + log log(2n))^6 is about 4500, and the product is about 4.0 * 10^7.
With p^{-1} = 2 included, the theorem promises fewer than
exp(8 * 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.

(M1) and (M2) in plain terms.

  (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.

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.

Proof of the deficit budget, in two parts.

Continuous bound. Since J_diamond is a polynomial, log|J_diamond| is
subharmonic on the unit disk. Evaluating at the interior point z_0 (push
the disk radius just past exp(-a) and let it return; capping at 100E can
only lower an average) and using the pivot bound gives

    Poisson average of log( 2 / |J_diamond| )  <=  log 2 + 1.0006 E,

which is at most 1.003 E for E >= 300. Under theta = a x, the Poisson
kernel at z_0 has density at least (1 - 2a)/(pi (1 + (x - v_0^*)^2)),
so the continuous integral dominates the Riemann sum that produces
G_cal + Q.

Discretization. By layer cake, f(x) = int_0^{infinity} 1_{f(x) > c} dc
for f(x) = min(100, E^{-1} log(2/|J_diamond(exp(-a + i a x))|)). A
polynomial with H_* index terms meets a fixed positive modulus at most
O(H_*) distinct values of x. So each level set f > c is a union of
O(H_*) intervals. Replacing each interval's length by the grid sum of
mesh h adds O(H_* h) = o(1) uniformly in n.

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.21, 7.23, 9.31, 13.25. Each exceeds the
corresponding maximum (4.96, 5.79, 6.94, 8.26, 9.71).

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)

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.


========================================================================
12. RADIAL COMPACTIFICATION
========================================================================

12.1 Applying the expansion.

Each fixed-omega discrepancy contribution, after the expansion, is a
merged sum with new outer factor

    exp( -t Y_0 + i (omega + theta_1 + ... + theta_r) Y_0 )

and inner phases exp(-i theta_s g_s), plus K(g, omega). Its modulus is
at most 2 Z_{j+1}/(1 - exp(-t)) <= exp(10^{-8} B), apart from the outer
factor exp(1.203 B), because every merged loose edge has coverage and
(11.1) applies. Latest offsets Y_0 < 0 cancel exactly, leaving Y_0 >= 0
and hence a geometric sum from the outer damping.

Where |omega + theta_1 + ... + theta_r| > 5t, the Gaussian factor is at
most

    exp(1.2031 B) * 2 exp( -(4.99)^2 (1 - zeta) 2.4 B / 8 )
       = o( exp(-6B) ),

uniformly over categories, after multiplying by P = exp(o(B)). So only
the range |omega + theta_1 + ... + theta_r| <= 5t matters.

12.2 Making the mask compact.

We must replace K(g, omega), for the retained omega, by an admissible
K^c(g, omega) with

    | K(g) - K^c(g) | <= exp(-20B) exp( -epsilon t sum_s g_s ).    (12.1)

Everything is uniform on |omega| <= 5.01 t and measurable in omega.

Fix omega. Put R = R_j; recall R/t = (j+1)^2 M_j/L_j >= 1. Choose a
large absolute constant H_0 and set

    W    = H_0 B / t
    tau^2 = W / R
    U(x) = (1/(tau sqrt(2 pi)))
           * int_{-W}^{W} exp( -(x - v)^2/(2 tau^2) ) dv.

So U is a smoothed indicator of the interval [-W, W].

CLAIM 12.1.

    | K(g) - K(g) prod_{s<r} U(g_s) |
       <= (1/4) exp(-20B) exp( -epsilon t sum_s g_s )    for g_s >= 0.

Proof. If all g_s <= W/2, use

    1 - U(g_s) <= 2 exp( -W^2/(8 tau^2) ) = 2 exp(-R W/8)
               = 2 exp(-H_0 B/8),

and then (11.1). If some g_s > W/2, use Q(g) >= W^2/(4r), by Cauchy
together with g_r = 0. That makes the Gaussian factor in the derivation
of (11.1) at most exp(-c W^2/sigma^2) = exp(-c H_0^2 B/2.4); and
0 <= U <= 1.

Now form the entire function of d' = r-1 variables

    J_cal(g) = exp( R sum_s g_s ) K(g) prod_{s<r} U(g_s),

and its inverse transform

    J_cal_hat(lambda) = (2 pi)^{-d'}
        int_{R^{d'}} J_cal(g) exp(-i lambda . g) dg.

CLAIM 12.2. There are absolute constants C_2, C_3, allowed to depend on
H_0, with

    int |J_cal_hat(lambda)| dlambda <= exp( C_2 R W ),            (12.2a)

    int_{|lambda|_inf > C_3 R} |J_cal_hat(lambda)| dlambda
       <= (1/4) exp(-20B).                                        (12.2b)

Proof. Continue the integral formula for K to complex
g = x - i y, with x and y real vectors in R^{d'} and the last coordinate
set to zero. Using Re Q(x - i y) = Q(x) - Q(y), one gets

    |K(x - i y)| <= exp( C_4 B + C_4 t sum_{s<r} (|x_s| + |y_s|)
                         + C_4 |y|_2^2/sigma^2 ).

Also

    |U(x_s - i y_s)| <= (2W/(tau sqrt(2 pi)))
        * exp( ( y_s^2 - max(|x_s| - W, 0)^2 ) / (2 tau^2) ).

Shifting g = x - i y in J_cal_hat by y, which is legal because of
Gaussian tail decay in every real direction, and using

    max(|x_s| - W, 0)^2 >= x_s^2/2 - W^2

to complete squares, gives

    |J_cal_hat(lambda)| <= tau^{d'}
        exp( C_5 R W + C_5 |y|_2^2 / tau^2 - lambda . y ).

Take y = tau^2 lambda/(2 C_5). This leaves a Gaussian envelope in
lambda. Integrating in lambda proves (12.2a) and (12.2b).

Finally define

    K^c(g) = exp( -R sum_s g_s )
        int_{|lambda|_inf <= C_3 R} exp( i lambda . g )
        J_cal_hat(lambda) dlambda.

Fourier inversion, valid because both J_cal and its transform are
integrable, together with Claim 12.1 and (12.2b), gives (12.1).
Combining (11.1) and (12.1) shows K^c also obeys the quadrant decay
|K^c| <= (some absolute constant) exp(-epsilon t sum g_s), and (12.2b)
gives the complex-measure representation with |xi|_inf <= C R_j and
total variation at most exp(C (R_j/t_j) B_j), after enlarging C once.

This completes the construction of stage j+1 from stage j.


========================================================================
13. STOPPING: FROM MASKS TO PHASE PARAMETERS
========================================================================

We now have a_j n <= .001, j <= b, and the pivot point
z_0 = exp(-a_j + i theta).

13.1 One frequency choice.

Expand all retained masks by their compact representations (M2). The
total variation of the product measure is at most

    exp( C sum_{h<j} 4^{j-h-1} (R_h/t_h) B_h )
       <= exp( O( ell_* 8^b (b+1)^6 ) ).                           (13.1)

Check. We have R_h/t_h = (h+1)^2 2^{b-h} and
B_h = 10^100 ell_* 8^h (b+1)^3, so

    4^{j-h-1} (R_h/t_h) B_h
       <= 4^{j-1} 2^b (h+1)^2 10^100 ell_* (b+1)^3,

and summing over h <= j <= b gives 8^b (b+1)^3 times a polynomial in b.

At every set of frequencies in the product integral, the difference
contribution is still a finite sum. So by the pivot magnitude and
(13.1), some frequency choice has an unintegrated discrepancy of
magnitude at least

    exp( -O( ell_* 8^b (b+1)^6 ) ).

Fix it.

13.2 Reading off a pure product probe.

A mask now contributes, for each of its differences g_s, a number

    exp( -R_h g_s + i (beta_s + xi_s) g_s ),

with |xi_s| <= C R_h and |beta_s| <= 10^4 a_h.

Assign these factors edgewise. Each g_s is a sum of consecutive gap
distances in the merged tuple, so each factor becomes a product of
per-edge factors. Drop factors on exact edges; they are scalars of
modulus at most 1 and nonzero. What remains is a pure product probe:
bit predicates, exact constraints, and, on each loose edge e, a complex
number

    z_e^* = exp( -eta_e + i phi_e ).

CLAIM 13.1. There is an absolute c_0 in (0,1] with

    kappa(z_e^*) := Re ( (1 + z_e^*)/(1 - z_e^*) )
                 = (1 - |z_e^*|^2) / |1 - z_e^*|^2
                 >= c_0 q.                                       (13.2)

Proof. Let the sum over covering differences be counted with
multiplicity. Because every loose edge is covered, eta_e, the sum of the
relevant R_h, is strictly positive. The frequency perturbations sum to
at most C eta_e in absolute value, and the beta's satisfy

    (sum beta_s)^2 <= eta_e sum_h (10^4 a_h)^2/R_h
                    <= C_6 (eta_e/q) sum_h 1/(h+1)^2
                    =  O( eta_e/q ),

using at most 3 * 4^{j-h-1} <= M_h coverings per level h. When
eta_e <= 1, we have

    1 - |z_e|^2 = 1 - exp(-2 eta_e) >= (1 - e^{-2}) eta_e

and |1 - z_e|^2 <= C_7 eta_e/q, because

    |1 - z|^2 = 1 - 2 exp(-eta) cos phi + exp(-2 eta)
              ~ eta^2 + phi^2 <= C (eta^2 + eta/q).

When eta_e > 1, use |1 - z_e| <= 2 and 1 - |z_e|^2 >= 1 - e^{-2}.

For the root itself, kappa(z_0) >= 4/p. Because |1 - z_0| <= 6 a_j and
1 - |z_0|^2 >= a_j, and a_j <= .001/n, we get

    1/(36 a_j) >= 1000 n/36 >= 4 n >= 4/p,

since p >= 1/n.


========================================================================
14. INTERPOLATION AND THE TRACE STATISTIC
========================================================================

14.1 A one-parameter family of edge arguments.

Put

    S = 4/(p c_0),

and let varsigma be a complex variable in the strip
1/2 <= Re varsigma <= S.

For each loose edge e, write

    alpha_e + i delta_e = (1 + z_e^*)/(1 - z_e^*).

By (13.2), alpha_e = kappa(z_e^*) >= c_0 q. For each exact edge, set the
target z_e^* = 0, i.e. alpha_e = 1 and delta_e = 0. Define

    z_e(varsigma) = (alpha_e varsigma + i delta_e - 1)
                    / (alpha_e varsigma + i delta_e + 1).

Then z_e(1) = z_e^* and

    kappa( z_e(varsigma) ) = alpha_e Re varsigma.

At Re varsigma = S this is at least c_0 q S = 4 q/p.

14.2 The interpolating function.

Let K_ex = k - A be the number of exact edges. Define G(varsigma) as the
difference sum at z_0 with the following substitutions:

    each exact-edge selector on source distance Delta becomes
    z_e(varsigma)^{Delta - 1}, so Delta = 1 gives 1, as before;

    each loose edge at source distance Delta uses z_e(varsigma)^Delta;

    the whole sum is normalized by n^{-A} (1 + varsigma)^{-2 K_ex};

    all negative latest-offset terms are cancelled.

So at an exact edge the weight relaxes from a strict constraint to a
geometric sum over possible source distances, and the value 1 of that
sum at varsigma = 1 recovers the original constraint.

CLAIM 14.1.

  (a) |G(varsigma)| <= 2 on the entire strip.

  (b) At varsigma = 1,

      log |G(1)| >= -A log n - 2 K_ex log 2
                    - O( ell_* 8^b (b+1)^6 )
                  = -O( ell_* 8^b (b+1)^6 ).                     (14.1)

Proof of (a). For an exact edge, the sum over Delta >= 1 of
|z_e|^{Delta - 1} equals 1/(1 - |z_e|). Since alpha = 1 and delta = 0,

    1 - |z_e|^2 = 4 Re varsigma / |1 + varsigma|^2,

so this costs at most 1 against the normalization |1 + varsigma|^{-2},
per exact edge. For the other edges, use at most n choices per gap and
modulus at most 1.

Proof of (b). Use k <= 2 P_0 4^b = O(8^b), A <= 4^b, and
log n <= ell_* 2^b.

14.3 Hadamard three-lines.

Apply the three-lines theorem on the vertical strip
1/2 <= Re varsigma <= S, with bounds |G| <= 2 on both vertical edges.
The right-line weight of the point Re varsigma = 1 is
(1/2)/(S - 1/2). So there is a point varsigma of the strip with
Re varsigma = S and

    |G(varsigma)| >= exp( -O( p^{-1} ell_* 8^b (b+1)^6 ) ).        (14.2)

14.4 The disk of radius p centred at q.

The following algebraic fact is the bridge from analytic parameters to
actual traces.

  FACT (product ratio). For u and v in the open unit disk with
  kappa(u), kappa(v) > 0,

      |1 - u v|^2 / (1 - |u v|^2)
         <= |1 - u|^2/(1 - |u|^2) + |1 - v|^2/(1 - |v|^2).        (14.3)

  Proof: 1 - uv = (1 - u) + u(1 - v), then the weighted
  Cauchy-Schwarz inequality with weights (1 - |u|^2) and
  |u|^2 (1 - |v|^2), whose sum is exactly 1 - |uv|^2.

Apply (14.3) with u = z_0, which has kappa >= 4/p, and v =
z_e(varsigma), which has kappa >= 4q/p. The right side is at most

    p/4 + p/(4q) <= p/(2q).

So, for y_e = z_0 z_e(varsigma),

    |1 - y_e|^2 / (1 - |y_e|^2) <= p/(2q).

Use the identity

    p^2 - |y - q|^2 = p(1 - |y|^2) - q |1 - y|^2.

Combining with the previous line gives, for y = y_e,

    p^2 - |y_e - q|^2 >= (p/2)(1 - |y_e|^2)
                       >= (p/2)(1 - |z_e|^2).

The same computation applies to z_0 itself.

So z_0 and y_1,...,y_{k-1} all lie in the closed disk of centre q and
radius p. Define

    lambda_e = (y_e - q)/p,       lambda_0 = (z_0 - q)/p.

All have modulus at most 1. Moreover

    1 - |lambda_e| >= (1 - |z_e|^2)/(4p).                         (14.4)

Proof: p - |y_e - q| >= (p^2 - |y_e - q|^2)/(2p) >= (1 - |z_e|^2)/4,
then divide by p.

14.5 The trace statistic.

Given a trace U of the extended input, define

    T(U) = n^{-A} (1 + varsigma)^{-2 K_ex}
           * sum over 0 <= s_1 < ... < s_k < |U|
             of  lambda_0^{s_1}
                 * prod_{e=1}^{k-1} lambda_e^{ s_{e+1} - s_e - 1 }
                 * prod_{e=1}^{k} chi_e( U_{s_e} ).

CLAIM 14.2. |T(U)| <= 1 for every trace.

Proof. Expand the sum over the exact-edge gaps. For an exact edge the
geometric weight over possible source distances Delta is
sum_{Delta >= 1} |lambda_e|^{Delta - 1} = 1/(1 - |lambda_e|). By (14.4)
together with 1 - |z_e|^2 = 4 S/|1 + varsigma|^2,

    1 - |lambda_e| >= (1 - |z_e|^2)/(4p) = S/(p |1 + varsigma|^2).

So the geometric sum is at most p |1 + varsigma|^2/S. After its share
|1 + varsigma|^{-2} of the normalization, each exact edge contributes at
most

    p/S = p^2 c_0/4 <= 1.

The loose edges, of which there are A-1, contribute at most n each,
because there are at most n possible gap lengths. The first index s_1
ranges over at most n values. Together with the normalization n^{-A},
the choices of the A block starts cost at most 1.

The bit predicates are bounded by 1.

14.6 The expectation gap.

Condition on the survival of the source symbols at the tuple positions.
Each survives independently with probability p. If all k chosen positions
survive, the symbols observed there are exactly x's bits at those
positions. The counts of survivors before s_1 and between consecutive
s_e have independent binomial laws.

By linearity of expectation, and the identity

    q + p lambda_0 = z_0,
    q + p lambda_e = z_0 z_e        for each loose edge e,

together with the fact that at an exact edge the exponent is zero (there
are no survivors strictly between two forced-consecutive positions, so
both sides contribute 1), we get

    E_x T - E_y T = p^k z_0^{d - k + 1}
                    * ( prod over loose e of z_e(varsigma)^{-1} )
                    * G(varsigma).

The powers of z_0 and z_e reorganize the inner sum into exactly
G(varsigma) times an overall factor, via the elementary identity

    z_0^{i_1} prod_{e=1}^{k-1} (z_0 z_e)^{i_{e+1}-i_e-1}
       = z_0^{i_k-k+1} prod_{e=1}^{k-1} z_e^{i_{e+1}-i_e-1},

which shows where each power of z_0, each power of a z_e, and each
factor of G(varsigma) comes from.

Modulus of the prefactor. We have

  |z_0|^{d - k + 1} >= exp(-a_j n) >= exp(-.001),
  because a_j n <= .001 and d - k + 1 <= n;

  |prod over loose e of z_e^{-1}| >= 1,
  because each |z_e| < 1 at Re varsigma = S, and |z_e| > 0 by (14.2);

  p^k >= exp(-k log(1/p)) >= exp(-O(8^b/p)),
  because k = O(8^b) and log(1/p) <= 3/p.

So, by (14.2),

    |E_x T - E_y T| >= exp( -O( p^{-1} ell_* 8^b (b+1)^6 ) ).     (14.5)


========================================================================
15. COUNTING TRACES
========================================================================

We now have a bounded statistic T with an expectation gap (14.5). That
is an information statement. To turn it into a reconstruction rule, use
Hoeffding and a union bound.

15.1 One pair.

Fix x and y. Let delta = |E_x T - E_y T|. By Claim 14.2, T takes
values in [0,1] after the harmless shift T -> (T + 1)/2, which does not
change the gap.

Hoeffding's inequality: with N independent traces, the empirical
average of T deviates from its mean by delta/2 with probability at most
2 exp(-N delta^2/2). So

    N = O( delta^{-2} log(1/eps) ).

With eps = e^{-2n} this is

    exp( O( p^{-1} ell_* 8^b (b+1)^6 ) + log(2n) ),

and since the first term grows like a fixed power of log n while
log(2n) does not, the whole exponent is O( p^{-1} ell_* 8^b (b+1)^6 ).
So this many traces distinguish x from y with error at most e^{-2n} per
pair, uniformly over pairs.

15.2 All pairs at once.

There are 2^n binary strings of the doubled length. By the union bound
over all 2^n challengers y for the true x,

    total error <= 2^n * e^{-2n} = (2/e^2)^n < 1/3       for n >= 1.

So that many traces suffice to recover x from its traces with
probability at least 2/3, simultaneously for every x.

15.3 The final bound.

Substitute the constants:

    N = exp( O( p^{-1} ell_* 8^b (b+1)^6 ) ).

Since 2^b <= 2X and ell_* X <= 1 + 2 log n, this is

    exp( O( p^{-1} (log n) X^2 (1 + log X)^6 ) ),

which is Theorem 2, and in turn

    exp( O( p^{-1} (log n)^3 (1 + log log(2n))^6 ) ),

which is Theorem 1.

For completeness, the steps in detail. We used the doubled length n
throughout. Switching back to the original length changes the exponent
by at most an absolute factor, since doubling n only changes log n by
an additive constant.

The exponent is p^{-1} ell_* 8^b (b+1)^6 times an absolute constant.
We have

    2^b <= 2 X              (by the definition of X),
    ell_* X <= 1 + 2 log n  (since ell_* = 1 + log(1/q) and q >= 1/n).

So

    p^{-1} ell_* 8^b (b+1)^6
       <= p^{-1} ell_* (2X)^3 (1 + log_2(1 + log n/ell_*) + 1)^6
       <= O( p^{-1} ell_* X^3 (1 + log X)^6 )
       <= O( p^{-1} (log n) X^2 (1 + log X)^6 ),

because ell_* X <= 3 log n for n >= 3. This is Theorem 2.

The cruder form, Theorem 1, follows because X <= 1 + log n, as
ell_* >= 1. Then X^2 (1 + log X)^6 <= (log n)^2 (1 + log log(2n))^6
times an absolute constant, and multiplying by log n gives
(log n)^3 (1 + log log(2n))^6.

The polynomial case, Theorem 3. Suppose q <= n^{-epsilon} for fixed
epsilon > 0. If q = 0, one trace suffices. Otherwise
ell_* >= 1 + epsilon log n, so X <= 1 + 1/epsilon while p >= 1/2 for
large n. Hence the exponent is at most O(log n), and polynomially many
traces suffice.

Two honest remarks about this counting.

  The bound counts how many traces suffice in principle. The rule that
  uses them is allowed to be infeasible to compute: it may need to
  search over all candidates, evaluate the test on each pair, and so
  on. That is consistent with the statement of the theorem, which only
  speaks about information.

  An equivalent, more elementary way to say the same thing: the
  expectation gap of a bounded real test implies, by total variation, a
  gap at least half as large for an optimal indicator test. The
  optimal indicator test for a binary pair is simply the one that
  counts how many times each trace appears as a subsequence of each
  candidate. Each trace word's probability under x is its embedding
  count as a subsequence of x, times a length-dependent factor that is
  the same for x and y. So the sample complexity already contains the
  subsequence-count test, and the bounded statistic is only an
  analytical device for proving the gap.


========================================================================
16. GENERAL ALPHABETS
========================================================================

Let Sigma be the alphabet, with |Sigma| <= n after the small-case
reduction.

16.1 A preliminary pass to see all symbols.

Collect ceil(p^{-1} log(12n)) preliminary traces. By a union bound over
the n source positions, with probability at least 11/12 every source
symbol appears somewhere in these traces. Indeed a given position
survives none of them with probability

    (1-p)^{ceil(p^{-1} log(12n))} <= exp(-log(12n)) = 1/(12n),

and there are n positions.

Let S be the set of observed symbols. Always |S| <= n.

16.2 Reducing to binary.

Condition on success.

On independent fresh traces, run the binary test for every ordered pair
of candidate strings in S^n. For a pair (a,b) whose first disagreement
is at position i, map symbols to bits by the indicator of the symbol
a_i, that is, send a_i to 0 and every other symbol to 1 (or do the
symmetric thing). Then a and b project to two distinct binary strings,
and the binary pair test distinguishes those projections, hence a from
b.

The Hoeffding sample size is multiplied by O(log n) so that the
per-pair error is at most e^{-3n log n}. Since there are at most
n^{2n} = e^{2n log n} ordered pairs, the union bound over all pairs gives
total error at most e^{-n log n} < 1/3. So a log n factor on top of the
same sample budget works.

Output a candidate that wins all its comparisons, or declare failure if
there is none.

The remaining cases are the bounded-size ones and the cases
min(p, 1-p) < 1/n; those are handled directly by the full-length-trace
strategy of section 6.1, and the increase does not affect the stated
asymptotic budget.

This extension concerns sample complexity only.


========================================================================
17. WHAT IS AND IS NOT PROVED
========================================================================

A short list, so that nothing here is oversold.

  Sample complexity, not computation. The theorem says: this many traces
  contain enough information. It does not give a fast algorithm. The
  union-bound recovery step explicitly has no assertion of efficiency
  and does not rely on computational checking or search as a proof.

  Not optimality. The true number of traces needed might be much
  smaller. No matching lower bound is claimed.

  Quasipolynomial, not polynomial, for fixed p in (0,1). The bound here
  does not reach polynomial order except when q is polynomially small
  in n. Getting polynomial for fixed p is open.

  Improvement relative to BVW's proved quantitative estimate, not to
  every possible optimization of their approach. The comparison of
  section 3 is against their Theorem 43 as stated, unwound through
  their Section 8.

  Absolute constants are enormous. 10^100 and 10^16 appear in the
  bookkeeping. They are chosen so that many inequalities are obviously
  true, not because anybody believes they are tight. A careful audit
  would shrink them. The interesting content is the shape of the
  exponent, not the size of the constants.

  Binary first, general alphabet second. The main proof is for binary
  strings with a known prefix; the general-alphabet extension reduces
  to the binary case after a preliminary sampling step. The surviving
  symbols are assumed to be observed exactly, with no substitution
  errors.


========================================================================
18. GLOSSARY
========================================================================

  trace              the survivors of a deletion channel, in order,
                     with their original positions unknown

  retention p        the probability that any given symbol survives

  trace reconstruction
                     the problem of recovering a string from traces

  sample complexity  the number of traces needed

  worst case         the bound must hold for every possible string

  quasipolynomial    of the form exp((log n)^C) for a fixed C

  Poisson kernel     the harmonic-measure weight on a circle

  Hadamard three-lines
                     log-convexity of the sup norm on vertical lines

  probe              a tuple of indices plus predicates and exact-gap
                     constraints plus masks

  mask               a factor with a bounded-modulus K and a compact
                     measure representation, attached to differences

  anchor             the last index of the tuple, taken relative to
                     the first mismatch of the two candidates

  Fine and Wilf      a word with periods p and q and length at least
                     p + q - gcd(p,q) has period gcd(p,q)

  Hoeffding          an average of N independent bounded variables
                     deviates from its mean with probability at most
                     2 exp(-2 N t^2)

  total variation    of a complex measure: the largest value of the sum
                     of absolute values over a partition

  vertex cover       a set of vertices hitting every edge

  fractional cover   the same with vertex numbers in [0,1] and
                     endpoint sums at least 1


========================================================================
19. A READING GUIDE
========================================================================

If you want the headline result only. Read sections 1, 2, and 17. Skip
everything else.

If you want the shape of the argument. Read sections 5, 6.4, 6.5, 6.6,
7, and 8. That gives the induction, the base case, and the meaning of
the probe.

If you want to understand the phase test. Read section 9 to the end of
9.3, including the table of numerical constants. The table is the only
place where concrete numbers are checked; everything else is
bookkeeping.

If you want to understand the shuffle and the spectral compactness.
Read sections 10, 11, and 12. The key sentence is (M2): a mask is an
average of pure phases with bounded frequencies. Once you believe that,
the induction is just bookkeeping.

If you want to understand the final counting. Read section 14.4, which
says z_0 and all the y_e lie in the disk of centre q and radius p, and
then 14.5 and 15. The disk is what makes the trace statistic bounded;
the expectation gap and Hoeffding do the rest.

If you want to check the arithmetic. Read the table in 9.4 and the
size bound in 8.5. Those are the two places where an actual number
needs to come out right.

What to skip on a first reading. Everything in 12.2 and 12.3. Those
sections make the mask admissible; they do not add intuition.


========================================================================
20. REFERENCES
========================================================================

  Burudgunte, Valiant, Wang. "A polynomial-time algorithm for the
  deletion channel." These are the BVW references in the text.

  Chase. "Separating words and trace reconstruction." The
  period-avoiding word base.

  Chen, De, Lee, Servedio, Sinha. The low-deletion setting, cited as
  CDLS.

  Fine and Wilf. "Uniqueness theorems for periodic functions." The
  period lemma.

  Hoeffding. "Probability inequalities for sums of bounded random
  variables." The concentration used in section 15.

  Rubinstein. "Explicit separating words." The period and separation
  ideas used in section 8.

  For later developments, see the recent work on trace reconstruction
  problem. arXiv:2207.11489.


========================================================================
END
========================================================================