Positive Univariate Polynomials:
SOS certificates, algorithms, bit complexity, and T-systems

Matías Bender *

Philipp di Dio †

Elias Tsigaridas ‡

December 27, 2025

Abstract

We consider certificates of positivity for univariate polynomials with rational coefficients that are positive over (an interval of) $\mathbb{R}$. Such certificates take the form of weighted sums of squares (SOS) of polynomials with rational coefficients.

We build on the algorithm of Chevillard, Harrison, Joldeş, and Lauter [6], and we introduce a variant that we refer to as uSOS. Given a polynomial of degree $d$ with maximum coefficient bitsize $\tau$, we show that uSOS computes a rational weighted SOS representation in $\mathcal{O}_B(d^3+d^2\tau)$ bit operations; the resulting certificate of posivitity involves rationals of bitsize $\tilde{\mathcal{O}}(d^2\tau)$. This improves the best-known complexity bounds by a factor of $d$ and completes previous analyses. We also extend these results to certificates of positivity over arbitrary rational intervals, via a simple transformation. In this case as well, our techniques yield a factor-$d$ improvement in the complexity bounds.

Along the same line, for univariate polynomials with rational coefficients, we introduce a new class of certificates, which we call perturbed SOS certificates. They consist of a sum of two rational squares that approximates the input polynomial closely enough so that nonnegativity of the approximation implies the nonnegativity of the original polynomial. This computation has the same bit complexity and yields certificates of the same bitsize as in the weighted SOS case.

We further investigate structural properties of these SOS decompositions. Relying on the classical result that any nonnegative univariate real polynomial is the sum of two squares of real polynomials, we show that the summands form an interlacing pair. Consequently, their real roots correspond to the Karlin points of the original polynomial on $\mathbb{R}$, establishing a new connection with the T-systems studied by Karlin [17]. This connection enables us to compute such decompositions explicitly. Previously, only existential results were known for T-systems. We obtain analogous results for positivity over $(0, \infty)$, and hence over arbitrary real intervals.

Finally, we present our open-source Maple implementation of the uSOS algorithm, together with experiments on various data sets demonstrating the efficiency of our approach.

<footnote>*Inria Saclay and CMAP, École Polytechnique, IPP, France</footnote>

<footnote>†University of Konstanz, Germany</footnote>

<footnote>‡Inria Paris and IMJ-PRG, Sorbonne Université, France</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>3</td></tr><tr><td>2 The uSOS algorithm</td><td>7</td></tr><tr><td>2.1 A detailed presentation of the uSOS algorithm</td><td>7</td></tr><tr><td>2.2 The bit complexity of the various steps</td><td>9</td></tr><tr><td>2.3 Overall complexity estimates</td><td>12</td></tr><tr><td>2.3.1 Dropping the assumptions</td><td>13</td></tr><tr><td>2.3.2 Witness point of non-nonnegativity</td><td>14</td></tr><tr><td>2.4 Positivity over an interval [a, b]</td><td>14</td></tr><tr><td>2.4.1 The separation bound of A<sub>φ</sub></td><td>14</td></tr><tr><td>2.4.2 The perturbation ε for A<sub>φ</sub></td><td>15</td></tr><tr><td>2.4.3 Overall complexity</td><td>16</td></tr><tr><td>3 Perturbed SOS certificate</td><td>16</td></tr><tr><td>4 Interlacing polynomials and T-systems</td><td>18</td></tr><tr><td>4.1 Strict positivity and Karlin points over ℝ</td><td>18</td></tr><tr><td>4.2 Strict positivity and Karlin points over [0, ∞)</td><td>19</td></tr><tr><td>5 Implementation and experiments</td><td>21</td></tr><tr><td>5.1 Some insights on the implementation</td><td>21</td></tr><tr><td>5.2 Experiments with positive polynomials</td><td>22</td></tr><tr><td>6 Conclusion</td><td>24</td></tr><tr><td>A Alternative proof for interlacing</td><td>28</td></tr><tr><td>B Useful algorithms and complexity bounds</td><td>29</td></tr><tr><td>B.1 Preliminaries on root separation and approximation</td><td>29</td></tr><tr><td>B.2 Bounds on the mimimum and the Fan-in algorithm</td><td>31</td></tr><tr><td>B.3 About T-systems</td><td>32</td></tr></table>
# 1 Introduction

A univariate polynomial $A$, with real coefficients that takes only nonnegative values over $\mathbb{R}$, admits a decomposition as a sum of (two) squares of real polynomials; this representation certifies the nonnegativity of $A$. Following Powers [28], a certificate of positivity is an algebraic identity<fnref n="1" /> that straightforwardly implies the nonnegativity of a (univariate) polynomial. Our focus is mainly on univariate polynomials with rational coefficients. Hence, in the analysis of the certificates, in addition to the involved polynomials, we also take into account the maximum number of bits we need to represent their coefficients. Along the same lines, the complexity of the corresponding algorithms refers to the number of bit operations.

The certificate(s) of positivity raises the following mathematical, algorithmic, and complexity-related questions:

(Q₀) What is the mathematical framework, usually a theorem, that implies an algebraic identity (or more than one) corresponding to the certificate?

(Q₁) Is there an algorithm to compute the certificate and what is its (bit) complexity?

(Q₂) What is the (bit) size of the certificate?

(Q₃) What is the (bit) complexity of verifying the algebraic identity induced by the certificate, that is, to verify that the certificate is correct?

(Q₄) If there is no certificate, that is, if $A$ could be negative, then, can we compute a witness point such that the evaluation of $A$ at this point is negative? What is the cost of computing the witness point and what is its (bit) size?

In the case of a polynomial $A$ with real coefficients of even degree, say $d = 2m$, it is well known that $A$ is nonnegative over $\mathbb{R}$ if and only if it is a sum of two squares of polynomials; this answers (Q₀). This is a special version of a more general mathematical foundation based on Karlin's description of nonnegative polynomials in T-systems [8]; we will exploit this connection further in the sequel. In particular, we have the following equivalence.

$$A(x) \geq 0 \text{ for all } x \in \mathbb{R} \Leftrightarrow A = P^2 + Q^2, \quad (1)$$

where $A, P, Q \in \mathbb{R}[x]$. The equivalence in (1) serves as a certificate of positivity for $A$. Regarding its size, we notice that the polynomials $P$ and $Q$ are of degree at most $m = d/2$; the degree of $A$ must be even. Because (1) involves polynomials with real coefficients, it is not relevant to discuss bounds on the bitsize of the certificate. To compute the polynomials $P$ and $Q$, we need to compute and manipulate the roots of $A$, e.g., [28, Chapter 8]. Hence, the complexity of the decomposition is dominated by the root-finding algorithm we employ for this task; this requires $\tilde{O}(d)$ arithmetic operations; see e.g., [24] and references therein. We can verify the algebraic identity of the certificate in (1), $A = P^2 + Q^2$, and so we can answer (Q₃), either deterministically with direct computations or probabilistically, by evaluating the left- and right-hand side polynomials at random numbers, as in polynomial identity testing. Finally, for (Q₄) it suffices to return a number, say $t \in \mathbb{R}$, such that $A(t) < 0$. If $A$ can be negative, then it has real roots of odd multiplicity, or it is everywhere negative. For the former case, we can choose as a witness point a number $t$ lying to the left or to the right of a real root. If $A$ is globally negative, then any number $t \in \mathbb{R}$ suffices.

If $A$ has rational coefficients and we opt for an SOS decomposition with polynomials having rational coefficients, then things are somewhat more complicated and there is a (slightly) different certificate. Specifically, there is a weighted SOS representation, that is, a representation as sum of squares of polynomials with rational coefficients multiplied by positive rational numbers. Namely, the certificate is:

$$A(x) \geq 0 \text{ for all } x \in \mathbb{R} \Leftrightarrow A = \sum_{j=1}^{\nu} w_j \cdot s_j^2, \quad (2)$$

<footnote n="1">¹There are cases where the certificate might consists of several algebraic identities. This is not the case for the problem we consider, so we do not explore further this direction.</footnote>
where $A, s_j \in \mathbb{Q}[x], w_j \in \mathbb{Q}_{\ge 0}$ for all $j \in [\nu]$, and some $\nu \in \mathbb{N}$.

Pourcet [27], improving a previous result of Landau [21], proved that only five or less squares are needed, $\nu \le 5$; but his proof is not constructive. We refer the reader to [19] for recent advances in this direction.

At the cost of having more summands, there is the constructive approach by Schweighofer [33] that, roughly speaking, successively subtracts positive quadratic polynomials from $A$. There is also the algorithm by Chevillard, Harrison, Joldeş, and Lauter [6] that computes a weighted SOS decomposition at the expense of involving $\nu = d + 3$ summands, where $d$ is the degree of the polynomial. This is the algorithm that we focus on; we call it uSOS. The crux of the algorithm is that it performs a sufficiently small perturbation to $A$ and then approximates the complex roots of the perturbed polynomial. We introduce a slight variant of this algorithm, that we also call uSOS and we refer to Section 2.1 for a detailed presentation.

If we assume that $A$ has integer coefficients, $A \in \mathbb{Z}[x]$, and the maximum coefficient bitsize is $\tau$, then Magron, Safey El Din, and Schweighofer [22] studied the bit complexity of uSOS. They demonstrated that uSOS computes the certificate in (2) in $\tilde{O}_B(d^4 + d^3\tau)$ bit operations ($Q_1$). They estimated the bitsize of the certificate, based on uSOS, to be $\tilde{O}(d^3 + d^2\tau)$ ($Q_2$). In particular, they show that the certificate involves $O(d)$ coefficients and their bitsize is $\tilde{O}(d^2 + d\tau)$. Finally, they show how to verify the certificate in $\tilde{O}_B(d^4 + d^3\tau)$ bit operations ($Q_3$).

Along the same lines, Krick, Mourrain, and Szanto [20] studied a more general (semi-algebraic) problem. That is, they presented an algorithm to support the following equivalence: $A \in \mathbb{Z}[x]$ is nonnegative on all real roots of $B \in \mathbb{Z}[x]$ iff $A$ is SOS (of polynomials with rational coefficients) modulo $B$. The authors do not provide a bit complexity analysis, however, we refer the reader to [1] for recent improvements and generalizations.

We emphasize that the certificates of positivity imposed by (1) or (2) are not the only ones. There are also certificates based on Bernstein basis [4], (dual) certificates based on the dual cone of weighted sums of square polynomials [7] or based on sums of circuits [10], just to mention a few alternatives; we do not proceed further in these directions.

Finally, let us also mention the relation of nonnegative polynomials to T-systems. The theory of T-systems goes back a long time and it is highly developed. We refer the interested reader to [18, 8] for a detailed exposition. For our purposes, Karlin's work [17], especially the following consequence [17, Cor. 1], is of utmost importance: Let $A \in \mathbb{R}[x]$ with $\deg(A) = 2m$. Then, the following are equivalent:

(i) $A(x) > 0$ for all $x \in \mathbb{R}$.

(ii) There exist unique constants $\alpha, \beta > 0$ and unique points $x_1, \dots, x_m, y_1, \dots, y_{m-1} \in \mathbb{R}$, that we call Karlin points of $A$ over $\mathbb{R}$, with

$$x_1 < y_1 < \dots < y_{m-1} < x_m, \qquad (3)$$

such that

$$A(x) = \alpha \cdot (x - x_1)^2 \cdots (x - x_m)^2 + \beta \cdot (x - y_1)^2 \cdots (x - y_{m-1})^2.$$

Clearly, $\alpha$ is the leading coefficient of $A$. There are variants for $A(x) > 0$ or $A(x) \ge 0$ for all $x$ in $[a, b]$ or $[0, \infty)$; we refer the reader to [18, 8] for a comprehensive treatment. We also refer to Section B.3 for a brief introduction and some additional details. Unfortunately, the non-constructive approach of Karlin does not provide us with an algorithm to compute the unique Karlin points $x_i$ and $y_j$. Hence, T-systems were mainly of theoretical interests. We present a constructive approach in Section 4.

**Our contribution.** We revisit and slightly modify the algorithm by Chevillard, Harrison, Joldeş, and Lauter [6] for decomposing a polynomial $A \in \mathbb{Z}[x]$, of degree $d$ and bitsize $\tau$, that is positive over $\mathbb{R}$, as an SOS of at most $d + 3$ polynomials with rational coefficients; we call the variant uSOS. We complete and improve by a factor of $d$ the complexity analysis by Magron, Safey El Din, and Schweighofer [22].

The main idea of uSOS consists in perturbing the input polynomial $A$. The perturbed polynomial is $A_{\epsilon} := A - \epsilon M$, where $M = \sum_{k=0}^{d/2} x^{2k}$ is a positive polynomial and $\epsilon$ is the perturbation. Our choice of a small enough $\epsilon = 2^{-b}$ is such that two things happen: (i) $A_{\epsilon}$ is positive, this is the main requirement of
[6] (Lem. 2.1), and (ii) the roots of $A_ε$ are close to the roots of A (Lem. 2.2); this is a new requirement. The second requirement guarantees that the separation bounds, that is the minimum distance between two distinct roots, of $A_ε$ and A are almost the same, even though their bitsizes are different (Lem. 2.3). Also, the asymptotic values of $ε$ for both requirements are the same. Notice that the appropriate value of $ε$ does not depend on the actual polynomial A but only on its degree, $d$ and bitsize $τ$.

This choice of $ε$ allows to deduce that the bit complexity of approximating, in sufficient precision, the roots of $A_ε$ is asymptotically the same as the complexity of approximating the roots of A; this saves us a factor of $d$ for this step of the algorithm (Cor. 2.4). Then, we exploit the fan-in algorithm from approximate multipoint evaluation [25] to compute good approximate SOS decomposition for $A_ε$ and thus for A (Cor. 2.5 and Lem. 2.6). In this way, we obtain a bound $\tilde{O}_B(d^2τ)$ for the complexity of computing a certificate of positivity for A (Theorem 2.7). This improves the previously known bound by a factor of $d$ and answers (Q₁). Our complexity bound matches the bitsize of the certificate and hence we improve the complexity of the algorithms supporting the answers to (Q₂) (Lem. 2.8) and (Q₃) (Theorem 2.9).

The following theorem summarizes our results on representing a positive polynomial as a weighted sum of squares of polynomials.

**Theorem (uSOS and weighted SOS representation).** Let $A \in \mathbb{Z}[x]$ be a square-free polynomial of degree $d = 2m$ and maximum coefficient bitsize $τ$. If A is positive over $\mathbb{R}$, then there is an algorithm to compute a weighted SOS representation of A as

$$A(x) = \sum_{j=1}^{\nu} w_j s_j^2(x),$$

where $w_j \in \mathbb{Q}_{>0}$, $s_j \in \mathbb{Q}[x]$, and $\nu \le d+3$, at the cost of $\tilde{O}_B(d^3 + d^2τ)$ bit operations. The bitsize of $w_j$'s and the coefficients of $s_j$'s is at most $\tilde{O}(dτ)$.

Based on the previous theorem and by applying a transformation from Chevillard et al [6], we also provide bit complexity estimates for certificates and algorithms for the positivity of a polynomial over any interval (Theorem 2.14). If the bitsize of the endpoints of the interval is $σ$, then the bit complexity of the algorithm is $\tilde{O}_B(d^3 + d^2τ + d^2σ)$. To achieve this bound, it is not enough to apply the algorithm supported by the previous theorem directly; this will give us an extra $d$ factor in the complexity bound. We save this factor by studying how the separation bound of A changes after the transformation, see Sec. 2.4 for details.

To demonstrate the efficiency of uSOS we present an open-source prototype implementation in MAPLE and experiments on various data sets (Section 5). The experiments (Section 5.2) verify the bounds on the bitsize of the certificates and demonstrate the efficiency of the algorithm.

An important and challenging ingredient of an implementation of uSOS is the computation of the perturbation $ε := 2^{-b}$. Other approaches compute $ε$ by relying on the minimum of A or by repeatedly dividing by 2. These choices do not result in efficient implementations because they force us to compute with numbers having the worst case bitsize right from the beginning and/or the number of steps they perform depends on the bitsize of $ε$, that is $\tilde{O}(b) = \tilde{O}(dτ)$; the latter leads to an overall complexity that is quadratic in the bitsize of the input, that is $\tilde{O}_B(τ^2)$. A natural question to ask is if we can come up with an efficient implementation to compute $ε$ that does not dominate or alter the overall worst case bit complexity of the algorithm. We compute $ε$ by performing exponential binary search to obtain b; this requires $O(\lg(b)) = O(\lg(dτ))$ steps. In this way, we obtain enormous speedups in the running times of our implementation of uSOS, while at the same time guarantee the best theoretical worst case bit complexity bounds. We refer to Section 5.1 for further details and an experimental evaluation.

**Perturbed SOS certificates.** Although the uSOS algorithm provides efficient nonnegativity certificates, it requires to approximate the roots of the auxiliary polynomial $A_ε$ with precision $O(dτ)$ bits and produces certificates consisting of $O(d)$ SOS summands. To overcome these limitations, for polynomials with rational coefficients, we introduce perturbed SOS certificates, which certify nonnegativity of a polynomial $A \in \mathbb{Z}[x]$ by constructing a rational SOS approximation B that is sufficiently close to A in the sup-norm. Specifically, if $||A - B||_\infty < 2^{-b}$ with $b > \tilde{O}(dτ)$, then the nonnegativity of B implies the nonnegativity of A (Thm. 3.1). For the case of square-free positive polynomials of even degree, the certificate has the simple form $B := \tilde{P}^2 + \tilde{Q}^2$,
where the polynomials $\tilde{P}, \tilde{Q} \in \mathbb{Q}[x]$ are obtained from rational approximations of the roots of $A$. The bit complexity of computing such a certificate is $\tilde{O}_B(d^3 + d^2\tau)$, while the bitsize of $\tilde{P}, \tilde{Q}$ is bounded by $\tilde{O}(d\tau)$ (Thm. 3.5). This establishes perturbed SOS certificates as a refinement of uSOS that yields minimal-size representations while retaining fully explicit complexity and precision guarantees.

