How to split a tera-polynomial

This article raises the state of the art of polynomial splitting to a degree about one trillion (i.e. $2^{40}\sim 10^{12}$), called a tera-polynomial²²2As a trinomial denotes a polynomial expression with 3 terms, we should call it a teranomial. We decided against it as the single letter difference with the prefix tetra (that denotes 4) could lead to confusion.. We do not claim full generality and focus on the families $p_{n}$ and $q_{\ell,n}$ defined above.

To achieve this result without brute force, the main idea is a new splitting algorithm. Roughly speaking, the classical Newton’s method consists in choosing [HSS01] a mesh of initial points that encircles the roots and to iterate the Newton map until the trajectories reach the roots; this method results in a dense mesh of most of the disk (see Figure 2).

We are releasing an exhaustive and proven list [Mand.DB] of all the roots of $p_{n}$ (called hyperbolic centers) for $n$ up to period $41$ and all the roots of $q_{\ell,n}$ (called pre-periodic or Misiurewicz-Thurston points) of order $\ell+n\leq 35$. By proven, we mean that each point $z_{j}$ in our database comes with a mathematical statement that guaranties two things: an exact root $\tilde{z}_{j}$ lies within a certain maximal tolerance from our numerical value, and a standard refinement technique (Newton’s method) applied to $z_{j}$ converges to $\tilde{z}_{j}$. To ensure a strict control of the computational errors that is robust enough to withstand intensive iterations, we have revisited and extended the theory of disk arithmetic.

As a companion to this article, we also release an implementation of our algorithm as a standard C library [Mandel]. This library can be used to check and explore our database [Mand.DB] and contains all the tools necessary to perform extractions, generate images, compute arbitrary precision refinements, etc. This library is optimized on all fronts: it can run on a single core for the smallest periods or can take advantage of the massively parallel architectures of modern HPC centers for the higher ones. For example, the splitting of $p_{21}$, which is reported as requiring at least 18.8 hours-core in [SS17] takes only 10 seconds with our algorithm, on a single core of a consumer-grade personal computer. Similarly, we can split $p_{29}$ in only 1 hour-core and $p_{33}$ in about 1 day-single-core. The computations automatically switch to high-precision, using the [MPFR] library, but only when necessary. For example, when $n\geq 30$ the roots of $p_{n}$ near $-2$ are so close that they cannot be resolved in the standard 64 bits floating point arithmetic. Up to $n\leq 41$, we use numbers with 128 significant bits. Throughout the library, custom data structures ensure near optimal memory usage, data access and computation times.

The case of the polynomials $q_{\ell,n}$ is interesting because, contrary to the $p_{n}$, they have multiple roots, which slow down drastically the dynamics of the Newton map. However, our algorithm remains robust in this case. We stopped at the order $\ell+n\leq 35$. The volume of data generated by the roots of all $(\ell,n)$ combinations is of the same order of magnitude to that of the roots of all $p_{n}$ for $n\leq 39$ (see Section 2.4).

Our database and its companion library constitute a powerful numerical microscope whose resolution extends far beyond the current state of the art and that can assist the mathematical community in the exploration of the properties of $\mathcal{M}$. While an in-depth analysis will take time, we already present in this first article a few properties, like the minimal distance between the roots of a given type (see Section 6.3.2).

Our article is organized as follows. Section §2 contains a brief self-contained introduction to the Mandelbrot set. Section §3 presents a quick review of the most common splitting (or at least root finding) algorithms for polynomials. A detailed presentation of our new algorithm is done in Section §4, along with a general discussion of its algorithmic and bit complexity. As we do not claim full generality, the specifics of working with polynomials tailored to the Mandelbrot set is discussed too. Readers interested in polynomials of high degrees may appreciate, on the side, our improvement of the standard evaluation scheme [AMV22], along with its companion library [FPE].

Section §5 is focused on the question of providing certifiable results. The backbone of this section is a theory of a disk arithmetic, which, contrary to interval arithmetic, tolerates a high number of iterations.

Our implementation [Mandel] is presented in Section §6. We discuss the practical cost of splitting the tera-polynomial $p_{41}$ and confirm, in practice, the complexity claims of Section §4. A listing of the tasks is given in Appendix D.

