Historical context and method

Metric \(k\)-median has been a testing ground for LP rounding, primal-dual algorithms, and local search. Charikar, Guha, Tardos and Shmoys gave the first constant-factor approximation through LP rounding (Charikar et al. 2002). Jain and Vazirani connected the cardinality constraint with metric facility location through Lagrangian relaxation and the primal-dual method (Jain and Vazirani 2001). The combinatorial improvements of Charikar and Guha (Charikar and Guha 1999) and the greedy dual-fitting approach of Jain et al. (Jain et al. 2003) developed this connection further. A different route, the bounded-swap local search of Arya et al. (Arya et al. 2004), achieves factor \(3+\varepsilon\).

Li and Svensson (Li and Svensson 2016) made constant-additive pseudo-approximations a route past the factor-three barrier: a constant number of extra facilities can be removed at an arbitrarily small additional approximation loss. Subsequent bi-point rounding, including the corrected \(2.675+\varepsilon\) guarantee of Byrka et al. (Byrka et al. 2017), refined this approach. Gowda et al. (Gowda et al. 2023) obtained \(2.613\) in their SODA 2023 version. Cohen-Addad et al. (Cohen-Addad et al. 2026) reached \(2+\varepsilon\), first announced in 2025. Byrka et al. (Byrka et al. 2026) subsequently developed a graph-based iterative rounding framework that also gives \(2+\varepsilon\) for metric \(k\)-median and extends to other clustering objectives. Theorem 1 reaches every fixed factor above the \(1+2/e\) hardness boundary in the general candidate-facility model.

Our construction follows the iterative graph framework, with the cost analysis replaced by a directed comparison inequality and a potential that retains some surviving assignments after an opening. The earlier framework’s general transfer theorem uses a parameter at least \(2\) in the metric \(k\)-median case (Byrka et al. 2026, Theorem 6.1); it does not supply the smaller factor proved here. We establish the grid preparation, physical updates, comparison, and deterministic implementation together.

A separate companion combines anchor recovery with bounded-price payment accounting to give an independent randomized approximation below factor \(2\) while respecting the facility budget (OpenAI 2026, Theorem 1.2).

Two other methodological antecedents are useful to distinguish. Splitting fractional facilities and organizing assignments around separated representatives and bundles are standard LP-rounding techniques; Charikar and Li (Charikar and Li 2012) provide a close predecessor. The pair updates used in our preparation belong to the dependent-rounding approach of Srinivasan (Srinivasan 2001) and Gandhi et al. (Gandhi et al. 2006), with weighted variants studied by Byrka et al. (Byrka et al. 2017). We prove the particular moment bounds and coupled supplier construction needed here. Likewise, reducing the support of a distribution while matching moments is a standard Carathéodory argument; see Bayer and Teichmann (Bayer and Teichmann 2005) for the finite-cubature perspective. Our implementation specifies which moments suffice throughout the physical process and proves the required rational encoding bounds.

Proof overview

The proof separates the preparation of fractional openings, their iterative rounding, and the removal of randomness and surplus. The intermediate guarantee is measured against the standard assignment LP: its variable \(y_i\) is the amount of opening at facility \(i\), and each client receives a unit of this opening mass. The LP has total opening mass at most \(k\) and cost at most \(\mathop{\mathrm{OPT}}_k\).

The fractional solution

We use the standard relaxation \[ \begin{aligned} \text{minimize}\quad &\sum_{j\in J}\sum_{i\in F}d(j,i)x_{ij}\\ \text{subject to}\quad& \sum_{i\in F}x_{ij}=1 &&(j\in J),\\ &0\le x_{ij}\le y_i\le1 &&(i\in F,\ j\in J),\\ &\sum_{i\in F}y_i\le k. \end{aligned} \tag{1}\] Its rational optimum can be found in polynomial time. For fixed openings \(y_i\), assign each client to its nearest unit of opening mass, using a fixed order to break ties. If \(h_j(t)\) is the distance needed to accumulate mass \(t\), put \[V_j=\int_0^1h_j(t)\,dt.\] Then \(L=\sum_jV_j\le\mathop{\mathrm{OPT}}_k\) for an optimal solution of (1).

As in LP-rounding constructions such as Charikar–Li (Charikar and Li 2012, sec. 2), it is convenient to split an opening into co-located pieces so that each client’s nearest-unit assignment uses whole pieces. There are only polynomially many needed split points: each client uses at most one partial piece in a fixed distance order. Refine each facility simultaneously at all the resulting cut points. Denote the whole-piece assignment of client \(j\) by \(F_j\). For a set \(U\) of pieces, write \(y(U)=\sum_{i\in U}y_i\); thus \(y(F_j)=1\) and \(V_j=\sum_{i\in F_j}y_id(j,i)\).

Internal pieces at one location are allowed to be distinct, so their distance can be zero. All geometric arguments use only nonnegativity, symmetry, and the triangle inequality. Whenever a procedure returns open pieces, merge them by their original facility identity. Merging preserves every client’s connection cost and cannot increase the number of facilities.

Constant-additive pseudo-approximation

We use the following established result.

The theorem permits separate client and candidate sets and applies the pseudo-algorithm to restricted candidate instances. Our construction works on every such instance. We will obtain a deterministic pseudo-algorithm before applying Theorem 4; the external reduction is not being used to convert an expected-cost statement into determinism.

For each fixed accuracy we permit all grid denominators, numbers of rounds, and additive allowances to be constants of any size. They may depend on that accuracy, but not on \(k\), the number of input points, their distances, or their bit lengths. Thus it suffices to prove, for each fixed \(\tau>0\), a deterministic \(c_\tau\)-additive \((1+2/e+\tau)\)-approximation. Our analysis in fact bounds the pseudo-solution against \(L\), the initial fractional cost. This does not assert the same factor for rounding the original relaxation to exactly \(k\) facilities: the additive allowance and the subsequent reduction are essential.

Probability and bookkeeping

We use \(\mathbb E\), \(\mathbb P\), and \(\mathop{\mathrm{Var}}\) for expectation, probability, and variance. Every distribution used below has finite support. Randomized procedures first provide the estimates; Section 6 replaces their choices by explicit polynomial-size lists.

Two elementary facts will be used repeatedly. First, for a nonnegative random cost \(C\) and an event \(E\) of probability at least \(1-\rho>0\), \[ \mathbb E[C\mid E]\le\frac{\mathbb EC}{1-\rho}. \tag{2}\] No cost–event independence is required. Second, a finite distribution contains an outcome whose cost is at most its expected cost. These facts allow us to discard budget failures and select a single output by total cost; we never require all individual client expectation bounds to hold simultaneously at that output.

The notation \(O_\Delta(\cdot)\) means that its implicit constant may depend on the fixed grid denominator \(\Delta\), and on no input parameter. All unqualified constants in the geometric and analytic estimates are absolute, unless an explicit previously fixed parameter is indicated.

The graph and the physical process

The source-opening and neighbor-erasure operation follows the framework of Byrka et al. (Byrka et al. 2026, Definitions 3.2 and 5.1). We specify the graph, taper, and batches here because their precise properties will be used in the sub-two drift analysis.