**T-systems and effective Karlin points.** Finally, we study structural properties of the SOS representation (Sec. 4). We establish a mathematical and algorithmic connection of the certificates of positivity of univariate polynomials with the T-systems, introduced by Karlin [18]. First, we consider the case a real polynomial, $A \in \mathbb{R}[x]$, that is positive over $\mathbb{R}$ (Sec. 4.1); hence it admits a representation as $A = P^2+Q^2$. The polynomials $P$ and $Q$ are interlacing (Lem. 4.1), as a direct consequence of Hermite-Bielher theorem (Thm. A.1) and thus their roots are the unique points, that we call Karlin points, Eq. (3), appearing in the decomposition induced by the T-systems of Karlin (Cor. 4.3). Even though, this seems to be a fundamental property, we were not able to find it in the literature. Besides on relying on Hermite-Bielher theorem [30] for the interlacing property, we give an alternative proof (Sec.A) based on algebraic manipulations of the roots of the polynomials. We refer to Fisk's survey [14] for a thorough study of real rooted and interlacing polynomials. The following theorem summarizes our results:

**Theorem (Positivity over $\mathbb{R}$, interlacing, and Karlin points).** If $A = \sum_{k=0}^d a_k x^k \in \mathbb{R}[x]$, of degree $d = 2m$, is square-free and positive over $\mathbb{R}$, then $A(x) = a_d P(x)^2 + a_d Q(x)^2$, where $P, Q \in \mathbb{R}[x]$ are interlacing of degrees $m$ and $m-1$, respectively. Moreover, the real roots of $P$ and $Q$ are the Karlin points of $A$ over $\mathbb{R}$.

We also present variants of the previous theorem for positivity over $(0, \infty)$, or any interval of $\mathbb{R}$, we discover the corresponding interlacing polynomials and we show that their real roots are the Karlin points of $A$ in the interval of interest; see Sec. 4.2.

**Organization** First, we present a detailed description of the various steps of the uSOS algorithm by Chevillard, Harrison, Joldeş, and Lauter [6] (Section 2.1). Then, Section 2.2 studies the bit complexity of the steps and the bitsize of the various quantities involved in the computations. In Section 2.3 we present the overall complexity of the algorithm and the certificate, and in Section 2.4 we consider nonnegativity over an interval. Section 5 presents our implementation and experiments. In Section 4 we establish the connection of positive real polynomials, SOS decompositions, T-systems, and interlacing polynomials. For a brief background on T-systems we refer to Section B.3 and [17, 8].

Finally, in the Appendix, we give and alternative proof of Lemma 4.1 (Sec. A) and we present auxiliary results on separation bounds and root approximation of univariate polynomials (Sec. B.1), bounds on the minimum of a polynomial and an approximation variant of the fan-in algorithm from multipoint evaluation (Sec. B.2). We also present a bird's eye view of T-systems (Section B.3).

**Notation** We denote by $\mathcal{O}$, resp. $\mathcal{O}_B$, the arithmetic, respectively bit, complexity; we also use $\tilde{\mathcal{O}}$, resp. $\tilde{\mathcal{O}}_B$, to ignore (poly-)logarithmic factors. For a polynomial $A = \sum_{k=0}^d a_k x^k \in \mathbb{C}[x]$ of degree $d$ we denote by $\|A\|_1$ resp. $\|A\|_\infty$, the one resp. the infinity, norm of the vector $(a_0, \dots, a_d)$. We denote by $\text{lct}(A) = a_d$, resp. $\text{tc}(A) = a_0$, the leading, resp. tailing, coefficient of $A$. If $A \in \mathbb{Q}[x]$, then the bitsize of $A$ is the maximum bitsize of its coefficients, including a bit for the sign. For $a \in \mathbb{Q}$, its bitsize is the maximum bitsize of the numerator and the denominator. If $A \in \mathbb{Q}[x]$ has degree $d$ and bitsize $\tau$, then we also say that $A$ has size $(d, \tau)$. We write $\Delta_\alpha(A)$ or just $\Delta_\alpha$ to denote the minimum distance between a root $\alpha$ of $A$ and any other root. We call this quantity *local separation bound*. We also write $\Delta_i$ instead of $\Delta_{\alpha_i}$. Also $\Delta(A) = \min_\alpha \Delta_\alpha(A)$ or just $\Delta$ denotes the *separation bound*, that is the minimum distance between all the roots of $A$. Finally, let $D(c, r) = \{x \in \mathbb{C} : |x - c| \le r\}$. Given a complex number $z \in \mathbb{C}$ such that $|z| < 2^\tau$, we say that $\tilde{z}$ is an approximation up to an absolute precision $\lambda \in \mathbb{N}$, if $|z - \tilde{z}| < 2^{-\lambda}$. Then, the bitsize of the approximation is at most $\tau + 2\lambda + 2$ and we can represent it as a dyadic fraction of the form $a2^{-b}$, for $a \in \mathbb{Z}$ and $b \in \mathbb{N}$.

We should note the constants in the various bounds we present are not the best possible. A more detailed analysis can improve them. We decided to present them, even in their rough form, to demonstrate that there are no hidden non-constant factors in the $\mathcal{O}$ notation of the bounds.
## 2 The uSOS algorithm

Consider the following polynomial that has even degree $d = 2m$ and is positive over $\mathbb{R}$:

$$A(x) = \sum_{k=0}^{d} a_k \cdot x^k = a_d \prod_{i=1}^{d} (x - \alpha_i).$$

Our goal is to provide a representation of $A$ as a weighted sum of squares of polynomials.

We assume that $A$ is square-free and has no real roots, that is $\alpha_i \notin \mathbb{R}$, for all $i \in [d]$. These assumptions are without loss of generality; we refer to Sec. 2.3.1 for a detailed discussion.

To simplify various calculations in the sequel, we also need to require the leading coefficient of $A$ to be such that

$$\frac{1}{2} \leq \mathrm{lc}(A) = a_d \leq 1,$$

and that all the other coefficients of $A$ are rational numbers of the same denominator and of bitsize bounded by $\tau$. Hence, when the input is a polynomial $A$ with integer coefficients, first, we multiply $A$ with a rational number in the interval $[1/(2a_d), 1/a_d]$ to ensure this condition. This operation (or requirement) does not change neither the positivity of $A$ nor the complexity bound, hence we will assume it in our analysis.

We present in detail (Sec. 2.1) the various steps of the uSOS algorithm by Chevillard, Harrison, Joldeş, and Lauter [6] that decomposes $A$ as a weighted sum of squares of polynomials with rational coefficients. The presentation leads to the precise bit complexity analysis of Sec. 2.3.

### 2.1 A detailed presentation of the uSOS algorithm

**Input:** A polynomial $A = \sum_{k=0}^{d} a_k x^k = a_d \prod_{i=1}^{d} (x - \alpha_i) \in \mathbb{Q}[x]$ of even degree, $d = 2m$.

**Assumptions:** (i) $A$ is square-free, (ii) $A$ is positive over $\mathbb{R}$, and (iii) The leading coefficient of $A$ is in $[\frac{1}{2}, 1]$, while the other coefficients are rationals of bitsize at most $\tau$, having a common denominator.

**Output:** A weighted SOS decomposition of $A$, that is a representation of the form

$$A(x) = \sum_{j=1}^{\nu} w_j s_j^2(x), \qquad (4)$$

where $w_j \in \mathbb{Q}_{>0}$ and $s_j \in \mathbb{Q}[x]$. It holds $\nu \leq d + 3$.

**[Step 1] Rewrite $A(x)$ using $M(x)$ and $\varepsilon$ [6, Sec. 5.2.2]**

Consider the polynomial

$$M(x) = \sum_{j=0}^{m} x^{2j},$$

that is the sum of even powers $x^{2j}$ less than or equal to $d$. Notice that $M(t) > 0$ for all $t \in \mathbb{R}$. Write $A$ as

$$A(x) = A(x) - \varepsilon M(x) + \varepsilon M(x) = A_{\varepsilon}(x) + \varepsilon M(x),$$

where $\varepsilon > 0$. In particular, $\varepsilon$ should be small enough so that $A_{\varepsilon} := A - \varepsilon M$ is strictly positive over $\mathbb{R}$. In Sec. 2.2 we estimate a precise value for $\varepsilon$.

**[Step 2] Approximate the roots of $A_{\varepsilon}$ and compute $\tilde{P}_{\varepsilon}$ and $\tilde{Q}_{\varepsilon}$ [6, Sec. 5.2.3]**

The polynomial $A_{\varepsilon}$ is strictly positive and has no real roots; let its factorization to linear factors be

$$A_{\varepsilon}(x) = \mathrm{lc}(A_{\varepsilon}) \cdot \prod_{i=1}^{d} (x - \alpha_{\varepsilon,i}), \quad \text{where } \alpha_{\varepsilon,i} \in \mathbb{C} \text{ and } a_{\varepsilon,d} := \mathrm{lc}(A_{\varepsilon}) > 0.$$
We can approximate the roots of $A_ε$, $α_{ε,i}$, with rationals, up to any desired accuracy, say $2^{-κ}$, for some positive integer $κ$. Let the approximations be

$$\tilde{\alpha}_{\varepsilon,j}^{\pm} = \tilde{\gamma}_{\varepsilon,j} \pm i \tilde{\delta}_{\varepsilon,j},$$

where $\tilde{\gamma}_{\varepsilon,j}, \tilde{\delta}_{\varepsilon,j} \in \mathbb{Q}$ and $j \in [m]$. Then, it holds

$$|\alpha_{\varepsilon,i}^{\pm} - \tilde{\alpha}_{\varepsilon,i}^{\pm}| \le 2^{-\kappa}.$$

In turn, the rational approximations of the roots lead to a polynomial $\tilde{A}_ε$ with rational coefficients, that is

$$\tilde{A}_{\varepsilon}(x) = a_{\varepsilon,d} \prod_{j=1}^{m} (x - \tilde{\alpha}_{\varepsilon,j}^{+})(x - \tilde{\alpha}_{\varepsilon,j}^{-}) = a_{\varepsilon,d} \prod_{j=1}^{m} (x - \tilde{\gamma}_{\varepsilon,j} + i \tilde{\delta}_{\varepsilon,j}) \prod_{j=1}^{m} (x - \tilde{\gamma}_{\varepsilon,j} - i \tilde{\delta}_{\varepsilon,j}),$$

where we can additionally assume that $\tilde{\delta}_{\varepsilon,j} \ge 0$, for all $j \in [m]$. Moreover,

$$\prod_{j=1}^{m} (x - \tilde{\gamma}_{\varepsilon,j} + i \tilde{\delta}_{\varepsilon,j}) = \tilde{P}_{\varepsilon}(x) + i \tilde{Q}_{\varepsilon}(x) \quad \text{and} \quad \prod_{j=1}^{m} (x - \tilde{\gamma}_{\varepsilon,j} - i \tilde{\delta}_{\varepsilon,j}) = \tilde{P}_{\varepsilon}(x) - i \tilde{Q}_{\varepsilon}(x),$$

which implies to the following representation of $\tilde{A}_ε$:

$$\tilde{A}_{\varepsilon}(x) = a_{\varepsilon,d} (\tilde{P}_{\varepsilon}(x) + i \tilde{Q}_{\varepsilon}(x)) (\tilde{P}_{\varepsilon}(x) - i \tilde{Q}_{\varepsilon}(x)) = a_{\varepsilon,d} (\tilde{P}_{\varepsilon}(x)^2 + \tilde{Q}_{\varepsilon}(x)^2).$$

Therefore, as the numbers $\tilde{\alpha}_{\varepsilon,j}^{\pm} = \tilde{\gamma}_{\varepsilon,j} \pm i \tilde{\delta}_{\varepsilon,j}$ approximate the roots $α_{ε,i}$, we also deduce that

$$\tilde{A}_{\varepsilon} = a_{\varepsilon,d}(\tilde{P}_{\varepsilon}^2 + \tilde{Q}_{\varepsilon}^2)$$

approximates the polynomial $A_ε$. The two polynomials $A_ε$ and $\tilde{A}_ε$ have the same leading coefficient. Let their difference be

$$A_{\varepsilon}(x) - \tilde{A}_{\varepsilon}(x) := B(x) = \sum_{k=0}^{d-1} b_k x^k.$$

In this way, we obtain the following relation for $A(x)$:

$$A(x) = A_{\varepsilon}(x) + \varepsilon M(x) = \tilde{A}_{\varepsilon}(x) + B(x) + \varepsilon M(x) = a_{\varepsilon,d} \tilde{P}_{\varepsilon}(x)^2 + a_{\varepsilon,d} \tilde{Q}_{\varepsilon}(x)^2 + B(x) + \varepsilon M(x). \quad (5)$$

### [Step 3] Write $B(x) + \varepsilon M(x)$ as SOS [6, Sec.5.2.5]

As the first two summands of (5) are weighted sum of squares, to represent $A$ as a weighted SOS, we should express $B(x) + \varepsilon M(x)$ as a weighted SOS. For this, we exploit the identities

$$x = (x + \frac{1}{2})^2 - (x^2 + \frac{1}{4}) \quad \text{and} \quad -x = (x - \frac{1}{2})^2 - (x^2 + \frac{1}{4}).$$

In this way, for any $c > 0$, we have

$$\pm c x^{2k+1} = c \left(x^{k+1} \pm \frac{x^k}{2}\right)^2 - c \left(x^{2k+2} + \frac{x^{2k}}{4}\right).$$

Using these identities, the odd-degree terms of $B(x) = \sum_{k=0}^{d-1} b_k x^k = A_ε(x) - \tilde{A}_ε(x)$ become

$$b_{2k+1} x^{2k+1} = |b_{2k+1}| \left(x^{k+1} + \operatorname{sgn}(b_{2k+1}) \frac{x^k}{2}\right)^2 - |b_{2k+1}| \left(x^{2k+2} + \frac{x^{2k}}{4}\right)$$
and consequently

$$B(x) + \varepsilon M(x) = \sum_{k=0}^{m-1} \underbrace{|b_{2k+1}| \left(x^{k+1} + \operatorname{sgn}(b_{2k+1}) \frac{x^k}{2}\right)^2}_{w_k} + \sum_{k=0}^{m} \underbrace{\left(\varepsilon + b_{2k} - |b_{2k-1}| - \frac{1}{4}|b_{2k+1}|\right)}_{w_{m+k}} x^{2k}, \quad (6)$$

where by convention $b_{-1} = b_{2m+1} = 0$. If, for every $k$, it holds

$$\varepsilon \ge \frac{1}{4} |b_{2k+1}| - b_{2k} + |b_{2k-1}|, \quad (7)$$

then (6) is a weighted SOS representation. For the inequalities in (7) to hold, we should approximate the roots of $A_{\varepsilon}(x) = A(x) - \varepsilon M(x)$ with enough precision, say $\kappa$, so that the polynomial $\tilde{A}_{\varepsilon}$ is close to $A_{\varepsilon}$, thus their different $B := A_{\varepsilon} - \tilde{A}_{\varepsilon}$ is small and so the coefficients $b_k$ are small (compared to $\varepsilon$).

## 2.2 The bit complexity of the various steps

We estimate the (bit)size of the various quantities appearing in the process of the uSOS and the complexity of the various operations. Along the way we estimate the value of $\varepsilon$ that suffices to perturb the original polynomial $A$ and the precision, $\kappa$, that we need to approximate the roots of the perturbed polynomial $A_{\varepsilon}$. We express both as a function of $d$ and $\tau$. Recall, that we assume that $A$ is a square-free polynomial of degree $d$, positive over $\mathbb{R}$, and its leading coefficient is in $[\frac{1}{2}, 1]$, while the other coefficients are rationals of bitsize at most $\tau$ having a common denominator.

### [Step 1] Rewrite $A(x)$ using $M(x)$ and $\varepsilon$

We estimate a suitable small value for $\varepsilon = 2^{-b}$ to ensure that the polynomial $A_{\varepsilon}(x) := A(x) - \varepsilon M(x)$ is positive for every $x \in \mathbb{R}$. We assume

$$0 < \varepsilon \le \frac{1}{8} \le \frac{a_d}{4}.$$

Then, the leading coefficient of $A_{\varepsilon}$ is $\operatorname{lc}(A_{\varepsilon}) = a_d - \varepsilon > 0$. Our analysis of this step follows closely [22].

An upper bound on the magnitude of the (real) roots of $A_{\varepsilon}$ [11, Theorem 1] is

$$1 \le R := 2 \frac{\|A_{\varepsilon}\|_{\infty}}{|\operatorname{lc}(A_{\varepsilon})|} \le 2^{\tau+5}. \quad (8)$$

Then, for any $x \in \mathbb{R}$ such that $|x| \ge R$, we have $A_{\varepsilon}(x) > 0$ (and $A(x) > 0$) as the leading coefficient of $A_{\varepsilon}$ is positive. If we choose $\varepsilon$ such that

$$\varepsilon = 2^{-b} \le \frac{\min_{|x| \le R} A(x)}{\max_{|x| \le R} M(x)},$$