A proof of our Theorem 2 stated in Section §2 on the factorization of $q_{\ell,n}$ is given in Appendix A. A few background auxiliary mathematical results used in section §5 are recalled very briefly in Appendix B. For the sake of convenience, the notations introduced throughout the article are listed in Appendix C.

We are grateful to Romeo, the HPC center of the University of Reims Champagne-Ardenne, for hosting our project graciously even though its scale (near 50TB of data and 1 century-core / 870 kh of CPU usage for both the development and production phase) is quite substantial for this regional structure.

It is our pleasure to thank our families who supported us in the long (5+ years) and, at times stressful, process that gave birth to the Mandelbrot library. We also thank Sarah Novak for a thorough proof-reading of this article.

For $n\geq 1$, the roots of $p_{n}$ are parameters $c\in\mathcal{M}$, called hyperbolic centers, whose critical orbits are periodic of period $n$. The roots of $p_{n}$ are simple [DH82, Buf18] and the polynomial $p_{n}(z)$ is divisible by $p_{k}(z)$ for any divisor $k$ of $n$. A. Douady and J.H. Hubbard [DH82, DH85] have shown that the hyperbolic centers are interior points of $\mathcal{M}$, with at most one hyperbolic center per connected component of the interior. The set of all hyperbolic centers is also dense in the boundary of $\mathcal{M}$ in the sense that its closure in $\mathbb{C}$ contains $\partial\mathcal{M}$.

Let us define the set of hyperbolic centers of order $n\geq 1$ as the subset of $p_{n}^{-1}(0)$ whose minimal (or fundamental) period is exactly $n$; in other terms:

$$\operatorname{Hyp}(n)=\left\{z\in p_{n}^{-1}(0)\,\big{|}\,\forall k\in\operatorname{Div}(n)^{\ast},\>p_{k}(z)\neq 0\right\}$$

(4)

where $\operatorname{Div}(n)=\{k\in\mathbb{N}^{\ast}\,;\,\exists k^{\prime}\in\mathbb{N},\enspace\enspace kk^{\prime}=n\}$ and $\operatorname{Div}(n)^{\ast}=\operatorname{Div}(n)\backslash\{n\}$ is the set of strict divisors of $n$. For example, $\operatorname{Hyp}(1)=\{0\}$ and $\operatorname{Hyp}(2)=\{-1\}$. The reduced hyperbolic polynomial, also known as Gleason’s polynomial, is:

$$h_{n}(z)=\prod\limits_{r\in\operatorname{Hyp}(n)}(z-r).$$

(5)

The polynomials $p_{n}$ and $h_{n}$ have integer coefficients. While expected, the irreducibility of $h_{n}$ over $\mathbb{Z}[z]$ remains conjectural (see [HT15, last remark of §3], [SS17, p.155]).

For integers $n\geq 1$ and $\ell\geq 2$, the roots of $q_{\ell,n}$ are parameters $c\in\partial\mathcal{M}$, called pre-periodic centers or Misiurewicz-Thurston points, whose critical orbits are pre-periodic, that is it becomes periodic of period $n$ after the first $\ell$ steps. For $\ell\in\{0,1\}$, one can check immediately that

$$q_{0,n}(z)=p_{n}(z)\qquad\text{and}\qquad q_{1,n}(z)=p_{n}^{2}(z).$$

(8)

The dynamical reason for the second identity in (8) is that $0$ is the only pre-image of $c$ under $f_{c}$ so a pre-periodic orbit of type $(1,n)$ starting at zero actually loops back to zero, so is periodic of period $n$. The polynomial $q_{\ell,n}(z)$ is divisible by $q_{\ell,k}(z)$ for any divisor $k$ of $n$ and by $q_{\ell^{\prime},n}(z)$ for any $\ell^{\prime}\leq\ell$. Some of the roots of $q_{\ell,n}(z)$ have multiplicity; however, the polynomial

$$s_{\ell,n}(z)=\frac{q_{\ell,n}(z)}{q_{\ell-1,n}(z)}\in\mathbb{Z}[z]$$

(9)

has simple roots [Buf18], [HT15, Lemma 3.1]. Douady-Hubbard [DH82, DH85] have shown that the set of all pre-periodic points is dense in the boundary of $\mathcal{M}$. Visually, those points are either branch tips, centers of spirals or points where branches meet (see Figure 3).

By analogy with the hyperbolic case, let us define the set of Misiurewicz points of type $(\ell,n)$ as the subset of $q_{\ell,n}^{-1}(0)$ whose dynamical parameters are exactly $\ell$ and $n$; in other terms:

$$\operatorname{Mis}(\ell,n)=\left\{z\in q_{\ell,n}^{-1}(0)\left|\begin{array}[]{l}q_{\ell-1,n}(z)\neq 0,\\[3.0pt] \forall k\in\operatorname{Div}(n)^{\ast},\>q_{\ell,k}(z)\neq 0\end{array}\right.\right\}.$$

(10)

We call $\ell+n$ the order of $\operatorname{Mis}(\ell,n)$ because $q_{\ell,n}$ is a polynomial of degree $2^{\ell+n-1}$. The reduced polynomial whose roots are exactly $\operatorname{Mis}(\ell,n)$ is denoted by

$$m_{\ell,n}(z)=\prod\limits_{r\in\operatorname{Mis}(\ell,n)}(z-r)\in\mathbb{Z}[z].$$

(11)

Note that $\operatorname{Mis}(0,n)=\operatorname{Hyp}(n)$ and $\operatorname{Mis}(1,n)=\emptyset$ because of (8). The count of the Misiurewicz points has been established by B. Hutz and A. Towsley [HT15, Cor. 3.3]:

$$\left|\operatorname{Mis}(\ell,n)\right|=\Phi(\ell,n)\left|\operatorname{Hyp}(n)\right|$$

(12)

where

$$\Phi(\ell,n)=\begin{cases}1&\text{if }\ell=0,\\ 2^{\ell-1}-1&\text{if }\ell\neq 0\text{ and }n|\ell-1,\\ 2^{\ell-1}&\text{otherwise.}\end{cases}$$

To get a complete factorization of $q_{\ell,n}(z)$ in terms of (5) and (11), one needs to understand the multiplicity of the hyperbolic factors. We claim the following result.

The sets of all hyperbolic-centers and of all Misiurewicz-Thurston points are respectively denoted by:

$$\operatorname{Hyp}=\bigcup_{n\in\mathbb{N}^{\ast}}\operatorname{Hyp}(n)\subset\overset{\circ}{\mathcal{M}}\qquad\text{and}\qquad\operatorname{Mis}=\bigcup_{\begin{subarray}{c}\ell\geq 2\\[1.0pt] n\in\mathbb{N}^{\ast}\end{subarray}}\operatorname{Mis}(\ell,n)\subset\partial\mathcal{M}.$$

As mentioned above, $\overline{\operatorname{Mis}}=\partial\mathcal{M}\subset\overline{\operatorname{Hyp}}$. The density measures of $\operatorname{Hyp}$ and $\operatorname{Mis}$ (i.e. the limit of the uniform counting measure on finite subsets of increasing order) both coincide with the harmonic measure for $\partial(\mathcal{M}^{c})$, which is the image of the uniform measure on the unit circle under the Riemann map $D(0,1)^{c}\to\mathcal{M}^{c}$. Intuitively, the harmonic measure expresses the asymptotic density of the external rays as they land on the boundary of the Mandelbrot set. The support of the harmonic measure is the only subset of $\partial\mathcal{M}$ that is “visible” to a random brownian motion starting at $\infty$ and is of Hausdorff dimension 1. For further details, see [Mak85], [Lev90], [BS05], [FG15], [GV17], [GV19] and Remark 7 below.

Using the library [Mandel] associated with this article, we have completely determined the sets $\operatorname{Hyp}(n)$ for all $n\leq 41$ and $\operatorname{Mis}(\ell,n)$ for $\ell+n\leq 35$. Both problems consists in solving polynomial equations, namely (1) or (2), of respective degrees $2^{n-1}$ or $2^{\ell+n-1}$. Overall, we have sifted through about $3.3\times 10^{12}$ polynomial roots and identified $2\,725\,023\,424\,662$ distinct parameters of the Mandelbrot set. About 80% of them are hyperbolic centers.

The tera-polynomial $p_{41}$ has degree $1\,099\,511\,627\,776$ and, as $41$ is a prime number,

$$p_{41}^{-1}(0)=\operatorname{Hyp}(41)\cup\{0\}.$$

More generally, the deflation of $\operatorname{Hyp}(n)$ due to strict divisors of $n$ constitutes an asymptotically negligible subset of the roots of $p_{n}$.

