Model and main result

Let \[b=\max\left\{1,\left\lceil \log_2\bigl(1+\max\{t,a_1,\ldots,a_n\}\bigr) \right\rceil\right\}\] be the maximum bit length of the input integers and target. The machine is a word RAM with word length \[ w=\left\lceil4\bigl(n+b+\log_2(n+2)\bigr)\right\rceil. \tag{1}\] Each input integer occupies one word in read-only random-access memory. Unit-cost operations are word reads and writes, indirect addressing, comparisons, addition, subtraction, multiplication, integer quotient and remainder with a nonzero divisor, bitwise Boolean operations, and shifts. Each operation has a constant number of \(w\)-bit inputs and outputs. Arithmetic exceeding these limits is implemented by multiple word operations and charged accordingly. A fresh independent uniform \(w\)-bit random word costs one operation. Randomness is one-way: an earlier random word can be reused only if it has been retained in writable memory. There is no random oracle, uncharged external storage, or separate workspace bound.

The algorithm must be a single finite program for all \(n\) and \(b\), without advice or externally supplied tables depending on \(n\) or the input. All preprocessing, table construction, input accesses, and random sampling are included in its running time.

The probability in the theorem is over the algorithm’s internal random words, for each fixed input. In particular, no distributional assumption is made on the input integers. The running-time bound applies even to executions on which the answer is incorrect. Our intermediate bound has the form \(2^{0.489995n}\operatorname{poly}(n+b+2)\); the strict exponent slack absorbs the polynomial factor separately for every fixed \(c\).

Related approaches

Two other directions explain both the progress on Subset Sum and the difficulty of a worst-case improvement in \(n\). Pseudopolynomial algorithms measure their complexity in the numerical target \(t\). Bellman’s dynamic-programming approach (Bellman 1957) yields \(O(nt)\) time. Koiliaris and Xu (Koiliaris and Xu 2017) improved this bound deterministically, and Bringmann (Bringmann 2017) obtained a randomized \(\widetilde O(n+t)\) algorithm. Chan (Chan 2026) subsequently obtained a deterministic near-linear pseudopolynomial algorithm. These bounds are particularly effective for small targets; polynomial bit length, however, permits much larger targets. Here and in historical comparisons, \(\widetilde O\) suppresses polylogarithmic factors.

Representation methods provide a different kind of saving. The method introduced by Howgrave-Graham and Joux (Howgrave-Graham and Joux 2010), and developed by Becker, Coron, and Joux (Becker et al. 2011), represents one solution by many pairs of partial objects and then searches a filtered fraction of those representations. Their improvements concern random instances, not a guarantee for every input. Austrin, Kaski, Koivisto, and Nederlof (Austrin et al. 2016) developed algorithms exploiting input blocks with many distinct sums through representations, or blocks with few distinct sums through faster enumeration of those sums. Nederlof and Węgrzycki (Nederlof and Węgrzycki 2021) used related ideas to improve the space bound below the Schroeppel–Shamir exponent.

Related work concerns Subset Balancing: given integers \(v_i\) and a fixed finite coefficient set \(C\subset\mathbb Z\), find a nonzero \(x\in C^n\) satisfying \(\sum_i v_ix_i=0\). Randolph and Węgrzycki (Randolph and Węgrzycki 2026) give worst-case exponent savings for several coefficient sets, though not for ordinary Subset Sum. Their compatibility certificates recover valid representation pairs without inspecting all invalid pairs. Our alias extraction uses the representation and distinct-sum perspective for a different task: listing modular false positives so that their contributions can be removed from an estimated checksum.

Proof overview

The proof separates two complementary tasks: estimating a modular target bin and identifying the false positives created by the modulus. First, random tie-breaking weights transform the original problem into polynomially many targets for new positive integer weights \(c_1,\ldots,c_n\). The isolation lemma of Mulmuley, Vazirani, and Vazirani (Mulmuley et al. 1987) ensures that, with high constant probability, some transformed target has exactly one solution whenever the original instance is positive. A preliminary meet-in-the-middle test handles inputs whose partial subset sums have sufficiently small support. Both reductions are developed in Section 2.