then we ensure that $A_{\varepsilon}$ is positive for all $|x| \le R$.

Regarding $M(x)$, 0 is the only real root of its derivative, thus $\max_{|x| \le R} M(x) = \max\{1, M(R)\} = M(R)$. It holds

$$\max_{|x| \le R} M(x) = M(R) = \sum_{k=0}^{d/2} R^k \le (\frac{d}{2} + 1) R^{d/2} \le (\frac{d}{2} + 1) 2^{d(\tau+5)/2}.$$

The global minimum of $A$ is reached at a critical value, that is the evaluation of $A$ at a root of its derivative, and, for obtaining worst case bounds, we can assume that is in the interval $(-R, R)$. So, Lemma B.3 implies

$$2^{-4d\tau-16d\lg d} \le \min_{|x| \le R} A(x) \le 2^{2d\tau+8d\lg d}.$$

Overall, to ensure that $A_{\varepsilon}$ is positive, we choose $\varepsilon = 2^{-b}$, where

$$b \ge 4d\tau + 8d \lg d. \quad (9)$$

**Lemma 2.1.** If $\varepsilon = 2^{-b}$, with $b \ge 4d\tau + 8d \lg d$, then $A_{\varepsilon}(x) := A(x) - \varepsilon M(x) > 0$ for all $x \in \mathbb{R}$.
[Step 2] Approximate the roots of $A_\epsilon$ and compute $\tilde{P}_\epsilon$ and $\tilde{Q}_\epsilon$

We approximate the roots of $A_\epsilon$ and then, using the approximations, we construct the polynomials $\tilde{P}_\epsilon$ and $\tilde{Q}_\epsilon$. To approximate the roots of $A_\epsilon$ with rationals up to any desired precision, say $2^{-\kappa}$, for a positive integer $\kappa$, first, we need to isolate them and then approximate them to any desired precision. We use well known algorithms for this task, the main ingredients of which are the splitting circle method and (variants of) the Newton operator, e.g., [24, 23] and references therein. The complexity of the rootfinding algorithms (mainly) depends on the (aggregate) separation bound of the roots of $A_\epsilon$; that is the minimum distance between the roots, e.g., [11] and references therein. We choose an $\epsilon := 2^{-\kappa}$ small enough so that the separation bounds of $A$ and $A_\epsilon$ are similar; this allows to bound the complexity of approximating the roots of $A_\epsilon$ in terms of the separation of $A$ and save a factor of $d$ in the overall complexity and the bitsize of the certificate.

We proceed as follows: first we compute the suitable value for $\epsilon$, then we bound the complexity of approximating the roots of $A_\epsilon$ up to precision $\kappa$, and, finally, we estimate the complexity of computing $\tilde{P}_\epsilon$ and $\tilde{Q}_\epsilon$ (as function of $d$, $\tau$, and $\kappa$).

**The separation bound of $A_\epsilon$.** We relate the separation bound of $A_\epsilon$ in terms of the separation bound of $A$, using Lemma B.1. In our case, $A$ plays the role of $p$ and $A_\epsilon$ plays the role of $\tilde{p}$. Notice that the leading coefficients of $A$ and $A_\epsilon$ are different, they are $a_d$ and $a_d - \epsilon$, respectively. Thus, to apply Theorem B.1 we need to consider the one norm of the following difference

$$ 
\begin{aligned}
\| \frac{a_d - \epsilon}{a_d} A - A_\epsilon \|_1 &\le \| \frac{a_d - \epsilon}{a_d} A - A + \epsilon M \|_1 & (A_\epsilon = A - \epsilon M) \\
&\le \epsilon \| \frac{1}{a_d} A - M \|_1 & (\|M\|_1 \le \frac{d}{2} + 1) \\
&= 16\epsilon d \frac{1}{8} \|A\|_1 & (\epsilon \le \frac{a_d}{4}) \\
&\le 2^{-b+\lg d+4} \cdot \| \frac{a_d - \epsilon}{a_d} A \|_1 & (\epsilon = 2^{-b})
\end{aligned}
$$

Notice that the polynomials $A$ and $\frac{a_d - \epsilon}{a_d} A$ have the same root and separation bounds. Then, Lem. B.1, for sufficiently small $\epsilon$ or equivalently for sufficiently big $b$, implies that the separation bounds of $A$ and $A_\epsilon$ are related with small constant depending on the degree, see (23). In particular, we should choose a positive integer $b$ that satisfies the following three conditions that correspond to eqs. (20) to (22) of Lem. B.1:

(i) $b \ge \lg d + 4 + \max\{8d, d \lg d\}$ and $b$ is a power of two.

(ii) $2^{(-b+\lg d+4)/2} \le \frac{\Delta_i(A)}{2d}$, or equivalently, $b \ge -2\lg \Delta_i(A) + 3\lg(d) + 6$. Using [11, Theorem 1] we can bound $\Delta_i(A)$ to obtain $b \ge 2d\tau + 12d \lg d + 8$.

(iii) $2^{-b/2} \le \frac{\prod_{j \ne i} (\alpha_i - \alpha_j)}{16 \cdot (d+1) \cdot 2^{\tau_A} \cdot M(\alpha_i)^d}$. It holds $M(\alpha_i) = \max\{1, |\alpha_i|\} \le 2^{\tau+3}$, using an upper bound on the roots of $A$, e.g., [11, Theorem 1]. Also using Lem. B.3, we obtain

$$ P_i := \prod_{j \ne i} |\alpha_i - \alpha_j| = \frac{1}{d a_d} |A'(\alpha_i)| \ge 2^{-3d\tau - 3d \lg d - \lg d}. \quad (10) $$

These bounds lead to the inequality $b \ge 5d\tau + 9d \lg d + 12$.

Therefore, if we choose $b$ such that

$$ b \ge 5d\tau + 9d \lg d + 12, \quad (11) $$

then (i), (ii), and (iii) are simultaneously satisfied. In this case, by Lem. B.1, Eq. (23), for any $i \in [d]$,

$$ (1 - \frac{1}{d}) \Delta_i(A) \le \Delta_i(A_\epsilon) \le (1 + \frac{1}{d}) \Delta_i(A). \quad (12) $$

The previous discussion leads to the following lemma:
**Lemma 2.2.** If $\varepsilon = 2^{-b}$, with $b \ge 5d\tau + 9d \lg d + 12$, then (12) holds.

By combining Lemmata 2.1 and 2.2 and considering all the roots of $A_\varepsilon$, we have the following

**Lemma 2.3** (The value of $\varepsilon$). If $\varepsilon = 2^{-b}$, with $b \ge 5d\tau + 9d \lg d + 12$, then $A_\varepsilon(x) = A(x) - \varepsilon M(x) > 0$ for all $x \in \mathbb{R}$, and

$$ (1 - \frac{1}{d})^d \prod_{i=1}^{d} \Delta_i(A) \le \prod_{i=1}^{d} \Delta_i(A_\varepsilon) \le (1 + \frac{1}{d})^d \prod_{i=1}^{d} \Delta_i(A). \quad (13) $$

**Bit complexity of approximating the roots of $A_\varepsilon$.** To approximate the roots of $A_\varepsilon$, up to precision $2^{-\kappa}$, for a positive integer $\kappa$, we employ the algorithm supported by Theorem B.2 that computes complex numbers $\tilde{\alpha}_{\varepsilon,j} \in \mathbb{Q}[i]$, where $j \in [m]$, such that

$$ \|A_\varepsilon - a_{\varepsilon,d} \prod_{j=1}^{m} (x - \tilde{\alpha}_{\varepsilon,j}^{+})(x - \tilde{\alpha}_{\varepsilon,j}^{-})\| \le 2^{-b} \|A_\varepsilon\|. $$

We let

$$ \tilde{A}_\varepsilon(x) := a_{\varepsilon,d} \prod_{j=1}^{m} (x - \tilde{\alpha}_{\varepsilon,j}^{+})(x - \tilde{\alpha}_{\varepsilon,j}^{-}). $$

The rootfinding algorithm returns the real and imaginary part of the $\tilde{\alpha}_{\varepsilon,j}^{\pm} = \tilde{\gamma}_{\varepsilon,j} \pm i\tilde{\delta}_{\varepsilon,j}$ as dyadic fractions of the form $a 2^{-b}$, where $a \in \mathbb{Z}$ and $b \in \mathbb{N}$; all fractions have the same denominator. The next lemma bounds the bit complexity of approximating the roots of $A_\varepsilon$, both the real and imaginary parts, up to precision $2^{-\kappa}$, as a function of $d$, $\tau$, and $\kappa$. We estimate the value of $\kappa$ in the next subsection.

**Corollary 2.4.** If $\varepsilon = 2^{-b}$, with $b \ge 5d\tau + 9d \lg d + 12$ (that is as in Lemma 2.3), then one can compute rational approximations (of the real and imaginary part) of the roots of $A_\varepsilon$, up to precision $2^{-\kappa}$, for a positive integer $\kappa$, at the cost of $\tilde{O}_B(d^3 + d^2\tau + d\kappa)$ bit operations.

*Proof.* We bound the various quantities appearing in the complexity bounds of Theorem B.2.

First, we estimate the bound on the coefficients, $\tau_{A_\varepsilon} \le \lfloor \frac{\|A_\varepsilon\|_\infty}{1c(A_\varepsilon)} \rfloor \le \tau + 5$.

To bound the (aggregate) separation bound of the roots of $A_\varepsilon$, we employ Eq. (13). Then,

$$ - \lg \prod_{i=1}^{d} \Delta_i(A_\varepsilon) \le - \lg(1 - \frac{1}{d})^d \prod_{i=1}^{d} \Delta_i(A) \le 2 - \lg \prod_{i=1}^{d} \Delta_i(A) = \tilde{O}(d\tau), $$

where for the last equality we refer to [11, Theorem 1].

It remains to bound $\tilde{P}_i := \prod_{j \ne i} |\tilde{\alpha}_{\varepsilon,i} - \tilde{\alpha}_{\varepsilon,j}| = \frac{1}{d a_{\varepsilon,d}} |A'_\varepsilon(\tilde{\alpha}_{\varepsilon,i})|$. It holds $\frac{P_i}{2} \le \tilde{P}_i \le 2P_i$ [23, proof of Theorem 4] and so

$$ - \lg \tilde{P}_i \le - \lg \frac{1}{2} - \lg P_i \le 3d\tau + 3d \lg d + \lg d = \tilde{O}(d\tau), $$

where the last inequality is due to (10). By combining all the bounds, we conclude the proof. $\square$

**The cost of constructing the polynomials $\tilde{P}_\varepsilon$ and $\tilde{Q}_\varepsilon$.** It remains to actually compute the polynomials $\tilde{A}_\varepsilon$, $\tilde{P}_\varepsilon$, and $\tilde{Q}_\varepsilon$ from the approximations $\tilde{\alpha}_{\varepsilon,j}$'s, based on Lemma B.4. It holds

$$ |\tilde{\alpha}_{\varepsilon,j} - \alpha_{\varepsilon,i}| \le 2^{-\kappa}. $$

The fan-in algorithm, supported by Lemma B.4, computes the polynomial

$$ \tilde{A}_{\varepsilon}(x) = 1c(A_{\varepsilon}) \prod_{j=1}^{m} (x - \tilde{\alpha}_{\varepsilon,j}^{+})(x - \tilde{\alpha}_{\varepsilon,j}^{-}), $$
and the polynomial $\tilde{P}_{\varepsilon}$ and $\tilde{Q}_{\varepsilon}$ (and $\tilde{A}_{\varepsilon}$) from the product

$$ \prod_{i=1}^{m} (x - \tilde{\gamma}_{\varepsilon,j} + i \tilde{\delta}_{\varepsilon,j}) = \tilde{P}_{\varepsilon}(x) + i \tilde{Q}_{\varepsilon}(x). $$

Recall, that $A_{\varepsilon}$ is a positive polynomial, so it admits a representation $A_{\varepsilon} = P_{\varepsilon}^2 + Q_{\varepsilon}^2$. As we approximate the roots of $A_{\varepsilon}$, we compute approximations $\tilde{P}_{\varepsilon}$ and $\tilde{Q}_{\varepsilon}$, so that, after simplifications,

$$ \|P_{\varepsilon} - \tilde{P}_{\varepsilon}\|_{\infty} \le 2^{-\kappa+(4m-4)(\tau+5)+32m-(\lg m+5)^2-7} \le 2^{-\kappa+2d\tau+26d-\lg^2 d-30}, $$

and similarly for $Q_{\varepsilon}$ and $\tilde{Q}_{\varepsilon}$. The cost is $\tilde{O}_B(d(\kappa + d\tau))$ bit operations, it holds $\|\tilde{P}_{\varepsilon}\|_{\infty}, \|\tilde{Q}_{\varepsilon}\|_{\infty} \le 2^{d\tau+4d}$, and the bitsize of all three polynomials is at most $\mathcal{O}(\kappa + d\tau)$.

**Corollary 2.5.** Assume $\varepsilon = 2^{-b}$, with $b \ge 5d\tau + 9d \lg d + 12$ (that is as in Lemma 2.3). If we are given a rational approximation of the roots of $A_{\varepsilon}$, with precision $2^{-\kappa}$, for a positive integer $\kappa$, then we can compute the polynomials $\tilde{P}_{\varepsilon}$ and $\tilde{Q}_{\varepsilon}$ at the cost of $\tilde{O}_B(d^2\tau + d\kappa)$ bit operations. The bitsize of the polynomials is $\tilde{O}(\kappa + d\tau)$.

### [Step 3] Write $B(x) + \varepsilon M(x)$ as SOS

Now, we have computed rational approximations of the roots of $A_{\varepsilon}$ such that $|\alpha_{\varepsilon,i} - \tilde{\alpha}_{\varepsilon,i}| \le 2^{-\kappa}$, that correspond to the polynomial $\tilde{A}_{\varepsilon}$. Based on Lem. B.4 we deduce that

$$ \|A_{\varepsilon} - \tilde{A}_{\varepsilon}\|_{\infty} \le 2^{-\kappa+2d\tau+26d-\lg^2 d-30}, $$

and, since $B := A_{\varepsilon} - \tilde{A}_{\varepsilon}$, it holds

$$ \|B\|_{\infty} \le 2^{-\kappa+2d\tau+26d-\lg^2 d-30}. $$

To satisfy the inequality (7), that is $\varepsilon \ge |b_{2k+1}|/4 - b_{2k} + |b_{2k-1}|$, the following inequality needs to hold

$$ \varepsilon = 2^{-b} \ge 2^{-\kappa+2d\tau+26d-\lg^2 d-30} \Rightarrow \kappa \ge 5d\tau + 40d \lg d, \quad (14) $$

where we also use the bound on $b$ from Lemma 2.3. Consequently, $\varepsilon$ and all the coefficients $b_k$ of $B$ have bitsize $\tilde{O}(d\tau)$. This leads to the following lemma:

**Lemma 2.6.** Let $\varepsilon = 2^{-b}$, with $b \ge 5d\tau + 9d \lg d + 12$ (Lem. 2.3). Then, if $\kappa \ge 5d\tau + 40d \lg d$, that is $\kappa = \tilde{O}(d\tau)$, then (7) holds for all $k$. The bitsize of the positive rationals $w_j$ in (6) is also $\tilde{O}(d\tau)$.

## 2.3 Overall complexity estimates

The previous two sections imply that for a given polynomial $A \in \mathbb{Z}[x]$, of degree $d = 2m$ and bitsize $\tau$, that is positive over $\mathbb{R}$, uSOS computes a representation of $A$ as a weighted SOS of polynomials with rational coefficients, as in (2). In particular, it represents $A$ as

$$ A(x) = a_{\varepsilon,d} \tilde{P}_{\varepsilon}(x)^2 + a_{\varepsilon,d} \tilde{Q}_{\varepsilon}(x)^2 + \sum_{k=0}^{m-1} w_k \left(x^{k+1} \pm \frac{x^k}{2}\right)^2 + \sum_{k=0}^{m} w_{m+k} x^{2k}, \quad (15) $$

where $a_{\varepsilon,d} = a_d - \varepsilon \in \mathbb{Q}_{\ge 0}$, $\varepsilon \in (0, \frac{1}{8})$, $\tilde{P}_{\varepsilon}, \tilde{Q}_{\varepsilon} \in \mathbb{Q}[x]$, and $w_j \in \mathbb{Q}_{\ge 0}$, for $0 \le j \le d$.

**Theorem 2.7 (Bit complexity of uSOS).** Let $A \in \mathbb{Z}[x]$ be a square-free polynomial of degree $d = 2m$ and maximum coefficient bitsize $\tau$. If $A$ is positive over $\mathbb{R}$, then uSOS computes a weighted SOS representation of $A$ as in (15), see also (2), at the cost of $\tilde{O}_B(d^3 + d^2\tau)$ bit operations.
*Proof.* The first step of the uSOS involves computing the polynomial $A_ε$. This requires $O(d)$ additions of numbers of bitsize $\tilde{O}(d\tau)$ (Lemma 2.3). So, the bit complexity is $\tilde{O}_B(d^2\tau)$.

The second step requires to approximate the roots of $A_ε$ up to precision $κ$. As $κ = O(d\tau)$ (Lemma 2.6), this costs $\tilde{O}_B(d^3 + d^2\tau)$ bit operations (Cor. 2.4). This also includes the cost of computing the polynomials $\tilde{P}_ε$ and $\tilde{Q}_ε$ (Cor. 2.5).

The last step requires $O(d)$ additions of numbers of bitsize $\tilde{O}(d\tau)$ to construct the coefficients $w_j$.

