Fast Evaluation of Real and Complex Polynomials

We propose an algorithm for quickly evaluating polynomials. It pre-conditions a complex polynomial $P$ of degree $d$ in time $O(d\log d)$, with a low multiplicative constant independent of the precision. Subsequent evaluations of $P$ computed with a fixed precision of $p$ bits are performed in average arithmetic complexity $O\big(\sqrt{d(p+\log d)}\big)$ and memory $O(dp)$. The average complexity is computed with respect to points $z \in \mathbb{C}$, weighted by the spherical area of $\overline{\mathbb{C}}$. The worst case does not exceed the complexity of H{\"o}rner's scheme. In particular, our algorithm performs asymptotically as $O(\sqrt{d\log d})$ per evaluation. For many classes of polynomials, in particular those with random coefficients in a bounded region of $\mathbb{C}$, or for sparse polynomials, our algorithm performs much better than this upper bound, without any modification or parameterization.The article contains a detailed analysis of the complexity and a full error analysis, which guarantees that the algorithm performs as well as H\''orner's scheme, only faster. Our algorithm is implemented in a companion library, written in standard C and released as an open-source project [MV22].Our claims regarding complexity and accuracy are confirmed in practice by a set of comprehensive benchmarks.


Introduction
The study of polynomials has sparked the interest of many generations of mathematicians and inspired major theoretical developments. In modern algebra, the notion of group stemmed from the impossibility of solving polynomials with radicals; abstract rings generalize the properties of Z and Z [X]. Modern number theory is indissociable from polynomials and algebraic curves.
Modern analysis too evolved from the prototype of a function space given by polynomials. A few obvious testimonies to this heritage are Descartes's notation of x, y,. . . for the variables of functions, the fact that successive approximations of a real number in base b are polynomials in b −1 , or the fact that polynomials in e ix 1 , . . . , e ixn (i.e. trigonometric polynomials) are the archetype of periodic functions over R n . Smooth functions can be approximated locally by Taylor's polynomial expansions or globally thanks to the Weierstrass approximation theorem.
Polynomials are also ubiquitous due to their practical interest. Greeks and Babylonians used quadratic equations circa 2000 BC to compute the boundaries of their agricultural fields in order to define fair taxes and trade rules. About 4000 years later, we still handle polynomials and solve polynomial equations, not just on school benches, but also in real life to find the natural modes of oscillations of engineering structures or the rate of spread of a virus. Polynomials are at the heart of numerical analysis and appear in particular as approximations and interpolations of other functions in finiteelement methods or through quadrature formulas [BM92], or as a unifying frame for Fast-Fourier transforms [Nus82]. Polynomials are found at the crossroads of science: computer-aided design relies heavily on geometric splines, polynomials arise naturally in finance [Ack17], in biology [MY20], etc. It is actually easy to find more than 50 different families of polynomials named after mathematicians and that play a central role in various applications.
In most of the applications, having the fastest evaluation algorithm for a given level of accuracy is of the utmost importance.
In this article, we propose a novel approach to evaluating complex polynomials in the case of fixed precision floating-point arithmetic. Our algorithm is designed for better speed without compromising precision, not directly for improving the precision of the results (though, for a given cost of computations, it may be used to achieve a higher precision than can be reached with the current, more costly, algorithms). The algorithm is designed for repeated single-point evaluations. A typical application is Newton's method to find one single root, where the sequence of evaluation points is not initially known. However, the algorithm can also be used as an embarrassingly parallel multipoint evaluator, which makes it highly versatile.

Existing evaluation schemes
Let us review briefly the state of the art regarding polynomial evaluation.

Definition 1 (Complexity).
• Let us denote by V d the arithmetic complexity (number of arithmetic operations, with the convention that 1 operation is a multiplication followed by an addition * ) of evaluating a polynomial of degree d. We will denote by V d (k) the arithmetic complexity of simultaneously computing k values of a polynomial of degree d.
• When all the computations are performed with a fixed precision of p bits * * , we denote by V d (k, p) the corresponding bit complexity * * * (number of bit operations).
For example, M (128) 10 3 bit instructions for a Toom-Cook multiplication. The choice between these methods is usually driven by the competition between the value of p and the size of the hidden prefactor. The acceptable level of technicality in the code can also be taken into consideration.

Single point evaluation
Hörner's method ensures that, in general, V d = d and V d (1, p) = M (p)d when the polynomial is defined by its coefficients. Since Ostrowski [Ost54] and Pan [Pan66], it has been well known that for evaluating a complex polynomial P ∈ C[X] of degree d at a given point z ∈ C, d multiplications and d additions are both necessary and sufficient. That is, Hörner's scheme is optimal for one point-wise evaluation of a general polynomial.
Some classical evaluation schemes offer a similar order of complexity with more balanced and better parallelizable intermediary computations to improve numerical stability and take advantage of modern hardware, like Estrin's divide and conquer method [Est60], [Mor13], which is e.g. implemented in the Flint library [HJP13]. As a side note, Hörner's method is at the heart of a beautiful graphical construction for finding the real roots of a polynomial, known as Lill's method [Kal08].
For iteratively defined polynomials, i.e. a family P n+1 (z) = Q(P n (z)), evaluation is obviously more efficient: in this case, one gets V d = O(r log r d) where r = deg Q. Similarly, any intermediary power (z k ) 0≤k≤d of z can be computed recursively in O(log 2 k). In particular, sparse polynomials that contain only σ non-zero coefficients can be evaluated in this fashion in O(σ log 2 d) operations.

Multi-point evaluation
If the same polynomial has to be evaluated repeatedly, one obviously seeks to obtain V d (k) k V d and there are better strategies to reduce the average computing time. Knuth [Knu62] proposed a preprocesing based on finding all the zeros of the odd part of the polynomial (with Eve's variant [Eve64] in the general case) that gains a factor of 2 for the number of multiplications. It then brings down the cost of subsequent evaluations to V d = [ 1 2 (d + 4)]. A common case where simultaneous evaluation brings a substantial benefit is the evaluation of trigonometric polynomials along a regular mesh on the unit circle. Computing the values P Ä e 2iπk/(d+1) ä for 0 ≤ k ≤ d where P (z) = a j z j can be performed by Fast-Fourier Transform (FFT) algorithms [DL42], [CT65], [Roc00] in V d (d) = O(d log d) operations by taking advantage of (i.e. factoring) the matrix structure with w = e 2iπ/(d+1) . As an evaluation algorithm, this method brings down the average cost per computed value to O(log d), the price to pay being that one single evaluation cannot be (efficiently) performed alone. The fact that the FFT is numerically well behaved [PST02] and essentially involutive brings evaluation and interpolation to an equal footing and is the key that unlocks most of its applications. Note that there are variants of the FFT method for computing approximations of the values of P along a non-uniform mesh [PST01]. The evaluation points are however constrained to remain on a circle. Anticipating on our algorithm (see Section 4), let us point out that there is indeed a benefit in sorting the evaluation points according to the size of |z| and that we are able to discard the geometric restriction of cocyclicity.
Fast multipoint evaluation on a general set of points is possible and relies on a few standard tricks in polynomial arithmetic, which we first recall briefly. Polynomial multiplication can be computed in O(d 1.585 ) with Karatsuba's algorithm. It is based on the identity (P z k + R)(Qz k + M ) = RM + ((P + R)(Q + M ) − RM − P Q) z k + P Qz 2k , which boils down to 3 multiplications of smaller polynomials (recursively optimized with k d/2) and coefficient shifts. Over C[X] or more generally if the field admits a discrete Fourier transform, one can use the FFT to conjugate the multiplication of polynomials to a pointwise multiplication of enough interpolation points on the unit circle, with an overall cost of O(d log d). Next, fast division is based on the reversal of the order of the coefficients of the polynomial P , i.e. ‹ P (z) = z deg P P (1/z) and the identity The quotient M (z) can then be computed as ‹ P (z)· ‹ Q(z) −1 mod z deg P −deg Q+1 . The series expansion of the inverse is computed recursively with Newton's method in C[[X]]: which ensures that J n (z) = ‹ Q(z) −1 mod z 2 n . Ultimately, this algorithm brings the overall cost of computing the division of a polynomial P to O(d log d) where d = deg P .
The fast multipoint evaluation algorithm allows us to compute simultaneously k values with a cost of V d (k) = O((d + k) log d log kd). The central idea is a divide and conquer recursion. One splits the evaluation points in two families Z 1 , Z 2 of size k/2. With the previous fast division scheme, one computes remainders modulo polynomials that vanish either on Z 1 or on Z 2 . The problem is thus reduced to the evaluation of the two remainders on sets of points that are half-sized: Conversely, one can reuse the structure of the intermediary computations to interpolate with a similar total cost, i.e. compute the coefficients of P from the values P (z j ) at k = d + 1 distinct points (z j ) j=0,...,d . The method can be refined [Pan95], [Rei99] to improve the poly-logarithmic factor in the arithmetic complexity.
For finite precision arithmetic, i.e. approximations of order 2 −p , advanced algorithms reach a theoretical bit complexity of V d (d, p) = O(d(d + p + ω) log F ) to compute the values of P (z) at d complex points. In this formula, one takes ω ≥ 0 such that |z| + max |a j | < 2 ω and log F denotes logarithmic factors; see [Sch82], [KS16] for more details.
Of course, the comparison with the standard multi-point evaluation is not clear-cut because bit complexity depends on the size of the data while arithmetic complexity does not. Theoretically, these variants thrive when d + p + ω M (p) i.e., roughly speaking, when for example d p 1.585 . Even though the multiplicative constants and logarithmic factor are either large or hard to track, one can expect these algorithms to be competitive for moderately large degrees (d 10 4 ) and a substantial fixed precision (128 ≤ p ≤ 300).

Practical considerations beyond arithmetic complexity
While the fast multipoint methods optimize the overall cost of multiple evaluations to V d (d) = O(d(log d) 2 ), the number of operations that are involved in the computation of one single value (or one single coefficient in the case of interpolation) exceed the number of naive operations that would be required to compute that value alone. Mechanically, one can thus expect a loss of precision.
A very detailed analysis of the numerical instabilities [KZ08] points out the Wilkinsontype [Wil84] expansion of j=1 (z − j) as the main culprit, which leads to ill-conditioned input to the subsequent polynomial divisions. On the other hand [KS16] exploit the fact that these same divisors are monic to ensure stable divisions when the required precision (in number of bits) dominates the degree.
In many situations, in order to guarantee the numerical stability of the the multipoint method, the precision of the numbers has to be much larger than d. For example, S. Köhler and M. Ziegler [KZ08] conclude that, for a specified level of accuracy, it is usually necessary to increase the precision of intermediary computations to the point where the benefit over Hörner's V d (d) = O(d 2 ) naive scheme is not significant. The conclusion of A. Kobel and M. Sagraloff [KS16] is more nuanced, as they insist on the fact that the extra precision is only required for intermediary computations. However, as one may need more than O(d log d) bits of memory per coefficient, it quickly becomes impractical as d increases (see [KS16,Corollary 8]).
On modern computing machines, even on super-computers, the workload is often dominated by data movements (disk, memory and cache access) and not by the computing power (usually measured in floating-point operations per second, or FLOPS). When implementing the multipoint evaluation algorithms, the memory size limitations impose rather tight bounds on the degree (say d 10 6 ). This raises the question of how to compute efficiently the values of a giga-polynomial. Our algorithm (see Section 4) does not present this limitation.

Alternatives
Let us close this tour of the literature by mentioning briefly some less common methods of evaluation, which have their own niche of applications.
If working with extended precision is not an option, various methods based on a compensation of Hörner's Algorithm [SW05], [LGL06], [Sut07] can improve the precision of the standard evaluation scheme for a moderate increase in complexity.
Choosing another basis instead of the canonical monomial basis of C[X] may provide better numerical stability. Evaluation algorithms on Newton's (interpolation) basis have similar complexity to the ones exposed above [BS05], sometimes even with better constants. We refer to [Far12], [Far08] for an in depth review of the benefits of the Bernstein basis and its industrial applications. The complexity of evaluating the Bernstein basis functions has recently been improved in [CW21] and is now O(d), which is still huge compared to evaluating z d but makes it a viable option for e.g. d 10 3 . In this article, we will not investigate further the question of generalizing our algorithm to other bases.
For the sake of completeness, let us also mention that better performances as low as O( √ d) can be achieved for non-scalar evaluations [MP73], [Fas19], i.e. if one computes polynomials of matrices where the complexity of scalar operations is simply discarded. However, these algorithms do not bring any improvement upon Hörner's method, when they are applied to the evaluation of scalar polynomials.

New evaluation scheme
In this article we propose a simple algorithm and its practical implementation as a C library [MV22] that brings down the average cost for the repeated evaluation of all polynomials and never exceeds Hörner's complexity in general. Figure 1: The average speed gain of our FPE (Fast Polynomial Evaluator ) algorithm over Hörner's method for computations with p = 53 bits (using MPFR numbers; [MPFR]). The solid curve corresponds to bound (57) that follows from Theorem 4. The dashed one corresponds to the example at the end of Section 4.3 where one evaluation point has the same complexity as Hörner, but the average complexity is favorable. Note that lacunary polynomials may lead to even higher gain factors. The data points are actual benchmarks (see Section 8.3). The case of the red squares (half-circle) is studied in detail in Section 5.3.
More precisely, our algorithm, called FPE or Fast Polynomial Evaluator (see Section 4.3) pre-conditions P ∈ C[X] in time O(d log d) with a low multiplicative constant that does not depend on the precision. Subsequent evaluations of P (z) with a fixed precision of p bits are performed in average arithmetic complexity The constant is small and explicit and is given by (40) below. The memory requirement is O(dp). The average of the complexity used in (2) is taken with respect to points z ∈ C weighted by the spherical area of C. A similar estimate holds for real polynomials and a uniformly distributed evaluation point along the circle R ∪ {∞}.
As illustrated in Sections 5.3 and 8.3, for many particular classes of polynomials, in particular for sparse polynomials or those with random coefficients confined in a bounded region of C, our algorithm performs much better than the upper bound (2).
The FPE algorithm has many interesting features. One can guarantee that the result is as precise as Hörner's method. Pointwise, the complexity of FPE does not exceed that of Hörner (i.e. V d ≤ Cd). In case of equality for some z 0 ∈ C (see Remark 15), one has which is even more advantageous than (2). This radical difference between the pointwise and the average complexity is a strong incentive in favor of studying averages. The FPE algorithm is embarrassingly parallel and can be implemented on any set of evaluation points, without constraints of size or of geometric structure. New evaluation points can be added on the fly. These properties make the FPE algorithm particularly well suited for a root finding scheme with Newton's method.
For low-precision computations the theoretical speed factor of FPE over Hörner's method is illustrated in Figure 1. Benchmarks of our implementation will be presented in Section 8 and, in particular, the analysis of the influence of the preprocessing phase over the global cost (it remains minimal).
The cornerstone idea at the foundation of the FPE algorithm is lazy polynomial evaluation (see Section 4.1): adding two finite precision numbers is only necessary if their orders of magnitude are close enough that their bits will interact. Monomials tend to have extremely diverse orders of magnitude, which means that the value of a polynomial at a given point is dictated by only a small subset of its monomials. The second ingredient is a geometric selection principle, i.e. the ability to identify this parsimonious representation with geometric tools, which, in practice, boil down to the computation of the concave cover of a dataset (see Figure 7). This article features a detailed analysis of the complexity (see Theorem 4) and a precise error analysis (see Theorem 3 and Figure 22) of the FPE algorithm. Both are put to the test in systematic benchmarks presented in Section 8.3. The preview offered in Figure 1 illustrates the extent to which we have explored the theoretical and practical envelope of this new algorithm.

Structure of the article
The structure of the text is the following.
In Section 2 we detail the fundamentals of finite precision arithmetic for a general audience and discuss the specificities of complex numbers. Subsection 2.4 is dedicated to the various notions of closeness in finite precision (equivalence, adjacency and similarity modulo phase-shift), which play a central role in the proof of the correctness of our algorithm.
In Section 3 we briefly introduce some geometric tools that will be needed to state the FPE algorithm and prove its complexity. This section contains only definitions and statements; the proofs of the corresponding theorems can be found in Appendix A.
In Section 4 we describe and analyze the FPE algorithm and state our two main results. The correctness of the algorithm and the associated error analysis is Theorem 3. The result regarding complexity is Theorem 4. The proofs are done in the next two sections: Theorem 4 is proved in Section 5 and Theorem 3 in Section 6.
Section 7 explores a few examples of possible applications of the FPE algorithm and should be of general interest. Section 8 is dedicated to presenting our implementation in the C language, which we are publishing [MV22] as an open source project. Our implementation uses both machine floating-point numbers and MPFR arbitrary precision numbers (see [MPFR]). Appendix C contains a listing and description of the tasks that can be performed with it. Extensive numerical benchmarks are presented in subsection 8.3 and confirm the theoretical predictions regarding the complexity and error analysis of the FPE algorithm.
In order to keep this article accessible to the widest possible audience, we provide comprehensive definitions of all notions and fully detailed proofs. We also tried to keep the sections as independent as possible. Overall, the key ideas are of a geometric nature. There are strong similarities between the geometrical reasonings of Section 3 and Appendix A and the presentation of the algorithm in Section 4 (for example, compare Figures 6 and 7, or 8 and 24). It is our belief that one may enlighten the other.
However, a reader interested in understanding the algorithm, but not the proof of its correctness and complexity, may safely read only the beginning of Section 2 and skip Subsection 2.4, then read Section 3 before proceeding to Sections 4, 7 and 8. The more theoretical Sections 5, 6 and Appendix A may be skipped.
Finally, as a convenience for all readers, Appendix B recapitulates the notations used throughout the article.

A bit of finite precision arithmetic
In this section we introduce notations and tools that will be useful in the statement of the algorithm and in all subsequent analysis. We refer to, e.g., [Ma18], [Gol91] or [754] for further details on finite precision arithmetic.

General considerations
Let us start with a word of caution: the exact evaluation of complex polynomials at arbitrary points is not possible in practice. Attempting to evaluate an explicit polynomial (that is, given by its coefficients) near one of its roots will produce a large cancelation of the digits. For example, P (z) = z − z 0 evaluated at z such that |z − z 0 | < 2 −n |z 0 | produces a result of order at most 2 −n |z 0 |, effectively losing n leading bits from the precision that was used to express z 0 and z. More generally, when P (z) has many terms, cancelations can occur not only at the roots but among all polynomial subexpressions of P (see Figure 2). To guarantee that the result has any number of significant exact digits, unbounded precision would have to be used for intermediary computations.
As this is not practical, we will focus only on computations done with some fixed precision p. Our algorithm does not attempt to produce results that are more accurate than those of Hörner's method: we want to produce results of similar accuracy, only faster. The error analysis of Hörner's scheme is classical and we refer the reader to [Oli79], [M83]. For a more general analysis with recursive basis functions, see [BJS13].
When using fixed precision numbers, additions and subtractions are the main sources of errors because they can produce a cancelation of the most significant digits. After such an occurence, the relative uncertainty is multiplied by 2 s , where s is the number of canceled bits. Fortunately, having s bits canceled is conditioned by the fact that the numbers have the same scale and then cancelation only occurs with a probability 2 −s .
Our algorithm exploits the limitations of finite precision arithmetic to discard unnecessary computations and thus obtain significant gains on the computation time.
The name stands for unit in the last place. To ensure a unique representation of all real numbers we always assume a rounding to the nearest representable number and choose a rounding away from zero at the tie.
A key observation for additions in finite precision p is that where n and m are the respective exponents of ξ and η. This means that by simply reading the values of n and m and comparing them to p, we can avoid costly operations, especially if computing one of the terms ξ or η requires additional steps as is the case when they are monomials a k z k .

Scale of a complex number
Let us define the scale of a number z ∈ C * by where α is the floor of α. By convention, s(0) = −∞. The scale is a logarithmic representation of the order of magnitude of z. For example, with the notations of (3), one has s(ξ) = n, i.e. for floating-point numbers, the scale coincides with the exponent; moreover ulp(ξ) = 2 s(ξ)−p−1 .
Lemma 2. For any z, z ∈ C and n ∈ Z, one has : s(z ± z ) ≤ max{s(z), s(z )} + 1 (13) and, for a family z 1 , . . . , z N ∈ C : Due to the cancelation of significant bits, the scale of a sum may be much smaller (even −∞) than the largest scale of the terms involved.

Equivalence and adjacency modulo finite precision
The topology induced by finite precision arithmetic is surprisingly subtle. In this subsection, we introduce three distinct binary relations that express the proximity between a floating-point number and a real or complex number. The complex case is even more subtle and is dealt with last.
Definition 3. For x, y ∈ R we use the notation x = p y to say that x and y have the same representation as floating-point numbers with precision p ∈ N * . It is an equivalence relation; the equivalence class of a floating-point number ξ > 0 is In general, the equivalence class (15) of a floating-point ξ > 0 is the interval An exception occurs at scale turnover where the class is asymmetric. For example, the consecutive numbers a = 2 n (1 − 2 −p−1 ) = 2 n × 0.111 . . . 11, b = 2 n = 2 n+1 × 0.100 . . . 00 and c = 2 n (1 + 2 −p ) = 2 n+1 × 0.100 . . . 01 satisfy ulp(b) = c − b and ulp(a Figure 3. When ξ < 0, the usual convention rounding away from zero at the tie implies that the interval (15) is flipped over and the ± signs must be reversed.
We say that equivalence classes for = p are adjacent if their respective closures in R have a non-empty intersection. The corresponding floating-point numbers are called adjacent too.
Definition 4. For x, y ∈ R, we denote by x p y if the floating-point representations of x and y with p bits are identical or adjacent.
Though not transitive because of the obvious overlap between the sets of real numbers that are adjacent to consecutive floating-point numbers, this relation is symmetric and simplifies the handling of scale turnover. For example, one can always find real points that are arbitrarily close to one another but whose floating-point representations are distinct; these points are however adjacent. See Figure 4. When x is a real number and ξ is its floating p-bit representation, it satisfies |x−ξ| ≤ 1 2 ulp(ξ) and s(x) ≤ s(ξ) ≤ s(x)+1. Moreover, s(x) = s(ξ) occurs only at scale turn-over, when 2 n (1 − 2 −p−2 ) ≤ x < ξ = 2 n for some n ∈ Z. Conversely, if ξ is a floating-point number such that |x − ξ| ≤ 1 4 ulp(ξ) then ξ is the p-bit representation of x. Along the real line, the characterization of adjacency in terms of scale is the following.
Lemma 5. Let x, y ∈ R and p ∈ N * . One has and, conversely, If x and y are p-bit floating-point numbers, then Proof. If x and y are not adjacent real numbers, then x and y are separated by a whole p-bit equivalence class of some floating-point number ξ (see Figure 4), i.e.
For complex numbers, the situation is more complicated because of phase shifts. In the real case, the only phase shift possible is a sign change, which does not affect precision; in particular, there are only 2 p+1 finite precision numbers for a given scale. The direct extension of the adjacency relation z p z as the conjunction of Re z p Re z and Im z p Im z can lead to an extreme scale imbalance between the real and imaginary parts, which is not compatible with phase shifts (i.e. complex rotations). Complex numbers whose argument is close to kπ/2 (with k ∈ Z, i.e. near the axes) have one component artificially over-resolved compared to the other one; consequently, there are infinitely many * finite precision numbers with a given scale (see Figure 5).
A similar instance of the same issue occurs if one choses two complex numbers close to the diagonal whose real and imaginary parts are p-bit adjacent, e.g. : . Once we rotate them to bring them along the real axis (or, equivalently, if we multiply by (1 + i) −1 = 1 2 (1 − i), which * In practice, the number of finite precision complex numbers at a given scale is limited by the extreme negative exponent value authorized in the implementation.
is an exact 0-bit number), their imaginary parts will not be adjacent anymore and will instead be separated by infinitely many scales: 1 + i2 −p−r p 1.
To address this issue, we introduce a looser version of p , which is based on the scale and inspired by the property (16) above.
Definition 6. For a pair of complex numbers z, z , we note z ≈ p z if and only if (20) We say that z and z have similar p-bit floating-point representations, modulo approximate phase-shift invariance.
Let us point out that, because s is an increasing function, the following criterion holds : Conversely, according to (8) Thanks to Lemma 5, the relation z ≈ p z implies adjacency when z, z ∈ R. If ζ is a finite precision number close to the x-axis, then {z floating-point such that z ≈ p ζ} ⊂ {x + iy ; x p Re ζ and y p−δ Im ζ} with δ = s(Re ζ) − s(Im ζ) ∈ N. Indeed, denoting ζ = ξ + iη and z = x + iy, (9) ensures in this case s(ζ) ≤ s(ξ) + 1 and s( This configuration is illustrated in Figure 5.
Remark 7. To build a finite arithmetic theory that is truly rotation invariant, one should use ball arithmetic. Its superiority is demonstrated in [MV], to provide estimates that remain significant after many iterations of a conformal map.
Conversely, if (22) holds, then In particular, if y = p −x, all bits cancel out and s(x + y) = −∞. The inequality (22) expresses that at least q + 1 leading bits, including the implicit leading ξ 0 = 1 in (3), cancel each other in the addition of x and y. It is therefore possible to estimate the loss of precision by comparing the scales of the operands with that of the result: when this value is the exact number of leading bits lost in the operation.

A bit of geometry
In this brief section we introduce some geometric tools that will be useful in the proof of Theorem 4. For any real valued map g defined on [0, 1], let us define the horizontal strip of height δ > 0 under the graph of g as follows: The set S(g, δ) is the intersection of the subgraph of g with the epigraph of g − δ.
where the segment These definitions are illustrated in Figure 6. One has L(f, δ, θ) ≤ 1 cos θ because |A B | = L(f, δ, θ) cos θ ≤ 1. The two following statements are key for estimating the complexity of the algorithm that is presented in Section 4.
Moreover, there exists a function f 0 ∈ C that satisfies a lower bound > 1.1128 √ δ for δ small enough. In the same conditions, one has We use the following variant for computing averages on the Riemann sphere.
Subsequently, we will apply these results to a function f that is a renormalized concave cover of the scales of the coefficients and a gap δ that depends on the precision of the computations (see equation (34) and Figure 10). We postpone the proof of Theorems 1 and 2 to Appendix A.

The FPE algorithm
We now focus on a new Fast Polynomial Evaluator algorithm, or FPE for short.

Key idea: lazy polynomial evaluation
Consider a polynomial P ∈ C[X], which we identify to the entire function The general idea of the FPE algorithm is to perform as many lazy additions (5) as possible. We cut down the cost by not computing the monomials that will have no influence on the final result. More precisely, with minimal overhead (preprocessing), we identify the favorable cases where it will be safe to perform lazy additions (this is the non-trivial and novel point, as the value P (z) is not yet known), and in the remaining cases, we apply a variant of Hörner's method. All unnecessary monomials are thus left out.
Example 9. The simplest case study of the FPE algorithm is the following. Let P (z) = 1 + z and assume we perform computations with precision p. If s(z) < −p − 1 then P (z) p 1 and if s(z) > p + 1 then P (z) p z. So in these two cases, we get the result for free. An actual computation is only required in the remaining case, i.e. when |s(z)| ≤ p + 1. If z is uniformly distributed on the Riemann sphere C, then on average (see Section 5.2), the result is computed in operations instead of 1. In this simplest case, the gain is negligible.
It turns out that the lazy evaluation method performs steadily better in terms of arithmetic complexity (and speed) as the degree of P gets higher. Compared to Hörner's scheme, the gain is substantial (e.g. for computations in machine precision) and holds for every polynomial once we average out over all the possible scales of the evaluation point. Before presenting the general case, we illustrate our algorithm and this phenomenon below, on a polynomial of degree 10 (see Figure 7).
Let us point out that the FPE algorithm thrives in the case of multiple evaluations, because the initial analysis needs not be repeated. We illustrate this fact on the example below too, by showing how the analysis at the points z = ±1 can easily be transposed to an arbitrary value of z (see Section 4.2.2). By construction, subsequent evaluations cannot exceed the complexity of Hörner's scheme and will, on average, be much better.
Contrary to the FFT or the Fast Multipoint Algorithm, the FPE algorithm is local in the sense that the evaluation at a given point is independent of the precise computations that are needed to evaluate P at another point. One can thus expect the algorithm to have the same numerical stability as Hörner's method. The memory requirements for each evaluation are also minimal because the overhead storage is negligible in comparison to that of the coefficients. For example, our algorithm will thrive in implementations of Newton's method to find a single root of a polynomial of large degree because the list of evaluation points is, obviously, not known in advance.
Finally, without increasing the arithmetic complexity, it is possible to complement the result P (z) with a confidence estimator that indicates how many of the p bits may have suffered from cancelations (for details, see Remarks 8 and 12). This feature is part of our implementation [MV22] (see Section 8, in particular Figure 23). It may help if finding the proper value of p is part of the problem, e.g. in the implementation of Newton's method with a dynamically adjusted precision. Note that changing the value of p will require a new preprocessing of the polynomial (see Remark 17).
In preparation for the general case, let us introduce the following notation.
Definition 10. The scales of the coefficients are modeled by the function We will denote by Ù E P the concave cover of E P , that is the minimal real concave function

A simple example detailed
We analyze an example depicted in Figure 7, which represents E P and Ù E P for a particular polynomial P of degree 10 with non-zero coefficients. The scales of the coefficients are readable on the graphic; the actual values of the phase of each coefficient are irrelevant to the discussion. To keep this example simple, we compute P (x) with a fixed precision p = 6. The strip S( Ù E P , p) defined by (23) is a polygonal band of vertical thickness p. For λ = tan θ, let us also denote by L −λ the longest segment of slope λ contained in S( Ù E P , p), that is with L defined by (24); see Section 3 if necessary. The reason for the sign convention will appear in Section 4.2.2.
The horizontal segment L 0 isolates the (three) coefficients of Q 0 (z), which is the suitable reduction of P (z) when the evaluation occurs on the unit circle.  Figure 7 by the affine map A λ (defined below) for λ = −3. The image of the concave cover Ù E P is Ù E P + λId. The horizontal segment A λ L λ isolates the (five) coefficients of Q λ (z), which is the reduction of P (z) when |z| = 2 λ .

Evaluation on the unit circle
For the previous example, let us first consider the case where the evaluation point z satisfies |z| = 1. For k ∈ 0, 10 , we have s(a k z k ) = s(a k ). To select the largest monomials in P (z), let us consider the segment L 0 on Figure 7, which is the horizontal segment situated at a scale p below the scale of the monomial with maximal scale. The monomials we keep are the ones above L 0 . In this example, we get Q 0 (z) = a 6 z 6 + a 8 z 8 + a 9 z 9 .
Let us now check that the floating-point value of P (z) coincides indeed with Q 0 (z). If k / ∈ {6, 8, 9}, then s(a k z k ) ≤ 18 and the inequality (14) with N = 8 implies On the other hand, unless there is an exceptional cancelation of the leading term while computing Q 0 (z), one has s(Q 0 (z)) ≥ 30. In particular, one has s(P as long as the leading terms of Q 0 (z) do not cancel each other.
Let us now investigate the possible cancelations within Q 0 . As computations are restricted to p = 6 bits, the last bit of both a 8 z 8 and a 9 z 9 represents a rounding interval of the real line of radius 1 2 × 2 30−6−1 = 2 22 . This is the best error bound for a 8 z 8 and a 9 z 9 that we can hope for. Thus, even if a cancelation of the most significant bits occurs, the error when computing a 8 z 8 + a 9 z 9 cannot be bounded to less than 2 23 . Consequently, as the bound for the error can only increase with more terms in the sum, the uncertainty on the value of Q 0 (z) will exceed 2 23 . On the other hand, we have checked that |P (z) − Q 0 (z)| < 2 22 , so this difference is always smaller than the error bound in the computation of both Q 0 (z) and P (z). Therefore, one can claim that, when |z| = 1, the value Q 0 (z) is always a good floating-point substitute for P (z), i.e. within the error bounds for 6 bits of precision throughout the computation, even if significant bits are canceled out.
To get a geometric feeling in this instance of our algorithm, observe that in the computation of Q 0 (z) we have only used terms that are above the line L 0 , which is is the longest horizontal segment contained in S( Ù E P , p). This line lies at the scale level max k s(a k z k ) − p. Note also that we did not use a precision larger than p for any intermediary result.

Towards the general case
Let us continue the analysis of the previous example for a general evaluation point, i.e. when z ∈ C. As P (0) = a 0 is immediately available, let us assume z = 0 and consider The product formulas (11) and (12) imply that, for all k ∈ 0, d , we have As before, we seek a simpler polynomial Q λ such that it is sufficient to evaluate Q λ at z instead of P when the computations are done with p = 6 bits. In order to visualize the polynomial Q λ , let us consider the image of Figure 7 by the following affine map of R 2 By definition, A λ maps lines of slope −λ to the horizontal. For λ = −3, we consider the reduced polynomial In other words, we select the powers k ∈ 0, d such that s(a k )+kλ is above the horizontal line A λ L λ .
Coefficient a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 The precision loss (4 bits out of 6) may seem significant in this example. However, in general, the loss is capped by s(d) + 3, which means that an offset on the thickness of the S( Ù E P , ·) strip will be enough to deal with the general case.
Remark 11. In the general statement of the FPE algorithm (see (44) in Section 4.3), we will use the scale threshold Figures 7 and 8) to prevent interaction between Q λ and P − Q λ and to secure upper bounds of s(a k z k ).
We may now link the statement of Theorems 1 and 2 to our algorithm. As the sum of two concave maps is concave, the map κ → Ù E P (κ) + λκ is concave. Therefore, there are at most |A λ L λ | + 1 terms in the reduced polynomial Q λ , where |L| is the length of a segment L. Observe that where θ ∈ − π/2, π/2 satisfies λ = tan θ. To estimate the average reduction in complexity of our algorithm over Hörner's scheme, we are interested in averaging the number of monomials of P (z) that are ultimately evaluated. We will therefore compute the average value of L(f, δ, θ) cos θ where f ∈ C is a renormalized version of Ù E P , defined for x ∈ [0, 1] by Depending on how the values z are chosen in C, various weights for θ ∈ (−π/2, π/2) are used (see Section 5.2).

Statement of the algorithm
We are given a precision p ≥ 1 and a polynomial expression We will also assume that a 0 = 0, otherwise we reduce the problem to a lower degree polynomial z −k P (z) for some k ≥ 1.
We can formalize our evaluation algorithm FPE p as follows. Each non-trivial operation has its time (bit) complexity marked as a comment on the right. The 3 for the definition of the scale functions s and the end of Section 4.1 for that of the concave cover Ù E P of s(a j ).

Statement of the main results
In the following subsections, we describe in detail each step of the algorithm. Subsequently, we prove its correctness, i.e. the statement of Theorem 3, and we compute the time complexity of FPE p as stated in Theorem 4 and equation (49).
Theorem 3. Given P as above and a precision p ≥ 1, for each z ∈ C * and λ = log 2 |z|, there exists a polynomial subexpression Q λ of P such that (see Section 2.4): where the number of canceled bits c ≥ 0 is defined by and M = max where ‹ a k is the p-bit floating-point representation of a k and where , r and G p are computed by the algorithm FPE p described above. The number of monomials of Q λ satisfies where the average is taken with respect to the uniform distribution of z ∈ C on the Riemann sphere. Moreover, for any d and p, there exists a polynomial P for which G p = 0, d and such that avg Remark 12. The proof of Theorem 3 ensures that The precision claimed by (35) is thus equivalent to that of Hörner's scheme, if all coefficients (limited to p bits) had been kept. Note also that, even though the number of cancelled bits c is defined by (36) and thus depends on the exact value of P (z), it is possible to give a precise upper-bound of c by using Remark 8 for each addition that occurs in the computation of (37).
Regarding complexity, the main result is as follows (see Figure 1).
Theorem 4. Given a polynomial P ∈ C[X] of degree d ≥ 1 and a bit precision p ∈ N * , the preconditioning phase of algorithm FPE p is performed on P in time d + 2d log 2 d and requires O(dp) in memory. Subsequently, for z uniformly distributed on the Riemann sphere C, the average evaluation time of P (z) by the algorithm FPE p is less than where M (p), recalled in (1), denotes the time of one multiplication followed by an addition of two floating-point numbers with precision p. If P ∈ R d [X] and x is uniformly distributed on the circle R = R ∪ {∞}, the average evaluation time of P (x) by the algorithm FPE p is less than In both cases, the bit complexity of the evaluator never exceeds 2 + 3 log 2 d + M (p)d.
Remark 14. Let us point out that the uniform average over C (or R) is unfavorable to our algorithm. Near the poles z = 0 and z = ∞, our algorithm will drop most terms and will thus be very quick. However, a uniform average does not favor those regions: the area of the region |z| > 10 (or |z| < 1/10) represents about 1% of the area of the sphere, which is of the same order of magnitude as that of the annulus ||z| − 1| < 10 −2 . Using the techniques exposed in Section 5, one can compute the average complexity for any particular distribution of evaluation points; for example, the case of a uniform distribution on D(0, 1) is treated in Remark 19 below, estimate (52). It is also possible to refine the estimate if the distribution of the coefficients of P is known (see e.g. Figure 12 for Chebyshev polynomials).
A point is worth underlining: if the algorithm FPE p encounters one "bad" case where one evaluation has the same complexity as Hörner, then, on average, it will perform much better than (40). More precisely, let us assume that one particular choice of z 0 with log 2 |z 0 | = λ 0 leads our algorithm to evaluate all the monomials of P (z 0 ), which is the worst case possible. Of course, for such a polynomial, our algorithm would not outperform Hörner if we were to evaluate only on z in an annulus |z| λ 0 . However, from this shortcoming, we learn that the graph of the concave cover of s(a k ) + kλ 0 rescaled to [0, 1] (see (34) and Figure 8) is comprised between two horizontal lines c and c + δ, i.e. the modulus of the coefficients of P are, roughly speaking, varying exponentially. If we briefly anticipate the computations of Section 5.2, the average number of terms when z is uniformly distributed on C can be estimated with a simple weight (see Figure 10) : The computation (75) from Example 1 then provides an explicit bound: This means that, if our evaluator performs once as poorly as Hörner, then it will, on average, perform as O(pM (p) log d) if the evaluation points are chosen uniformly over the Riemann sphere C and log 2 d ≤ p d or instead as O(M (p) log 2 d) if p ≤ log 2 d. This is a much better behavior than the one claimed by Theorem 4 in general and it is the best that we have observed in practice (see Figure 1 and, for details, Section 8.3).
Remark 15. More generally, if the coefficients of P are a union of a few long geometric progressions (even possibly intertwined), the graph of Ù E P will be composed of only a few piecewise straight lines, say N d. Each straight line will only be visible on a finite range of values of |z| and will contribute a logarithmic complexity bounded by (42). The overall average complexity of the FPE p evaluator will then be bounded by if the evaluation points are chosen uniformly over the Riemann sphere C (or R in the real case).
For further details and the construction of an example that saturates the upper bound (40), see Section 5.3.

Complexity analysis and proof of Theorem 4
In this section, we describe the details of the algorithm FPE p and prove Theorem 4 regarding complexity.

Analysis of the preconditioning phase
We describe briefly the computation of the concave hull Ù E P . The first step is standard and consists in obtaining an enumeration (k n ) of 0, d to sort the values s n = s(a kn ) in decreasing order, i.e. such that for all n ∈ 0, d − 1 , In case of equality, k n is chosen in increasing order (i.e. k n ≤ k n+1 if s(a kn ) = s(a k n+1 )). This step can be performed in d log 2 d operations. Observe that k 0 = argmax Ù E P .
Lemma 16. For n ∈ 0, d , we construct a sequence of concave maps E n : [ n , r n ] → R such that, for all n, [ n , r n ] is the convex hull of {k 0 , . . . , k n } and ∀k ∈ n , r n , s(a k ) ≤ E n (k) ≤ Ù E P (k) .
Constructing E n+1 knowing E n is performed in log 2 n operations.
Observe in particular that E d = Ù E P and that it is obtained in less than d log 2 d steps once (k n ) is known. An example of this construction is given on Figure 9. Proof. It is immediate that 0 = r 0 = k 0 , E 0 (k 0 ) = s 0 and that the graph of E 1 is the segment [(k 0 , s 0 ), (k 1 , s 1 )]. Assume now that E n is constructed. Then s n+1 ≤ min E n and k n+1 / ∈ { n , r n }. If k n+1 ∈ ( n , r n ) we let E n+1 = E n (on Figure 9, this is the case for E 5 = E 4 and E 6 = E 5 ). Otherwise, by symmetry, we may assume that k n+1 < n and we build E n+1 as a left extension of E n (alternatively, if k n+1 > r n , then we would build E n+1 as a right extension of E n in a similar way). Let ‹ E n+1 be the extension of E n to [k n+1 , r n ] with the property that ‹ E n+1 is affine on the interval k n+1 , n + 1/2 . If s n+1 ≤ ‹ E n+1 (k n+1 ) (on Figure 9, this is the case for n + 1 = 4), we define E n+1 to coincide with E n on [ n , r n ], E n+1 (k n+1 ) = s n+1 and be affine on [k n+1 , n ].
If s n+1 > ‹ E n+1 (k n+1 ) (on Figure 9, this is the case for n + 1 = 7), the previous construction would not be concave. Therefore, we search to restrict E n to some interval [ n , r n ], n > n and to add a new segment that includes s n+1 , such that E n+1 is concave. We remove the segments whose slopes are smaller than the slope of the new segment (on Figure 9, we drop the segment of E 6 that is below the dashed line in E 7 ).
Precisely, let Σ n 0 , . . . , Σ n s be the lines containing a segment in the graph of E n , enumerated from left to right. As E n is concave, the y-coordinate of Σ n j ∩ {x = k n+1 }, denoted by y n j , is increasing in j. A binary search finds j such that s n+1 ∈ (y n j−1 , y n j ] in time log 2 n. The value n is the leftmost abscissa x such that (x, E n (x)) ∈ Σ n j . Then E n+1 is defined as affine on [k n+1 , n ] with E n+1 (k n+1 ) = s n+1 and E n coincides with E n on [ n , r n ].
As (s n ) n∈ 0,d is decreasing, zero coefficients are sorted last and are only treated when the graph is already complete. Indeed, as we have assumed that a 0 = 0 and that a d = 0, if s n+1 = −∞ then k n+1 ∈ ( n , r n ) and in that case we set E d = E n .
A final parsing of the list of s(a k ) is performed to mark the indices k such thats We denote by G p ⊆ 0, d the set of these good indices and by B p = 0, d \ G p the set of bad indices. This step has a linear time complexity. In subsequent evaluations, only the monomials (a k z k ) k∈Gp are kept. In what follows, we show that those associated with B p cannot influence the first p bits of the result. The set G p will be thinned even more during the evaluation phase, depending on |z|.
Let us emphazise that the complexity of the preconditioning does not depend on the precision p. If the coefficients a k are provided in machine floating-point numbers, obtaining s(a k ) is performed in constant time using hardware-accelerated functions. In the case of an arbitrary precision p, the value s(a k ) is already computed and stored in the number format and there is nothing to do. All computations for the preconditioning phase can thus be performed with machine floating-point numbers.
Remark 17. Let us mention a slight variant of our algorithm, which is based on the fact that the lines 2 and 3 of the algorithm FPE p are independent of the value of p. For certain applications, one could split the preconditioning in two parts. The computation of the concave map Ù E P could be done during the compilation (if P is known in advance) or at early runtime without any knowledge of p (if a low-precision version of P is availlable). Once the precision p is known, one will finish the preconditioning (i.e. determine the set G p , line 4 of FPE p ) in time O(d). Subsequent evaluations of P will be performed as before, using only lines 5-11 of FPE p .

Analysis of the evaluation phase
To compute k λ = argmax Ù E P + λId , observe that Ù E P + λId is concave. That is, its derivative (in our case the slope of the segments from some point (k, s(a k ) + λk) to the next one (k , s(a k ) + λk )) is decreasing. Therefore, a binary search finds k λ in log 2 d operations. The maximum value is Next, as Ù E P + λId has at most two monotone branches (separated by k λ ) we can perform a binary search on each of them to find respectively the two indices < k λ and r > k λ such that , r is the largest integer interval that satisfies Each of these searches costs at most log 2 d operations. Therefore lines 5-10 of FPE p cost 2 + 3 log 2 d operations, which is the first part of (40) in Theorem 4.
Let us now focus on the complexity analysis of the last step (line 11) of the algorithm FPE p . Formula (33) reads L(f −(tan θ) Id, δ, 0) = L(f, δ, θ) cos θ; joined with (46), it implies where f and δ are defined by (34) and λ = − tan θ (the minus sign reflects that positive slopes correspond to evaluation points z such that |z| < 1). The metric on the Riemann sphere C that is associated with a uniform probability measure is given by The corresponding volume element is For a radial function and r 2 = x 2 + y 2 , the volume element becomes Therefore, the average number of monomials that are required to evaluate a polynomial of degree d with our algorithm, when the point z = x + iy is chosen uniformly on the Riemann sphere C, is bounded from above by with (note that ω is even): Therefore, Theorem 2 implies with C ω < 1.9046 and whose exact numerical value is given by (28). Remark 18. Note that using the looser estimate ω(θ) < 0.46 cos θ (see Figure 10) and the numerical constant of Theorem 1 overshoots the value of C ω by 42%.
To conclude the evaluation of P (z), we compute z in 2 log 2 steps and then use Hörner's method to compute Q λ (z). The average arithmetic complexity of line 11 of the algorithm FPE p is thus bounded by while the bit/time complexity is M (p) times larger. Putting (48) and (49) together gives the last part of (40) in Theorem 4.
In the case of a real polynomial evaluated along the real line, the uniform measure on the circle R = R ∪ {∞} obtained by stereographic projection is With the change of variable log 2 |x| = tan θ, one gets Using (28) again provides The average complexity in the real-valued case is given by (49) with C ω instead of C ω .
Remark 19. One can easily adapt the computation to the case of a complex polynomial evaluated at z uniformly distributed over the unit disk D(0, 1). The weight becomes The asymmetry of ω D implies that one must restrict the integral of Lemma 23 to y ≥ 0, i.e. x 2 ≥ x 1 . Theorem 2 remains valid with The numerical value satisfies C ω D < 2.3141 and one can claim

Example that (almost) saturates the upper bound on complexity.
Because of the fast decay of ω(θ) as θ → ±π/2, it is not possible to reuse directly the lower bound obtained in Theorem 1 for the weight cos 2 θ. However, the examples of Section A.4 can be adapted easily to saturate the complexity of the algorithm FPE. Inspired by the second example, let us consider a polynomial P whose coefficients have a scale profile that follows a half-circle, for example: Reasoning as in Section 5.2 and using the maximality of , r in (46), we get r − + 2 d ≥ L(f, δ, θ) cos θ for any z ∈ C * such that log 2 |z| = tan θ and f , δ defined by (34). The average arithmetic complexity of the evaluation of P when z ∈ C (resp. z ∈ R) is bounded from below by or, respectively, the same integral with ω in place of ω. One can check easily that E P defined by (29) satisfies As d → ∞, the graph of f converges uniformly to that of x(1 − x), which is concave.
With the notations of Example 2 of Section A.4, the average complexity is thus asymptotically bounded from below by The leading coefficient C 3 (δ) is maximal at δ → 0 i.e. d → ∞. The asymptotic value is C 3 (0) 1.32178. In the real case, the average complexity is asymptotically bounded from below by C 4 (δ) d(p + s(d) + 3) with C 4 (0) 1.04074.
The theoretical predictions of this section have been confirmed, in practice: polynomials (53) are the slowest to evaluate (see Section 8.3 and, in particular, Figure 21). Figure 11: Comparison between a polynomial (53) whose coefficients obey the half-circle law (above) and a polynomial of the same degree whose coefficients obey a normal law (below). The complex plots (left) illustrate phase (in color), level lines (gray) and roots (white). The rate of cancelation of monomials is computed along the real axis (right) and illustrates that half-circle polynomials have extreme cancelations along the negative real axis, which extend far beyond the immediate vicinity of their roots.
Evaluation benchmarks with polynomials (53) whose coefficients obey the half-circle law lead to an interesting observation. As the degree increases, these polynomials appear to be extremely difficult to evaluate precisely along the real line. At degree 1000, only half of the computations with 600 bits along the real line are fully trustworthy; about 20% of the computations lead to at least 200 bits being identified as uncertain by our FPE algorithm. At degree 33 113, half of the 600-bit computations report that no bit is trustworthy. This exceptional situation piqued our interest because the roots of these polynomials appear to concentrate mostly along a sub-arc of the unit circle, which means that evaluations along the real line are usually not in the direct vicinity of a root.
A deeper analysis (see Figure 11) suggests that polynomials in this family have an extremely high cancelation rate of the monomials, which at a given z ∈ C is defined as the proportion of the monomials (a n z n ) n∈ 0,d such that |a n z n | > |P (z)|. Of course, along the positive half of the real axis, no cancelations can occur because all the coefficients are positive. In comparison, the cancelation rates for other families of polynomials seem to spike in much narrower regions of the complex plane. This observation is consistent with our statements on the complexity of the FPE algorithm and draws a parallel between slow FPE evaluations and precision loss.

Error analysis and proof of Theorem 3
In this section, we prove the correctness of the algorithm FPE p , i.e. Theorem 3. We adopt the notations from Section 4.3 and show that where c ≥ 0 is defined by (36). In the light of the property (21), it is enough to show instead that Let us assume first that |P (z)| ≥ |a k λ z k λ | where k λ = argmax Ù E P +λId , i.e. (roughly speaking) that there is no cancelation of leading bits. Using (31) and (6), one gets: with N λ defined by (45). Thanks respectively to the definitions (46) and (44), one has In both cases, we get s(a k ) i.e. for each dropout monomial. Using (31) a second time, we get any k ∈ D p (λ): We may now estimate R(z) = P (z) − Q λ (z) using #D p (λ) ≤ d < 2 s(d) or, equivalently, using (14): In this case, (54) holds with c = 0 (with a margin of 1 bit) and P (z) ≈ p Q λ (z). The lazy algorithm is therefore essentially exact when no leading bits get canceled.
In the case where |P (z)| < |a k λ z k λ |, some of the most significant bits cancel each other. More precisely, let us define c ∈ N by According to Remark 8, exactly c leading bits have been canceled while computing P (z). In this case, as we carry all computations with a fixed precision of p bits, only the first p−c bits of the result are meaningful (plus one implicit leader). One still has On the other hand, using (32), one has now and (54) holds in this case too.

Applications
In this section, we expose a few possible applications of the FPE algorithm, at both the theoretical and practical levels.

Parsimonious representation of polynomials
At a theoretical level, Theorem 3 states the existence of a parsimonious representation of any polynomial. This reduction can be computed algorithmically, is valid on any given annulus of C and guarantees a fixed arbitrary bound on the relative error.
For example, let us consider the Chebyshev polynomials T n (cos x) = cos(nx). They are the archetype of evaluations with extreme cancelations because each T n maps the interval [−1, 1] onto itself while the coefficients (a n,j ) 0≤j≤n of T n grow exponentially (namely max j |a n,j | 2 1.26n , as indicated by the maximum point of Figure 12, left). The scale profile of the coefficients of T n renormalized with (34), i.e. s(a n,j )/n appears to converge towards a fixed profile (red curve on Figure 12). Taking this fact for granted, Theorem 3 predicts the degree q(n, α) such that the reduced polynomial Q n,α (x) = T n (x) mod x q(n,α) provides an accurate approximation of T n (x) over the interval [−α, α]. For example, for α = 0.3 2 −1.74 , Figure 12 shows that the maximum of |a n,k 0.3 k | is achieved for k n 0.285n and that s(a n,kn ) 0.9n. For n = 200, max |a n,k 0.3 k | 2 200(0.9−1.73×0.285) 2 81 , which means that c 81 leading bits will be lost in the computation of T 200 (x) when x 0.3. Theorem 3 with p = 85 ensures that T 200 (0.3) can be computed with at least 3 significant bits if we keep the coefficients above the dashed line on Figure 12 (offset δ 96 bits), i.e. if we drop the last 25% of the coefficients. In general, this proportion is independent of n and we can claim that q(n, α)/n too is asymptotically independent of n. A direct proof of this result (without Theorem 3) does not seem obvious.  The Jacobi polynomials P (α,β) n and, in particular, the Legendre polynomials P n enjoy a similar property (Figure 13), which may be of interest for mathematical physics.
The engineering pressure towards better onboard electronics, using microcontrolers and field-programmable gate arrays, requires that some non-linear functions be computed quickly, often in reduced precision (e.g. 32, 16 and even 8 bits), with hardwarespecific optimizations. This problem has revived interest * in the Remez algorithm on the polynomial approximation of an arbitrary function that minimizes the L ∞ -error, i.e. minimax approximation [Rem34], [Hoc20]. For example, when dealing with periodic functions, engineers are interested in bypassing a costly reduction mod π/2 if a suitable interpolator provides accurate values on the natural range of angles for their problem.
For a given range of evaluation points, the algorithm FPE will either provide a further reduction of the number of coefficients needed at a given precision, or conversely, it will show that no further reduction is possible (see e.g. Figure 12). In both cases, such a result provides theoretical backing for the implementation choices. The -analyse task in our implementation [MV22] (see Section 8) provides a rudimentary tool to perform this analysis.
In practice, the level of parsimony achieved by the FPE algorithm can be remarkably high. For example, Figure 14 illustrates the proportion of monomials that are kept in Q λ (z) and therefore lead the value of P (z).

Application to root finding with Newton's method
On the practical side, Theorem 4 ensures the following two benefits.
Firstly, for a given allotment of computation time, one can perform k H evaluations of a certain polynomial with Hörner's method, or k FPE evaluations with the FPE algorithm. For a given precision p and large d, the asymptotic ratio can be extracted from (40), provided that the set of evaluation points is statistically diverse. One gets: The corresponding asymptotic gain factor is illustrated on Figure 1.
Secondly, using Remark 8, it is also very easy to detect cancelations of leading bits, which means that the FPE algorithm allows not only faster computations, but also provides a hint at runtime on the precision that should be used to achieve a certain level of accuracy (typically, the desired accuracy plus the number of canceled bits). Running error bounds (i.e. estimates of the absolute error committed during the evaluation process) are also available for Hörner [Hig02]; however they do not directly indicate the number of leading bits that where canceled, contrary to FPE, which can compare the scale of the largest monomial to the final result at no extra cost (see Figure 23). A typical application that takes advantage of these two properties is finding roots with Newton's method. Given a starting point z 0 ∈ C, one computes the sequence Almost surely, the sequence will converge towards a root of P ; divergence occurs when z 0 is in the Julia set of N P (which is of Hausdorff dimension < 2, see [Mil90], [CG93]). Costly excursions near ∞ occur also if the sequence visits a small neighborhood of a critical point. Using enough starting points (see the algorithm described in [HSS01]), one can compute all the roots of P . We refer the reader to our work [MV] for a refinement of [HSS01] that allowed us to split a tera-polynomial, i.e. deg P = 2 40 10 12 using a set of carefully chosen starting points for Newton's method. Here, we focus on the simpler task of showing the benefits of applying FPE to compute (58) instead of Hörner's scheme.
The first benefit is that the preconditioning of P and P can be done simultaneously. Indeed, if a j are the coefficients of P , then those of the derivatives satisfy: In a first approximation, the maps Ù E P and Ù E P are thus simply offset from one another by the concave map j → log 2 j. In practice however, Ù E P may have more segments than Ù E P . If speed is of the essence, one can choose to keep a low-resolution profile Ù E P and increase the safety margin δ (our implementation choice for the Newton demonstrator). Alternatively, one could perform a separate preprocessing for P and P .
The second and main key point is that the computation of P (z)/P (z) can be largely improved if one takes into account the cancelation of valuation * induced by FPE. Precisely, if one assumes that Q 1 (z) = k∈ ,r ∩Gp a j z j and Q 2 (z) = k∈ ,r ∩G p ja j z j−1 are the respective p-bit reductions of P (z) and P (z), then where m = min{ , }. Taking the simplification of z m into account improves both the speed and the accuracy. It is especially important in the early phase of Newton's sequence, where P (z) and P (z) may still be huge, which would cause a substantial loss of precision in the computation of the increment, or even an overflow. In [MV], we encounter examples where neither P (z) nor P (z) can be represented accurately with the precision chosen, but where P (z)/P (z) can be computed flawlessly. * The valuation of a polynomial is the lowest degree of its non-zero monomials, i.e. the multiplicity of zero as a root.
For example, with P (z) = z 64 + 1 and p = 24 bits, the evaluation of P (10) and P (10) with our FP32 implementation produces inf because of the obvious overflow. However, we can compute the correct Newton increment 10 − N P (10) 0.15625 with FP32 hardware arithmetic, which is actually an accurate value up to 2 × 10 −65 .
The third point in favor of the FPE algorithm occurs when the sequence z n eventually approaches a root of P , as it should; the computation of P (z n ) then leads to an increasing number of cancelations. Using Remark 8, we can easily issue a warning when it is time to switch the computations to a higher precision.
The last point is that FPE is embarrassingly parallel, which means that multiple roots can be searched for simultaneously on different cores using the method of [HSS01]. Also, contrary to more global algorithms that can be influenced negatively if some of the evaluation points lead to overflow values (e.g. if z n is a near miss of a root of P ), each computation with FPE is carried out independently of the others, even on a single core.
An example of root-finding using our implementation is illustrated in Figure 15.

Perspectives
Quadrature methods are at the heart of numerical analysis [BM92], [SSD04]. Using Gaussian quadrature, one may use n evaluations to compute exactly the integral of a polynomial of degree 2n − 1 over a given interval. The evaluation points (and the weights of the linear combination) are determined by orthogonal polynomials. The FPE algorithm can be used to speed up the evaluations without compromising precision in the case of high-degree polynomials (typically n 100).
Clenshaw's algorithm [Cle55] generalizes Hörner's method in order to evaluate recursively linear combinations of a polynomial basis, which is itself defined by a three-term recurrence relation. The principal of lazy addition at the heart of the FPE algorithm could be used in this general context to reduce finite precision computations to a parsimonious summation. The practical condition is the ability to compute easily the scale of the basis functions at a given point (like log 2 z k = k log 2 z).
Extending the FPE algorithm to the multivariate case would be a welcome generalization because the number of terms increases drastically. There are (d+n−1)! d!(n−1)! monomials of total degree d in n variables, i.e. O(d n−1 ) if d n. For example, a polynomial of degree 68 in 4 variables contains more than a million monomials, which is an instance of the well known curse of the dimension. For a recent study of the error estimates in the multivariate Hörner algorithm, we refer the reader to [PS00].
The key idea of the FPE algorithm (namely the lazy addition) is independent of the dimension. The transfer of the analysis of the dominant coefficients to an arbitrary evaluation point (x 1 , . . . , x n ) ∈ C n remains similar to the 1D case (32): log(|a k 1 ,...,kn x k 1 1 . . . x kn n |) = log |a k 1 ,...,kn | + n j=1 k j log |x j | .
The main question will be to estimate the average complexity of the FPE algorithm, which is essentially equivalent to the question of computing the average area of the horizontal projection of the largest hyperplane wafer that can be sandwiched between two copies, vertically offset by δ, of the graph of a concave function. A preliminary numerical exploration with a half-sphere function, i.e. f (x 1 , . . . , x n ) = 1 − |x| 2 , confirms that the area does scale as δ n/2 for n = 1, 2, 3 when δ → 0, which is encouraging.
Finally, let us mention that the FPE algorithm is of interest when evaluating polynomials or analytic functions on a disk. Remark 19 gives the appropriate weight to compute the average complexity when z is chosen at random uniformly on a disk. In Section 8.3, this case is benchmarked, along with the Riemann sphere and the real line.

Implementation and benchmarks
We have implemented our algorithm in the C language and we release the implementation as an open-source project [MV22]. Our implementation aims for the highest versatility and user-friendliness, without compromising performance. As a general rule, special cases that can lead to a substantial optimization are automatically recognized and dealt with.

General considerations
The main function, FastPolyEval, is called at the command line. Polynomials are specified as CSV files (passed as arguments) in which each coefficient, starting with a 0 , is written as a pair of its real and imaginary part in decimal form. For example, the polynomial P (z) = 2 + (3 − 5i)z is represented by the listing 2, 0 3, -5 Similarly, the set of evaluation points is specified as a CSV file that obeys the same format.
The first argument is systematically the precision at which the result of the operation is desired. If the requested precision is at most 24, 53 or 64 (depending on compile time options), FastPolyEval uses machine floating numbers, respectively FP32, FP64 or FP80. Otherwise, arbitrary-precision MPFR floating numbers [MPFR] are used. It is therefore possible to store the values of a polynomial with a high precision in a file and only use machine precision in a first set of low-precision evaluations, without worrying about a performance loss. On the contrary, if the precision requested exceeds that of the input, the input is considered exact (in decimal form) and padded with zero trailing bits if necessary.
FastPolyEval automatically identifies the case of real polynomials (all imaginary parts of coefficients are identically zero) because one can preprocess this case faster. Similarly, evaluations along the real line are also silently optimized by the evaluator.
Computing the scale of a real number is indeed about twice as fast as computing the scale of a complex number. In all cases, when using high-precision numbers, the scale, which is integer valued, is computed efficiently using only machine precision.
Our implementation of the FPE algorithm is complemented by a set of tasks that can generate polynomials (interpolation from a given set of roots, four common orthogonal families, the family of polynomials associated with the hyperbolic centers of the Mandelbrot set) and to manipulate them (sum, products, derivatives). Rescaling can be done by evaluating λz on the coefficients. We also provide a comprehensive set of tools to build and operate on sets of complex numbers. See Appendix C.
The tasks -eval, -evalD and -evalN can be used directly in production cases to evaluate a polynomial, its derivative or one Newton step with the FPE algorithm. An optional argument can be passed to generate a report on the number of bits that can reasonably be trusted in each evaluation, in accordance with Remark 8. Additional arguments enable the benchmark mode (timing, comparison with Hörner). One important optional argument is the errorsFile specification, that generate a complementary report on the estimated quality of the evaluation at the given precision (see Remarks 8 and 12). For each evaluation point, it contains an upper bound for the evaluation errors (in bits), a conservative estimate on the number of correct bits of the result, and the number of terms that where kept by the FPE algorithm.
The -iterN task is for the convenience of the user and provides a reasonably optimized stopping criterion for Newton's method. For best results, we recommend multiple runs, each with a limited number of iterations, and where the precision is gradually increased. The choice of the starting point and the pruning of duplicates is left to the end-user; see [HSS01] for guidance. For a complete implementation of a splitting algorithm, we refer the reader to [MV].
The -analyse task computes the concave cover Ù E P , the strip G p , and the intervals of |z| for which the evaluation strategy (i.e. the reduced polynomial Q λ (z)) changes. It is intended mostly for an illustrative purpose on low degrees, when the internals of the FPE algorithm can still be checked by hand. However, the intervals where a parsimonious representation is valid may also be of practical use; see Section 7.1.
The question of parallelization is a legitimate one if one wishes to get the most out of modern hardware. If the number of evaluation points is high compared to the core count, the algorithm FPE is embarrassingly parallel. Further optimization could be achieved by performing evaluations at points of similar size on the same core. To avoid an excessive complexity of the code that may only be of use in some specialized application, we chose to only implement a single-core version of FPE.

Implementation notes
The fact that E P , defined by (29), is discrete valued helps build a concave cover Ù E P with few segments (see Section 5.1), which in turn speeds up the binary searches for k λ , and r in the evaluation phase. Note that even if log 2 (|a k |) is concave, the scale function is integer valued, which, in practice, may prevent E P from being concave. See Figure 16.
In the course of sorting the values s(a k ), we could check whether the profile is concave and, if it is indeed concave, we could identify the maximum in an overall of 2d operations and reduce the complexity of the preconditioning to only O(d) operations. However, in general, it induces a loss in the evaluator (more segments in Ù E P ) and it is not worth the trouble. Similarly, using a non-integer scale would be ill advised. Figure 16: Even if log 2 |a k | is concave, the scale function E P is "pixelated" and not necessarily concave. It is a good thing because Ù E P contains fewer segments, which optimizes the evaluator.
At the end of each evaluation, one needs to compute the valuation monomial, i.e. z . The canonical method consists in writing = 2 α + β with β < 2 α . As z 2 α can be computed with α successive squaring that can be kept in memory to compute z β , the evaluation costs at most 2 log 2 multiplications.
In the preconditioning phase, we mark a set of indices B p ⊆ 0, d that never need to be computed for the given precision p (see below (44)). If their density is close to one in , r , which is the case for sparse polynomials, then, to evaluate z , we pre-compute other powers of z than the canonical z 2 j with j ≤ log 2 . For example, the gaps can be filled more efficiently and be dynamically optimized for the interval , r , with a negligible overhead. This remark is implemented in FastPolyEval, which ensures an equal treatment of all possible types of lacunarity, be it regular or not.

Benchmarks
We have tested the correctness and efficiency of our implementation on a few classes of polynomials that are either of large interest, hard to handle in general, or both.
The Chebychev, Legendre, Laguerre and Hermite polynomials are classic. The hyperbolic polynomials play a central role in the study of the Mandelbrot set and are defined recursively by p 1 (z) = z and p n+1 (z) = p n (z) 2 + z. Normal polynomials are of the form P ω (z) = a n (ω)z n where a n (ω) ∼ N (0, 1) are either real-or complex-valued random variables following normal law. The so-called half-circle family is defined by (53); the coefficients are real valued and form a half-circle when drawn in logarithmic coordinates. The complex version of the half-circle family is obtained by multiplying the coefficients of the previous family by a random phase uniformly distributed on the unit circle. As explained in Section 5.3, these polynomials are remarkably hard to evaluate accurately.
Systematic benchmarks were performed on Romeo, in the HPC center of the University of Reims. The overall benchmark time depends on the family of polynomials and, obviously, on the degrees and precisions used; in our case, it took a total CPU time of 70-90h per family. Multiple identical runs (typically two consecutive runs of the Hörner algorithm, and ten runs of the FPE algorithm) ensured that the average time is not biased by the loading time of a library or by fluctuations in the ambient load of the server. We performed complementary benchmarks on our personal computers to confirm the data points for the smaller degrees.
As explained in Section 5.2, we compute the average complexity when the evaluation points are chosen uniformly on either C, R or the unit disk. Before going further in our analysis, let us comment on the number of evaluation points used for the benchmarks.
We used 10 084 points uniformly distributed on the Riemann sphere (modulus ranging from 8 × 10 −3 to 2 × 10 2 ), 5 000 points uniformly distributed on the unit disk (modulus ranging from 4×10 −3 to 1) and 5 000 on the real line (ranging from ±6×10 −5 to ±2×10 3 ). Using more points does not change the average time significantly (see Figure 17), however it can dramatically and unnecessarily extend the CPU time. An order of magnitude of the computation time is given in Figure 18. On a modern laptop * , the evaluation of a polynomial of degree 1024 with a precision of 100 bits with Hörner's method takes in the ballpark of 143µs ±23% * * . With the FPE algorithm, the computation time depends significantly on the shape of E P . It can affect the preprocessing time negatively if Ù E P has many segments. Conversely, it can affect the evaluation time positively if Ù E P has few segments (ideally, with radically different slopes) or if G p contains a small number of terms (see Section 4.3 for the definition of G p ). The order of magnitude of the preprocessing step is 84µs ±32% and subsequent evaluations with FPE boil down to 46µs ±62%.
For a one time evaluation, the balance tilts slightly in favor of FPE, which is an interesting practical update on the optimality of the Hörner scheme. Note that our experiment does not contradict the theoretical result of Ostrowski [Ost54] and Pan [Pan66] because our advantage holds on average and only for computations with a fixed precision.
For a one time evaluation, FPE outperforms Hörner when M (p) log d. In the preprocessing phase, we only need to read the exponents of the coefficients, which remains a small amount of data to handle (O(d log d) with a small fixed constant). The evaluator then performs a minimal number of costly high-precision operations. Hörner on the other hand, has O(dM (p)) bit-operations to perform and may end up being slower. The advantage is especially pronounced in the complex case, where each numerical product costs 4 real multiplications.     If the polynomial is evaluated repeatedly (which is what FPE is designed for), the preprocessing overhead becomes negligible (see Figure 19) and the asymptotic gain obtained by FPE becomes substantial. Speedup in excess of ×10 occurs for some families like Hermite or Laguerre for degrees as low as 1000 (see Figure 18). In accordance with Theorem 4, if the degree is high enough, the speed gain is bounded from below by O( d/ log 2 d), which is observed in practice (see Figures 20 and 21).
In the best cases (Laguerre and Hermite; see Figure 21), the complexity of the FPE evaluation scales, in practice, as O(d 0.26 ) along the real line and O(d 0.28 ) on the Riemann sphere. In our data range, this complexity is consistent with O(log 2 d)±10%, which is the theoretical bound suggested by the last example of Section 4.4. Note that Figure 20 also hints that, in general, the exponent of the scaling law of the complexity does not depend on the precision used for the computations. Finally, let us point out that Chebychev, Legendre and Hermite polynomials are mildly lacunary (alternately odd or even); the others are not.
Sorting the polynomial families by increasing asymptotic gain as in Figure 18 and 21 is effectively a way of measuring the complexity in the variability of the scales of the coefficients. As explained in Section 5.3, the slowest case is that of the half-circle family. Similarly, the hyperbolic polynomials are slow to evaluate from their coefficients because of systematic compensations among monomials on the Mandelbrot set, which represents a substantial part of the Riemann sphere (about 29%). On the contrary, if Ù E P is composed of only a few segments, there will be very few values of |z| for which massive compensations among monomials can occur; in this case, the FPE algorithm produces a very parsimonious representation of the polynomial (see Section 7.1), which in turn is responsible for extreme speed gains.
Most polynomial families behave qualitatively the same on C, R and on the unit disk. The only substantial anomaly in this classification occurs with the normal family (both real and complex), which is asymptotically evaluated significantly faster on R than on C: for a polynomial of degree 30 000 and a precision of 53 bits MPFR, FPE evaluations are asymptotically 100 times faster than Hörner's on the real line but only 30 times faster on the Riemann sphere (see Figure 21). A reasonable explanation beyond the fact that real powers are easier to compute than complex ones, is the fact that the roots of normal polynomials accumulate uniformly along the unit circle (Hammersley's theorem [Ham56], [SZ03]); therefore, one may expect fewer cancelations along the real line than for other families. However, the anisotropic example at the end of Section 5.3 suggests caution and further studies would be required to confirm this explanation. In particular, the reason why the evaluation time of the normal family on the Riemann sphere fails to obey a power law contrary to all other families is not clear.
Our implementation [MV22] handles both hardware number formats FP32, FP64 or FP80 and arbitrary-precision MPFR floating point numbers [MPFR]. The main limitation of hardware formats is the short range of exponents: roughly speaking, one can only represent numbers whose absolute value lies between 10 ±38 with FP32 numbers and between 10 ±308 with FP64 numbers. Concretely, this means it is simply impossible to compute a monomial z n in FP64 when |z| ≥ 2 and n > 1024. The FP80 format provides a slightly more comfortable range between 10 −4 951 and 10 4 932 , but it is still not enough to handle polynomials of degree 50 000 as in our benchmark.
In the range of exponents where a comparison was possible, hardware numbers behave about 4 times faster than 53-bit MPFR numbers; however, the gain factor of FPE over Hörner obeys the same scaling law as for MPFR. In practice, the sweet spot for using the FPE algorithm with FP80 numbers is for polynomials of degree 1 000 to 5 000 and |z| < 10. When using MPFR numbers, this range is extended to essentially any degree above 100 with almost no practical limitation on |z|.
The last crucial part in our benchmarks is the question of the accuracy of the FPE algorithm, which is guaranteed by Theorem 3. To put it to the test, we systematically computed a 600-bit evaluation of our polynomials with a Hörner scheme, which served as a reference value. For each benchmarked precision (up to 304 bits), the outputs of both Hörner and FPE algorithms with the current precision were compared to the reference value to identify the absolute computation error.
The most significant data that can be extracted from this computation is the accuracy bias, defined as the difference of the number of exact bits between the two algorithms, which is presented in Figure 22. The practical conclusion is that the values computed either by FPE or by Hörner are essentially identical, up to 1 exact bit. Of course, when cancelations occur, the displayed result may differ by many bits, but the divergence only affects the non significant bits at the end. Based on this extensive benchmark, we can now confirm, in practice, that the FPE algorithm holds the promise of Theorems 3 and 4 and performs as accurately as Hörner, only faster.
Let us conclude this section by pointing out that our implementation of FPE provides additional tools for analyzing polynomial evaluations like the proportion of leading monomials at a given evaluation point (see Figure 14). Similarly, the localization of cancelations in the evaluation process can easily be deduced from the output files (see Figure 23), which may guide practical decisions to ensure the precision of subsequent computations. We hope that the ideas presented in this article will inspire future developments, either theoretical or applied. We also thank the reader for reaching this point.

A Proof of the geometric statements
In this appendix, we prove the geometric results stated in Section 3.
The proof of Theorem 1 is based on two complementary geometric constructions and is split into several lemmas. In Lemma 20 we show the existence of a segment [A(τ )B(τ )] in S(f, δ) of maximal length and of slope τ = tan θ. This segment touches the graph of f − δ in a "tangent" way (in the convex sense, i.e. as a subderivative). In Lemma 21 we express the main integral from Theorem 1 in terms of x R − x L where x L and x R are the respective abscissae of A(τ ) and B(τ ). Next, we make an alternative geometric construction of x R and x L based on area computations. This second construction is the key to Lemma 22 where we estimate the diagonal of a square built upon the graph of f . A change of variable in the plane (Lemma 23) allows us to collect all the prior estimates and leads to the proof of the upper bound in Theorems 1 and 2. Finally, we obtain the lower bounds by constructing explicit examples.

A.1 First geometric construction based on the graph of f
Let us start by presenting the construction in a simple case.
A typical example. For now, we suppose f ∈ C 2 ([0, 1]) with f (x) < 0 for all x ∈ (0, 1) and we single out x 0 ∈ (0, 1). The line of slope f (x 0 ) and passing through the point (x 0 , f (x 0 ) − δ) is tangent to the graph of f − δ. This line intersects the graph of f in at most two points A and B, one on each side of x 0 . Indeed, the points of intersection with the graph of f correspond to the solutions of the equation The function F is of class C 2 ([0, 1]) and F (x) = f (x) < 0 for all x ∈ (0, 1), thus F is concave. Therefore Consequently, there exist at most two points, one on each side of x 0 , such that F (x) = −δ. When they exist, we denote them by 0 ≤ x L < x 0 < x R ≤ 1. When they do not, we simply take respectively x L = 0 or x R = 1. In Figure 6, the common abscissa of the points A, A is x L while that of B, B is x R .
General case. Let us now come back to the general setting where f is concave, but not necessarily of class C 2 ([0, 1]). The next statement extends the simpler case presented in the previous paragraph. It is essentially an elementary version of F. Riesz's rising sun lemma [Rie32] in a concave setting.
Lemma 20. For any real number τ ∈ R, there exists a unique segment [A(τ )B(τ )] of maximal length, of slope τ = tan θ, contained in the strip S(f, δ) and that touches the graph of f − δ in a tangent way in the convex sense.
We denote by g(x 0 ± 0) or g(x ± 0 ) the sided limits of a function g: A segment [AB] of slope tan θ is said to be tangent in the convex sense to the graph of g ∈ C if, at any contact point (x 0 , g(x 0 )) ∈ [AB], one has g (x 0 + 0) ≤ tan θ ≤ g (x 0 − 0). If g is smooth, then g (x 0 ) = tan θ. If an endpoint x 0 ∈ {0, 1} is a contact point, then the requirement is lightened respectively to tan θ ≥ g (0 + ) or tan θ ≤ g (1 − ).
The abscissa of the endpoints of the maximal segment [A(τ )B(τ )] will be respectively denoted by x L (τ ) for the left side and x R (τ ) for the right side. Figure 24: The affine map (x, y) → (x, y − τ x) rearranges Figure 6 in a so-called "rising sun" configuration [Rie32]. Recall that τ = tan θ. Two non-optimal segments illustrate the last part of the proof: the longest segment is the one "tangent" to the graph of f (x)−τ x−δ.
Proof. As f is a concave function, it is continuous, it is differentiable almost everywhere, and its derivative is decreasing. Even when the derivative is not continuous at a point, it necessarily has left and right limits. Therefore it has at most a countable number of jump points. Let us first construct the segment [A(τ )B(τ )]. For any τ ∈ R, the concave function f (x) − τ x presents a maximum m τ in [0, 1], which is reached on some non-empty compact sub-interval J τ ⊂ [0, 1] (usually a singleton). The function f (x) − τ x is monotone on each of the connected components of [0, 1]\J τ . As a consequence, the set is an interval increasing in δ that contains J τ = I(τ, 0). Let us define x L (τ ) = inf I(τ, δ) , x R (τ ) = sup I(τ, δ) .
For any x 0 ∈ J τ , one has m τ = f (x 0 ) − τ x 0 and the line of equation y = y τ (x) with does not depend on the actual choice of x 0 within J τ . Let us define By definition (60), for all x ∈ [x L (τ ), x R (τ )] = I(τ, δ), one has f (x)+δ ≥ m τ +τ x ≥ f (x), thus the segment [A(τ )B(τ )] is of slope τ and is included in S(f, δ).
Conversely, any segment [AB] of slope τ included in S(f, δ) is supported by a line of equation y = τ x + β and must satisfy (denoting by x A , x B the abscissa of A and B) Similarly, one has x B ≤ x R (τ ). In other words, one has If β < m τ − δ, the constraint (62) cannot be satisfied for x ∈ J τ because The concavity of f (x) − τ x implies that one can increase x B − x A by shifting the interval towards J τ . More precisely, let us assume for example that [x A , x B ] is on the right side of J τ and that x R (τ ) < 1 (otherwise nothing needs to be proved). The inclusion [AB] ⊂ S(f, δ) implies The function τ − f is defined almost everywhere and is positive and increasing on the right-hand side of J τ . The inequality (63) can thus be rephrased Similarly, for any x 0 ∈ J τ , one has f (x 0 ) − τ x 0 = m τ and x R (τ ) < 1 implies : This proves that the maximal length is obtained for β = m τ − δ. As the values of the projection of A, B on the x-axis are unique, the equations (60)-(61) ensure that the segment [A(τ )B(τ )] is unique.
The quantity that interests us for Theorem 1 is obviously related to this first geometric construction.
Lemma 21. With the notations of Theorem 1, we have We also have and for any positive measurable weight ω on − π 2 , π 2 : Proof. The definition (61) ensures that the length of the projection on the x-axis of the segment [A(τ )B(τ )] is L(f, δ, θ) cos θ = x R (tan θ) − x L (tan θ). The identities are then obtained by the change of variable y = tan θ.

A.2 A second geometric construction based on the graph of f
We are now going to provide a second geometric construction of x L (τ ), x R (τ ). We consider the graph of f and we complete it to a continuous curve in the following way. At a jump point we add the vertical segment that joins the left and the right limits.
We thus obtain a continuous curve, which we denote by γ f , that contains the graph of f , whose projection on the x-axis contains ]0, 1[ and is included in [0, 1] and whose projection on the y-axis is R.

Let us introduce
The sets Γ ± f are the closed subsets of the strip [0, 1] × R that are respectively above and below γ f . For every y 0 ∈ R we consider is the intersection between γ f and the horizontal line y = y 0 ; it reduces to a point when x + 0 = x − 0 . We now build a family of "triangles" whose hypotenuse rests on γ f and that collapse on the segment [x − 0 , x + 0 ] × {y 0 } (see Figure 25). For x ∈ [0, 1], let The continuity of γ f implies that the area |T f (y 0 ; x)| of this triangle is a continuous function of x and the monotonicity of f implies that the area vanishes along [x − 0 , x + 0 ] and is respectively strictly decreasing on [0, x − 0 ] and strictly increasing on [x + 0 , 1]. We are interested in the two points where either the area of the triangle equals δ or the triangle hits the edge of the strip: Let us prove that x L (y 0 ) = x L (y 0 ) and x R (y 0 ) = x R (y 0 ), i.e. they are the same values as the ones defined by (60). Using elementary calculus, we know that the area of the triangle is The conditions defining x L (y 0 ) and x R (y 0 ) thus boil down to S(f, δ) and x → y 0 x + m y0 − δ Figure 25: Images of the curves γ f in black, γ R in red and γ L in blue. The triangles T f (y 0 ; x L ) and T f (y 0 ; x R ) are grayed out; their area is δ. Inset: the profiles of f and f − δ used to generate the figure; the dashed segment [A(y 0 )B(y 0 )] has a slope τ = y 0 . Note that, contrary to the profile presented in Figure 6, this one has an unbounded slope near zero.
We now have two equivalent constructions of x L (y) = x L (y) and x R (y) = x R (y) for every y ∈ R. The second construction ensures that the maps y → x R (y) and y → x L (y) are continuous and decreasing on R. They are strictly decreasing respectively when x L (y) > 0 and x R (y) < 1. We denote by γ L (δ) = {(x L (y), y) ; y ∈ R} and γ R (δ) = {(x R (y), y) ; y ∈ R} the two curves that are "offset" from γ f by a triangular area of δ (see Figure 25).
Lemma 22. For every y ∈ R let us construct the unique square with an upper-right corner at (x R (y), y) ∈ γ R (δ) and a lower-left corner on γ L (δ). The area of this square is smaller than 2δ and therefore its diagonal is smaller than 2 √ δ.
Proof. The construction of the square is obvious. The curve γ L is below the curve γ f , which is itself below γ R . As y → x L (y) is decreasing, the curve γ L intersects a line of slope π/4 passing through (x R (y), y) in a unique point whose coordinates are, by definition, of the form (x L (y ), y ) for some y < y. The two points (x R (y), y) and (x L (y ), y ) are the opposite corners of a square, which we will denote by Q(y , y) in the rest of this proof; see Figure 26 (left). Notice that this square is always included in the strip between the curves γ L and γ R .
γ L γ f γ R γ L + 2 √ δ #» e Q(y , y) (x R (y), y) (x L (y ), y ) Q(y , y) Figure 26: Left: An example of squares Q(y , y) with opposite corners (x L (y ), y ) ∈ γ L and (x R (y), y) ∈ γ R . According to Lemma 22, their diagonal is bounded by 2 √ δ: the curve γ R is below the offset γ L + 2 √ δ #» e (dashed blue) where #» e = (1, 1)/ √ 2 is the unit vector along the diagonal. In this example, note how tight the estimate is near the point where f (x) −1. Right: The four possible configurations corresponding to how γ f can enter or exit Q(y , y). The area of Q(y , y) complemented by the two highlighted triangles is, by construction, exactly 2δ.
For θ ∈ [−θ 0 , θ 0 ] the longest segment [AB] ⊂ S(f, δ) of angle θ has both of its ends on the graph of f and is tangent to the graph of f − δ. Using the symmetry of the graph, one can assume that θ ≥ 0. Let (x θ , f (x θ ) − δ) denote the point where [AB] is tangent to the graph of f − δ. As tan θ = f (x θ ) = , the equation of [AB] gives 4x 2 θ − 4x θ + cos 2 θ = 0.
The functions C 1 and C 2 are continuous and strictly decreasing, with where K and Γ are classical special functions (respectively the complete elliptic integral of the first kind and the Gamma function).
See Section 5.3 for an adaptation of this example to polynomials that saturate the upper-bound on the complexity of the FPE algorithm, both theoretically and in practice.

Example 3 :
We study f (x) = −a(x − 1 2 ) 2 for x ∈ [0, 1] with a > 0 and with δ < a/4. Notice the symmetry with respect to the line x = 1/2 and that the maximum of f − δ is −δ and is superior to f (0) = −a/4. We denote by θ 0 ∈ (0, π 2 ) the angle between the x-axis and the line tangent to the graph of f − δ that passes through (0, −a/4) and by      G p , B p : list of a-priori good (resp. ignored) coefficients for a given precision p.
, r : left/right edges to further reduce G p for a given λ. Q λ (z) : reduced polynomial produced by the FPE algorithm.
FPE p : new algorithm proposed in this article, for computations with a fixed precision p. avg C : average operator for z uniformly distributed over C. avg R : average operator for z uniformly distributed over R. avg D(0,1) : average operator for z uniformly distributed over the unit complex disk.
C Listing of tasks implemented in [MV22] In our implementation [MV22], the tasks listed in this section are called in the command line with FastPolyEval -task [arguments]. The first argument is systematically the precision of the computation, in bits. Use -task -help for more detailed informations.
Tools for generating and handling polynomials -sum computes the sum of two polynomials and writes the result to a CSV file -diff computes the difference of two polynomials -prod computes the product of two polynomials -der computes the derivative of a polynomial -roots computes the polynomial with a given list of roots -Chebyshev writes the coefficients of the Chebyshev polynomial -Legendre writes the coefficients of the Legendre polynomial -Hermite writes the coefficients of the Hermite polynomial -Laguerre writes the coefficients of the Laguerre polynomial -hyperbolic writes the coefficients of the hyperbolic polynomial Tools for generating and handling sets of complex numbers -cat concatenates two CSV files containing complex numbers -re writes the real part of the list of complex numbers -im writes the imaginary part of the list of complex numbers -conj writes the conjugates of the list of complex numbers -join joins the real part of two sequences into one sequence of complex numbers -tensor computes the tensorial product of the two lists of numbers (c i = a i × b i ) -grid computes the set product of the real parts of two sequences -exp computes the complex exponential of a list of points -rot maps complex numbers (a, b) to a * exp(ib) -unif writes real numbers in arithmetic progression -rand writes real random numbers uniformly distributed in an interval -normal writes real random numbers with Gaussian distribution -sphere writes polar coordinates approximating a uniform distribution on the sphere -polar computes the points given by polar coordinates on the sphere -comp compares two lists of points Fast Polynomial Evaluator algorithm for production use and benchmarking -eval quickly evaluates a polynomial on a set of points -evalD quickly evaluates the derivative of a polynomial on a set of points -evalN quickly evaluates one Newton step of a polynomial on a set of points -iterN quickly iterates the Newton method (partial search of roots of the polynomial) -analyse computes the concave cover and the intervals of |z| for which the evaluation strategy changes