Represent each opening by active copies of mass \(g\), grouped by facility. Each group has at most \(\Delta\) active copies. Co-located facilities introduced by earlier splitting may remain distinct groups throughout the process; they can be merged in the output. A group always means one such labeled facility group, rather than an entire zero-distance class. An integral group has exactly \(\Delta\) copies.

Fix a small time step \(\ell>0\) and a large time horizon \(\Lambda\). There are \(T=\lceil\Lambda/\ell\rceil\) rounds. At the start of round \(r\in\{0,\ldots,T-1\}\) let \(\xi=r\ell\) and \(x=e^{-\xi}\). For that round choose a fixed rational number \(\widehat P\) satisfying \[ 0<\widehat P\le\tfrac12, \qquad |\widehat P-P(x)|\le\zeta, \qquad 0<\zeta<\tfrac14. \tag{43}\] The finitely many rational approximations may be fixed in advance. The analytic value \(x\) is used only in the analysis.

For each active target copy \(i\), independently of the other targets conditional on the current state, choose a set \(\mathcal H_i\) of exactly \(\Delta\) active source copies. If the colocated active mass is below one, choose \(N_i>0\) minimally so that \[ q_{hi}=\min\left\{1,\frac{(1-d(i,h)/N_i)_+}{\widehat P}\right\}, \qquad \sum_h gq_{hi}=1, \tag{44}\] and give \(\mathcal H_i\) these marginals. Such a radius exists: the sum is continuous and nondecreasing in the positive radius, starts below one, and reaches the total active mass at a finite radius. This last observation also handles total active mass exactly one. The equation is piecewise linear in \(1/N_i\). Fixed-sum rounding realizes the marginals since their unscaled sum is the integer \(\Delta\).

If the colocated active mass is already at least one, set \(N_i=0\) and choose \(\Delta\) colocated copies, including every active copy in the target’s own group. Any fixed such choice suffices in this case; denote its marginals by \(q_{hi}\) as well. In either case every own-group copy belongs to \(\mathcal H_i\) with certainty, and \[ \begin{aligned} N_i(1-\widehat Pq_{hi})&\le d(i,h)&&\text{if }q_{hi}<1,\\ d(i,h)&\le N_i(1-\widehat Pq_{hi})&&\text{if }q_{hi}>0. \end{aligned} \tag{45}\] For \(N_i>0\) these statements follow from the three cases of the clipped linear formula; for \(N_i=0\) all sources of positive marginal are colocated. Define the out-set of source \(h\) by \[E_h=\{i:h\in\mathcal H_i\}.\] A full update at \(h\) erases \(E_h\) and opens the facility of \(h\).

We use balanced and biased unbalanced batches, with forced openings for large updates, as in the approach of Byrka et al. (Byrka et al. 2026, Algorithms 4–5 and Section 5.4). The exact batch and chunk rules below are part of our construction. To implement many updates while controlling the total opening mass, put \[u_h=1-g|E_h|, \qquad \mathcal Z=\{h:u_h=0\},\quad \mathcal P=\{h:u_h>0\},\quad \mathcal N=\{h:u_h<0\}.\] Every target has \(\Delta\) incoming sources, so \(\sum_hu_h=0\). Write \[ A_* =\sum_{h\in\mathcal P}u_h =\sum_{h\in\mathcal N}|u_h|, \qquad |\mathcal P\cup\mathcal N|\le 2A_*/g. \tag{46}\] The last inequality uses that each nonzero \(u_h\) is a nonzero multiple of \(g\). Choose a further small constant \(b_*>0\), to be fixed last, and put \[B_* =\max\{1,b_*A_*\}.\] A source is heavy if \(g|E_h|>B_*\). Force open the facilities of all heavy sources, whether or not any random update is triggered.

All random choices and erasure sets in a round refer to its pre-round graph and active copies. Set \(p_*=2\ell g\), and fix a small bias \(0<\phi\le1\). Flip an independent fair coin to choose one of the following two types of round.

1. In a balanced round, propose each \(h\in\mathcal Z\) independently with probability \(p_*\). Select a proposed source only if no other proposed source has an out-set intersecting its out-set.

2. In an unbalanced round, select sources of \(\mathcal P\) independently with probability \(p_*\), and nonheavy sources of \(\mathcal N\) independently with probability \(p_*(1+\phi)\). For each heavy source, partition its out-set into chunks of mass at most \(B_*+1\) and trigger each chunk independently with probability \(p_*(1+\phi)\). A triggered chunk erases just that chunk; its source facility is already included among the forced openings.

All probabilities are at most one for the parameter choices below. Take the union of the selected full out-sets and triggered chunks and erase it. Then reset each forced or selected facility to \(\Delta\) active copies. Copies removed by either erasure or reset are distinct from the regenerated copies.

Two invariants are useful. First, an integral target group’s own \(\Delta\) copies exhaust each of its incoming sets. Therefore an erasure of an integral group is caused by one of its own sources; that source either is selected or is heavy. The group is reset in either case. Thus integral facilities persist, including when other facilities are colocated with them. Second, every erasure has an associated selected or forced opening. Hence total active mass always stays at least one. The construction of the next round’s graph is consequently well defined.

The opening-count bound

We first bound the number of facilities independently of connection cost. The next lemma controls the total active mass with high probability; completion will later turn that mass bound into a facility-count bound.

For later use, the mean and variance of \(Z_*\) depend only on indicator moments through degree four, since \(Z_*\) is a quadratic polynomial. The probability bounds for balanced selections below need only first and second moments of their proposal indicators. These are the moment requirements used in Section 6.

A persistent backup and final completion

We next establish a persistent bound on each client’s eventual connection distance, including trajectories on which none of its assigned copies opens. This is the backup-and-completion role in iterative rounding (Byrka et al. 2026, sec. 3.3 and Lemmas 6.3 and 6.5); the following argument proves the bound for our tapered graph and reset rule.

Fix a client for the remainder of the analysis. Let \(S_0\) be the \(\Delta\) copies in its initial nearest unit, chosen using any fixed tie-breaking rule. Write \(d_i\) for its distance to copy \(i\), \[V_0=\sum_{i\in S_0}gd_i, \qquad R=\max_{i\in S_0}d_i.\] In particular \(R\le\Delta V_0\).

At any intermediate time define \(b\) to be the distance from the fixed client to its nearest integral facility, capped at \(5R\); when no integral facility exists use \(b=5R\). By persistence and Lemma 11, the eventual completed solution has connection distance at most \(b\) at every time on the trajectory. The cap need not already be attained by an open facility. It is a valid bound on the eventual connection because completion supplies \(5R\) on every trajectory. The quantity \(b\) never increases.

The potential and its single-action bounds