Hence, the overall bit complexity is $\tilde{O}_B(d^3 + d^2\tau)$. □

**Lemma 2.8** (Bitsize of the certificate). *The representation in (15) involves at most $d+3$ summands. The bitsize of the rationals involved in the representation is at most $\tilde{O}(d\tau)$, while their total bitsize is $\tilde{O}(d^2\tau)$.*

*Proof.* The number of summands follows directly from the algorithm [6] and (15).

We notice that $ε, a_{ε,d}$ (Lemma 2.3), the coefficients of $\tilde{P}_ε$ and $\tilde{Q}_ε$ (Cor. 2.5), and the rationals $w_j$ (Lemma 2.6) all have bitsize $\tilde{O}(d\tau)$.

The polynomials $\tilde{P}_ε$ and $\tilde{Q}_ε$ have at most $d+2$ coefficients, so the total number of rational numbers in the right hand side of (15) is at most $2d+4$. Thus, the bitsize of all the rational appearing in the right hand side of (15) is $\tilde{O}(d^2\tau)$. □

**Theorem 2.9** (Bit complexity of verifying the certificate). *We can verify the positivity certificate in (2), that is the identity in (15) using $\tilde{O}(d^2\tau)$ bit operations.*

*Proof.* The right hand side of (15) involves the squaring of two polynomials, that is, $\tilde{P}_ε$ and $\tilde{Q}_ε$, that have degree $d/2$ and bitsize $\tilde{O}(d\tau)$. Each squaring corresponds to one polynomial multiplication that costs $\tilde{O}_B(d^2\tau)$ bit operations. Then, it suffices to compare the coefficients of the left and right hand sides, which we can do in linear time. □

### 2.3.1 Dropping the assumptions

To certify the nonnegativity of a univariate polynomial using a weighted SOS representation, it suffices to provide such a representation for square-free polynomials with no real roots.

To justify this, assume that $A$ is not square-free and consider its square-free factorization:

$$A = \prod_{\mu} A_{\mu}^{2\mu} \prod_{\nu} A_{\nu}^{2\nu+1}.$$

Each factor raised to an even power, say $A_\mu^{2\mu}$, satisfies $A_\mu(x)^{2\mu} \ge 0$ for all $x \in \mathbb{R}$, and thus does not affect the sign of $A$. Consequently, we may disregard such factors when certifying nonnegativity. For factors raised to odd powers, $A_\nu^{2\nu+1}$, it suffices to analyze the contribution of $A_\nu$ to the sign of $A$. Therefore, we may restrict attention to square-free polynomials.

Moreover, if $A$ has real roots and is nonnegative over $\mathbb{R}$, then all real roots must have even multiplicity. Hence, the real roots are roots of the even-powered factors $A_\mu$.

To construct a weighted SOS representation for a polynomial $A$ that is not square-free, we consider again its square-free factorization:

$$A = \prod_{\mu} A_{\mu}^{2\mu} \prod_{\nu} A_{\nu}^{2\nu+1} = \underbrace{\prod_{\mu} (A_{\mu}^{\mu})^2 \prod_{\nu} (A_{\nu}^{\nu})^2}_{S^2} \prod_{\nu} A_{\nu} = S^2 \prod_{\nu} A_{\nu}.$$

If $A$ is nonnegative over $\mathbb{R}$, then each $A_\nu$ is a positive, square-free polynomial with no real roots. We may apply the uSOS algorithm to each $A_\nu$ to obtain a representation of the form $A_\nu = \sum_{j_\nu} w_{j_\nu} s_{j_\nu}^2$. Multiplying by $S^2$ yields a weighted SOS representation for $A$.

If $A$ is not nonnegative, i.e., if there exists $t \in \mathbb{R}$ such that $A(t) < 0$, then at least one of the $A_\nu(t)$ must be negative, and in fact, an odd number of them must satisfy $A_\nu(t) < 0$. Such a point $t$ necessarily lies between two real roots of the product $\prod_\nu A_\nu$.
### 2.3.2 Witness point of non-nonnegativity

What if A is not nonnegative? In this case, A has at least one real root. By isolating the real roots of A (that is, by computing intervals with rational endpoints, each containing exactly one real root) we can compute rational points between successive real roots. We refer to these as intermediate points.

If A is not nonnegative, then there exists at least one intermediate point $t \in \mathbb{Q}$ such that $A(t) < 0$. The cost of computing these intermediate points is asymptotically the same as that of isolating the real roots of A, which is $\tilde{O}_B(d^3 + d^2\tau)$ [24]. The bitsize of such a t is $\tilde{O}(d\tau)$, e.g., [11, Theorem 1]. This matches the bitsize of the root separation bound of A, since t lies between two distinct real roots.

The evaluation of A at t can be performed in $\tilde{O}_B(d^2\tau)$ bit operations [2, 15].

**Lemma 2.10 (Witness point).** Let A be a univariate polynomial of size $(d, \tau)$. If A is not nonnegative, then there exists a rational number $t \in \mathbb{Q}$ of bitsize $\tilde{O}(d\tau)$ such that $A(t) < 0$. We can compute $t$ in $\tilde{O}_B(d^3 + d^2\tau)$ bit operations and we can verify the inequality $A(t) < 0$ in $\tilde{O}_B(d^2\tau)$.

## 2.4 Positivity over an interval $[a, b]$

We study certificates of positivity over an interval $[a, b]$, where $a, b \in \mathbb{Q}_{\ge 0}$, for the polynomial

$$A = \sum_{k=0}^{d} a_k x^k = a_d \prod_{i=1}^{d} (x - \alpha_i) \in \mathbb{Z}[x].$$

Let A has size $(d, \tau)$ and the bitsize of a and b be bounded by $\sigma$. Following Chevillard et al. [6, Sec. 5.2.5], we consider the transformation,

$$\phi : x \mapsto \frac{a + by^2}{1 + y^2},$$

that in turn induces the following transformation for A:

$$A_{\phi} := (1 + y^2)^d \phi(A) = (1 + y^2)^d A(\frac{a+by^2}{1+y^2}) \in \mathbb{Z}[y],$$

Now, $A_{\phi}$ is nonnegative over $\mathbb{R}$ if and only if A is nonnegative over $[a, b]$. Thus, we can use the results of Sec. 2.1 to certify that $A_{\phi}$ is positive over $\mathbb{R}$, instead of certifying directly that A is positive over an interval.

Notice that if the bitsize of a and b is at most $\sigma$, then the bitsize of $A_{\phi} \in \mathbb{Z}[y]$ is $\tilde{O}(\tau + d\sigma)$. Thus, if we straightforwardly apply the complexity bounds of the previous section, then we end up with a bit complexity bound of $\tilde{O}_B(d^3\sigma + d^2\tau)$. However, we can save a factor of d from the term involving $\sigma$, if we study the effect of $\phi$ on the separation bound, that is the minimum distance between the roots, of A.

A close look in the complexity analysis of uSOS reveals that the two important quantities are the value of the perturbation, $\varepsilon := 2^{-b}$ and the separation bound of the input polynomial A. We study both them for the transformed polynomial $A_{\phi}$. First we consider the separation bound.

### 2.4.1 The separation bound of $A_{\phi}$

If we apply $\phi$ to A and we clear denominators, then the resulting polynomial is

$$A_{\phi}(y) = (1 + y^2)^d \cdot A(\phi(x)) = a_d \cdot \prod_{i=1}^{d} \left( (b - \alpha_i)y^2 + (a - \alpha_i) \right) \in \mathbb{Q}[y],$$

and its roots, for $i \in [d]$, are

$$\zeta_i^{\pm} = \pm \sqrt{\frac{\alpha_i - a}{b - \alpha_i}}.$$

We need to (lower) bound the separation bound for $A_{\phi}$, that is the quantity

$$\Delta_i(A_{\phi}) = |\zeta_i^{\pm} - \zeta_j^{\pm}|,$$
where $\zeta_j^{\pm}$ is the closest root to $\zeta_i^{\pm}$.

We consider the polynomial

$$F(\delta) = (\delta - \zeta_i^+ - \zeta_j^+)(\delta - \zeta_i^+ - \zeta_j^-)(\delta - \zeta_i^- - \zeta_j^+)(\delta - \zeta_i^- - \zeta_j^-) \\ = \delta^4 + 2 \left( \frac{a - \alpha_i}{b - \alpha_i} + \frac{a - \alpha_j}{b - \alpha_j} \right) \delta^2 + \frac{(\alpha_i - \alpha_j)^2 (a - b)^2}{(b - \alpha_i)^2 (b - \alpha_j)^2}.$$

Notice that among the roots of $F \in \mathbb{C}[\delta]$ is the separation bound of $A_{\phi}$, that is the difference $\zeta_i^{\pm} - \zeta_j^{\pm}$. Hence, if we compute a lower bound for the roots of $F$, we also obtain a lower bound for $\Delta_i(A_{\phi})$.

Following [11, Theorem 1], we have that

$$\Delta(A_{\phi}) \ge 2 \frac{\mathrm{tc}(F)}{\|F\|_{\infty}} \ge \frac{\Delta_i(A)^2 (b-a)}{(b-\alpha_i)(b-\alpha_j)[2(ab+\alpha_i\alpha_j) - (\alpha_i+\alpha_j)(a+b)]} \ge 2^{-6\sigma-4\tau-8} \Delta_i(A)^2,$$

where we use the inequalities

$$|\alpha_i| \le 2^{\tau+2}, \quad |b-a| \ge 2^{-2\sigma}, \quad \text{and} \quad |b-\alpha_i| \le 2^{\sigma+\tau+2}.$$

Consequently,

$$\prod_i \Delta_i(A_{\phi}) \ge \prod_i \Delta_i(A) 2^{-6\sigma-4\tau-8} = 2^{-\tilde{O}(d\sigma+d\tau)}.$$

This leads to the following lemma

**Lemma 2.11.** Assume $A \in \mathbb{Z}[x]$ of size $(d, \tau)$ and $a, b \in \mathbb{Q}_{\ge 0}$ of bitsize $\sigma$. Then, $-\lg \prod_i \Delta_i(A_{\phi}) = \tilde{O}(d\sigma + d\tau)$.

### 2.4.2 The perturbation $\varepsilon$ for $A_{\phi}$

The bound on $\varepsilon := 2^{-b}$ depends on a (lower) bound of the minimum of a polynomial $A$, in our case, on a minimum of $A_{\phi}$. If $\xi_i$ are the roots of the derivative of $A_{\phi}$, then $\varepsilon$ depends on a lower bound on $|A_{\phi}(\xi_i)|$. To obtain this lower bound we will not rely on Lem. B.3 as we did in Sec. 2.2. Instead we will rely on the following lemma from this will save us a factor of $d$ in the complexity.

**Lemma 2.12.** [26, Lemma 3] Consider a square-free $A \in \mathbb{R}[x]$ of degree $d$, and let its real roots be $\alpha_i$. Let $x_0 \in \mathbb{R}$ be such that $|x_0 - \alpha_i| \ge \Delta_i/c$ for all real $\alpha_i$ such that $i \ne 1$ and $c \ge 2$. Then

$$|A(x_0)| > |\mathrm{tc}(A)| |x_0 - \alpha_1| c^{1-d} \mathcal{M}(A)^{-1} 2^{\lg \prod_i \Delta_i(A)-1},$$

where $\mathcal{M}(A)$ is the Mahler measure of $A$.

Let $\xi \in \mathbb{R}$ be the root of $A'_{\phi}$ where the minimum of $A_{\phi}$ is attained. Based on Dimitrov [9, Theorem 1] $|\zeta_i - \xi| \ge \Delta_i(A_{\phi})/d$, for all $i$. Also $\mathcal{M}(A_{\phi}) \le \|A_{\phi}\|_2$. In addition $\prod_i \Delta_i(A_{\phi}) \ge 2^{-\tilde{O}(d\sigma+d\tau)}$ (Lem. 2.11).

Thus,

$$|A_{\phi}(\beta)| \ge |\mathrm{tc}(A_{\phi})| |\beta - \alpha_1| d^{1-d} \|A_{\phi}\|_2 2^{\lg \prod_i \Delta_i(A_{\phi})-1},$$

which results in an $\varepsilon = 2^{-b}$, where $b = \tilde{O}(d\sigma + d\tau)$.

**Lemma 2.13** (The value of $\varepsilon$ for $A_{\phi}$). If $\varepsilon = 2^{-b}$, with $b = \tilde{O}(d\sigma + d\tau)$, then $A_{\phi,\varepsilon}(x) = A_{\phi}(x) - \varepsilon M(x) > 0$ for all $x \in \mathbb{R}$, and $-\lg \prod_i \Delta(A_{\phi}) = -\lg \prod_i \Delta(A_{\phi,\varepsilon}) = \tilde{O}(d\sigma + d\tau)$.
### 2.4.3 Overall complexity

**Theorem 2.14 (uSOS at an interval).** Let $A \in \mathbb{Z}[x]$ of size $(d, \tau)$. The algorithm uSOS provides a certificate of positivity of $A$ over an interval $[a, b]$, where $a$ and $b$ are rationals of bitsize $\sigma$, in $\tilde{O}_B(d^3 + d^2\sigma + d^2\tau)$.

*Proof.* The computation of $A_\phi$ consists of a series of Mobius transformations. The most computationally expensive is the (Taylor) shift, that is the transformation $x \mapsto x+a$, where $a$ is a rational of bitsize $\sigma$. This costs $\tilde{O}_B(d^2\sigma + d\tau)$ bit operations and results in a polynomial of bitsize $\tilde{O}(\tau + d\sigma)$, [36].

The input to uSOS is $A_\phi$. By choosing an $\varepsilon := 2^{-b}$, such that $b = \tilde{O}(d\sigma + d\tau)$, we guarantee that $A_{\phi,\varepsilon} = A_\phi - \varepsilon M$ is positive over $\mathbb{R}$ and that the separation bounds of $A_\phi$ and $A_{\phi,\varepsilon}$ are asymptotically the same; in our case $\tilde{O}(d\sigma+d\tau)$. The approximation of the roots of $A_{\phi,\varepsilon}$ up to precision $\kappa$ costs $\tilde{O}_B(d^3+d^2\tau+d^2\tau+d\kappa)$, where $\kappa = \tilde{O}(d\sigma + d\tau)$. This cost dominates the overall complexity of the algorithm. $\square$

**Remark 2.15.** If we target a certificate of positivity for the interval $(0, \infty)$, then it suffices to consider the map $x \mapsto y^2$, and so there is no dependency on $\sigma$.

To recover a weighted SOS representation of $A$ from a weighted SOS decomposition of $A_\phi$, we proceed as follows [6]: In the representation of $A_\phi$, we decompose each polynomial $s_j$ into terms of odd and even degree, that is

$$A_{\phi}(y) = \sum_{j=1}^{\nu} w_j s_j(y)^2 = \sum_j w_j s_{j,e}(y^2)^2 + w_j y^2 s_{j,o}(y^2)^2.$$

Then, we invert the change of variable using $y^2 \mapsto \frac{x-a}{b-x}$ and $1 + y^2 \mapsto \frac{b-a}{b-x}$. Consequently, after clearing denominators, we distinguish two cases. If $d$ is even, then

$$A(x) = \sum_j \frac{w_j}{(b-a)^d} \left[ (b-x)^{\frac{d}{2}} s_{j,e}\left(\frac{x-a}{b-x}\right) \right]^2 + (x-a)(b-x) \sum_j \frac{w_j}{(b-a)^d} \left[ (b-x)^{\frac{d}{2}-1} s_{j,o}\left(\frac{x-a}{b-x}\right) \right]^2.$$

If $d$ is odd, then

$$A(x) = (b-x) \sum_j \frac{w_j}{(b-a)^d} \left[ (b-x)^{\frac{d-1}{2}} s_{j,e}\left(\frac{x-a}{b-x}\right) \right]^2 + (x-a) \sum_j \frac{w_j}{(b-a)^d} \left[ (b-x)^{\frac{d-1}{2}} s_{j,o}\left(\frac{x-a}{b-x}\right) \right]^2.$$

In both cases, it is important to notice that the bitsize of the polynomials in the representation is $\tilde{O}(d\sigma + d\tau)$.

## 3 Perturbed SOS certificate

Even though uSOS algorithm is very efficient, it has the drawback that it forces us to compute (approximate) the roots of the polynomial $A_\varepsilon$, instead of the input polynomial $A$. This has the consequence, that we should work with precision $O(d\tau)$ bits right from the beginning of the algorithm. Even more, the positivity certificate it corresponds to, involves $O(d)$ sums of squares, while we know that it is theoretically possible to construct certificates with a smaller number of summands [27]. We try to leverage these weaknesses, by introducing an alternative, yet closely related, certificates that we call perturbed SOS certificate for rational univariate polynomials. These, certify the nonnegativity for carefully chosen approximations (or perturbations) of the input polynomials. We demonstrate that, if the approximation is closed-enough (under a norm), then the nonnegativity of the perturbed polynomial is equivalent to the nonnegativity of the original polynomial. In this way, the perturbed SOS certificate involves only two squares of polynomials with rational coefficients.

**Theorem 3.1.** Let $A \in \mathbb{Z}[x]$ be of even degree $d$ and maximum coefficient bitsize $\tau$. If there is a nonnegative polynomial $B \in \mathbb{Q}[x]$ such that $\|A - B\|_{\infty} < 2^{-b}$, with $b > 4d\tau + 16d \lg d$, then $A$ is nonnegative.
*Proof.* For all the (real) roots, say $\alpha$, of $A$ it holds $|\alpha| < 2^{\tau+2}$, e.g. [11, Theorem 1]. Moreover, as the leading term of $A$ is positive and $d$ is even, when $x$ goes to $\pm\infty$, then $A(x)$ goes to $+\infty$. Hence, $A(x) > 0$ for every $x$ such that $|x| > 2^{\tau+2}$.

