  This article presents a new algorithm to compute all the roots of two
families of polynomials that are of interest for the Mandelbrot set
$\mathcal{M}$ : the roots of those polynomials are respectively the parameters
$c\in\mathcal{M}$ associated with periodic critical dynamics for $f_c(z)=z^2+c$
(hyperbolic centers) or with pre-periodic dynamics (Misiurewicz-Thurston
parameters). The algorithm is based on the computation of discrete level lines
that provide excellent starting points for the Newton method. In practice, we
observe that these polynomials can be split in linear time of the degree.
  This article is paired with a code library \citelib{MLib} that implements
this algorithm. Using this library and about 723 000 core-hours on the HPC
center \textit{Rom\'eo} (Reims), we have successfully found all hyperbolic
centers of period $\leq 41$ and all Misiurewicz-Thurston parameters whose
period and pre-period sum to $\leq 35$. Concretely, this task involves
splitting a tera-polynomial, i.e. a polynomial of degree $\sim10^{12}$, which
is orders of magnitude ahead of the previous state of the art. It also involves
dealing with the certifiability of our numerical results, which is an issue
that we address in detail, both mathematically and along the production chain.
The certified database is available to the scientific community.
  For the smaller periods that can be represented using only hardware
arithmetic (floating points FP80), the implementation of our algorithm can
split the corresponding polynomials of degree $\sim10^{9}$ in less than one
day-core. We complement these benchmarks with a statistical analysis of the
separation of the roots, which confirms that no other polynomial in these
families can be split without using higher precision arithmetic.