For analysis only, track an alive set \(S\subseteq S_0\) of original copies. A physically removed copy leaves \(S\) permanently, even when its facility is reset with regenerated copies. We also permit the additional discards specified below. Put \[y=g|S|,\qquad V=\sum_{i\in S}gd_i.\] At elapsed time \(\xi\), with \(x=e^{-\xi}\), define \[ f=aV+(1-y)b, \qquad a=A(x)+B(x)(y-x)=1+m(x)+B(x)y. \tag{53}\] The coefficient \(a\) is the value at \(y\) of the tangent line to the concave function \(A\) at \(x\); \(x\) is a deterministic clock, while \(y\) is the realized alive mass. Its affine dependence on \(y\) makes the erasure correction quadratic in removed mass and cost, so the graph average of the ideal bound below requires only first and second edge moments. In particular \(a\ge1\). The initial potential is \[ f_0=A(1)V_0=(1+2/e)V_0. \tag{54}\] If \(C_j\) denotes the actual connection distance in the completed output, the preceding cap property gives \[ C_j-f_T\le y_Tb_T-a_TV_T\le5Ry_T. \tag{55}\] The potential and every single-action bound below are \(O(R)\) uniformly over time: \(V\le R\), \(0\le y\le1\), \(b\le5R\), and \(a\) is bounded.

We first freeze time and the graph. Write \(B=B(x)\) and \(b_h=\min\{b,d_h\}\). For a source \(h\notin S\), let \[\delta_h=g|E_h\cap S|, \qquad Y_h=\sum_{i\in E_h\cap S}gd_i.\] If only this full source update affects \(S\), its removal changes the coefficient to \(a-B\delta_h\) and makes the new cap at most \(b_h\). Its potential change is consequently at most \[\begin{align*} &(a-B\delta_h)(V-Y_h)+(1-y+\delta_h)b_h-f\\ &\quad=\delta_hb_h-aY_h-B\delta_h(V-Y_h) +(1-y)(b_h-b)\\ &\quad\le\delta_hb_h-aY_h-B\delta_h(V-Y_h). \end{align*}\] The same statement holds when \(h\) is not alive but its group contains alive copies: all those copies belong to \(E_h\) by the own-group edges, so its reset creates no unaccounted removal.