Now let $|x| \le 2^{\tau+2}$, that is $x \in [-2^{\tau+2}, 2^{\tau+2}] =: J$. We will bound the maximum value of function $|(A - B)(x)|$ in this interval. It holds

$$ \max_{x \in J} |(A - B)(x)| \le (2^{(\tau+2)})^d (d+1) \|A - B\|_{\infty} \le 2^{(\tau+2)d+\lg d+1-b} < 2^{-4d\tau-16d\lg d}. $$

As $B(x)$ is a positive function, we have that for any $x \in J$, it holds $0 \le B(x) \le A(x) + 2^{-4d\tau-16d\lg d}$. As $A(x) > 0$ when $x \notin J$, we deduce that

$$ 0 \le A(x) + 2^{-4d\tau-16d\lg d}, \quad \text{for all } x \in \mathbb{R}. \qquad (16) $$

To conclude, it remains to prove that the previous inequality, (16), implies the nonnegativity of $A(x)$. We argue by contradiction. Assume that $A(x)$ is not nonnegative. As $A(x) \ge 0$, $|x| > 2^{\tau+2}$, there must have a negative critical value. Assume that $x^*$ is a critical point realizing the biggest, strictly-negative, critical value. By Lem. B.3, we have that $A(x^*) < -2^{-4d\tau-16d\lg d}$, as $x^*$ is a root of $A'$ but not root of $A$. However, this implies that $A(x^*) + 2^{-4d\tau-16d\lg d} < 0$, which contradicts (16). $\square$

**Definition 3.2.** Let $A \in \mathbb{Q}$ be of size $(d, \tau)$, where $d$ is even, that is nonnegative over $\mathbb{R}$. Fix rational polynomials $s_1, \dots, s_r \in \mathbb{Q}[x]$ and positive rational constants $\lambda_1, \dots, \lambda_r \in \mathbb{Q}_{\ge 0}$, and let $B(x) := \sum_{i=1}^r \lambda_i s_i^2$. If $\|A - B\|_{\infty} < 2^{-b}$, with $b > 4d\tau + 16d\lg d$, then $B$ is a perturbed SOS certificate of nonnegativity of $A$.

**Remark 3.3.** Recall that, if $A(x)$ is square-free, then $A(x) \ge 0$ is equivalent to $A(x) > 0$. Moreover, we can certify that a polynomial is square-free by certifying that the greatest common divisor of $A(x)$ and $A'(x)$ is 1 via Bézout's identity. Therefore, in the square-free case, the perturbed SOS certificates of nonnegativity lead to certificates of positivity.

By combining all the previous ideas, we can exploit the SOS representation of (1) to obtain a perturbed SOS certificate for a square-free polynomial $A$. For this, we will approximate the roots of $A$ up to a precision that will guarantee that the induced polynomials $\tilde{P}$ and $\tilde{Q}$ form a perturbed SOS certificate of positivity for $A$. The following theorem gives the details of this approach.

**Theorem 3.4.** Let $A \in \mathbb{Z}[x]$ be a positive square-free of even degree $d = 2m$ and maximum coefficient bitsize $\tau$. Let $\alpha_j^+ = \gamma_j + i\delta_j$, where $\gamma_j \in \mathbb{R}$, $\delta_j \in \mathbb{R}_{\ge 0}$, for $j \in [m]$; that is $\alpha_j^+$ are the roots with positive imaginary part.

Let $\tilde{\gamma}_j \in \mathbb{Q}$, $\tilde{\delta}_j \in \mathbb{Q}_{>0}$, for $j \in [m]$, be rational approximations up to precision $2^{-\lambda}$, for $\lambda = 9d\tau + 60d\lg d$. That is

$$ |\alpha_i - (\tilde{\gamma}_j - i\tilde{\delta}_j)| \le 2^{-\lambda}. $$

Then, the polynomial $\tilde{P}^2 + \tilde{Q}^2$, where $\tilde{P}, \tilde{Q} \in \mathbb{Q}[x]$ are defined as

$$ \tilde{P}(x) + i \tilde{Q}(x) = \prod_{j=1}^{m} (x - \tilde{\gamma}_j - i \tilde{\delta}_j). $$

is a perturbed SOS certificate of positivity of $A$.

*Proof.* All the roots of $A$ have magnitude smaller than $2^{\tau+2}$, e.g. [11, Theorem 1]. Using this, and the fact that $\lambda > 9d\tau + 60d\lg d$, Lem. B.4, implies the following inequality

$$ \left\| A(x) - 1c(A) \prod_{j=1}^{m} (x - \tilde{\gamma}_j - i \tilde{\delta}_j)(x - \tilde{\gamma}_j + i \tilde{\delta}_j) \right\|_{\infty} < 2^{-\lambda+(4d(\tau+2)+32d-7)+1} < 2^{-5d\tau-19d\lg d}. $$

Moreover, following the arguments in Step 2 of Sec. 2.1

$$ \tilde{P}(x)^2 + \tilde{Q}(x)^2 = 1c(A) \prod_{j=1}^{m} (x - \tilde{\gamma}_j - i \tilde{\delta}_j)(x - \tilde{\gamma}_j + i \tilde{\delta}_j). $$

Therefore, from Def. 3.2 and Thm. 3.1, $\tilde{P}^2 + \tilde{Q}^2$ is a perturbed SOS certificate of nonnegativity for $A(x)$. As $A(x)$ is square-free, we conclude that $A(x)$ is positive. $\square$
Regarding the complexity of computing the perturbed certificate, we have the following theorem.

**Theorem 3.5.** Let $A \in \mathbb{Z}[x]$ be a square-free positive polynomial of over $\mathbb{R}$, of size $(d, \tau)$. We can compute a perturbed SOS certificate of positivity of $A$ in $\tilde{O}_B(d^3 + d^2\tau)$. The certificate consists of the polynomial $\tilde{P}^2 + \tilde{Q}^2$ that has bitsize $\tilde{O}(d\tau)$.

*Proof.* The computation of a perturbed SOS certificate for $A$, requires us to (efficiently) approximate the roots $A(x)$ (up to any desired precision), say $\lambda$. In our case, $\lambda = \tilde{O}(d\tau)$ (Thm. 3.4). The approximation of all the roots requires $\tilde{O}_B(d^3 + d^2\tau)$ bit operations (Thm. B.2). The computation of $\tilde{P}$ and $\tilde{Q}$ is based on Lem. B.4 and also costs $\tilde{O}_B(d^3 + d^2\tau)$. Lem. B.4 shows that the bitsizes of $\tilde{P}$ and $\tilde{Q}$ are $\tilde{O}(d\tau)$. $\square$

# 4 Interlacing polynomials and T-systems

For a given square-free polynomial $A \in \mathbb{R}[x]$, positive over $\mathbb{R}$, $(0, \infty)$, or any interval $[a, b] \subset \mathbb{R}$, we establish a connection, actually an equivalence, between an SOS representation of $A$, and the T-systems and Karlin points [8]. In this way, we make the theory of T-systems constructive we show that the Karlin points are the real roots of certain interlacing polynomials. We refer the reader to Section B.3 for the (very) basic definitions and to [8] for further details on T-systems and positive polynomials.

## 4.1 Strict positivity and Karlin points over $\mathbb{R}$

Consider the square-free polynomial $A(x) = \sum_{k=0}^d a_k x^k \in \mathbb{R}[x]$, of even degree $d$ and $a_d > 0$. We further assume that $A$ is strictly positive over $\mathbb{R}$, hence it does not have any real roots; its factorization to linear factors is

$$A(x) = a_d \prod_{i=1}^{d} (x - \alpha_i) = a_d \prod_{i=1}^{d/2} (x - \gamma_i + i \delta_i) \prod_{i=1}^{d/2} (x - \gamma_i - i \delta_i), \quad (17)$$

where we assume that $\delta_i > 0$ for all $i \in [d/2]$. We define the polynomials $P, Q \in \mathbb{R}[x]$ in such a way that the following equations hold,

$$\prod_{i=1}^{d/2} (x - \gamma_i + i \delta_i) = P(x) + iQ(x) \quad \text{and} \quad \prod_{i=1}^{d/2} (x - \gamma_i - i \delta_i) = P(x) - iQ(x). \quad (18)$$

**Lemma 4.1 (Positivity over $\mathbb{R}$ and interlacing).** If $A = \sum_{k=0}^d a_k x^k \in \mathbb{R}[x]$, of degree $d = 2m$ and $a_d > 0$, is strictly positive over $\mathbb{R}$, let $P, Q \in \mathbb{R}[x]$ be as defined above. Then, we have that $A(x) = a_d P(x)^2 + a_d Q(x)^2$ and the polynomials $P$ and $Q$ are interlacing of degrees $m$ and $m - 1$, respectively.

*Proof.* Equations (17) and (18) imply that $A$ is a sum of two squares. If there are no assumptions on $\delta_i$'s, then there are $2^{d/2-1}$ possible distinct pairs of $\{P(x) \pm iQ(x)\}$, e.g., [31], and consequently, $2^{d/2-1}$ inequivalent representations of $A$ as a sum of two squares. However, by the uniqueness of the decomposition in Theorem B.7, the only representation that results on $P$ and $Q$ being interlacing are the ones that we obtain when we impose the condition $\delta_i > 0$, or equivalently $\delta_i < 0$, for all $i \in [d/2]$. This is no coincidence, as the claim is a direct consequence of the Hermite-Biehler (HB) theorem (Theorem A.1), that states that $P(x) + iQ(x)$ has not roots in the upper half plane if and only if $P$ and $Q$ interlace. $\square$

In the Appendix, we present an alternative proof of Lemma 4.1, which also could be seen as a proof of HB theorem. It is not based on analytic properties of polynomials, as is the proof(s) of HB, e.g. [30, Thm. 6.3.4], but relies directly on the factorization of a univariate polynomial to linear factors. It seems, from our point of view at least, that our proof is conceptually simpler, and potentially of independent interest. For a restatement of the HB theorem in the modern language of real-stable polynomials and proper position we refer the reader to the work of Bränden, e.g. [5].
**Remark 4.2.** An interesting remark is that since P and Q are real rooted, their number of positive real roots is exactly the number of sign variations in their coefficients list. This follows from Descartes' rule of sign.

From the point of view of T-systems, Section B.3 and Theorem B.6, for a strictly positive polynomial A over $\mathbb{R}$ the polynomials $P^2$ and $Q^2$ of Lamma 4.1 play the role of $f^*$ and $f_*$. In our case $\mathcal{F} = \{1, x, x^2, \dots, x^d\}$. Also, the Karlin points, see Equation (3), are unique and interlacing. The latter is also a property of the real roots of P and Q. Hence, we have the following corollary that characterizes the Karlin points and allows us to compute them.

**Corollary 4.3 (Karlin points over $\mathbb{R}$).** Let $A \in \mathbb{R}[x]$ of degree $d = 2m$ be strictly positive over $\mathbb{R}$, and $A(x) = P(x)^2 + Q(x)^2$, where $P, Q \in \mathbb{R}[x]$, as in Lemma 4.1. The real roots of P and Q are the Karlin points.

## 4.2 Strict positivity and Karlin points over $[0, \infty)$

Consider $A(x) = \sum_{k=0}^{d} a_k x^k \in \mathbb{R}[x]$, where $a_d > 0$ and $a_0 \neq 0$. We assume that A is strictly positive over $[0, \infty)$; so, it does not have nonnegative real roots. In this case, Polýa and Szegö proved, the now classical result, that A has the representation $A(x) = f(x) + xg(x)$, where $f$ and $g$ are sum of squares of real polynomials and $\deg(f), \deg(xg) \le \deg(A)$, e.g., [29, Prop. 2].

We will show that if A is strictly positive over $[0, \infty)$, then A has a representation $A = \mathbb{P}^2 + x\mathbb{Q}^2$, where $\mathbb{P}$ and $\mathbb{Q}$ are interlacing polynomials and their real roots are the Karlin points of A in $[0, \infty)$. We use the notation $\mathbb{P}$ and $\mathbb{Q}$ to highlight that these polynomials are different from the polynomials P and Q of Lemma 4.1.

**Theorem 4.4 (Strict positivity over $[0, \infty)$).** Let $A = \sum_{k=0}^{d} a_k x^k \in \mathbb{R}[x]$, of degree $d$ and $a_d > 0$, be strictly positive over $[0, \infty)$. Then, $A(x) = a_d \mathbb{P}(x)^2 + a_d x \mathbb{Q}(x)^2$, where the polynomials $\mathbb{P}, \mathbb{Q} \in \mathbb{R}[x]$ are interlacing and have only positive real roots. If A has even degree, $\deg(A) = d = 2m$, then $\deg(P) = \frac{d}{2}$ and $\deg(Q) = \frac{d}{2} - 1$. If A has odd degree, $\deg(A) = d = 2m + 1$, then $\deg(P) = \frac{d-1}{2} = \deg(Q)$.

To simplify the calculations, we assume that A is monic, that is $a_d = 1$.

*Proof.* Assume that $A(x)$ is positive in the interval $(0, \infty)$. Then, $A(x)$ admits the following factorization

$$A(x) = \prod_{i=1}^{r} (x + \gamma_i) \prod_{j=1}^{t} (x - \alpha_j - i\beta_j)(x - \alpha_j + i\beta_j),$$

where $d = r + 2t$, $\gamma_i, \beta_i > 0$ and $\alpha_i \in \mathbb{R}$.

Consider now the polynomial $A(x^2)$. The factorization of this polynomial is

$$A(x^2) = \prod_{i=1}^{r} (x - i\sqrt{\gamma_i})(x + i\sqrt{\gamma_i}) \prod_{j=1}^{t} (x - \eta_j^+ + i\eta_j^-)(x + \eta_j^+ - i\eta_j^-)(x - \eta_j^+ - i\eta_j^-)(x + \eta_j^+ + i\eta_j^-),$$

where $\eta_j^{\pm} = \sqrt{\frac{\alpha_j^2 + \beta_j^2 \pm \alpha_j}{2}}$, $\gamma_i, \beta_i > 0$ and $\alpha_i \in \mathbb{R}$.

Notice that $A(x^2)$ is strictly positive in $\mathbb{R}$, so by Lemma 4.1, there are two real polynomials P and Q with only reals interleaving roots, such that $A(x^2) = P^2 + Q^2$ and,

$$\begin{aligned} P + iQ &= \prod_{i=1}^{r} (x - i\sqrt{\gamma_i}) \prod_{j=1}^{t} (x + \eta_j^+ - i\eta_j^-)(x - \eta_j^+ - i\eta_j^-) \\ &= (i)^r \prod_{i=1}^{r} (\sqrt{\gamma_i} - i x) \prod_{j=1}^{t} (x^2 - (\eta_j^+)^2 - (\eta_j^-)^2 - i 2\eta_j^- x). \end{aligned}$$
We introduce a new variable $\omega$ and we set $\omega = i x$. Then, there are real constants $\{c_{a,b}\}_{a,b}$ such that,

$$ P + iQ = (i)^r \prod_{i=1}^{r} (\sqrt{\gamma_i} - \omega) \prod_{j=1}^{t} (x^2 - (\eta_j^+)^2 - (\eta_j^-)^2 - 2\eta_j^- \omega) = (i)^r \sum_{a,b} c_{a,b}(x^2)^a \omega^b \\ = (i)^r \left( \sum_{a,k} c_{a,2k}(x^2)^a \omega^{2k} + \sum_{a,k} c_{a,2k+1}(x^2)^a \omega^{2k+1} \right). $$

By combining all the identities, we get that

$$ P + iQ = (i)^r \left( \sum_{a,k} (-1)^k c_{a,2k}(x^2)^{a+k} + i x \sum_{a,k} (-1)^k c_{a,2k+1}(x^2)^{a+k} \right) =: (i)^r (\mathcal{P}(x^2) + i x \mathcal{Q}(x^2)). $$

Therefore, it holds

$$ A(x^2) = P^2 + Q^2 = (i)^r (-i)^r (\mathcal{P}(x^2) + i x \mathcal{Q}(x^2))(\mathcal{P}(x^2) - i x \mathcal{Q}(x^2)) = \mathcal{P}(x^2)^2 + x^2 \mathcal{Q}(x^2)^2, $$

and so $A(x) = \mathcal{P}(x)^2 + x \mathcal{Q}(x)^2$. To see that the roots of $\mathcal{P}$ and $\mathcal{Q}$ interlace, with no loss of generality, assume $r$ is even. Then, $P(x) = (-1)^{r/2}\mathcal{P}(x^2)$ and $Q(x) = (-1)^{r/2}x\mathcal{Q}(x^2)$. By Lemma 4.1, $P$ and $Q$ are real rooted and interlacing. Hence, $\mathcal{P}$ and $\mathcal{Q}$ have only real positive roots. Also they are interlacing, as the square root of every root of $\mathcal{P}$, respectively $\mathcal{Q}$, is a root of $P$, respectively $Q$. $\square$

As in the case of $\mathbb{R}$, the next corollary follows straightforwardly if we combine Lemma 4.4 with the uniqueness property of Karlin points.