The situation is not exactly the same for $\operatorname{Mis}(\ell,n)$. There are $33$ non-hyperbolic pre-periodic polynomials $q_{\ell,n}$ of order $\ell+n=35$, each of degree $17\,179\,869\,184$. However, the exact count $271\,590\,481\,057$ of $\bigcup\operatorname{Mis}(\ell,n)$ for all types of order 35 represents only $47.9\%$ of the total number of roots of those polynomials. More generally, (12) ensures that, for large $N$:

$$\sum_{\begin{subarray}{c}\ell+n=N\\ \ell\neq 0\end{subarray}}|\operatorname{Mis}(\ell,n)|\simeq\frac{N-2}{2}\>|\operatorname{Hyp}(N)|.$$

The deflation of $\operatorname{Mis}(\ell,n)$ due to divisors is therefore, asymptotically, 50% as $\ell+n\to\infty$.

To convey a sense of the orders of magnitude involved for storing and processing the roots of polynomials of giga- and tera-scale, let us underline that a computation that takes, on average, only 1 milli-second per root requires about $12$ days-core for a polynomial of degree $10^{9}$. This requirement jumps to $32$ years-core ($\sim 280$k hours-core) for a polynomial of degree $10^{12}$.

Our actual computation times are discussed in Section §6. The success of our entreprise is due to the efficiency of our new splitting algorithm for $p_{n}$ and $q_{\ell,n}$ (see Section 4), which is bounded by $O(d\log d)$ and runs as $O(d)$ in practice.

A naive splitting algorithm consists in finding a root of $P$ by any available method, and then reduce the degree of $P$ by performing a euclidian division.

In practice, this method is plagued by numerical instabilities, illustrated on Wilkinson’s polynomial [Wil84]

$$W_{20}(x)=\prod_{n=1}^{20}(x-n)$$

whose largest root is $20^{19}$ times less stable than its smallest root when the coefficients are perturbed in the canonical basis (though other basis may behave better). In general, determining the roots of a polynomial from the list of its coefficients is a highly ill-conditioned problem.

Consequently, for very large polynomials, computing the coefficients is not always advisable. For example, as $p_{3}(z)=z^{4}+2z^{3}+z^{2}+z$, the largest coefficient of $p_{n}(z)$ for $n\geq 3$ exceeds $2^{2^{n-3}}$ while the smallest one remains $1$. Storing the coefficients of $p_{41}(z)$ exactly is, at best, impractical and any reasonable approximation of those coefficients would not be sufficiently precise to compute the roots accurately.

Bracketing (or root isolation) consists in identifying an interval of the real line or a domain of the complex plane where the number of roots of $P$ is prescribed.

Methods based on this principle are usually developed on the real line where they constitute a textbook application of the intermediary value theorem (bisection method, or the more involved ITP method).

In the complex plane, the splitting circle method and the Lehmer–Schur algorithm are based on the residue theorem, which implies that, for all $n\in\mathbb{N}$

$$\sum_{\begin{subarray}{c}z\in\Omega\\[2.0pt] P(z)=0\end{subarray}}m_{P}(z)z^{n}=\frac{1}{2i\pi}\int_{\partial\Omega}\frac{P^{\prime}(z)}{P(z)}z^{n}dz$$

for any smooth domain $\Omega$ whose boundary avoids the zeros of $P$ and where $m_{P}(z)$ denotes the multiplicity of $z$ as a root of $P$. Successive refinements of $\Omega$ will either isolate single roots or clusters of close roots that can be identified using Newton’s identities. Aside from the cost of computing the integral to a sufficient precision, those methods are plagued by the geometrical problem of finding a “good” circle, i.e. one which is proper to initiate a sequence of divide and conquer iterations.

Of similar nature is Graeffe’s method, where the sequence

$$P_{n+1}(z^{2})=(-1)^{d}P_{n}(z)P_{n}(-z)$$

is computed from $P_{1}=P$. The roots of $P_{n}$ are $\{z^{2^{n}}\,;\,P(z)=0\}$ and are thus, generically, exponentially separated from each other. In that case Vieta’s formulas can easily be inverted approximately, which gives the roots of $P_{n}$ and ultimately, those of $P$.

Root isolation methods excel in finding roots of a given polynomial, but they are not a splitting algorithms per se. We will not insist on comparing their complexity.