For a prime \(P\), write \(e_P(x)=\exp(2\pi\mathrm i x/P)\). With independent uniform residues \(\phi_1,\ldots,\phi_n\) modulo \(P\), assign a subset \(I\) the phase \[\chi(I)=e_P\left(\sum_{i\in I}\phi_i\right).\] For each transformed target \(t'\), a subset is an alias if its sum is \(t'\) modulo \(P\) but not over the integers. The identity \[ \sum_{\sum_{i\in I}c_i=t'}\chi(I) =\sum_{\sum_{i\in I}c_i\equiv t'\pmod P}\chi(I) -\sum_{\substack{I\text{ an alias}\\\text{for }t'}}\chi(I) \tag{2}\] is the organizing principle. On a negative instance the left side is zero for every transformed target; an isolated positive target leaves a single phase of modulus one. The two sums on the right can be much larger, so their difference must be recovered to small absolute error. A relative approximation to the modular sum would not suffice. We choose a random prime modulus larger than the desired running time, so neither scanning all residues nor evaluating every Fourier frequency is affordable. We estimate the modular sum and subtract the phases of a complete list of aliases.

Isolating one exact target

The following application of the Isolation Lemma of Mulmuley, Vazirani, and Vazirani (Mulmuley et al. 1987) allows a complex-valued checksum to detect a solution without an anti-concentration argument. Choose a power of two \(N_0\) with \(12n\le N_0<24n\), and draw independent uniform integers \(u_i\in\{1,\ldots,N_0\}\). Define \[ C=1+nN_0,\qquad c_i=Ca_i+u_i,\qquad t_s=Ct+s\quad(0\le s\le nN_0). \tag{3}\] The transformed weights and targets have \(O(b+\log n)\) bits, and there are \(nN_0+1=O(n^2)\) targets.

Deficiency and the preliminary search

For the rest of this section, fix the transformed weights \(c_i\). For \(V\subseteq[n]\) define its subset-sum support and its deficiency by \[ \Sigma(V)=\{c(I):I\subseteq V\},\qquad d(V)=|V|-\log|\Sigma(V)|. \tag{4}\] Thus \(d(V)\) is the exponent saved by storing distinct subset sums instead of all \(2^{|V|}\) subsets: \(|\Sigma(V)|=2^{|V|-d(V)}\). Large deficiency means that many subsets have the same sum. These definitions make sense for arbitrary integer weights, not only the positive transformed weights.

We now use distinct sums in the meet-in-the-middle algorithm of Horowitz and Sahni (Horowitz and Sahni 1974). Related distinctions between small subset-sum supports and the complementary regime appear in earlier Subset Sum algorithms (Austrin et al. 2016; Nederlof and Węgrzycki 2021). The procedure can stop without deciding the instance; that outcome means only that the later algorithm must be run.

The expectation in Proposition 5 is used only in the analysis; the algorithm does not evaluate it. For each fixed choice of the tie weights, either the preliminary search decides with probability \(1-O(1/n)\), or \[ \mathbb E_{\substack{V\subseteq[n]\text{ uniform}\\ |V|=\lfloor n/2\rfloor}}d(V) \le0.01004n. \tag{5}\] The alias-recovery argument will use this latter inequality. A CONTINUE outcome is not a certificate that it holds.

A collision bound for random primes

Reduction modulo a random prime makes any fixed false integer equality unlikely. The following elementary bound is the only number-theoretic estimate needed for that purpose.

Frequency patterns and block tables

Character orthogonality gives the exact identity \[ V_\rho=\frac1P\sum_{k=0}^{P-1}J_\rho(k),\qquad J_\rho(k)=e_P(-k\rho)\prod_{i=1}^n(1+e_P(x_i(k))),\qquad x_i(k)=\phi_i+kb_i\pmod P. \tag{15}\] Indeed, expanding the product leaves one term for every subset, and its average character is one exactly when its sum is \(\rho\) modulo \(P\). In particular, \[ |J_\rho(k)|^2=\prod_{i=1}^n M(x_i(k))^2, \tag{16}\] where residues are identified with their angles in the finite circle. All functions in this section are the finite-circle functions of Lemma 8, with their finite normalizations.

Partition \([n]\) deterministically into \(q=\Theta(\sqrt n)\) blocks \(B_1,\ldots,B_q\) of sizes \(d_j=\Theta(\sqrt n)\). A frequency’s pattern records, for each block, rational approximations \(D_j,A_j,E_j\) to the block averages of \(D,A,E\) at its coordinates \(x_i(k)\), each with error at most \(\epsilon\). We choose these approximations so that \[ D_j\ge0,\qquad A_j\le D_j,\qquad E_j\le D_j. \tag{17}\] Here and throughout, the approximation rule is fixed and deterministic. For example, approximate each average to error \(\epsilon/8\), round its upper endpoint upward for \(D_j\), and round its lower endpoint downward for \(A_j,E_j\), on a rational grid of step at most \(\epsilon/4\). The exact inequalities \(D\ge0,D\ge A,D\ge E\) then imply (17). Uniform boundedness of these functions puts all the rounded values in a fixed finite grid. Lemma 8 computes these approximations in polynomial time per frequency, without summing over \(P\) points. Fix likewise deterministic approximations to the individual scores \(\log f_*,\log g_*\) with error at most \(\epsilon\).

Let \(\mathcal P\) be the collection of all blockwise grid triples satisfying (17), including patterns that no frequency attains. Its size is \[ Q_n:=|\mathcal P|=2^{O(\sqrt n)}. \tag{18}\] We enumerate this collection, not the \(P\) frequency assignments. The collection and its grid indices are fixed independently of the phases; only a frequency’s membership in a pattern depends on them. Fix \(\pi\in\mathcal P\) and put \[(D_\pi,A_\pi,E_\pi) =\frac1n\sum_{j=1}^q d_j(D_j,A_j,E_j),\qquad U_\pi=A_\pi+E_\pi-D_\pi,\qquad W_\pi=2D_\pi-A_\pi-E_\pi.\] For the moment let \(0<\alpha\le1\) and \(0\le H_\pi\le\eta\) be rational parameters. The parameter \(\alpha\) is the probability of testing a coordinate, and \(H_\pi\) determines the relative lengths of the two digit lists. After defining the tests and their entries, we choose them to keep both lists short. Choose the integer base \(B_\pi=2^{\lceil(\eta-H_\pi)n\rceil}\) and use the digit ranges \[0\le h<\lceil P/B_\pi\rceil, \qquad 0\le\ell<B_\pi.\] Their vectors have coordinates \[y_i(h)=\phi_i+hB_\pi b_i\pmod P, \qquad z_i(\ell)=-\ell b_i\pmod P.\] The pair \((h,\ell)\) represents \(k=hB_\pi+\ell\) uniquely, and \(y_i(h)-z_i(\ell)=x_i(k)\). Pairs representing \(k\ge P\) will be discarded only when contributions are evaluated. The high and low digit ranges have sizes \(2^{H_\pi n+O(1)}\) and \(2^{(\eta-H_\pi)n+O(1)}\), respectively; together they represent \(2^{\eta n+O(1)}\) pairs.

We test each digit vector separately against random block tables. This organization is an asymmetric filtering construction; related filter frameworks also distinguish the two marginal acceptance probabilities from their joint acceptance probability (Christiani 2017). The estimates required here are proved below. For each block \(B_j\), generate \[ R_j=2^{\lceil(\alpha D_j+4\epsilon)d_j\rceil} \tag{19}\] independent tables. At each coordinate of a table independently, its label is \(*\) with probability \(1-\alpha\) and otherwise a uniform \(\theta\in\mathbb Z/P\mathbb Z\) with probability \(\alpha\). Tables in distinct blocks and patterns are independent, and all are independent of the phases. Distinct table indices are retained even when their labels coincide.

A label \(*\) gives score zero on either side. A non-\(*\) label \(\theta\) gives the high and low scores \[\log f_*(y_i(h)-\theta),\qquad \log g_*(z_i(\ell)-\theta),\] respectively. Use the fixed approximations to these scores when testing a vector. It passes a block table if its average approximate score is at least \(\alpha A_j-3\epsilon\) on the high side, or at least \(\alpha E_j-3\epsilon\) on the low side. A whole table is a tuple of one table index from each block. A vector passes it when it passes all its block tables. A table entry is a digit index and a whole table it passes. A match is a high entry and a low entry with the same whole table. Write \(L_\pi\) for the number of matches. These definitions do not restrict the represented frequency to have pattern \(\pi\).

Coverage and list sizes

The parameters must ensure both that each frequency has a match and that its two lists of entries can be constructed within the time bound. The estimates below give the leading exponents for the expected entry counts as \[H_\pi+\alpha(D_\pi-A_\pi),\qquad \eta-H_\pi+\alpha(D_\pi-E_\pi),\] apart from \(O(\epsilon)\) rounding slack and subexponential factors. The digit ranges contribute \(H_\pi\) and \(\eta-H_\pi\); the number of tables and the score thresholds contribute the remaining terms. Their sum is \(\eta+\alpha W_\pi\), so both exponents can be at most \(\tau\) exactly when \(\alpha W_\pi\le2\tau-\eta=\Gamma\). Choose the largest \(\alpha\in(0,1]\) satisfying this inequality, and make the high-side exponent equal to \(\tau\): \[ \alpha_\pi= \begin{cases}1,&W_\pi\le\Gamma,\\ \Gamma/W_\pi,&W_\pi>\Gamma, \end{cases} \qquad H_\pi=\tau-\alpha_\pi(D_\pi-A_\pi). \tag{20}\] Use these parameters in the construction, abbreviating \(\alpha=\alpha_\pi\). Since \(D_\pi-A_\pi,D_\pi-E_\pi\ge0\) and \(\alpha W_\pi\le\Gamma\), \[ \Delta\le H_\pi\le\tau, \qquad \eta-H_\pi+\alpha(D_\pi-E_\pi) =\Delta+\alpha W_\pi\le\tau. \tag{21}\] In particular both digit ranges have size at most \(2^{\tau n+O(1)}\). We now prove the coverage and size estimates for this choice.

Constructing and sampling the matches

The preceding bounds would not be useful if finding the table entries required visiting every whole table. We instead construct the lists block by block. For each digit vector, test it against the \(R_j\) tables of every block and record its passing index sets. If any set is empty, the vector produces no entries. Otherwise, emit the Cartesian product of these sets, with the digit index attached to each tuple. Since all \(D_j\) are uniformly bounded, \[\sum_jR_j=2^{O(\sqrt n)}\operatorname{poly}(n).\] Consequently the work per digit, apart from emitted entries, is \(2^{O(\sqrt n)}\operatorname{poly}(n)\); a nested counter emits each tuple with polynomial overhead. Coincident vectors still retain their different digit indices.

Sort both entry lists lexicographically by table-index tuple. If the high and low group lengths for a tuple are \(a_\nu,b_\nu\), then \[ L_\pi=\sum_\nu a_\nu b_\nu. \tag{27}\] Prefix sums of these products allow a uniformly random match to be chosen without enumerating them: choose a uniform integer below \(L_\pi\), locate its group by binary search, and use quotient and remainder to select one entry on each side. A zero product contributes no interval. The same groups also allow sequential match enumeration.

For any represented frequency \(k=hB_\pi+\ell\), define \(\mu_\pi(k)\) to be its number of matches in this pattern procedure. Its digits are unique, so \[ \mu_\pi(k)=\prod_{j=1}^q \bigl|\{r\in[R_j]:\text{block table }r \text{ accepts both }y(h)\text{ and }z(\ell)\}\bigr|. \tag{28}\] We compute this integer directly by retesting the stored block tables, in \(2^{O(\sqrt n)}\operatorname{poly}(n)\) time. Repeated table labels are counted with their distinct indices, exactly as in the lists. The very same deterministic score approximations are used in list construction and multiplicity computation. Thus the identity is exact, not approximate. Whenever a match is sampled or enumerated, its multiplicity is positive, even on outcomes where coverage fails.

These operations never scan the frequency range or the set of whole tables. They also permit an immediate cap: abort before either entry list exceeds \[ C_n=2^{\lceil(\tau+20\epsilon)n\rceil} \tag{29}\] entries. Once the lists exist, \(L_\pi\le C_n^2\) on every execution, so the exact match count and its prefix sums have \(O(n)\) bits. Large match counts therefore do not force large enumerations.

Direct summation and importance sampling

Let \(\mathcal G\) denote the event in Lemma 10(iii). We first analyze exact arithmetic and ideal discrete uniform draws. Entry and enumeration caps will be imposed after this analysis; on \(\mathcal G\) they do not interrupt it.

For each pattern define its exact contribution \[v_{\pi,\rho}=\frac1P \sum_{\substack{k<P\\\text{pattern}(k)=\pi}}J_\rho(k).\] If \(\alpha U_\pi\ge\Delta\), enumerate every match and give it contribution zero unless \(k<P\) and its deterministic pattern is \(\pi\). For a retained match use \[ \frac{J_\rho(k)}{P\mu_\pi(k)}. \tag{30}\] On \(\mathcal G\), every frequency of this pattern is present, and its \(\mu_\pi(k)\) contributions sum to \(J_\rho(k)/P\). Also (23) gives \(L_\pi\le2^{(\tau+18\epsilon)n}\) in this case. We cap this enumeration at \(C_n\), aborting before visiting a further match.

If \(\alpha U_\pi<\Delta\), take \(N=2^{\lceil\tau n\rceil}\) independent uniform matches with replacement. When \(L_\pi=0\), return zero for this pattern. Otherwise, average the sample values \[ X_{\pi,\rho}= \begin{cases} \displaystyle\frac{L_\pi}{P\mu_\pi(k)}J_\rho(k), &k<P\text{ and }\text{pattern}(k)=\pi,\\[2mm] 0,&\text{otherwise}. \end{cases} \tag{31}\] The inverse-multiplicity correction is the usual principle of counting each object once across overlapping representations (Karp et al. 1989, sec. 5). Here the summands are complex and we need an absolute error bound, so we give the relevant calculation. Conditional on phases and tables in \(\mathcal G\), a frequency of this pattern is sampled with probability \(\mu_\pi(k)/L_\pi\). Hence the sample mean is unbiased for \(v_{\pi,\rho}\) and its mean squared error is at most \[ \frac1N\,\frac{L_\pi}{P^2} \sum_{\substack{k<P\\\text{pattern}(k)=\pi}}|J_\rho(k)|^2. \tag{32}\] Indeed, its second moment before averaging is the same sum with an additional factor \(1/\mu_\pi(k)\le1\), and independent complex sample errors have zero cross terms. If \(L_\pi=0\), coverage implies that there are no frequencies of this pattern, so the zero estimate is exact.

It remains to show that (32) is small. The phases have been conditioned on in that formula; they must not now be treated as independent under the conditional law of \(\mathcal G\). Instead, we bound the formula times \(\mathbf1_{\mathcal G}\) pointwise and only then take expectation under the original phase law.

For a frequency of pattern \(\pi\), let \(\bar u,\bar v\) be the coordinate averages of \(u,v\) at \(x_i(k)\). The pattern errors imply \[ |U_\pi-\bar u|\le3\epsilon, \qquad |W_\pi-\bar v|\le4\epsilon. \tag{33}\] There are two cases within the sampling branch. When \(W_\pi\le\Gamma\), \(\alpha=1\) and \(U_\pi<\Delta\). Thus \[ -\alpha U_\pi \le-\Delta+8.5(\Delta-\bar u)+30\epsilon. \tag{34}\] For example, writing \(\delta=\Delta-U_\pi>0\), the right side minus the left side is at least \(7.5\delta+(30-25.5)\epsilon\). When \(W_\pi>\Gamma\), we have \(\alpha aW_\pi=\Delta\), and the sampling condition implies \(aW_\pi-U_\pi>0\). Since \(\alpha\le1<8.3\), \[\begin{align*} -\alpha U_\pi &=-\Delta+\alpha(aW_\pi-U_\pi)\\ &\le-\Delta+8.3(a\bar v-\bar u)+30\epsilon. \tag{35}\end{align*}\] Here the rounding loss is at most \(8.3(4a+3)\epsilon<30\epsilon\).

Set \[R_1(x)=8.5(\Delta-u(x)),\qquad R_2(x)=8.3(av(x)-u(x)).\] For a fixed sampled pattern, let \(r=1\) or \(2\) according to the preceding case. Combining (23), (16), and (34) or (35) gives, for every \(k<P\), \[ \mathbf1_{\mathcal G}\mathbf1_{\{\text{pattern}(k)=\pi\}} L_\pi|J_\rho(k)|^2 \le2^{(\eta-\Delta+48\epsilon)n} \prod_{i=1}^n\bigl(M(x_i(k))^2\,2^{R_r(x_i(k))}\bigr). \tag{36}\] The nonnegative right side no longer restricts the pattern or the good event. At each fixed \(k\), under the original law the \(x_i(k)\) are independent uniform residues, regardless of the values of \(b_i\). Lemma 8 therefore bounds the expected product by \(2^{0.5169n}\).

Let \(E_{\pi,\rho}\) be the error of this pattern’s ideal estimate. Taking conditional expectation over its samples first in (32), and then using (36), yields \[\begin{align*} \mathbb E\bigl[\mathbf1_{\mathcal G}|E_{\pi,\rho}|^2\bigr] &\le\frac{2^{(\eta-\Delta+48\epsilon+0.5169)n}}{NP}\\ &=2^{(-\tau-\Delta+0.5169+48\epsilon)n+O(1)} \le2^{-0.0007n} \tag{37}\end{align*}\] for sufficiently large \(n\). The exponent before the \(O(1)\) term is \(-0.00081952\). Directly summed patterns have zero error on \(\mathcal G\). Cauchy–Schwarz gives \[\mathbb E\!\left[\mathbf1_{\mathcal G} \left|\sum_{\pi\in\mathcal P}E_{\pi,\rho}\right|^2\right] \le Q_n^2\,2^{-0.0007n}.\] Markov’s inequality and a union bound over \(\mathcal T\) now show that, with probability \(1-o(1)\), \(\mathcal G\) holds and every total ideal estimate has absolute error at most \(1/16\). The subexponential number of patterns and the polynomial number of targets are both absorbed by the negative linear exponent. This conclusion uses no independence between pattern errors or between target estimates.

Resource and precision bounds

We finish the proof of Theorem 9 by making the preceding procedure bounded on all histories. Generate block tables and digit entries as described, aborting at the entry cap (29). For a direct pattern, abort at the same cap on enumerated matches. Sampled patterns use exactly \(N\) match draws, or return zero when there is no match. The event \(\mathcal G\) implies that none of these caps is reached. On any other history, the caps still bound the work.

Each grid value has a fixed rational denominator. Thus the weighted averages and \(\alpha\) have numerators and denominators of \(O(\log n)\) bits, even for unattained patterns. The ceilings in the digit and table counts are computed exactly by rational arithmetic. For each pattern, digit processing costs \(2^{\tau n+O(\sqrt n)}\operatorname{poly}(n)\) before the emitted entries. Sorting, prefix sums, and random indexing have polynomial cost per entry or sample. Each valid contribution requires at most \(2^{O(\sqrt n)}\operatorname{poly}(n)\) work for its multiplicity, and a polynomial factor for the target list. All tuple keys, match counts, and multiplicities have \(O(n)\) bits: for multiplicities use (25), and for match counts on retained lists use \(L_\pi\le C_n^2\). Multiplication by the \(Q_n\) patterns gives (14).

The coarse fixed approximations in pattern and acceptance tests define the lists themselves; they need not be refined later. The complex contributions require greater accuracy. On every nonaborted history there are only \(2^{O(n)}\operatorname{poly}(|\mathcal T|)\) of them, each formed from at most \(n+1\) factors of modulus at most two. The rational weights in (30) and (31) have \(O(n)\)-bit numerators and denominators and magnitude at most \(2^{O(n)}\). Section 6.2, under Complex evaluations and summation, implements precisely such sums by evaluating each phase to error \(2^{-2n^2}\) and rounding each contribution to \(n^2\) fractional bits on a common dyadic grid. Its total error here is \(2^{-n^2+O(n)}\operatorname{poly}(|\mathcal T|)=o(1)\), hence less than \(1/16\) for sufficiently large \(n\), with polynomial work per contribution.

The rational probabilities in the label choices, finite uniform draws, and phase evaluations admit bounded implementations on the specified word RAM; Section 6 gives these implementations and their coupling to the ideal draws. Their additional abort probability is \(o(1)\), and their cost is polynomial per operation just counted. The retained phases and block labels are ordinary writable data; the procedure requires no access to unstored randomness. The ideal error bound \(1/16\), the arithmetic error \(1/16\), and the vanishing probabilities of cap or sampling failure prove Theorem 9.

Representations and filtered lists

Set \[ m=2\left\lfloor\frac{.0567n}{2}\right\rfloor, \qquad \gamma=.48925. \tag{40}\] Thus \(m\) is even. All arguments below concern sufficiently large \(n\), so \(m>0\) and \(17m/2\le\lfloor n/2\rfloor\). Write \(h(x)=-x\log x-(1-x)\log(1-x)\) for binary entropy, with \(h(0)=h(1)=0\). The numerical bounds \(h(.0567)+.0567<.40\) and \(h(1/4)<.81128\) used below follow, for example, from the logarithm series (66), after scaling its argument into \([1,2]\).

Every alias \(I\) of \(t_s\) has a carry \(q\in\{0,\ldots,n\}\) such that \[ \sum_{i\in I}b_i=T_{s,q}, \qquad T_{s,q}=(t_s\bmod P)+qP. \tag{41}\] For each \(s\) and every \(q=0,\ldots,n\), perform \(n^{10}\) independent trials, conditional on \(P\) and the transformed input. A trial chooses a uniform \(m\)-element set \(M\subseteq[n]\) and an independent uniform subset \(F\subseteq[n]\). These two draws have the same law for every \(P\) and use randomness independent of it. The block \(M\) supplies multiple representations: selected positions in it will be split between two masks. For any fixed subset \(I\), the flipped set \(J=I\mathbin\triangle F\) selects exactly half of this block with inverse-polynomial probability, independently of the original cardinality of \(I\). Introduce signed weights and translated targets \[ \begin{aligned} b_i^F&=(-1)^{\mathbf1_{i\in F}}b_i, &T^F&=T_{s,q}-\sum_{i\in F}b_i,\\ c_i^F&=(-1)^{\mathbf1_{i\in F}}c_i, &t_s^F&=t_s-\sum_{i\in F}c_i. \end{aligned} \tag{42}\] For every subset \(I\), these definitions give \[\sum_{i\in J}b_i^F=\sum_{i\in I}b_i-\sum_{i\in F}b_i, \qquad \sum_{i\in J}c_i^F=\sum_{i\in I}c_i-\sum_{i\in F}c_i.\] Thus flipping is a bijection, not a relaxation of the subset problem.

Within the trial, consider every \(k=0,\ldots,m/2\). Partition the positions outside \(M\) into ordinary left and right sets \(L,R\), of sizes \(l,r\), where \(l+r=n-m\). Choose \(l\) so that \[ N_L=2^l\binom mk, \qquad N_R=2^r\binom m{m/2-k} \tag{43}\] are within a factor of two. Such a choice follows by rounding \[l=\frac{n-m+\log\binom m{m/2-k}-\log\binom mk}{2}\] to the nearest integer. The unrounded value lies between \((n-2m)/2\) and \(n/2\), so the required partition exists. The choice can be made by exact comparisons of integers, without logarithmic arithmetic. For definiteness, take the first \(l\) positions of \([n]\setminus M\) as \(L\).

A left record is a pair \((A,X)\) with \(A\subseteq L\), \(X\subseteq M\), \(|X|=k\); a right record is a pair \((B,Y)\) with \(B\subseteq R\), \(Y\subseteq M\), \(|Y|=m/2-k\). Each record stores both its signed reduced weight, obtained using \(b_i^F\), and its signed unreduced weight, obtained using \(c_i^F\). At this point \(X\) and \(Y\) need not be disjoint. If they are disjoint, their union with \(A,B\) represents a subset with exactly \(m/2\) selected positions in \(M\).

For this \(k\), independently draw a prime \[P'\in[2^{\lfloor\gamma m\rfloor},2^{1+\lfloor\gamma m\rfloor})\] uniformly and then draw a uniform residue \(j'\bmod P'\). Keep left records whose signed reduced weight is \(j'\bmod P'\), and right records whose signed reduced weight is \(T^F-j'\bmod P'\). Only these filtered records are formed explicitly.

To generate a filtered list, split its ordinary coordinate set into two parts as evenly as possible. Enumerate subsets of one part together with all its permitted masks in \(M\), and enumerate subsets of the other part. Sorting these two lists by their weights modulo \(P'\) permits their matching pairs to be emitted in time polynomial in \(n+b+2\) times the sum of their lengths and the number of emitted records. Indeed, for each entry in one list, a binary search finds the interval of complementary residues in the other, after which only output records are visited. The setup work for both sides is at most \[ \operatorname{poly}(n+b+2) 2^{m+\lceil\max(l,r)/2\rceil} \le \operatorname{poly}(n+b+2)2^{.3067n+O(1)}. \tag{44}\] Here \(\max(l,r)\le n/2+O(1)\) follows from the balancing rule.

Why the visited lists and pairs are small

The next estimates concern the uncapped construction. They hold without the low-deficiency assumption and average over the original random choice of \(P\). In particular, we do not restrict \(P\) to primes for which the later coverage argument succeeds.

Conditional on all choices except \(j'\), any fixed left or right record survives with probability \(1/P'\). Balancing (43), the entropy bound for binomial coefficients, and concavity of \(h\) give \[\begin{align*} \log\mathbb E[\text{one filtered list length}] &\le \frac{n-m}{2} +\frac m2\left(h(k/m)+h(1/2-k/m)\right) -\gamma m+O(1)\\ &\le \frac n2+\left(h(1/4)-\frac12-\gamma\right)m+O(1) <.48992n \tag{45}\end{align*}\] for sufficiently large \(n\). The last inequality uses \(h(1/4)<.81128\); the coefficient obtained from this upper bound is \(.489909101\).

To bound visited pairs, fix \(s,q,k\), the block \(M\), the flip \(F\), and two raw records. The indices \(s,q,k\) range over predetermined integers; in particular \(q\) is fixed while \(P\) is averaged. These records and their ordinary-coordinate partition can be fixed independently of \(P\). Let \(\mu_i\in\{0,1,2\}\) be their combined multiplicity at position \(i\), and define \[ D=\sum_i(-1)^{\mathbf1_{i\in F}}\mu_i c_i -t_s+\sum_{i\in F}c_i. \tag{46}\] This is a fixed integer independent of \(P\); in particular neither the reduced weights nor the carry target occurs in its definition. The exact-integer exclusion requires \(D\ne0\). The reduced-weight equality implies \(P\mid D\), because \(T_{s,q}\equiv t_s\pmod P\) and \(b_i\equiv c_i\pmod P\). The integer \(D\) has \(O(b+\log n)\) bits, even when masks overlap. Lemma 6 therefore gives \[ \Pr_P[\text{both integer tests hold}] \le\operatorname{poly}(n)(b+\log n)2^{-\eta n}. \tag{47}\] For \(D=0\) the pair is always excluded, so its contribution is zero.

Conditional on \(P,P',M,F\) and the raw records, the smaller-modulus filters can both hold for at most one value of \(j'\). Their probability is therefore at most \(1/P'\le 2^{1-\gamma m}\), regardless of whether the reduced-weight equality holds. Combining this bound with (47) and summing over the \(N_LN_R\) raw pairs yields \[\begin{align*} \mathbb E[\text{visited pair count}] &\le \operatorname{poly}(n)(b+\log n) 2^{n+(2h(1/4)-1-\gamma)m-\eta n}\\ &<2^{.48988n} \tag{48}\end{align*}\] uniformly under (38), for sufficiently large \(n\). Indeed, using \(h(1/4)<.81128\) and the guard, the coefficient of \(n\) before polynomial factors is at most \(.489848677\).

This is the reason to exclude exact equalities before testing mask disjointness. If exact equalities for overlapping pairs were visited, they could have \(D=0\) and would not satisfy the nonzero-divisor estimate. The sorted-interval exclusion avoids their cost entirely.

Cap each filtered list and each visited-pair count at \[ C_{\mathrm{alias}}=2^{\lceil\tau n\rceil}. \tag{49}\] The procedure returns failure before emitting a record that would exceed a cap. There are only polynomially many choices of \(s,q\), trial, and \(k\). Markov’s inequality and (45)–(48) show that a cap is exceeded with probability \(o(1)\), uniformly over guarded inputs. On every execution, including one that fails, list generation, interval searches, mask checks, storage, and deduplication take at most \(2^{\tau n}\operatorname{poly}(n+b+2)\) bit operations, by (44) and the caps.

Why every alias is recovered

It remains to prove that the constructed representations cover every alias. The low-deficiency condition will provide many distinct mask weights on a randomly selected block. We first show that, with high probability, reduction modulo \(P\) preserves all distinctions of this kind.

In particular, fix any block of at most \(m\) positions and any assignment of signs on that block. Distinct integer subset sums of the signed \(c_i\) remain distinct modulo \(P\) on the event of Lemma 12. Differences of two such sums have coefficients in \(\{-1,0,1\}\) on that block. Consequently they also give distinct integer subset sums of the corresponding signed \(b_i\). The event covers all blocks and all sign assignments at once, including those chosen by the random trials.

Assume now the low-deficiency condition (5), and fix a prime satisfying Lemma 12. We continue in the uncapped experiment. Fix a target \(s\) and an alias \(I\in\mathcal A_s(P)\), and use its carry \(q\) in (41). The constants in the following probability bounds do not depend on \(I,s\), or the chosen prime.

Subtracting aliases

Choose a fixed sufficiently large integer \(n_*\), whose requirements will be specified below. If \(n<n_*\) or the bit-length guard (38) fails, use exact meet-in-the-middle on the original instance. Otherwise choose the tie weights and form \(c_i,t_s\) as in (3). Run the preliminary procedure of Proposition 5. Return its answer if it decides the instance; its outcome Continue means only that we proceed to the following steps.

Run the alias procedure of Theorem 11 to choose the large prime \(P\) and obtain its alias lists. Write \(\mathcal A_s\) for its deduplicated set of aliases for target \(t_s\). Every retained subset is checked using the original transformed integers, so even on an unsuccessful extraction run the set \(\mathcal A_s\) contains only genuine aliases. On a successful run it contains all of them. Next draw and retain the phases \(\phi_1,\ldots,\phi_n\) used in Theorem 9, and compute its estimates \(\widehat V_{t_s\bmod P}\) for all target residues. Repeated residues can share an estimate. For each \(s\), form the finite-precision quantity \[ \widehat B_s =\widehat V_{t_s\bmod P} -\sum_{I\in\mathcal A_s} e_P\left(\sum_{i\in I}\phi_i\right). \tag{55}\] Each phase in this subtraction is evaluated from the retained integer residues \(\phi_i\); it is not resampled. Deduplication is by the binary subset itself, separately for each \(s\), not by its sum or phase.

Output Yes if \(|\widehat B_s|\ge 1/2\) for some \(s\), and No otherwise. Since the computed real and imaginary parts are dyadic rationals, this test is implemented exactly by comparing the sum of their squares with \(1/4\). Any failure declared by a capped main procedure or a main-procedure random sampler results in No. This convention does not apply to the preliminary outcome Continue.

The cancellation behind this decision rule is exact before numerical rounding. When extraction is complete, \[ V_{t_s\bmod P} -\sum_{I\in\mathcal A_s}e_P\left(\sum_{i\in I}\phi_i\right) = B_s, \qquad B_s:=\sum_{\substack{I\subseteq[n]\\\sum_{i\in I}c_i=t_s}} e_P\left(\sum_{i\in I}\phi_i\right). \tag{56}\] In particular, an empty exact bin has \(B_s=0\), and a bin containing one subset has \(|B_s|=1\), for every choice of phases. Thus isolation, rather than an anti-concentration estimate for a many-term sum, is what turns checksum estimation into a decision algorithm.

Finite arithmetic and bounded sampling

The preceding sections describe finite combinatorial constructions, but their analyses use exact uniform draws and exact complex values. We now implement both without unbounded running times or real-arithmetic instructions. All constants below are fixed independently of the input and of the exponent \(c\) in the polynomial-bit convention.

Correctness and the uniform running-time bound

Reducing the scores to two convolutions

Use the continuous functions of Section 3. Define \[ \begin{split} \mathcal P(x)&=\int F(x+w)G(w)\,dw,\\ \mathcal Y(x)&=\int F(x+w)G(w)(r(x+w)+s(w))\,dw,\\ L(x)&=\log Q(x),\qquad B(x)=\frac{\mathcal Y(x)}{\mathcal P(x)}-\log(Z_fZ_g),\\ Z_{f,\kappa}&=\int F^\kappa,\qquad Z_{g,\kappa}=\int G^\kappa,\\ H&=\log Z_{f,\kappa}+\log Z_{g,\kappa} -\kappa\log(Z_fZ_g). \end{split} \tag{62}\] Here \(Q=\mathcal P/(Z_fZ_g)\), and direct substitution in the score definitions gives \[ D=B-L,\qquad u=L+(\kappa-1)B-H,\qquad v=(2-\kappa)B-2L+H. \tag{63}\] Indeed, \(B\) is the \(K_x\)-mean of \(\log(f(x+\cdot)g)\), whereas \(A+E=\kappa B-H\). These identities reduce the computation to the two convolutions in (62).

Write \(x_j=2\pi j/64\) and \(T_{64}(h)=64^{-1}\sum_{j=0}^{63}h(x_j)\). A hat means that each integral defining a quantity has been replaced by \(T_{64}\); nonlinear operations are performed after these replacements. In particular \[\widehat{\mathcal P}_j=\frac1{64}\sum_{i=0}^{63}F(x_{i+j})G(x_i), \quad \widehat{\mathcal Y}_j=\frac1{64}\sum_{i=0}^{63} F(x_{i+j})G(x_i)(r(x_{i+j})+s(x_i)),\] where indices are taken modulo 64. Let \(\widehat X_j\) and \(\widehat X'_j\) be the values of the first and second integrands in (11), respectively, obtained from \(\widehat L,\widehat B,\widehat H\) using (63).

Finite arithmetic certificate

Put \(\beta=0.19\), and let \(r_\beta,s_\beta\) be obtained from \(r,s\) by multiplying their \(j\)th cosine coefficients by \(\cosh(\beta j)\). The following enclosures will suffice. In the first four rows the error is at most \(5\cdot10^{-7}\), and in the last two it is at most \(5\cdot10^{-6}\).

Table 1 gives the remaining node enclosures. The errors in its four numerical columns are at most \(4\cdot10^{-6}\), \(2\cdot10^{-5}\), \(0.0004\), and \(0.0004\), respectively. All columns are invariant under \(j\mapsto64-j\). Consequently, averaging the last two columns with weight one at \(j=0,32\) and weight two otherwise, including their stated errors, gives \[ \frac1{64}\sum_{j=0}^{63}\widehat X_j<1.4295, \qquad \frac1{64}\sum_{j=0}^{63}\widehat X'_j<1.4295. \tag{64}\]

We specify elementary rational operations that reproduce these enclosures. The numbers \(\cos x_j\), for \(0\le j\le16\), lie within \(3\cdot10^{-10}\) of the following values, in order: \[\begin{array}{rrrrr} 1&0.9951847267&0.9807852804&0.9569403357&0.9238795325\\ 0.8819212643&0.8314696123&0.7730104534&0.7071067812&0.6343932842\\ 0.5555702330&0.4713967368&0.3826834324&0.2902846773&0.1950903220\\ 0.0980171403&0 \end{array}\] For example, these intervals follow by using the cosine series through degree 24, with remainder at most \(|t|^{26}/26!\), and \[3.1415926535897<\pi<3.1415926535899.\] The latter enclosure follows from \(\pi/4=\arctan(1/2)+\arctan(1/3)\) and the alternating arctangent series through degree 159. Reflection and a half-turn give the remaining cosine intervals. Taking the dot products with the coefficients in (10), with indices reduced modulo 64, gives rational intervals for every \(r(x_j),s(x_j)\).

For all the exponential evaluations one may use \[ \exp(t)\in \sum_{i=0}^{45}\frac{t^i}{i!} +\left[-R_{\exp}(t),R_{\exp}(t)\right],\qquad R_{\exp}(t)=\frac{|t|^{46}}{46!\,(1-|t|/47)}. \tag{65}\] All arguments needed for the finite tables have \(|t|\le8\). For logarithms put \(d=(y-1)/(y+1)\) and use \[ \ln y\in 2\sum_{i=0}^{15}\frac{d^{2i+1}}{2i+1} +\left[-R_{\log}(y),R_{\log}(y)\right],\qquad R_{\log}(y)=\frac{2|d|^{33}}{33(1-d^2)}. \tag{66}\] These follow by bounding the remaining absolutely convergent tails; in particular they give \(0.693147180559<\ln2<0.693147180561\). Compute powers of two as \(\exp(t\ln2)\), and compute \(\cosh t=(\exp t+\exp(-t))/2\). Apply interval addition, multiplication, and division to the displayed finite sums and then to \[\begin{split} \widehat L_j&=\log\widehat{\mathcal P}_j -\log\widehat Z_f-\log\widehat Z_g,\\ \widehat B_j&=\widehat{\mathcal Y}_j/\widehat{\mathcal P}_j -\log\widehat Z_f-\log\widehat Z_g,\\ \widehat H&=\log\widehat Z_{f,\kappa}+\log\widehat Z_{g,\kappa} -\kappa(\log\widehat Z_f+\log\widehat Z_g). \end{split}\] These operations give the stated enclosures. The intermediate intervals may be kept at their rational endpoints; the displayed decimal columns are final enclosures, not instructions to round intermediate calculations. Thus the tables and the two explicit series remainder bounds constitute a finite arithmetic certificate.

Controlling the inner quadratures

We use the following form of the analytic trapezoid estimate (Trefethen and Weideman 2014, Theorem 3.2). If a \(2\pi\)-periodic function \(h\) is analytic on a neighborhood of the closed strip \(|\operatorname{Im}z|\le\sigma\) and \(|h(z)|\le B_0\) there, then \[ \left|T_{64}(h)-\int h\right| \le\frac{2B_0}{e^{64\sigma}-1}. \tag{67}\] Indeed, shifting the contour for its \(k\)th Fourier coefficient toward the appropriate boundary gives \(|\widehat h(k)|\le B_0e^{-|k|\sigma}\). Averaging at the 64 nodes retains precisely the coefficients indexed by multiples of 64; the nonzero multiples have total absolute value at most the right side of (67).

For the inner integrals defining the hatted values, take \(\sigma=0.45\). The inequality \(|\cos(jz)|\le\cosh(j\sigma)\) and the displayed coefficients give \[|r(z)|<6.50,\qquad |s(z)|<8.71, \qquad |2r_\beta(z)|<18.4,\qquad |2s_\beta(z)|<22.4.\] For instance these are obtained by summing the absolute coefficients times \(\cosh(0.45j)\); the two last bounds also include the factor \(2\cosh(0.19j)\). Formula (67) gives the following errors, distinct from the finite-table rounding errors: \[ \begin{array}{c|c} \text{Quantity}&\text{Error caused by quadrature}\\\hline Z_f,Z_g,Z_{f,\kappa},Z_{g,\kappa}&4\cdot10^{-10}\\ \mathcal P(x_j)&3\cdot10^{-8}\\ \mathcal Y(x_j)&5\cdot10^{-7}\\ \int 2^{2r_\beta},\ \int 2^{2s_\beta}&5\cdot10^{-6} \end{array} \tag{68}\] In detail the respective strip bounds can be taken as \(2^{\kappa\,8.71}\), \(2^{6.50+8.71}\), \((6.50+8.71)2^{6.50+8.71}\), and \(2^{22.4}\). Substituting them in (67) proves every row of (68) by elementary rational bounds on the exponential.

Using the ranges in Table 1 to propagate these errors through (62) gives, at every node, \[ |L-\widehat L|<10^{-7},\qquad |B-\widehat B|<2\cdot10^{-6},\qquad |H-\widehat H|<5\cdot10^{-9}. \tag{69}\] For example, division is safe because \(\widehat{\mathcal P}_j>1.06789\), and \(|\log x-\log y|\le |x-y|/(\min(x,y)\ln2)\) controls the logarithms. Applying (63) to (69) shows that either true node integrand is at most \(1.00001\) times its hatted value. This also holds at \(x_{32}=\pi\), where both values are zero. We have therefore proved \[ T_{64}(M^2 2^{8.5(\Delta-u)})<1.4295(1.00001),\qquad T_{64}(M^2 2^{8.3(av-u)})<1.4295(1.00001). \tag{70}\] It remains to control the outer quadrature, whose integrands involve the logarithm and quotient in (62).

An analytic strip for the outer integrands

The first table and (68) give \[ \int|f(w+i\beta)|^2\,dw =Z_f^{-2}\int2^{2r_\beta(w)}\,dw<1.761, \qquad \int|g(w+i\beta)|^2\,dw<1.38. \tag{71}\] The same bounds hold for imaginary shifts of absolute value at most \(\beta\). To see this directly, the functions are real on the real circle, so Parseval’s identity pairs their positive and negative Fourier coefficients into nonnegative multiples of \(\cosh(2k y)\); their squared norms are thus nondecreasing in \(|y|\). Their shifted means remain one by contour shifting. Hence the corresponding squared norms for \(f-1\) and \(g-1\) are less than \(0.761\) and \(0.38\).

For \(z=x+iy\) with \(|y|\le0.38=2\beta\), periodic contour shifting in (6) yields \[Q(z)=\int f(x+w+iy/2)g(w-iy/2)\,dw.\] Subtracting the product of the shifted means and applying Cauchy–Schwarz gives \[ |Q(z)-1|\le\sqrt{0.761\cdot0.38}<0.54. \tag{72}\] In particular \(Q\) has no zeros on this strip, and \(L=\log Q=(\ln Q)/\ln2\) is analytic there with the branch real on the real axis. Moreover \[-\operatorname{Re}L\le-\log0.46<1.121.\] The numerator \(QB\) is also analytic. The same split contour gives, for every \(|y|\le2\beta\), \[\begin{align*} Q(z)B(z) =\int &f(x+w+iy/2)g(w-iy/2)\\ &\mathord{}\cdot \bigl(r(x+w+iy/2)+s(w-iy/2)-\log(Z_fZ_g)\bigr)\,dw. \end{align*}\] This uses the analytic expression \(r+s-\log(Z_fZ_g)\) for the logarithm of the product, not a new branch choice. Each factor is evaluated within the strip of half-width \(\beta\), where the coefficient sums give \(|r|<2.96\) and \(|s|<4.67\). Consequently, on the full strip of half-width \(2\beta\), \[ |Q B|\le\sqrt{1.761\cdot1.38} \bigl(2.96+4.67+\log Z_f+\log Z_g\bigr)<13, \qquad |B|<13/0.46. \tag{73}\] The normalizer enclosures also give \(H<0.016\).

Now (63), (72), and (73) imply \[8.5\operatorname{Re}(\Delta-u) \le8.5\left(0.02776+1.121+0.0315\frac{13}{0.46}+0.016\right) <17.6.\] For the second exponent use the exact identity \[av-u=-(1+2a)L+ \bigl(a(2-\kappa)-(\kappa-1)\bigr)B+(1+a)H.\] The same bounds give \[8.3\operatorname{Re}(av-u) \le8.3\left((1+2a)1.121+ |a(2-\kappa)-(\kappa-1)|\frac{13}{0.46} +(1+a)0.016\right)<17.6.\] Finally \(|2+2\cos z|\le2+2\cosh(0.38)<4.16\). Both outer integrands therefore extend analytically to the strip \(|\operatorname{Im}z|\le0.38\) and have absolute value below \(2^{20}\) there. (The strict bounds also give a neighborhood of the closed strip.) Formula (67) bounds either outer quadrature error by \[\frac{2^{21}}{e^{64\cdot0.38}-1}<0.00007.\] Combining this with (70), each integral in (11) is strictly less than \[1.4295(1.00001)+0.00007 =1.429584295<2^{0.5168}.\] The last comparison follows already from \(2^{0.5168}>1.4307781\), obtained from (65). This proves both filtering moments.