**Corollary 4.5** (Karlin points over $[0, \infty)$). Let $A \in \mathbb{R}[x]$, of degree $d$, be positive over $(0, \infty)$, and $A(x) = \mathcal{P}(x)^2 + x\mathcal{Q}(x)^2$, where $\mathcal{P}, \mathcal{Q} \in \mathbb{R}[x]$, as in Lemma 4.4. The real roots of $\mathcal{P}$ and $\mathcal{Q}$ are the Karlin points of $A$ in $(0, \infty)$.

What about strict positivity over an arbitrary interval $I \subseteq \mathbb{R}$? To study this case, we follow closely Powers and Reznick [29]. If a polynomial $A$ is strictly positive over an interval $[a, b] \subseteq \mathbb{R}$, then by the transformation

$$ B(x) = A\left(\frac{(b-a)x+(b+a)}{2}\right), $$

we obtain a polynomial $B$ that is strictly positive over $[-1, 1]$. The transformation involves homothecy and translation; both of them preserve interlacing. Hence, it suffices to consider the strict positivity of $A$ over $[-1, 1]$. Even more, strict positivity over $[0, \infty)$ and strict positivity over $[-1, 1]$ are closely related through the Goursat transform. The $d$-th degree Goursat transform is

$$ \mathcal{G}[A](x) = (1 + x)^d A\left(\frac{1-x}{1+x}\right). $$

If we apply $\mathcal{G}$ two times, then we notice that it is almost, up to a constant depending on the degree, its own inverse, that is

$$ \mathcal{G}[\mathcal{G}[A]](x) = (1 + x)^d \mathcal{G}[A]\left(\frac{1-x}{1+x}\right) = (1 + x)^d \left(1 + \frac{1-x}{1+x}\right)^d A\left(\frac{1-\frac{1-x}{1+x}}{1+\frac{1-x}{1+x}}\right) = 2^d A(x). $$

So, Goursat's lemma states that $A$ is strictly positive over $[-1, 1]$ if and only if $\mathcal{G}[A]$ is strictly positive over $[0, \infty)$ and $\deg(\mathcal{G}[A]) = d$ [29, Lemma 1]. Similarly, $A$ is strictly positive over $[0, \infty)$ if and only if $\mathcal{G}[A]$ is strictly positive over $[-1, 1]$ and $\deg(\mathcal{G}[A]) \le d$.

The application of $\mathcal{G}$ consists of a composition of homothecies, translations, and inversions. All three preserve interlacing, hence it suffices to study the positivity, and compute the Karlin points, on $(0, \infty)$.

Even though from a theoretical point of view $(0, \infty)$ and $[-1, 1]$ are almost equivalent, with respect to positivity certificates, this is not, exactly the case, from a practical point of view. Goursat's transform might increase the norm of the polynomial and hence we might be forced to compute with bigger coefficients. We do not exploit this direction further.
# 5 Implementation and experiments

We provide an open source implementation of uSOS in MAPLE, available in UPOS. We also perform some experiments to demonstrate the efficiency of the algorithm and our implementation. The software is easy to use and, besides MAPLE, does not rely on external software packages.

## 5.1 Some insights on the implementation

Regarding the actual implementation of uSOS, there are certain choices that play an important role in practice.

Computing the perturbation $\epsilon$. The first and most important decision concerns how to compute the perturbation $\epsilon$ in Step 1 of the algorithm. An obvious choice would be to use the a-priori worst case bound for the minimum of the input polynomial $A$, see Lem. B.3. This is a bad choice, as we are forcing the implementation to always operate with the worst case number of bits and slows down the running times dramatically. We do not consider this case at all. Our goal should be to compute with the actual minimum of $A$ and not rely on worst case bounds; which very often overestimate a lot the bitsize.

In this context, a straightforward strategy to compute eps consists in starting with an initial value, say $\frac{1}{2}$, and repeatedly dividing by 2, until the polynomial $A_\epsilon := A - \epsilon M$ has no real roots and positive leading coefficient. Hence, each iteration performs a call to a real root-finding algorithm. This is the strategy implemented in the [22]. Unfortunately, this approach also has significant theoretical and practical drawbacks. The worst case bound for $\epsilon$ is, roughly, $2^{-\tilde{O}(d\tau)}$ (Cor. 2.4, Lem. 2.3), thus, with this approach, we need to perform $\tilde{O}(d\tau)$ iterations, in the worst case, to compute $\epsilon$. This leads to a bit complexity bound $\tilde{O}_B(d^3\tau^2)$, which is not anymore linear in the bitsize of the input! From a practical point of view, this choice has the consequence that it forces the number of iteration to depend on the bitsize of $\epsilon$, which, if it is big, slows down significantly the running times, even when the worst case bitsize bounds are not attained.

A natural question to ask is if we can come up with an efficient implementation to compute $\epsilon$ that, at the same time, does not dominate or alter the overall worst case bit complexity of uSOS.

To overcome the previous, theoretical and practical, obstacles we rely on the important and simply observation that we can consider the dyadic representation of the perturbation, $\epsilon = 2^{-b}$; notice that $b = \tilde{O}(d\tau)$, roughly (Lem. 2.3). Then, instead computing $\epsilon$, our goal should be to compute $b$. We do so by performing exponential binary search starting from zero. To find $2^{-b}$, we start at exponent $e = 0$, i.e., $2^0 = 1$, and successively double it: 1, 2, 4, ... until $e$ exceeds $b$. This exponential phase yields an interval $(e/2, e]$ that must contain $b$. Then, we perform binary search on this interval to compute the exact exponent $b$. The exponential phase doubles the exponent until it exceeds $b$, requiring $O(\log b)$ steps. The subsequent binary search on the interval $(e/2, e]$ also takes $O(\lg b)$ steps. Thus, the overall complexity of exponential binary search on the exponent to find $\epsilon := 2^{-b}$ is $O(\log b) = \tilde{O}(\lg(d\tau))$ iterations. Recall that corresponds to number of times we call the real root-finding algorithm to test whether that polynomial $A - 2^{-e}M$ has real roots or not. This approach has the advantages of all the other approaches and none of their drawbacks. It does not alter the worst case complexity bound and from a practical point of view the gain is also enormous.

In Table 1 we present a comparison on various approaches for computing $\epsilon$, on some positive polynomials coming from [6]; see also the next section for further details. The first three columns are: the id of the polynomial, its degree and its input bitsize. The fourth columns is the computed bound on $\epsilon := 2^b$; actually its exponent. The last three columns present the number of steps that various algorithms perform to compute $\epsilon$. The column SubDiv is the subdivision method, starting from $\frac{1}{2}$ and repeatedly dividing by 2. The column Min+ExpBin corresponds to a hybrid algorithm that first approximates the minimum the polynomial and then finds the best dyadic bound using exponential search. Finally, the last column, ExpBin, is the (pure) exponential binary search algorithm that we described earlier. Evidently exponential binary search is the most efficient algorithm for computing $\epsilon$. For a demonstration of how the computation of $\epsilon$ affects the actual running time of uSOS, we refer to Table 4. In this set of experiments, we consider the three algorithms on modified Wilkinson polynomials, see Eq. 19 and also next section. Again, the exponential binary search approach is best approach. Our implementation of uSOS offers all three approaches for computing $\epsilon$.
<table><caption>Table 1: Comparison of the number of steps that various approaches perform to compute a bound on $ε$; the initial perturbation of uSOS. The polynomials are from [6].</caption><thead><tr><th>#</th><th>Degree</th><th>Input<br>bitsize</th><th>ε = 2<sup>-b</sup></th><th>SubDiv<br>#(iters)</th><th>Min+ExpBin<br>#(iters)</th><th>ExpBin<br>#(iters)</th></tr></thead><tbody><tr><td>1</td><td>13</td><td>359</td><td>60</td><td>61</td><td>15</td><td>13</td></tr><tr><td>3</td><td>32</td><td>439</td><td>116</td><td>117</td><td>27</td><td>15</td></tr><tr><td>4</td><td>22</td><td>492</td><td>103</td><td>104</td><td>24</td><td>15</td></tr><tr><td>5</td><td>34</td><td>775</td><td>235</td><td>236</td><td>27</td><td>17</td></tr><tr><td>6</td><td>17</td><td>190</td><td>69</td><td>70</td><td>30</td><td>15</td></tr><tr><td>7</td><td>43</td><td>371</td><td>70</td><td>71</td><td>27</td><td>15</td></tr><tr><td>8</td><td>22</td><td>275</td><td>74</td><td>75</td><td>24</td><td>15</td></tr><tr><td>9</td><td>20</td><td>353</td><td>45</td><td>46</td><td>27</td><td>13</td></tr><tr><td>10</td><td>25</td><td>312</td><td>67</td><td>68</td><td>24</td><td>15</td></tr></tbody></table>

*Deciding the working precision.* Another important choice concerns the (initial) precision we use to approximate the roots of the polynomial $A_ε$. The implementation of (complex) root finding should avoid working right from the beginning with the worst case theoretical bounds. It does not seem very likely that we will need to compute with (that) many bits of precision. We expect that the separation bound, that is the minimum distance between the roots of a polynomial not to be very small. Indeed this is the case for (a wide variety of) “random” polynomials [12]. Following this discussion, our implementation starts with some initial precision and if the corresponding inequalities are not satisfied, Eq. (7), then we double it. Our experiments suggest that a good practical heuristic is to consider a starting accuracy that depends on the bitsize of the perturbation $ε$. To approximate the complex roots of a polynomial we use the build-in *Isolate* function of *MAPLE*. It realizes the algorithm of Imbach and Moroz [16]. To refine roots, up to any desired precision, we use the build-in *MAPLE* function *hefroots:-refine* which, unfortunately, to the date, is without documentation.

*A disclaimer.* We should mention, that even though it is easier to provide mathematical software on top of well established computer algebra packages, like *MAPLE*, this convenience comes with certain limitations. For example, we cannot work exclusively with numbers that are powers of two; these are the only type of numbers that we need for the perturbation $ε$ and the approximation of the roots of $A_ε$. Such implementation tricks can speed the actual running times by several orders of magnitude, e.g. [35, 32]. In *MAPLE*, it is very difficult, if possible at all, to implement such tricks, while it is not that complicated in standard programming languages, like C/C++.

## 5.2 Experiments with positive polynomials

We perform various experiments to demonstrate the efficiency of our implementation and to study its practical behavior. All the experiments were performed using *MAPLE* 25, on MacBook Air with an Apple M2 cpu, having 24GB of memory, and running *Sonoma* 14.7.5. The running times presented, at the last column of every table, are the average of 10 runs.

The first set of experiments that we performed are on polynomials coming from [6]. This is a set of nine polynomials that we need to certify that they are positive in a (small) interval having rational coefficients. The results appear in Table 2. The first column is the index of the polynomial, the second its degree, the third its bitsize, the fourth is the maximum bitsize of the polynomials and rationals in the SOS decomposition, and the last one is the time needed by our implementation of uSOS to compute the rational SOS decomposition. Even though our implementation is quite efficient and the bitsize of the output is reasonable, it is difficult to draw general conclusion as the polynomials are varying difficulty, i.e., different bitsizes, minimum, separations bounds, etc.

Another set of experiments considers modified Wilkinson polynomials, inspired from [22]. These are
<table><caption>Table 2: Positive polynomials in an interval from [6].</caption><thead><tr><th>#</th><th>Degree</th><th>Input bitsize</th><th>Output bitsize</th><th>Time (ms)</th></tr></thead><tbody><tr><td>1</td><td>13</td><td>359</td><td>2 655</td><td>7</td></tr><tr><td>3</td><td>32</td><td>439</td><td>8 084</td><td>62</td></tr><tr><td>4</td><td>22</td><td>492</td><td>4 351</td><td>28</td></tr><tr><td>5</td><td>34</td><td>775</td><td>16 675</td><td>171</td></tr><tr><td>6</td><td>17</td><td>190</td><td>3 227</td><td>15</td></tr><tr><td>7</td><td>43</td><td>371</td><td>14 540</td><td>176</td></tr><tr><td>8</td><td>22</td><td>275</td><td>4 040</td><td>26</td></tr><tr><td>9</td><td>20</td><td>353</td><td>2 620</td><td>8</td></tr><tr><td>10</td><td>25</td><td>312</td><td>3 226</td><td>6</td></tr></tbody></table>

polynomials of the form

$$A(X) = \prod_{i=2}^{n} (x - i)^2 - x^2/11237 + 1. \qquad (19)$$

The output data of this set appears in Table 3. Again we observe a linear fit concerning the theoretical and the experimental bound on the bitsize of the SOS decomposition. The graph and the linear equation appears in Fig. 1. We notice a fluctuation in the output bitsize, which is probably due to the fact that these polynomials are not good representatives of “generic” positive polynomials. Nevertheless, our implementation computes a weighted SOS representation efficiently.

<table><caption>Table 3: Modified Wilkinson polynomials</caption><thead><tr><th>Degree</th><th>Input bitsize</th><th>Output bitsize</th><th>Time (ms)</th></tr></thead><tbody><tr><td>10</td><td>17</td><td>721</td><td>31</td></tr><tr><td>12</td><td>22</td><td>870</td><td>27</td></tr><tr><td>14</td><td>28</td><td>1 283</td><td>49</td></tr><tr><td>16</td><td>34</td><td>2 472</td><td>53</td></tr><tr><td>18</td><td>41</td><td>2 568</td><td>77</td></tr><tr><td>20</td><td>48</td><td>7 527</td><td>146</td></tr><tr><td>22</td><td>55</td><td>8 133</td><td>196</td></tr><tr><td>24</td><td>62</td><td>8 557</td><td>265</td></tr><tr><td>26</td><td>70</td><td>9 118</td><td>145</td></tr><tr><td>28</td><td>78</td><td>8 675</td><td>162</td></tr><tr><td>30</td><td>86</td><td>8 447</td><td>302</td></tr><tr><td>32</td><td>94</td><td>24 294</td><td>822</td></tr><tr><td>34</td><td>102</td><td>25 567</td><td>1 080</td></tr><tr><td>36</td><td>111</td><td>11 077</td><td>1 052</td></tr><tr><td>38</td><td>119</td><td>28 004</td><td>1 536</td></tr><tr><td>40</td><td>128</td><td>28 733</td><td>1 826</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>Figure 1: The linear fit of positive Wilkinson polys is Output Bitsize ≈ 4.8866(d · Input Bitsize) − 1085.69</figcaption></figure>

It is rather a difficult task to consider random positive univariate polynomials. Such polynomials would be useful to study the practical behaviour of uSOS. Unfortunately, this task requires us to sample uniformly from the convex cone of positive polynomials. This is computationally very expensive, even if we consider a polyhedral approximation of the corresponding cone [13]. A rather good compromise<fnref n="2" /> is to consider random

<footnote n="2">²We thank Alperen Ergür for this useful suggestion.</footnote>
<table><caption>Table 4: Running times of uSOS to compute the SOS representation of modified Wilkinson polynomials, Eq. (19), based on different methods to compute the initial perturbation $\epsilon$.</caption><thead><tr><th>Degree</th><th>Input<br>bitsize</th><th>Output<br>bitsize</th><th>SubDiv<br>Time (ms)</th><th>Min+ExpBin<br>Time (ms)</th><th>ExpBin<br>Time (ms)</th></tr></thead><tbody><tr><td>10</td><td>17</td><td>721</td><td>40</td><td>36</td><td>21</td></tr><tr><td>12</td><td>22</td><td>870</td><td>45</td><td>31</td><td>24</td></tr><tr><td>14</td><td>28</td><td>1283</td><td>101</td><td>54</td><td>47</td></tr><tr><td>16</td><td>34</td><td>2472</td><td>116</td><td>66</td><td>49</td></tr><tr><td>18</td><td>41</td><td>2568</td><td>179</td><td>139</td><td>75</td></tr><tr><td>20</td><td>48</td><td>7527</td><td>362</td><td>214</td><td>145</td></tr><tr><td>22</td><td>55</td><td>8133</td><td>635</td><td>323</td><td>207</td></tr><tr><td>24</td><td>62</td><td>8557</td><td>1075</td><td>450</td><td>255</td></tr><tr><td>26</td><td>70</td><td>9118</td><td>728</td><td>159</td><td>143</td></tr><tr><td>28</td><td>78</td><td>8675</td><td>879</td><td>164</td><td>152</td></tr><tr><td>30</td><td>86</td><td>8447</td><td>1316</td><td>525</td><td>291</td></tr><tr><td>32</td><td>94</td><td>24294</td><td>2425</td><td>1309</td><td>770</td></tr><tr><td>34</td><td>102</td><td>25567</td><td>4328</td><td>1690</td><td>983</td></tr><tr><td>36</td><td>111</td><td>11077</td><td>6295</td><td>1833</td><td>994</td></tr><tr><td>38</td><td>119</td><td>28004</td><td>9777</td><td>2837</td><td>1470</td></tr><tr><td>40</td><td>128</td><td>28733</td><td>12717</td><td>3396</td><td>1840</td></tr></tbody></table>

polynomials and the sum of their squares. Hence, the second set of experiments considers polynomials of the form $A(x) = \sum_{i=1}^{\nu} A_i(x)^2$, where we sample the integer coefficients of polynomials $A_i$ uniformly at random from the interval $[-2^{40}, 2^{40}]$; thus, the bitsize of $A$ is around 80. Tables 5, 6, and 7 present the results of the experiments for $\nu \in \{3, 11, 31\}$ SOS summands. Figures 2 and 3 present the graphs of bitsize of the polynomials in the decomposition vs the product of the degree of $A$ and its bitsize, for the two extreme cases 3 and 31. Recall, by Theorem 2.8, the bitsize of the decomposition is $\tilde{O}(d\tau)$. We notice that there is a precise linear fit, as predicted by the theory. In the figures we also mention the explicit linear equations.