Our algorithm (see §4 below) presents a strong familiarity with the splitting circle method: the level lines that we will use provide a tight enclosure of the roots of $P$ by a Jordan curve. Tighter level sets would enclose the roots in a finite union of disjointed connected sets (see Figures 2 and 6).

Splitting a polynomial $P$ can always be restated as the question of finding all the eigenvalues of the companion matrix $A_{P}\in\mathcal{M}_{d}(\mathbb{C})$, at least provided that the coefficients of $P$ are known.

$$A_{P}=\begin{pmatrix}0&1&0&\ldots&0\\ 0&0&1&\ddots&\vdots\\ \vdots&&\ddots&\ddots&0\\ 0&\cdots&\cdots&0&1\\ c_{0}&c_{1}&\cdots&\cdots&c_{n-1}\end{pmatrix}\qquad\text{if}\qquad P(x)=x^{n}-\sum_{k=0}^{n-1}c_{k}x^{k}.$$

The vector $(1,x,x^{2},\ldots,x^{n-1})$ is an eigenvector of $A_{P}$ if and only if $P(x)=0$.

The complexity of eigenvalue algorithms usually depends on the structure of the matrix under consideration. For example, the modern version of the QR algorithm for Hessenberg matrices⁵⁵5A Hessenberg matrix has zero entries above the first superdiagonal, which is the case for $A_{P}$. costs $O(d^{2})$ per iteration [BDG04], [CGXZ07]. The sequence starts with $A_{0}=A$, the matrix whose spectrum is being computed. The principle of each step consists in an orthogonal similarity transform $A_{k+1}=Q_{k}^{T}A_{k}Q_{k}$ where

$$A_{k}=Q_{k}R_{k}\quad\text{with}\quad\begin{cases}Q_{k}\text{ orthogonal}\\ R_{k}\text{ upper triangular}\end{cases}$$

As $Q_{k}^{T}Q_{k}=I$, it boils down to computing $Q_{k}$, $R_{k}$ and then forming $A_{k+1}=R_{k}Q_{k}$. The number of iterations used in practice until $A_{k}$ reaches the Schur form of $A$ (upper triangular) may be $O(d)$ or higher. Companion matrices belong to the Hessenberg class, so, as a splitting algorithm, the QR method costs at least $O(d^{3})$.

Note that the polynomial $P$ is not explicitly evaluated here, so the factor $\operatorname{V}_{d}$ is not included in the arithmetic complexity. However, it appears implicitly in each step as a cost of matrices operations.

The core of many iterative splitting algorithms are Newton’s iterations:

$$z_{n+1}=N_{P}(z_{n})\qquad\text{where}\qquad N_{P}(z)=z-\frac{P(z)}{P^{\prime}(z)}\cdotp$$

(15)

The fixed points of the dynamics of $N_{P}$ are the roots of $P$.

In the vicinity of each simple root, this method converges bi-exponentially. More precisely, if $z_{\ast}=\lim z_{n}$ with $P^{\prime}(z_{\ast})\neq 0$, Taylor’s formula for $P$ can be rewritten

$$z_{n+1}-z_{\ast}=\frac{(z_{n}-z_{\ast})^{2}}{P^{\prime}(z_{n})}\int_{0}^{1}(1-t)P^{\prime\prime}(tz_{\ast}+(1-t)z_{n})dt.$$

For $\rho>0$, one defines $C_{\rho}=\frac{1}{2}\sup|P^{\prime\prime}(z)|/\inf|P^{\prime}(z)|$, each extremum being computed for $z\in\operatorname{D}(z_{\ast},\rho)$. One can always assume that $\rho C_{\rho}<1$ because $C_{\rho}\to\frac{1}{2}\left|\frac{P^{\prime\prime}(z_{\ast})}{P^{\prime}(z_{\ast})}\right|\in\mathbb{R}_{+}$ as $\rho\to 0$, so $\rho C_{\rho}\to 0$. In that case, if $z_{n_{0}}\in\operatorname{D}(z_{\ast},\rho)$ then

$$\forall n\geq n_{0},\quad|z_{n+1}-z_{\ast}|\leq C_{\rho}|z_{n}-z_{\ast}|^{2}\quad\textit{i.e.}\leavevmode\nobreak\ \quad|z_{n}-z_{\ast}|\leq C_{\rho}^{-1}(\rho C_{\rho})^{2^{n-n_{0}}}.$$

(16)