If an alive source \(h\in S\) opens, retain among the physical survivors only copies satisfying \(ad_i<d_h\). If their mass and cost are \(y',V'\), their new coefficient is at most \(a\) and their new cap is at most \(d_h\). Therefore \[f'\le aV'+(1-y')d_h =d_h-\sum_{\substack{i\in S\\h\notin\mathcal H_i}} g(d_h-ad_i)_+.\] Copies satisfying equality can be discarded without changing this bound. Figure 1 separates physical erasure from these additional analysis discards. This continuation after an opening supplies the subtraction that will enter Theorem 7.

The ideal first-order bound, assigning rate \(g\) to each full source update, is thus \[ \begin{aligned} \mathcal B={}& \sum_{h\in S}g\left[ d_h-\sum_{\substack{i\in S\\h\notin\mathcal H_i}} g(d_h-ad_i)_+-f\right]\\ &+\sum_{h\notin S}g\left[ \delta_hb_h-aY_h-B\delta_h(V-Y_h)\right]. \end{aligned} \tag{56}\] At a fixed alive set and cap, the time derivative is \[ \dot f=-x(y-x)B'(x)V. \tag{57}\] Indeed the derivative of \(A(x)+B(x)(y-x)\) with respect to \(x\) is \(B'(x)(y-x)\).

From the physical batch to the ideal bound

For the actual analysis update, first account for forced openings by removing from \(S\) all copies whose groups are forced and updating the cap. Then apply the round’s erasures and resets. If a source that was alive after force removals but before these erasures was selected, choose one such source \(h\) and discard every remaining copy with \(ad_i\ge d_h\), using the pre-round value of \(a\). This is purely an analysis rule and need not be computed by the algorithm.

The order of conditioning in Lemma 12 is essential. Whether a group is forced depends on the graph. We have proved the drift comparison for every realized graph, in both force cases; we will now average that comparison over the original graph law, without conditioning the law on a forcing event.

Averaging the graph and applying the finite inequality

The physical batch has now been compared with the ideal generator for every realized graph. We average that generator using edge marginals and cross-target second moments, then verify the hypotheses of the abstract finite comparison inequality. The resulting bound is the cost estimate that will be summed over rounds.

Survival, accumulated error, and parameter choices

We finish the proof of Theorem 9. First run the process without rejecting count failures, and complete every output by Lemma 11. Average Lemma 12 over the graph and use Lemma 13. Since \(V\le V_0\), \(B\le1\), and \(R\le\Delta V_0\), summing over \(T\) rounds gives \[\mathbb Ef_T\le \left[1+2/e+(\Lambda+1) \{g+12\zeta+C_\Delta\Delta(\phi+\ell)\}\right]V_0\] for a constant \(C_\Delta\), provided \(\ell\le1\). Equations (55) and (71) consequently imply the unconditional estimate \[ \mathbb EC_j\le \left[1+2/e+5\Delta e^{-\Lambda/2} +(\Lambda+1)\{g+12\zeta+C_\Delta\Delta(\phi+\ell)\}\right]V_0. \tag{72}\] All these estimates also apply when \(V_0=0\): then \(R=0\), the completed connection cost is zero, and no normalization by \(V_0\) is used.

Here is an explicit order of choices that avoids circular dependence. It suffices to consider \(0<\tau\le1\).

1. Choose a sufficiently large constant integer \(\Delta\) and put \(g=1/\Delta\). For each such choice take, for example, \[\Lambda=\left\lceil2\log(40\Delta/\tau)\right\rceil+1.\] Then \(5\Delta e^{-\Lambda/2}\le\tau/8\). Taking \(\Delta\) sufficiently large also makes \((\Lambda+1)g\le\tau/8\), since \(g\log(1/g)\longrightarrow0\). Any additional upper bound on \(g\) required by the preparation stage may be imposed at the same time.

2. With \(\Delta,\Lambda\) fixed, choose a small positive rational \(\zeta\) so that \(12(\Lambda+1)\zeta\le\tau/8\). Next choose a small positive rational \(\phi\), and then a small positive rational \(\ell\), so that \[C_\Delta\Delta(\Lambda+1)(\phi+\ell)\le\tau/8.\] Require in addition every upper bound on \(\ell\) used in the count, batch, and survival lemmas, and \(p_*(1+\phi)\le1\). The batch constants are uniform in the remaining threshold \(b_*\), so these requirements can be imposed before choosing that threshold. This fixes the constant number of rounds \(T=\lceil\Lambda/\ell\rceil\).

3. Put \(\alpha=1+2/e\) and choose a fixed \(0<\eta\le\tau/(4(\alpha+\tau))\). Finally choose rational \(b_*>0\) as in Lemma 10 so that the final mass is at most \(M+C\) with probability at least \(1-\eta\).

The unconditional cost in (72) is then at most \((\alpha+\tau/2)V_0\) for every client. Repeat the entire process on failure of the final mass test. The accepted output is the unconditional output conditioned on that test, whose success probability is at least \(1-\eta\). Since connection costs are nonnegative, \[\mathbb E[C_j\mid\text{mass test succeeds}] \le\frac{\alpha+\tau/2}{1-\eta}V_0 \le(\alpha+\tau)V_0.\] Completion uses at most the final active mass, so an accepted output has at most \(M+C\) facilities. This proves the count and cost assertions of Theorem 9.

For fixed accuracy all parameters are instance-independent constants. The number of current facilities never exceeds the number in the split input, and every group has at most \(\Delta\) active copies. The graph, its marginals, a fixed-sum rounding, the action sets, their intersections, and the greedy completion can therefore be computed in polynomial time. The radius equations and all physical probabilities are rational; the approximations in (43) are fixed constants. Only the analysis uses the transcendental functions or alive-set discards. The expected number of repetitions is at most \(1/(1-\eta)\), giving expected polynomial running time. Section 6 replaces the local random choices by polynomial-size lists and obtains a deterministic output with the corresponding total-cost bound.

Compression of finite distributions

A distribution list is a finite list of states, each equipped with a nonnegative rational weight, whose weights sum to one. The state records the data needed to continue a construction; a feature vector records the statistics whose expectations must be retained.

The next lemma is the constructive finite-dimensional Carathéodory argument for preserving an expectation. Its connection with moment matching is classical (Bayer and Teichmann 2005); we include the exact elimination and its bit bound because both are needed by the algorithm.

Support reduction preserves every property satisfied by each state in the list, whether or not that property appears among the features. It does not preserve the entire joint distribution or arbitrary future adaptive behavior. We therefore specify the retained features at each use and verify the required estimates for the resulting process.

A distribution list for the preparation stage

Let \(M_\star\) denote the number of stars in the preparation construction, and write \(X_e\) for the current count of star \(e\). Each count remains within its original floor and ceiling. A weighted pair update within a strict sign preserves the weighted sum and has conditional mean equal to the current count vector. For every pair of distinct coordinates \(e,f\), it also satisfies \[ \mathbb E[X'_e\mid X]=X_e, \qquad \mathbb E[X'_eX'_f\mid X]\le X_eX_f. \tag{74}\] Indeed, unless both coordinates are the selected pair, the product is affine in the update parameter. If they are the selected pair, their increments have opposite signs and the product is concave. A one-coordinate rounding with the correct mean preserves every distinct pair-product expectation.

These assertions hold for every admissible state and every permissible choice of the next pair. Thus the next pair may be selected adaptively, for example as the first two fractional coordinates in a fixed order. After branching one update at every state, apply Lemma 15 with the features \[X_e \quad\text{and}\quad X_eX_f \qquad(e<f).\] Induction using (74) preserves the initial first moments and the upper bounds on distinct pair moments. Full independence between the positive, negative, and zero-coefficient classes is unnecessary for these conclusions.

The weighted sums and floor/ceiling restrictions are properties of each retained state. Each productive pair update makes at least one fractional count integral, and integral counts are never updated again. There are therefore only \(O(M_\star)\) updates, with identity updates padding states that finish earlier. Ordinary one-coordinate rounding also handles the zero coefficients and the final positive-side fraction. We leave the final negative-side fraction unrounded until these operations are complete. Each compressed list has \(O(M_\star^2)\) entries.

A distribution list for each physical round

Fix a physical state and a round. For each target copy, first construct a small list of source sets of cardinality \(\Delta\) with the prescribed edge marginals. Use ordinary fixed-sum pair rounding: after each pair update, branch each current state to its at most two endpoints and immediately compress the resulting list while preserving first moments. Thus no full tree of a target’s pair updates is ever enumerated, and every intermediate list has polynomial support. All sets in this list contain the target’s own group and respect the edges whose prescribed marginal is zero or one.

Combine the target lists by successive independent products. After each product, compress while retaining every edge first moment and every product of two edge indicators belonging to distinct targets. On adjoining a new target independently, its edge products with earlier targets have the required expectations; compression then preserves them. Induction gives a graph list with all the first and cross-target second moments of the independent-target graph construction. Its feature dimension and support size are polynomial in the number of active copies.

The ideal generator in (56) is linear in edge indicators except for the terms \(\delta_hY_h\). Expanding these terms gives products of edges from a common source to distinct targets, and diagonal terms reducible to first moments. The retained statistics therefore suffice for its graph average. The heavy-source rule and forced openings depend on more than these moments, but this introduces no extra averaging requirement: the batch drift bound (58) holds conditional on every realized graph, with that graph’s own forced openings and chunks. Only after that bound is applied do we average the ideal generator, obtaining (61) for the compressed graph distribution.

Rational implementation and bit complexity

We give the encoding details because support size alone does not imply a polynomial-time algorithm. Fix rational values of all algorithmic parameters. In each of the finitely many rounds, fix in advance a rational taper value \(\widehat P\) satisfying its stipulated accuracy. The analytic functions of elapsed time, the potential, and the analysis-only alive-set decisions need not be evaluated by the algorithm.

An optimal basic LP solution has polynomial encoding length. Splitting facilities to make assignments use whole pieces introduces only polynomially many pieces. Piece masses, bundle masses, and intersection masses are sums and differences of the LP coordinates. The extra star weights multiply a deficit by the fixed rational \(\gamma\), and conversion to grid counts multiplies by the integer \(\Delta\). Hence the initial star counts have a common denominator \(D_0\) of polynomial bit length.

For weighted pair rounding, write \(a_e=|u_e|\). If a pair update makes \(X_e\) integral and leaves \(X_f\) fractional, preservation of its weighted sum gives \[a_fX'_f=a_fX_f+a_eX_e-a_eX'_e.\] Throughout these updates, therefore, every weighted count lies in \(D_0^{-1}\mathbb Z\). Recovering \(X_f\) only divides by the polynomially bounded integer \(a_f\). Counts remain in their original integer bounds, so their numerators and their endpoint probabilities have polynomial bit length as well. Final single-coordinate rounding makes the coordinate integral; zero-coefficient counts retain their initial values until rounded directly. Clipping and rescaling give fresh rational inputs of polynomial bit length for the within-bundle rounds, where the same argument applies with unit coefficients.

The taper radii can also be computed exactly. On an interval where \(m\) source copies have saturated marginal one and \(r\) source copies have marginal strictly between zero and one, the radius equation is \[m+\frac{r}{\widehat P} -\frac{\sum_{h\text{ partial}}d(i,h)}{\widehat P N_i} =\Delta.\] When \(r>0\), an interior solution is \[ N_i= \frac{\sum_{h\text{ partial}}d(i,h)} {r+\widehat P(m-\Delta)}. \tag{76}\] The breakpoints are the rational numbers \(d(i,h)\) and \(d(i,h)/(1-\widehat P)\). Sorting these values and checking the successive intervals finds the smallest solution, including a breakpoint solution. If \(r=0\), the equation is constant on the interval, and the first point of a solution interval is a breakpoint. The separately prescribed case of sufficient co-located mass uses \(N_i=0\). All radii, marginals, and resulting branch probabilities have polynomial bit length. The physical states and their feature values are determined without reference to the list weights.

It remains to control the list weights across repeated branching and compression. Suppose a feature vector has \(d\) coordinates of bit length at most \(L\). In the affine-dependence computation, clearing denominators in each row of the matrix on at most \(d+2\) columns produces integer entries of bit length \(O(dL)\). A nonzero kernel vector can be chosen using minors, with bit length \(O(d^2L+d\log(d+1))\) by the determinant bound. These bounds are independent of the current list weights.

By (73), an elimination multiplies a common weight denominator by just one coefficient of such a dependence. An exact branch with probability denominator \(q\) multiplies a common denominator by at most \(q\). Branching a whole list adds the bit lengths of only polynomially many such denominators, and taking a product of lists adds their denominator bit lengths. There are polynomially many list operations, including the entire constant-depth physical history tree. Their common denominators and numerators consequently retain polynomial bit lengths. Exact linear algebra, comparisons, and cost evaluations therefore give a deterministic algorithm whose bit complexity is polynomial in the input encoding length for every fixed choice of accuracy parameters.

Loss allocation and the approximation theorem

We first record explicitly how conditioning and deterministic selection are combined. For a nonnegative cost \(C\) and an event \(G\) with \(\mathbb P(G)\ge1-r\), \[ \mathbb E[C\mid G] \le\frac{\mathbb EC}{1-r}. \tag{77}\] Thus a finite distribution list contains a successful outcome whose cost is no more than the right-hand side. No independence between cost and budget success is required.

Applying the reduction on rational inputs

Before applying Theorem 4, we record its interface with the deterministic algorithm just proved. Every call of that reduction uses the same clients and distances with a restricted candidate set \(F'\subseteq F\); our pseudo-algorithm uses no additional input promise or prescribed-open facilities. Its internal forced openings are decisions of that call, not constraints passed from another call. Copies are merged back to original candidates before returning its output. If \(0<|F'|<k\), opening every candidate gives the optimum for that restricted instance. If \(F'=\varnothing\) and clients remain, discard the branch. The distinguished branch retaining a padded optimal \(k\)-set has neither exception. When the additive allowance is zero, the pseudo-algorithm already satisfies the required budget.

The zero-optimum case can be recognized before the reduction. Every client location must be an available candidate location, and there can be at most \(k\) distinct client locations. Equivalently, every zero-distance class containing a client must contain a candidate, and there are at most \(k\) such client classes. These conditions are necessary and sufficient in the original metric; when they hold, open those locations. Thus the positive-optimum analysis of the reduction does not require division by zero. Co-located pieces occur only inside our pseudo-algorithm and do not change this original-metric test.

We also make explicit a rational implementation of the deletion threshold in the proof of Li–Svensson’s reduction (Li and Svensson 2016, sec. 2, Algorithm 2). That proof analyzes a threshold determined by the unknown optimum. Its numerical value need not be computed. Suppose a call returns a set \(T\) with \(k<|T|\le k+c\). For every set \(T'\subseteq T\) obtainable by at most \(c\) deletions with \(|T'|>k\), and every \(i\in T'\), compute \[\Delta(T',i)=\mathop{\mathrm{cost}}(T'\setminus\{i\})-\mathop{\mathrm{cost}}(T').\] These numbers are nonnegative rationals of polynomial bit length. There are \(n^{O(c)}\) such numbers. Run the reduction’s deletion phase for each distinct number as threshold, and for zero, with a fixed order for choosing deletions; follow each run by its prescribed enumerative completion and retain the cheapest feasible result.

To justify this replacement, let \(B^*\ge0\) be the threshold used in the reduction’s analysis. Choose the largest enumerated value at most \(B^*\), including equality, or zero if there is none. Every test of the form \(\Delta(T',i)\le B^*\) has the same answer at this representative threshold. Induction through the at most \(c\) deletions therefore reproduces the ideal deletion history, including its fixed tie choices. Its output is present among the enumerated outputs. The approximation analysis still uses \(B^*\); the algorithm uses only exact rational comparisons. This adds \(n^{O(c)}\) overhead to the stated reduction. All candidate enumeration, distance comparisons and cost evaluations likewise use polynomial-bit rational arithmetic.

Finally, the approximation parameter supplied to the reduction can be a fixed rational upper bound on our ratio, with its excess charged to the allowed accuracy. The auxiliary positive constants in its proof may then be chosen rational as well. These observations preserve both determinism and polynomial binary-input time; they require no oracle for the optimum.

Combining Theorem 1 with the hardness bound recalled in Section 1 proves Corollary 2: under \(P\ne NP\), the infimum of deterministic polynomial-time approximation factors in the specified candidate-facility model is \(1+2/e\). The upper bound holds for every fixed positive loss and does not require attainment of the endpoint.

Rational approximations and their errors

Use the following polynomials in \(X\), with coefficients displayed in ascending order: \[ \begin{aligned} p(X)={}&.26424112+.16053514X+.05751895X^2\\ &\hspace{12mm}+.01313705X^3+.00456774X^4,\\ v(X)={}&1.73575888-.52832269X-.16191093X^2\\ &\hspace{12mm}-.03442938X^3-.01109588X^4,\\ k(X)={}&.0767+.0564X+.0142X^2+.0246X^3. \end{aligned} \tag{79}\] We will prove \[ |p-P|\le\epsilon_P:=.0000075,\qquad |v-A|\le2\epsilon_P=.000015,\qquad k_0\le k. \tag{80}\]

For the first two claims, the integral Taylor series are \[P(x)=\sum_{n\ge0}\frac{(-x)^n}{n!(n+2)},\qquad A(x)=3-2\sum_{n\ge0}\frac{(-x)^n}{n!(n+1)}.\] On \(0\le x\le1\) the magnitudes of the summands decrease. Thus the partial sums through orders \(12\) and \(13\) enclose each defining integral. Evaluating these rational enclosures at \(X=0,1/4,1/2,3/4,1\) proves that the nodal errors of the displayed \(p,v\) are, respectively, less than \(2.343560\cdot10^{-9}\) and \(2.344324\cdot10^{-9}\). In particular, each is less than \(4\cdot10^{-9}\). This nodal calculation can be performed without evaluating any transcendental number: at each node substitute \(x=1-X\) into the two finite sums just given and compare both rational endpoints with the displayed polynomial value.

Let \(I_P,I_A\) be the exact degree-four interpolants at these five nodes. The sum of the absolute values of the Lagrange basis polynomials is at most \[\sum_{i=0}^4\frac{4^4}{i!(4-i)!}=\frac{512}{3},\] because every numerator factor has absolute value at most one. Differentiating with respect to \(X\) gives the fifth-derivative bounds \(1/7\) for \(P\) and \(1/3\) for \(A\). The absolute node product, with \(u=|X-1/2|\), is \[u(1/4-u^2)|1/16-u^2|.\] For \(u\le1/4\) it is at most \(1/256\). For \(1/4\le u\le1/2\), the product of the two quadratic factors is at most \(((1/4-1/16)/2)^2=9/1024\), so the entire expression is at most \(9/2048\). Both bounds are less than \(.0045\). The interpolation remainder formula therefore gives \[\begin{align*} |p-P|&\le\frac{512}{3}(4\cdot10^{-9})+\frac{.0045}{840} =\frac{31709}{5250000000}<\epsilon_P,\\ |v-A|&\le\frac{512}{3}(4\cdot10^{-9})+\frac{.0045}{360} =\frac{9887}{750000000}<2\epsilon_P. \end{align*}\]

The coefficient signs and endpoint values give \[ \begin{gathered} .26424112\le p\le.5,\qquad 1\le v\le1.73575888,\\ .0767\le k\le.1719,\qquad .014\le H_0\le.394. \end{gathered} \tag{81}\] For \(H_0\), use \(H_0'(X)=X(-.84+.12X)\le0\). The first row of Table 1 below proves \[ T_k:=4k-(p+\epsilon_P)(1+k)^2>0. \tag{82}\] Since \(0<k<1\) and \(t\mapsto4t/(1+t)^2\) is increasing on \((0,1)\), while \(4k_0/(1+k_0)^2=P\), inequality (82) and \(P\le p+\epsilon_P\) imply \(k_0\le k\). This proves all of (80).

Define polynomial versions of the quantities used in the following reductions by \[ \begin{gathered} c=1-sp,\qquad N=c-(1-s)k,\qquad l=c+1-v+spx,\\ d=d_s,\qquad h=1.6+.4d,\qquad u_0=d/h,\qquad u_1=(2-s)(1+d)/(2h). \end{gathered} \tag{83}\] The polynomial families in the reductions are \[\begin{align*} T_1&=1-p-1.12p(1+x)+.4p^2X(1+x),\\ T_2&=((1-p)^2-.38pX+px)(1+x)-2x(1-p),\\ Q_0(W)&=N-dv+pW(hW-2d),\\ Q_1&=N-(2-s)px-d(cp-l)-hp u_1^2,\\ Y&=N-p\bigl((2-s)x+(2-s)^2/6.4\bigr),\\ Y'&=N+p(1.6X^2-2)\qquad(s=0). \end{align*}\] The notation \(T(P,A)\) means replacement of \(p,v\) by \(P,A\), keeping \(X,k\) and the rational constants fixed.

Reduction of the first, second, and final scalar bounds

We first establish the first scalar inequality \[ m+b-.4z/r\le\frac{C(1-P)}{1-Pr}\qquad(0<r\le1). \tag{84}\] Put \(t=1/r-P>0\). By (78), the right side of (84) minus its left side equals \[\begin{align*} &C(1-P)-\frac{2x(1-P)}{1+x}+sPx-(2-s)Pw+.4P^2w^2\\ &\hspace{35mm}+P\left(\frac{C(1-P)}t+.4w^2t\right). \end{align*}\] Applying AM–GM to the last parentheses gives the lower bound \[ C(1-P)-\frac{2x(1-P)}{1+x}+sPx -(2-s-2\sqrt{.4C(1-P)})Pw+.4P^2w^2. \tag{85}\] For \(s=0\) the coefficient \(2-2\sqrt{.4(1-P)}\) is at most \(1.12\). The resulting quadratic decreases for \(0\le w\le X\): its derivative is at most \(-1.12P+.8P^2<0\). Substituting \(w=X\) into this lower bound gives \(X T_1(P)/(1+x)\), where \(T_1\) is defined above. For \(s=1\), the coefficient \(1-2\sqrt{.4}(1-P)\) is at most \(.38\). Dropping the nonnegative quadratic term and using \(w\le X\) gives \(T_2(P)/(1+x)\). The numerical coefficient comparisons also follow by squaring positive quantities: \(.2>.44^2\) and \(.4>.62^2\). Table 1, with the substitution budget proved below, gives \(T_1(P),T_2(P)>0\). At \(X=0\), the extra factor \(X\) in the \(s=0\) reduction is zero, which still proves the required non-strict inequality.

For \(r>0\), (84) implies \[L\ge Cr\left(1-\frac{1-P}{1-Pr}\right)\ge0.\] For \(r=0\), one has \(L=.4z\ge0\) directly.

To check \(K\ge d_sG\), write \(d=d_s\) and \(h=1.6+.4d\). The difference \(K-dG\) is affine in \(r\). With \(N_P=C-(1-s)k\), replacing \(k_s\) by \((1-s)k\) only decreases this difference. Put \(\underline K=N_P-br+1.6z\le K\). The endpoint expressions after this replacement are \[\begin{align*} (\underline K-dG)|_{r=0}&=N_P-dA+Pw(hw-2d), \tag{86}\\ (\underline K-dG)|_{r=1}&=Q_1(P,A)+hP(w-u_1)^2, \tag{87}\end{align*}\] where \[u_0=d/h,\qquad u_1=\frac{(2-s)(1+d)}{2h}\] and \(Q_1\) is defined above. For \(s=1\) we may use the unrestricted minimum \(w=u_0\) in (86). For \(s=0\), the constrained minimum is at \(w=X\) if \(X\le u_0\), and at \(w=u_0\) otherwise. Equation (87) is bounded below by \(Q_1(P,A)\) for both values of \(s\). These are precisely the \(Q_0,Q_1\) checks in Table 1.

We next reduce the final two scalar bounds. Recall \[\kappa_s(q)=C\frac{1-q}{1-Pq}-(1-q),\qquad 0\le q\le1.\] For \(s=1\), \(k_s=0\) and \(\kappa_s(q)\le0\). Thus both relevant left sides are bounded below by \(C+1.6z-b\). Its unrestricted minimum in \(w\) is the true-parameter version of \(Y\) above, so the rows \(Y-H_0\) and \(Y-.25\) suffice.

For \(s=0\) and \(X\ge.5\), use \(q\le1\) and \(\kappa_0(q)\le k_0\). Again a common lower bound is \(1-k_0+1.6z-b\), whose unrestricted minimum is the corresponding \(Y\). For \(s=0\) and \(X\le.5\), we have \(x\ge.5\) and \[P(x)\le P(.5)<.47.\] For example the inequality \(e^{-u}\le1-u+u^2/2\) yields \(P(.5)\le1/2-1/6+1/32=35/96<.47\). It follows that \(b\ge .52P/(1-P)\), because \(x+w\ge.5\) and \(1-P>.53\). Since \[\kappa_0(q)=\frac{Pq(1-q)}{1-Pq} \le\frac{P}{1-P}(1-q),\] we obtain \(bq+.52\kappa_0(q)\le b\). The fourth scalar left side is consequently at least \[1-.48k_0+1.6z-b.\] The quadratic \(1.6Pw^2-2P(x+w)\) decreases on \(0\le w\le X\le.5\), so its minimum is \(P(1.6X^2-2)\). Thus the fourth and fifth inequalities follow, respectively, from \(Y'+.52k-H_0\ge0\) and \(Y'-.25\ge0\) after substitution. In the former expression the net coefficient of \(k\) is \(-.48\), so replacing \(k_0\) by the upper bound \(k\) is conservative.

The last scalar bound also gives positivity of \(K\): since \(b\ge0\), \[K=C-k_s-br+1.6z\ge C-k_s-b+1.6z\ge.25>0.\]

The first coefficient table and its exact verification

Replacing \(p,v\) by \(P,A\) in these expressions changes their values by less than \[ \sigma=.0002. \tag{88}\] For completeness, the following derivative bounds hold on the line segment between \((p,v)\) and \((P,A)\), using \(0\le p,P\le.5\), \(1\le v,A\le2\), and \(0\le W\le1\): \[\begin{array}{c|rr|r} T&|\partial_pT|&|\partial_vT|& 10^6\,|T(p,v)-T(P,A)|\text{ upper bound}\\ \hline T_1&3.64&0&27.3\\ T_2&8.76&0&65.7\\ Q_0(W)&3.848&.52&36.66\\ Q_1&4.848&.52&44.16\\ Y&2.625&0&19.6875\\ Y'&2&0&15 \end{array}\] Here \(T_{1,p}=-1-1.12(1+x)+.8pX(1+x)\), with \(X(1+x)\le1\), and \(T_{2,p}=[-2(1-p)-.38X+x](1+x)+2x\). For \(Q_0\), use \(h\le1.808\) and \(2d\le1.04\). For \(Q_1\), \((cp-l)_p=1-2sp+sX\) has absolute value at most \(2\), \((cp-l)_v=1\), and \(u_1\le1\). The last two derivative bounds follow directly from their displayed linear dependence on \(p\). These estimates prove the claimed budget even on the full interval, rather than only the listed subintervals. Adding \(.52k\) or subtracting \(H_0,.25\) does not change this budget. The row for \(T_k\) already includes \(\epsilon_P\) explicitly and does not need this additional substitution estimate.

Here is the finite verification mechanism used in both coefficient tables. A row with scale \(10^a\) and integers \(\lambda_0,\ldots,\lambda_n\) asserts \[ T(X)\ge t(X):=10^{-a}\sum_{j=0}^n\lambda_jX^j \qquad(0\le X\le1). \tag{89}\] To verify it, expand \(T=\sum_jc_jX^j\) and check \[ 10^ac_j\ge\lambda_j\quad(0\le j<n),\qquad 10^a\sum_{j\ge n}\min(0,c_j)\ge\lambda_n. \tag{90}\] Positive tail terms may be discarded; for every negative tail term, \(c_jX^j\ge c_jX^n\). This proves (89).

All short minorants have degree at most six. On any listed interval \([\alpha,\beta]\), their degree-six Bernstein coefficients are \[ \mathfrak b_m(t;\alpha,\beta)=10^{-a} \sum_{h=0}^m\frac{\binom mh}{\binom6h}(\beta-\alpha)^h \sum_{j=h}^n\lambda_j\binom jh\alpha^{j-h}, \quad 0\le m\le6. \tag{91}\] Indeed, expand at \(X=\alpha+(\beta-\alpha)u\) and convert ordinary powers of \(u\) to degree-six Bernstein form. Thus if every coefficient in (91) is at least \(b_*\), then \(t(X)\ge b_*\) throughout the interval, since the Bernstein basis is nonnegative and sums to one. An endpoints entry with more than two numbers prescribes each interval between consecutive numbers.

Substitution of (79) into the definitions verifies (90) for all fourteen rows of Table 1. As additional explicit sign certificates, the minima of the finite sets in (91), over all prescribed subintervals for each row, are at least \(10^{-6}\) times the following integers, in row order: \[ \begin{gathered} 114, 20000, 10000, 10000, 390000, 40000, 33818,\\ 2743, 15000, 80000, 3625, 31250, 356, 21250. \end{gathered} \tag{92}\] Both (90) and these lower bounds are finite rational comparisons specified completely by the table and (91). For example, the first row has smallest Bernstein coefficient \(73/640000\), and the thirteenth row has smallest coefficient \(27341/76800000\). Consequently \(T_k>0\), and every other short minorant exceeds \(\sigma=.0002\) on its indicated intervals. After the substitution budget, this completes the first and second groups of scalar bounds and the final two scalar bounds.

The third scalar bound: perturbation estimate

It remains to prove \[ ((J_s-b)r+1.6z)M+(C-k_s-J_sr)L\ge0. \tag{93}\] Since \(L\ge0\) is already established, replacing \(k_s\) with \((1-s)k\) only decreases the left side. We next replace \(P,A\) by \(p,v\). This replacement changes the expression by at most \[ \sigma(r+w^2),\qquad \sigma=.0002. \tag{94}\] Here are full bounds for this assertion. For either endpoint choice \((\pi,\nu)=(P,A)\) or \((p,v)\), put \(c_\pi=1-s\pi\), \(H=H_s\), and write the two products as \(F_1F_2+F_3F_4\), where \[\begin{align*} F_1&=[c_\pi H-(2-s)\pi(x+w)]r+1.6\pi w^2,\\ F_2&=(1-r)(\nu+2\pi w)+c_\pi\pi r,\\ F_3&=c_\pi(1-Hr)-(1-s)k,\\ F_4&=[2-\nu-s\pi X-(2-s)\pi w]r+.4\pi w^2. \end{align*}\] Writing \(\rho=r+w^2\), the ranges already established imply \[ |F_1|\le1.1\rho,\qquad |F_2|\le3,\qquad |F_3|\le1,\qquad |F_4|\le2\rho. \tag{95}\] For \(F_1\), both nonnegative terms in its coefficient of \(r\) are at most one, while \(1.6\pi\le.8\). The expression \(F_2\) is a convex combination of two quantities between zero and three. For \(F_3\), the lower bounds are \(1-.394-.1719=.4341\) when \(s=0\) and \(.5(1-.25)=.375\) when \(s=1\), and the upper bound is one. Finally the coefficient of \(r\) in \(F_4\) belongs to \([-1,1]\): \(0\le2-\nu\le1\) and \(0\le s\pi X+(2-s)\pi w\le1\).

Let \(\delta F_j\) denote the change between the two choices. With \(\epsilon_P=.0000075\), direct subtraction gives \[\begin{align*} |\delta F_1|&\le2\epsilon_P\rho,& |\delta F_2|&\le4\epsilon_P,\\ |\delta F_3|&\le\epsilon_P,& |\delta F_4|&\le4\epsilon_P\rho. \end{align*}\] For the first estimate, the coefficient multiplying \(\epsilon_Pr\) is \(sH+(2-s)(x+w)\), at most \(2\) for \(s=0\) and at most \(1.25\) for \(s=1\); the remaining coefficient is \(1.6\). For the second, use the convex combination in \(r\), together with \(|\delta(\nu+2\pi w)|\le4\epsilon_P\) and \[|\delta(c_\pi\pi)| =|P-p|\,|1-s(P+p)|\le\epsilon_P.\] The third estimate is immediate from \(0\le1-Hr\le1\). For the fourth, use \(2+sX+(2-s)w\le4\) and \(.4\le4\). Combining these differences with (95) and the identity \(ab-a'b'=(a-a')b+a'(b-b')\) yields \[\begin{align*} |\delta(F_1F_2+F_3F_4)| &\le(3\cdot15+1.1\cdot30+2\cdot7.5+30)10^{-6}\rho\\ &=.000123\rho<\sigma\rho\qquad(\rho>0). \end{align*}\] If \(\rho=0\), both expressions vanish. Thus (94) holds throughout the domain. The factor bounds hold for both choices, so this telescoping estimate already accounts for products of perturbations.

Bernstein expansion and the second coefficient table

Use the polynomial quantities in (83), and set \[j=cH_s,\qquad b'_0=(2-s)px,\qquad b'_1=(2-s)pX, \qquad f=2pX,\qquad Z'=pX^2.\] In the polynomial version of the left side of (93), subtract \(\sigma(r+w^2)\) and put \(w=X\eta\), \(0\le\eta\le1\). Its degree-three Bernstein coefficients in \(\eta\), for \(i=0,1,2,3\), are \[ B_0+B_1r+B_2r^2, \tag{96}\] where, writing \(e=i/3\), \(e_2=i(i-1)/6\), and \(t=\mathbf{1}_{\{i=3\}}\), the exact formulas are \[ \begin{aligned} B_0={}&Z'\bigl((1.6v+.4N)e_2+1.6ft\bigr)-\sigma X^2e_2,\\ B_1={}&(j-b'_0)v+Nl\\ &+\bigl((j-b'_0)f-b'_1v-Nb'_1\bigr)e-b'_1fe_2\\ &-Z'\bigl((1.6(v-cp)+.4j)e_2+1.6ft\bigr)-\sigma,\\ B_2={}&(b'_0-j)(v-cp)-jl\\ &+\bigl((b'_0-j)f+b'_1(v-cp)+jb'_1\bigr)e+b'_1fe_2. \end{aligned} \tag{97}\] These formulas give all eight expansions, one for each pair \((s,i)\). They can also be derived directly by expanding \[\begin{align*} &\bigl([j-(2-s)p(x+w)]r+1.6pw^2\bigr) \bigl((1-r)(v+2pw)+cpr\bigr)\\ &\quad +(N-jr)\bigl([l-(2-s)pw]r+.4pw^2\bigr) -\sigma(r+w^2) \end{align*}\] in powers of \(\eta\) after \(w=X\eta\), and using the degree-three conversion rule: if \(a_h\) is the coefficient of \(\eta^h\), the \(i\)th Bernstein coefficient is \(\sum_{h=0}^i a_h\binom ih/\binom3h\).

In all eight cases \(B_0\ge0\). For \(i=0,1\), it is zero. Otherwise \[B_0=X^2\bigl(p(1.6v+.4N)-\sigma\bigr)e_2+1.6pX^2ft.\] Here \(N\ge.5\), \(v\ge1\), and \(p\ge.26424112\), so the bracket is at least \(.475434016>0\) and the other term is nonnegative.

Define \[D'=4B_0B_2-B_1^2,\qquad S'=B_0+B_1+B_2.\] For \(i\le1\), positivity of \(B_1\) and \(B_1+B_2\) proves nonnegativity of (96) on \(0\le r\le1\). For \(i\ge2\), any of the following is sufficient:

1. \(B_1,B_2\ge0\);

2. \(D'>0\);

3. \(S'>0\) and \(B_0-B_2>0\).

For the second criterion, \(B_0\ge0\) and \(D'>0\) force \(B_0,B_2>0\); completing the square proves positivity for every real \(r\). For the third, use \[B_0+B_1r+B_2r^2 =rS'+(1-r)(B_0-rB_2),\] where \(B_0-rB_2=(1-r)B_0+r(B_0-B_2)\ge0\).

Table 2 supplies these signs. It uses exactly the coefficient rule (90) and the degree-six Bernstein formula (91).

In the printed row order, the Bernstein coefficients of these short minorants, over all of each row’s subintervals, are at least \(10^{-6}\) times the following integers: \[ \begin{gathered} 10000, 10000, 90000, 3050, 78500, 171, 710,\\ 89250, 60, 50000, 80000, 40000, 40000, 1600,\\ 66000, 3560, 4590, 69000, 2518, 30000, 52656. \end{gathered} \tag{98}\] The bound \(90000\) in the third position applies to both indices in that row. These are all positive. The smallest exact Bernstein coefficient among the table’s checks is \(193875074979/3200000000000000\), in the ninth printed row. Thus all the required strict signs follow from the finite rational formula (91).

For clarity, the negative-tail calculations needed for the four discriminant rows are detailed next. Let \(b_{h,j}=[X^j]B_h\) and write \[c_j=[X^j]D' =\sum_u\bigl(4b_{0,u}b_{2,j-u}-b_{1,u}b_{1,j-u}\bigr),\] where coefficients outside the degree range are zero. For a short minorant of degree \(n\), group \(\min(0,c_j)\) in consecutive blocks of five starting at \(j=n\), truncating at degree \(28\). The five block sums are bounded below by \(.0001\) times the following integers: \[ \begin{array}{cc|c|c|rrrrr} s&i&n&\deg D'&1&2&3&4&5\\ \hline 0&2&6&20&-930&-11&-1&0&0\\ 0&3&6&22&-1070&-180&-2&-1&0\\ 1&2&4&28&-1312&-18&-1&-1&-1\\ 1&3&4&28&-8240&-1560&-29&-1&-1 \end{array} \tag{99}\] These twenty rational comparisons follow by substituting (97) and (79) into the displayed coefficient formula. Their row sums are, respectively, \(-.0942,-.1253,-.1333,-.9831\), which are at least the final scaled minorant coefficients \(-.1,-.14,-.135,-.995\). For all entries before the last, one checks simply \(10^a[X^j]T\ge\lambda_j\) as in (90). The same coefficient rule verifies the non-discriminant rows directly.

The tail cutoff is exact: the four discriminant degrees are \(20,22,28,28\), as recorded in (99), so no term above degree \(28\) has been omitted. More generally, from (97), the respective triples of degrees \((\deg B_0,\deg B_1,\deg B_2)\) are \[\begin{array}{c|rrrr} &i=0&i=1&i=2&i=3\\ \hline s=0&(-\infty,9,9)&(-\infty,10,10)&(10,10,10)&(11,11,9)\\ s=1&(-\infty,9,13)&(-\infty,10,13)&(10,14,13)&(11,14,12). \end{array}\] Here \(\deg0=-\infty\). These bounds also justify a uniform cutoff of \(28\) for every polynomial in the two coefficient tables.

For \((s,i)=(0,2)\), the first criterion covers \(X\in[0,.05]\) and the discriminant criterion covers \(X\in[.05,1]\). For \((s,i)=(0,3)\), the corresponding split is \(.015\), with the listed additional subdivisions used only for Bernstein positivity. For \((s,i)=(1,2)\), the split is \(.2\). For \((s,i)=(1,3)\), the first criterion covers \([0,.1]\), the discriminant criterion covers \([.1,.75]\), and the final criterion covers \([.75,1]\). For \(i=0,1\), the full interval is covered by \(B_1\) and \(B_1+B_2\). All endpoints belong to the stated intervals. Thus every coefficient (96) is nonnegative for \(0\le X,r\le1\). Nonnegativity and unit sum of the degree-three Bernstein basis in \(\eta\) prove the desired polynomial inequality for \(0\le w\le X\); when \(X=0\) one has \(w=0\) directly. Restoring the error budget (94) proves (93). Together with the preceding reductions, this completes the proof of Lemma 8.