# 6 Conclusion

We have presented improved complexity bounds for computing rational weighted sums of squares (SOS) certificates for univariate polynomials positive over $\mathbb{R}$ or rational intervals, refining previous analyses and reducing the bit complexity by a factor of $d$. Beyond the algorithmic advances, we uncovered new structural properties: the SOS summands form an interlacing pair, revealing a connection to Karlin points and T-systems. These insights contribute to our understanding of positivity certificates and their geometric interpretation. Our open-source MAPLE implementation confirms the practical efficiency of the algorithm. Future directions include considering sparse certificates for sparse univariate polynomials.

# Acknowledgements

We thank Mario Kummer for pointing out that Lemma 4.1 corresponds to Hermite-Biehler theorem. The authors thank John May and Jürgen Gerhard for their help with the undocumented Maple function `hefroots:-refine`. This research benefited from the support of the FMJH Program Gaspard Monge for optimization and operations research and their interactions with data science. MB and ET are partially supported by the PGMO grant SOAP, ANR JCJC PeACE (ANR-25-CE48-3760), and ANR PRC ZADyG (ANR-25-CE48-7058). PJdD is supported by the Deutsche Forschungsgemeinschaft DFG with the grant DI-2780/2-1 and his research fellowship at the Zukunftskolleg of the University of Konstanz, funded as part of the Excellence Strategy of the German Federal and State Government.
<table><caption>Table 5: Sum of 3 squares of random polynomials.</caption><thead><tr><th>Deg</th><th>Input bsz</th><th>Output bsz</th><th>Time (ms)</th></tr></thead><tbody><tr><td>20</td><td>82</td><td>1 161</td><td>2</td></tr><tr><td>40</td><td>81</td><td>2 248</td><td>2</td></tr><tr><td>60</td><td>83</td><td>3 296</td><td>4</td></tr><tr><td>80</td><td>83</td><td>4 379</td><td>5</td></tr><tr><td>100</td><td>83</td><td>5 440</td><td>6</td></tr><tr><td>120</td><td>83</td><td>6 520</td><td>7</td></tr><tr><td>140</td><td>83</td><td>7 577</td><td>15</td></tr><tr><td>160</td><td>84</td><td>8 671</td><td>17</td></tr><tr><td>180</td><td>84</td><td>9 754</td><td>39</td></tr><tr><td>200</td><td>84</td><td>10 777</td><td>49</td></tr><tr><td>220</td><td>84</td><td>11 871</td><td>91</td></tr><tr><td>240</td><td>84</td><td>12 980</td><td>248</td></tr><tr><td>260</td><td>84</td><td>14 001</td><td>1 715</td></tr></tbody></table>

<table><caption>Table 6: Sum of 11 squares of random polynomials.</caption><thead><tr><th>Deg</th><th>Input bsz</th><th>Output bsz</th><th>Time (ms)</th></tr></thead><tbody><tr><td>20</td><td>83</td><td>1 139</td><td>2</td></tr><tr><td>40</td><td>84</td><td>2 222</td><td>1</td></tr><tr><td>60</td><td>84</td><td>3 292</td><td>2</td></tr><tr><td>80</td><td>84</td><td>4 346</td><td>3</td></tr><tr><td>100</td><td>84</td><td>5 426</td><td>5</td></tr><tr><td>120</td><td>84</td><td>6 520</td><td>5</td></tr><tr><td>140</td><td>85</td><td>7 581</td><td>6</td></tr><tr><td>160</td><td>85</td><td>8 631</td><td>12</td></tr><tr><td>180</td><td>85</td><td>9 718</td><td>17</td></tr><tr><td>200</td><td>84</td><td>10 830</td><td>22</td></tr><tr><td>220</td><td>85</td><td>11 848</td><td>45</td></tr><tr><td>240</td><td>85</td><td>12 905</td><td>91</td></tr><tr><td>260</td><td>85</td><td>14 029</td><td>369</td></tr></tbody></table>

<table><caption>Table 7: Sum of 31 squares of random polynomials.</caption><thead><tr><th>Deg</th><th>Input bsz</th><th>Output bsz</th><th>Time (ms)</th></tr></thead><tbody><tr><td>20</td><td>84</td><td>1 143</td><td>1</td></tr><tr><td>40</td><td>84</td><td>2 214</td><td>1</td></tr><tr><td>60</td><td>85</td><td>3 281</td><td>3</td></tr><tr><td>80</td><td>85</td><td>4 349</td><td>3</td></tr><tr><td>100</td><td>85</td><td>5 406</td><td>4</td></tr><tr><td>120</td><td>85</td><td>6 532</td><td>4</td></tr><tr><td>140</td><td>85</td><td>7 595</td><td>6</td></tr><tr><td>160</td><td>85</td><td>8 662</td><td>9</td></tr><tr><td>180</td><td>85</td><td>9 704</td><td>16</td></tr><tr><td>200</td><td>86</td><td>10 789</td><td>26</td></tr><tr><td>220</td><td>86</td><td>11 841</td><td>54</td></tr><tr><td>240</td><td>85</td><td>12 935</td><td>88</td></tr><tr><td>260</td><td>86</td><td>14 019</td><td>483</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 2: The linear fit of 3 squares is Output Bsz ≃ 0.6342(d · Input Bsz) + 163.91</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: The linear fit of 31 squares is Output Bsz ≃ 0.6239(d · Input Bsz) + 120.82</figcaption></figure>

## References

[1] L. Baldi, T. Krick, and B. Mourrain. An effective positivstellensatz over the rational numbers for finite semialgebraic sets. *arXiv preprint arXiv:2410.04845*, 2024.

[2] M. Bodrato and A. Zanoni. Long integers and polynomial evaluation with estrin's scheme. In *2011 13th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing*, pages 39–46. IEEE, 2011.

[3] A. Borodin and R. Moenck. Fast modular transforms. *Journal of Computer and System Sciences*, 8(3):366–386, June 1974.

[4] F. Boudaoud, F. Caruso, and M.-F. Roy. Certificates of positivity in the bernstein basis. *Discrete & Computational Geometry*, 39(4):639–655, 2008.

[5] P. Brändén. Polynomials with the half-plane property and matroid theory. *Advances in Mathematics*, 216(1):302–320, 2007.

[6] S. Chevillard, J. Harrison, M. Joldeş, and C. Lauter. Efficient and accurate computation of upper bounds of approximation errors. *Theoretical Computer Science*, 412(16):1523–1543, 2011.
[7] M. M. Davis and D. Papp. Rational dual certificates for weighted sums-of-squares polynomials with boundable bit size. *Journal of Symbolic Computation*, 121:102254, Mar. 2024.

[8] P. J. di Dio. An introduction to T-Systems-with a special emphasis on sparse moment problems, sparse Positivstellensätze, and sparse Nichtnegativstellensätze. arXiv preprint arXiv:2403.04548, 2024.

[9] D. Dimitrov. A refinement of the gauss-lucas theorem. *Proceedings of the American Mathematical Society*, 126(7):2065-2070, 1998.

[10] M. Dressler, S. Iliman, and T. De Wolff. A positivstellensatz for sums of nonnegative circuit polynomials. *SIAM Journal on Applied Algebra and Geometry*, 1(1):536-555, 2017.

[11] I. Emiris, B. Mourrain, and E. Tsigaridas. Separation bounds for polynomial systems. *Journal of Symbolic Computation*, 101:128-151, 2020.

[12] A. Ergür, J. Tonelli-Cueto, and E. Tsigaridas. Beyond worst-case analysis for root isolation algorithms. In *Proc Intl. Symposium on Symbolic and Algebraic Computation (ISSAC)*, pages 139-148, 2022.

[13] A. A. Ergür. Approximating nonnegative polynomials via spectral sparsification. *SIAM Journal on Optimization*, 29(1):852-873, 2019.

[14] S. Fisk. Polynomials, roots, and interlacing. *arXiv preprint math/0612833*, 2006.

[15] W. Hart and A. Novocin. Practical divide-and-conquer algorithms for polynomial arithmetic. In *International Workshop on Computer Algebra in Scientific Computing*, pages 200-214. Springer, 2011.

[16] R. Imbach and G. Moroz. Fast evaluation and root finding for polynomials with floating-point coefficients. In *Proc Intl. Symposium on Symbolic and Algebraic Computation (ISSAC)*, pages 325-334, 2023.

[17] S. Karlin. Representation theorems for positive functions. *Journal of Mathematics and Mechanics*, 12(4):599-617, 1963.

[18] S. Karlin and W. J. Studden. *Tchebycheff systems: With applications in analysis and statistics*. Inter-science Publishers, New York, London, Sydney, 1966.

[19] P. Koprowski, V. Magron, and T. Vaccon. Pourchet's theorem in action: decomposing univariate nonnegative polynomials as sums of five squares. In *Proc 2023 International Symposium on Symbolic and Algebraic Computation (ISSAC)*, pages 425-433, 2023.

[20] T. Krick, B. Mourrain, and A. Szanto. Univariate rational sums of squares. *Revista de la Unión Matemática Argentina*, 64(2):215-237, 2022.

[21] E. Landau. Über die darstellung definiter funktionen durch quadrate. *Mathematische Annalen*, 62(2):272-285, 1906.

[22] V. Magron, M. Safey El Din, and M. Schweighofer. Algorithms for weighted sum of squares decomposition of non-negative univariate polynomials. *Journal of Symbolic Computation*, 93:200-220, 2019.

[23] K. Mehlhorn, M. Sagraloff, and P. Wang. From approximate factorization to root isolation with application to cylindrical algebraic decomposition. *Journal of Symbolic Computation*, 66:34-69, 2015.

[24] V. Y. Pan. Univariate polynomials: Nearly optimal algorithms for numerical factorization and root-finding. *Journal of Symbolic Computation*, 33(5):701-733, May 2002.

[25] V. Y. Pan and E. Tsigaridas. Nearly optimal computations with structured matrices. *Theoretical Computer Science*, 681:117-137, 2017.
[26] V. Y. Pan and E. P. Tsigaridas. Nearly optimal refinement of real roots of a univariate polynomial. *Journal of Symbolic Computation*, 74:181–204, May 2016.

[27] Y. Pourchet. Sur la représentation en somme de carrés des polynômes à une indéterminée sur un corps de nombres algébriques. *Acta Arithmetica*, 19(1):89–104, 1971.

[28] V. Powers. *Certificates of Positivity for Real Polynomials*. Springer, 2021.

[29] V. Powers and B. Reznick. Polynomials that are positive on an interval. *Transactions of the American Mathematical Society*, 352(10):4677–4692, 2000.

[30] Q. I. Rahman and G. Schmeisser. *Analytic theory of polynomials*. Number 26. Oxford University Press, 2002.

[31] B. Reznick. Some concrete aspects of hilbert’s 17th problem. *Contemporary mathematics*, 253(251-272), 2000.

[32] F. Rouillier and P. Zimmermann. Efficient isolation of polynomial’s real roots. *Journal of Computational and Applied Mathematics*, 162(1):33–50, 2004.

[33] M. Schweighofer. Algorithmische Beweise für Nichtnegativ- und Positivstellensätze. Diplomarbeit, Universität Passau, 1999.

[34] M. Sombra. The height of the mixed sparse resultant. *American Journal of Mathematics*, 126(6):1253–1260, 2004.

[35] E. Tsigaridas and I. Emiris. On the complexity of real root isolation using continued fractions. *Theoretical Computer Science*, 392(1-3):158–173, 2008.

[36] J. Von Zur Gathen and J. Gerhard. Fast algorithms for taylor shifts and certain difference equations. In *Proc. International Symposium on Symbolic and Algebraic Computation (ISSAC)*, pages 40–47, 1997.
# A Alternative proof for interlacing

We present an alternative proof of Lemma 4.1, which is also an alternative proof of Hermite-Bielher theorem (Thm. A.1) that is based on induction and direct manipulation of the roots of the polynomials. The HB theorem is as follows:

**Theorem A.1** (Hermite-Biehler [30, Thm. 6.3.4]). Let $F(x) = P(x) + iQ(x) \in \mathbb{C}[x]$, where $P, Q \in \mathbb{R}[x]$. The following statements are equivalent:

(i) All zeros of $F$ lie in the open upper half-plane $\{\zeta \in \mathbb{C} : \Im(\zeta) > 0\}$.

(ii) The real polynomials $P$ and $Q$: (a) are real rooted, (b) have roots that interlace, and (c) have a Wronskian $W(x) := P'(x)Q(x) - P(x)Q'(x)$ of constant nonzero sign on $\mathbb{R}$.

**Lemma** (Positivity over $\mathbb{R}$ and interlacing). If $A = \sum_{k=0}^d a_k x^k \in \mathbb{R}[x]$, of degree $d = 2m$ and $a_d > 0$, is strictly positive over $\mathbb{R}$, let $P, Q \in \mathbb{R}[x]$ be as defined above. Then, we have that $A(x) = a_d P(x)^2 + a_d Q(x)^2$ and the polynomials $P$ and $Q$ are interlacing of degrees $m$ and $m-1$, respectively.

*Proof.* First observe that the identity $A(x) = a_d P(x)^2 + a_d Q(x)^2$ follows from the fact that, by definition of $P(x)$ and $Q(x)$, $A(x) = a_d(P(x) + iQ(x))(P(x) - iQ(x))$.

Hence, we prove by induction on $m$ that $P(x)$ and $Q(x)$ have interlacing roots. Notice that $m$ also corresponds to the number of products required for $P$ and $Q$. To make this explicit, we write $P_m$ and $Q_m$. It holds $\deg(P_m) = m$ and $\deg(Q_m) = m-1$.

Induction start ($m=2$, i.e., $d=4$): Following (18),

$$ (x - \gamma_1 - i \delta_1)(x - \gamma_2 - i \delta_2) = P_2(x) - i Q_2(x), $$

where

$$ P_2(x) = x^2 - (\gamma_1 + \gamma_2)x + \gamma_1\gamma_2 - \delta_1\delta_2 \quad \text{and} \quad Q_2(x) = -(\delta_1 + \delta_2)x + \gamma_1\delta_2 + \gamma_2\delta_1. $$

Regarding $Q_2$, it is of degree 1 and has one real root

$$ \xi := \frac{\gamma_1 \delta_2 + \gamma_2 \delta_1}{\delta_1 + \delta_2} \in \mathbb{R}. $$

As for $P_2$, it has degree 2 and its discriminant is

$$ \mathrm{disc}(P_2) = (\gamma_1 - \gamma_2)^2 + 4\delta_1\delta_2 \geq 0. $$

Hence, $P_2$ has 2 real roots, $\zeta_{\pm} \in \mathbb{R}$,

$$ \zeta_{-} := \frac{\gamma_{1} + \gamma_{2}}{2} - \frac{\sqrt{(\gamma_{1} - \gamma_{2})^{2} + 4\delta_{1}\delta_{2}}}{2} < \frac{\gamma_{1} + \gamma_{2}}{2} + \frac{\sqrt{(\gamma_{1} - \gamma_{2})^{2} + 4\delta_{1}\delta_{2}}}{2} =: \zeta_{+}. $$

It is straightforward that $\zeta_{-} \leq \zeta_{+}$. Overall, both $P_2$ and $Q_2$ are real rooted.

It remains to show that the roots of $P_2$ and $Q_2$ interlace, that is, $\zeta_{-} < \xi < \zeta_{+}$. Regarding the right
inequality, we have

$$
\begin{align*}
\xi < \zeta_+ &\Leftrightarrow \frac{\gamma_1 \delta_2 + \gamma_2 \delta_1}{\delta_1 + \delta_2} < \frac{\gamma_1 + \gamma_2}{2} + \frac{\sqrt{(\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2}}{2} \\
&\Leftrightarrow 0 < \frac{(\delta_1 + \delta_2)(\gamma_1 + \gamma_2) + (\delta_1 + \delta_2)\sqrt{(\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2} - 2(\gamma_1 \delta_2 + \gamma_2 \delta_1)}{2(\delta_1 + \delta_2)} \\
&\Leftrightarrow 0 < (\delta_1 + \delta_2)\sqrt{(\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2} + (\gamma_1 - \gamma_2)(\delta_1 - \delta_2) \\
&\Leftrightarrow (\gamma_2 - \gamma_1)(\delta_1 - \delta_2) < (\delta_1 + \delta_2)\sqrt{(\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2} \\
&\qquad \text{(if the } \text{lhs is negative, then the inequality holds, so we assume it is positive)} \\
&\Leftrightarrow (\gamma_2 - \gamma_1)^2(\delta_1 - \delta_2)^2 < (\delta_1 + \delta_2)^2((\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2) \\
&\qquad \text{(the } \text{lhs is positive, so we square both sides)} \\
&\Leftrightarrow 0 < -(\gamma_2 - \gamma_1)^2(\delta_1 - \delta_2)^2 + (\delta_1 + \delta_2)^2((\gamma_1 - \gamma_2)^2 + 4\delta_1 \delta_2) \\
&\Leftrightarrow 0 < 4\delta_1 \delta_2((\delta_1 + \delta_2)^2 + (\gamma_1 - \gamma_2)^2)
\end{align*}
$$

Similarly we prove that $\zeta_- < \xi$.

Induction step ($m \to m+1$): Assume that the result is true for some $m$. That is, $P_m$ and $Q_m$ are interlacing and their degrees are $m$ and $m-1$, respectively. We obtain the polynomials $P_{m+1}$ and $Q_{m+1}$ as follows

$$
\begin{align*}
P_{m+1} - iQ_{m+1} &= \prod_{i=1}^{m+1} (x - \gamma_i - i \delta_i) = (x - \gamma_{m+1} - i \delta_{m+1}) \prod_{i=1}^{m} (x - \gamma_i - i \delta_i) \\
&= (x - \gamma_{m+1} - i \delta_{m+1})(P_m - iQ_m) \\
&= ((x - \gamma_{m+1})P_m - \delta_{m+1}Q_m) - i(\delta_{m+1}P_m + (x - \gamma_{m+1})Q_m).
\end{align*}
$$

If we write the computation of $P_{m+1}$ and $Q_{m+1}$ from $P_m$ and $Q_m$ in matrix form, then we have

$$
\begin{pmatrix} x - \gamma_{m+1} & -\delta_{m+1} \\ \delta_{m+1} & x - \gamma_{m+1} \end{pmatrix} \begin{pmatrix} P_m \\ Q_m \end{pmatrix} = \begin{pmatrix} P_{m+1} \\ Q_{m+1} \end{pmatrix}.
$$

Following Fisk [14, Cor. 3.54(4)], as $\delta_{m+1} > 0$, the matrix multiplying the vector $(P_m)$ preserves interlacing. Hence, if $P_m$ and $Q_m$ are interlacing, then so are the polynomials $P_{m+1}$ and $Q_{m+1}$. By inspecting the computations, we deduce that the degrees of resulting polynomials are $m+1$ and $m$, respectively. This concludes the proof. $\square$

# B Useful algorithms and complexity bounds

We present some known results that are useful in our analysis. In particular, we use results on the separation bounds of polynomials that are “close” with respect to the one norm, algorithms for multiplication of polynomials, for approximating their (complex) roots to any desired accuracy, and bounds on the minimum of a univariate polynomial.

## B.1 Preliminaries on root separation and approximation

We exploit the work of Mehlhorn, Sagraloff, and Wang [23] on root approximation and refinement for univariate polynomials. The lemmata that we present are simplified variants of the original ones, as we assume we are working with square-free polynomials. We refer the reader to [23] for the general versions and further details. The root isolation algorithm assumes that there is an oracle that is able to provide rational approximations of the coefficients of the input polynomial up to arbitrary precision.
Consider the square free polynomial

$$p(x) = \sum_{i=0}^{n} p_i x^i = p_n \prod_{i=1}^{n} (x - z_i) \in \mathbb{C}[x],$$

where $z_i \in \mathbb{C}$ are its roots. Also let

$$\tilde{p}(x) = p_n \prod_{i=1}^{n} (x - \tilde{z}_i).$$

We need the following notations:

* $M(x) := \max\{1, |x|\}$, for $x \in \mathbb{R}$,

* $\tau_p$ is the minimal nonnegative integer such that $\frac{|p_i|}{|p_n|} \le 2^{\tau_p}$ for all $i = 0, \dots, n-1$,

* $\Gamma_p := M(\max_i(\log M(z_i)))$ denotes the logarithmic root bound of $p$, and

* $\Delta_i = \min_{j \neq i} |z_i - z_j|$ is the local separation bound.

The following lemma relates the separation bound(s) of the roots of $p$ with the separation bound(s) of the roots of $\tilde{p}$, when the polynomials are "sufficiently" close.

**Lemma B.1** ([23, Lemma 3]). Consider $p = \sum_{i=0}^{n} p_i x^i \in \mathbb{C}[x]$. Also, let $\tilde{p} \in \mathbb{C}[x]$ be such that
$$\|p - \tilde{p}\|_1 \le 2^{-b} \|p\|_1.$$

If, for all $i \in [n]$,

$$b \ge \max\{8n, n \log(n)\}, \text{ and } b \text{ is a power of two,} \quad (20)$$

$$2^{-b/2} \le \frac{\Delta_i}{2n}, \quad \text{and} \quad (21)$$

$$2^{-b/2} \le \frac{\prod_{j \neq i} |z_i - z_j|}{16(n+1)2^{\tau_p} M(z_i)^n}, \quad (22)$$

then the disk $D(z_i, 2^{-b/2})$ contains exactly one root approximation. For $i \neq j$, let $\tilde{z}_i$ and $\tilde{z}_j$ be arbitrary approximations of $z_i$ and $z_j$ in the disks $D(z_i, 2^{-b/2})$ and $D(z_j, 2^{-b/2})$, respectively. Then,

$$(1 - \frac{1}{n}) \cdot |z_i - z_j| \le |\tilde{z}_i - \tilde{z}_j| \le (1 + \frac{1}{n}) \cdot |z_i - z_j|. \quad (23)$$

Based on the previous lemma and Pan's root approximation algorithm [24], Mehlhorn, Sagraloff, and Wang [23] developed an algorithm for isolating and approximating the roots of a univariate polynomial up to any desired precision. The algorithm has the additional capability to consider polynomials with bitstream coefficients. A simplified version of their theorem [23, Theorem 4] that considers square-free polynomials is as follows:

**Theorem B.2.** Consider a square-free polynomial $p(x) = \sum_{i=0}^{n} p_i x^i \in \mathbb{C}[x]$ such that $\frac{1}{4} \le p_n \le 1$. If $z_i$ are the roots of $p$, then $P_i := \prod_{j \neq i} |z_i - z_j|$, for $i \in [n]$. We can compute isolating discs $D(\tilde{z}_i, R_i)$ with radius $R_i < 2^{-\kappa}$, for the roots $z_i$ of $p$, in a number of bit operations upper bounded by

$$\tilde{O}_B(n^3 + n^2 \tau_p + n \sum_{i=1}^{n} \lg M(P_i^{-1}) + n \sum_{i=1}^{n} \lg M(\Delta_i^{-1}) + n\kappa).$$

For this bound, we need rational approximation of the coefficients of $p$ up to precision of $L$ bits, where

$$L = \tilde{O}\left(n\Gamma_p + \sum_{i=1}^{n} \lg M(P_i^{-1}) + \sum_{i=1}^{n} \lg M(\Delta_i^{-1}) + n\kappa\right).$$

The numbers $\tilde{z}_i$ are the approximations to the roots $z_i$ and $R_i \le \Delta_i/(64n)$.

The algorithm supported by Theorem B.2 and its complexity depend on the geometry of the roots, that is the (aggregate) separation bound. Moreover, they do not depend on the type and the size of the coefficients.
## B.2 Bounds on the mimimum and the Fan-in algorithm

Besides root approximations, we also need to bound the minimum of a univariate polynomial, the evaluation of this polynomial at the roots of its derivatives, and an algorithm to compute an approximation of this polynomial from approximations of its roots.

**Lemma B.3.** Consider the polynomial

$$A(x) = \sum_{k=0}^{n} a_k x^k = a_n \prod_{i=1}^{n} (x - \alpha_i) \in \mathbb{Q}[x],$$

such that $\frac{1}{2} \le a_n \le 1$ and the other coefficients are rationals with the same denominator of bitsize at most $\tau$. Let $A'$ be the derivative of $A$ with respect to $x$ and let $\beta_j, j \in [n-1]$, be its roots. Then, either $A'(\alpha_i) = 0$ (or $A(\beta_j) = 0$), or it holds

$$2^{-4n\tau-16n \lg n} \le |A'(\alpha_i)|, |A(\beta_j)| \le 2^{2n\tau+8n \lg n}$$

for all $i \in [n]$ and $j \in [n-1]$.

*Proof.* The proof is a direct application of the bounds of the resultant, appeared in [11]. For the first bound, we consider the resultant

$$H = \text{Res}(A(x), y - A'(x), x) = a_n^{n-1} \prod_{i=1}^{d} (y - A'(\alpha_i)) \in \mathbb{Q}[y],$$

that eliminates $x$; the last equality is due to the Poisson formula of the resultant. Then, $H$ is a univariate polynomial in $y$, and its roots are the evaluation of $A'$ at the roots of $A$. To bound the coefficients of $H$, proceed as follows. We notice that $H$ is a homogeneous polynomial of degree $n-1$ in the coefficients of $A$ and homogeneous of degree $n$ in the coefficients of $y - A'(x) \in (\mathbb{Q}[y])[x]$. Specifically, $H$ is of the form

$$H = \dots + \varrho \mathbf{a}_1^{n-1} \mathbf{a}_2^n + \dots,$$

where $\varrho \in \mathbb{Z}$, $\mathbf{a}_1^{n-1}$ denotes a monomial in the coefficients of $A$ of total degree $n$, and $\mathbf{a}_2^{n-1}$ denotes a monomial in the coefficients of $y - A'$ of total degree $n$.

We can bound $\varrho$ using [34], see also [11, Table 1 and Eq. (1)] as

$$|\varrho| \le (n+1)^{n-1} n^n \le n^{2n}.$$

Since the bitsize of $A$ is at most $\tau$, we can upper bound $\mathbf{a}_1^{n-1} \in \mathbb{Q}$ as

$$|\mathbf{a}_1^{n-1}| \le (2^\tau)^{n-1} = 2^{\tau(n-1)};$$

also the denominator of $\mathbf{a}_1^{n-1}$ is an integer at most $2^{\tau(n-1)}$.

To upper bound $\mathbf{a}_2^n$, we consider the worst case scenario that every coefficient of $y - A'$ is $y - n2^\tau$. So $\mathbf{a}_2^{n-1}$ is a polynomial in $y$ with coefficients rational number with numerator at most $\binom{n}{2} 2^{n\tau+n \lg n}$ and denominator at most $2^{n\tau}$.

Taking all these bounds into account, $H$ is a polynomial in $\mathbb{Q}[y]$ of degree at most $n$, its leading coefficient is in $[(\frac{1}{2})^{n-1}, 1]$, and the other coefficients are rationals having a numerator with magnitude at most $2^{2n\tau+7n \lg n}$, and a denominator of magnitude at most $2^{2\tau(n-1)}$.

Consequently, using [11, Theorem 1], we can bound the roots of $H$, and thus the evaluations $A'(\alpha_i)$, as follows

$$2^{-4n\tau-16n \lg n} \le \frac{\text{tc}(H)}{2\|H\|_{\infty}} \le |A'(\alpha_i)| \le 2 \frac{\|H\|_{\infty}}{\text{lc}(H)} \le 2^{2n\tau+8n \lg n}.$$

The same bounds hold for $A(\beta_j)$, where we use the resultant $\text{Res}(A'(x), y - A(x), x) \in \mathbb{Q}[y]$. $\square$
The following theorem supports an algorithm to compute an approximation of a polynomial A when we have approximations of its roots.

**Lemma B.4 ([25, Theorem 17]).** Assume that we are given n complex numbers $z_i$ known up to an absolute precision $\lambda$, that is, we know a rational $\tilde{z}_i$ such that

$$|z_i - \tilde{z}_i| \le 2^{-\lambda}.$$

Also assume that $|z_i| \le 2^{\tau}$, for a positive integer $\tau$.

Based on the Fan-in process of the Moenck-Borodin algorithm [3] we can approximate the (coefficients of the) polynomial with the (rational coefficients of the) polynomial

$$\tilde{m}(x) = \prod_{i} (x - \tilde{z}_i),$$

so that it holds

$$\|m - \tilde{m}\|_{\infty} \le 2^{-\ell+(4n-4)\tau+32n-(\lg n+5)^2-7},$$

at the cost of $\tilde{O}_B(n(\ell + n\tau))$ bit operations. Moreover, it holds

$$\lg \|m\|_{\infty} \le n\tau + 8n - 2 \lg n - 8.$$

## B.3 About T-systems

We borrow the following definition from [8], where we also refer the reader for further details. For the following definition and further properties of the T-systems see [8, Chapter 4].

**Definition B.5 (T-systems).** Let $n \in \mathbb{N}_0$, $\mathcal{X}$ be a set with $|\mathcal{X}| \ge n + 1$, and $\mathcal{F} = \{f_i\}_{i=0}^n$ be a family of real functions $f_i: \mathcal{X} \to \mathbb{R}$. In this setting, a polynomial the following linear combination

$$f = \sum_{i=0}^{n} a_i \cdot f_i \quad \in \mathrm{lin}(\mathcal{F}) := \{a_0 f_0 + \cdots + a_n f_n \mid a_0, \ldots, a_n \in \mathbb{R}\}. \qquad (24)$$

The family $\mathcal{F}$ on $\mathcal{X}$ is a Tchebycheff system (or short T-system) of order $n$ on $\mathcal{X}$ if every polynomial $f \in \mathrm{lin}(\mathcal{F})$ with $\sum_{i=0}^n a_i^2 > 0$ has at most $n$ zeros in $\mathcal{X}$.

If additionally $\mathcal{X}$ is a topological space and $\mathcal{F}$ is a family of continuous functions, then we call $\mathcal{F}$ a continuous T-system.

The following theorem, due to Karlin, see [8, Chapter 7, Theorem 7.1] and references therein, leads to Positivstellensatz for positive univariate polynomials.

**Theorem B.6.** Let $n \in \mathbb{N}_0$, $\mathcal{F} = \{f_i\}_{i=0}^n$ be a continuous T-system of order $n$ on $[a, b]$ with $a < b$, and let $f \in C([a, b], \mathbb{R})$ with $f > 0$ on $[a, b]$ be a strictly positive continuous function. The following hold:

(i) There exists a unique polynomial $f_* \in \mathrm{lin}(\mathcal{F})$ such that

(a) $f(x) \ge f_*(x) \ge 0$ for all $x \in [a, b]$,

(b) $f_*$ vanishes on a set with index $n$,

(c) the function $f - f_*$ vanishes at least once between each pair of adjacent zeros of $f_*$,

(d) the function $f - f_*$ vanishes at least once between the largest zero of $f_*$ and the end point $b$,

(e) $f_*(b) > 0$.

(ii) There exists a unique polynomial $f^* \in \mathrm{lin}(\mathcal{F})$ which satisfies the conditions (a) to (d) of (i) and

$$(e') \quad f^*(b) = 0.$$
As consequences of Theorem B.6 we get the following (much sharper versions of known) results, where we emphasize the uniqueness of $\alpha, \beta > 0$, $x_i$, and $y_i$.

**Corollary B.7.** Let $p \in \mathbb{R}[x]$ with even degree $\deg p = 2m$, $m \in \mathbb{N}$, and $p \ge 0$ on $\mathbb{R}$. Let $z_1, \dots, z_k$ be the zeros of $p$ in $\mathbb{R}$ with (even) multiplicities $m_1, \dots, m_k \in 2\mathbb{N}$. Then there exist unique $\alpha, \beta > 0$ and unique

$$x_1 < y_1 < x_2 < \dots < y_{l-1} < x_l$$

with $\deg p = m_1 + m_2 + \dots + m_k + 2l$ such that

$$p(x) = \prod_{i=1}^{k} (x - z_i)^{m_i} \cdot \left( \alpha \cdot \prod_{i=1}^{l} (x - x_i)^2 + \beta \cdot \prod_{i=1}^{l-1} (x - y_i)^2 \right).$$

In this setting, we consider polynomials in $\mathbb{R}[x]$, thus we can assume that we can assume that we have factored out their real zeros and we consider only polynomials without real roots. Hence, we state the following two consequences of Theorem B.6 only for $p > 0$.

**Corollary B.8.** Let $p \in \mathbb{R}[x]$ with $p > 0$ on $[0, \infty)$.

(i) If $\deg p = 2m$ with $m \in \mathbb{N}_0$, then there exist unique $\alpha, \beta > 0$ and unique

$$0 < x_1 < y_1 < \dots < y_{m-1} < x_m < \infty$$

such that

$$p(x) = \alpha \cdot \prod_{i=1}^{m} (x - x_i)^2 + \beta \cdot x \cdot \prod_{i=1}^{m-1} (x - y_i)^2.$$

(ii) If $\deg p = 2m + 1$ with $m \in \mathbb{N}_0$, then there exist unique $\alpha, \beta > 0$ and unique

$$0 < x_1 < y_1 < \dots < x_m < y_m < \infty$$

such that

$$p(x) = \alpha \cdot \prod_{i=1}^{m} (x - x_i)^2 + \beta \cdot x \cdot \prod_{i=1}^{m} (x - y_i)^2.$$

**Corollary B.9.** Let $p \in \mathbb{R}[x]$ with $p > 0$ on $[a, b]$ for some $a, b \in \mathbb{R}$ with $a < b$.

(i) If $\deg p = 2m$ for some $m \in \mathbb{N}_0$, then there exist unique $\alpha, \beta > 0$ and unique

$$a < x_1 < y_1 < \dots < y_{m-1} < x_m < b$$

such that

$$p(x) = \alpha \cdot \prod_{i=1}^{m} (x - x_i)^2 + \beta \cdot (x - a) \cdot (b - x) \cdot \prod_{i=1}^{m-1} (x - y_i)^2.$$

(ii) If $\deg p = 2m + 1$ for some $m \in \mathbb{N}_0$, then there exist unique $\alpha, \beta > 0$ and

$$a < x_1 < y_1 < \dots < x_m < y_m < b$$

such that

$$p(x) = \alpha \cdot (x - a) \cdot \prod_{i=1}^{m} (x - x_i)^2 + \beta \cdot (b - x) \cdot \prod_{i=1}^{m} (x - y_i)^2.$$

Even though we have presented the statements that if $p > 0$ (on $\mathbb{R}$, $[0, \infty)$, or $[a, b]$), then there exists a unique representation as a sum of squares (with additional factors $x$, $x - a$, and $b - x$), the opposite implications are straightforward. Thus, all the previous imply $p > 0$. Therefore, these equivalences are positivity certificates.