## Overview

This is a summary of the splitting circle method as implemented1 in Pari. The
overall idea is to recursively split the polynomial $p(x)$ into two equal-
degree factors $F(x),G(x)$ at each step until we reach linear factors which
are the roots.

We can divide the algorithm into two parts: determining a splitting circle
that contains roughly half of the roots on each side and estimating the factor
$F(x)$ consisting of the product of the roots inside the circle. In this part,
we will go in detail on using the splitting circle.

## Finding a splitting circle

### Overview

We start with a monic polynomial $P(x)=x^{n}+a_{1}x^{n-1}+...+a_{n}$ of degree
$n$ where the roots have moduli $\varrho_{1}(P),...,\varrho_{n}(P)$ arranged
in increasing order. We would like to find a circle that contains some roots
but not all and that the two collections of roots are as far apart as
possible. Intuitively, if we see that the ratio
$\frac{\varrho_{j+1}(P)}{\varrho_{j}(P)}$ is sufficiently large, then the
geometric mean $\sqrt{\varrho_{j+1}(P)\varrho_{j}(P)}$ would make for a good
radius for the splitting circle centered at the origin. Unfortunately, the
roots of $P$ can have the exact same moduli.

Our algorithm has two steps: finding an appropriate center for the splitting
circle and then calculating the radius.

### Finding the center

For our splitting circle, we care about the range of the moduli, as expressed
by the ratio $\frac{\varrho_{n}(P)}{\varrho_{1}(P)}$. We would like to pick a
center such that the range from the new center is sufficiently large.

First, we standardize the polynomial by shifting the coordinates to have the
center of mass of the roots be the origin. Using Vieta’s formulas, this is
equivalent to considering the new polynomial $P_{1}(x)=P(x-a_{1}/n)$. Next, we
would like to scale the coordinates such that the largest modulus is close to
1, which is the change of coordinates $P_{2}(x)=P_{1}(\varrho_{n}(P_{1})x)$.
Finally, for the shift, we consider four points outside the circle:
$v=(2,2i,-2i,-2)$ and consider the shifted polynomials $Q_{j}=P_{2}(x-v_{j})$.
We choose as center the point $v_{j}$ that maximizes the ratio
$\frac{\varrho_{n}(Q_{j})}{\varrho_{1}(Q_{j})}$. Let $Q$ be that maximizing
polynomial. As the center of mass was defined to be the origin in the first
step, we can easily show that $\frac{\varrho_{n}(Q)}{\varrho_{1}(Q)}$ is
bounded below by a constant, specifically $e^{0.3}$.

In the example below, we show the final shift of the center for a sample
centered polynomial $P$. Note that the ratio
$\frac{\varrho_{n}(Q)}{\varrho_{1}(Q)}$ substantially increases after choosing
$(-2,0)$ as the new center.

(0,2i)(0,-2i)(2,0)(-2,0)$Re$$Im$

Figure 1: Finalizing the choice of center of the splitting circle

### Finding the radius

After recentering, we now have a polynomial $P$ with ratio
$\frac{\varrho_{n}(P)}{\varrho_{1}(P)}>e^{\Delta}$ for $\Delta>0$. We would
like to find an annulus $\Gamma$ that contains no roots of $P$ and such that
$n/2$ roots are inside the interior circle.

We define $\Gamma$ by the interior radius $r$ and exterior radius $R$. We can
then initialize $r_{0}=\varrho_{1}(P)$ and $R=\varrho_{n}(P)$ and note that
there is $i_{0}=1$ root in the inner circle and $j_{0}=n-1$ roots in the outer
circle. Afterwards, we will iteratively shrink our annulus until $n/2$ roots
are located in the inner circle, with none in the annulus. At each step
$\ell$, we take the geometric mean $\rho=\sqrt{r_{\ell-1}R_{\ell-1}}$ and
choose our new annulus to be either $\Gamma_{\ell}=(r_{\ell-1},\rho)$ or
$\Gamma_{\ell}=(\rho,R_{\ell-1})$ using a subroutine that finds on the number
of roots $k$ within the disk $\\{|z|=\rho\\}$. If
$k<(i_{\ell-1}+j_{\ell-1})/2)$, we choose the former, otherwise we choose the
latter. We stop when $i_{\ell}=j_{\ell}$, that is, when the annulus contains
no roots. Then there must be an index $k$ such that
$\varrho_{k}(P)<r_{e}ll<R_{\ell}<\varrho_{k+1}(P)$. We conclude by finding
$m=\varrho_{k}(P),M=\varrho_{k+1}(P)$ and setting our splitting circle to the
geometric mean $\\{|z|=\sqrt{mM}\\}$.

We claim that the initial guarantee of
$\frac{\varrho_{n}(P)}{\varrho_{1}(P)}>e^{\Delta}$ implies that
$\frac{R_{\ell}}{r_{\ell}}>e^{\delta/(n-1)}$. Then scaling our splitting
circle to a unit circle by a change of coordinates $z\rightarrow z/R$ gives us
a root-free annulus of the form $(e^{-\delta},e^{\delta})$ where
$\delta=\frac{1}{2}\log(M/m)$.

### Graeffe’s method

To numerically evaluate the modules of the roots of a polynomial $P$, we will
use Graeffe’s method, which follows from the observation that if

$P(x)=(x-z_{1})(x-z_{2})...(x-z_{n})$

then

$P(x)P(-x)=(-1)^{n}(x^{2}-z_{1}^{2})(x^{2}-z_{2}^{2})...(x^{2}-z_{n}^{2})$

We note that the roots of $Q(x^{2}):=P(x)P(-x)$ are the squares of the roots
of $P(x)$, and we will denote this operation as $Q=Graeffe(P)$. If

$P(x)=\sum_{i=0}^{n}a_{i}x^{i}$

We can efficiently calculate $Graeffe(P)$ thanks to the formula

$Q=A^{2}-zB^{2},\qquad A=\sum_{i}a_{2i}z^{i},\qquad B=\sum_{i}a_{2i+1}z^{i}$

We iteratively define $P_{m+1}=Graeffe(P_{m})$ and note that

$P_{m}(x)=a_{0}^{(m)}+a_{1}^{(m)}x+...+a_{n}^{(m)}x^{n}$

Through Vieta’s formulas, we see that at the $m$th iteration,

$\frac{a_{k}^{(m)}}{a_{n}^{(m)}}\sum_{i_{1}<...<i_{n-k}}(z_{i_{1}}*...*z_{i_{n-k}})^{2^{m}}$

where $z_{1},...,z_{n}$ are the roots of $P$. If the moduli are ordered
$|z_{1}|<...<|z_{n}|$, then we get the approximation

$\frac{a_{k}^{(m)}}{a_{0}^{(m)}}=(z_{k+1}*...*z_{n})^{2^{m}}$

so that

$\varrho_{k}(P)=\lim_{m\rightarrow\infty}\left|\frac{a_{k}^{(m)}}{a_{0}^{(m)}}\right|^{2^{-m}}$

Of course, using the Graeffe iteration directly only guarantees convergence
when all roots have distinct moduli and there is no bound as to the rate of
convergence. We will use Graeffe for three distinct cases: calculating the
number of roots in a disk, calculating the kth modulus $\varrho_{k}(P)$, and
calculating the $n$th modulus $\varrho_{n}(P)$, where there are more efficient
algorithms for estimating the largest modulus than in the general case.

### Number of roots in a disk

Given a polynomial $P$, error parameter $\tau$ and a radius $R$, we would like
to find an index $k$ such that $\varrho_{k}(P)e^{-t}<R<\varrho_{k+1}(P)e^{t}$.
Doing so tells us that there are exactly $k$ roots in the disk of radius $R$.
Applying a change of coordinates, we can reduce this problem to the case where
$R=1$. We will get an initial error tolerance of $t$ using a result of
Schonhage and apply Graeffe’s method until the $m$ th iteration of the error
tolerance $t_{m}<\tau$.

We will use the following theorem2 from Schonhage without proof:

Theorem 1 Let $P(x)=a_{0}+a_{1}x+...+a_{n}x^{n}$ be a complex-valued
polynomial with $k$ an integer such that $1\leq k\leq n$. Given the existence
of $c,q>0$ such that

$|a_{k-m}|\leq cq^{m}|a_{k}|,\qquad m=1:k$

we have that

$\varrho_{k}(P)\leq(c+1)q(n-k+1)$

Corollary 1 In addition, if we have that

$|a_{k+m}|\leq cq^{m}|a_{k}|,\qquad m=1:n-k$

then

$\varrho_{k+1}(P)\geq\frac{1}{(c+1)q(n-k+1)}$

Proof. Consider the reciprocal polynomial $\tilde{P}(x)=x^{n}P(\frac{1}{x})$.
We see that $\varrho_{n-k}(\tilde{P})=\varrho_{k}(P)$ and we apply the
inequality from the previous theorem.

Moreover, we see that $c=1,q=1$ satisfies the condition of the theorem for
$k=\arg\max_{0\leq i\leq n}|a_{i}|$. Then applying the theorem and corollary
tells us that $\varrho_{k}(P)\leq 2n$ and $\varrho_{k+1}(P)\geq\frac{1}{2n}$.
Setting $t>\log(2n)$ gives us the desired bound. However, this error tolerance
is unworkably high, with many potentially valid values of $k$. As stated
earlier, we can tighten the error bound with Graeffe’s method.

Let $P_{m}(x)=b_{0}+b_{1}(x)+...+b_{n}x^{n}$ be the $m$th Graeffe iteration of
$P$ and $k=\arg\max_{0\leq i\leq n}|b_{i}|$. As before,
$\varrho_{k}(P_{m})\leq 2n$ and $\varrho_{k+1}(P_{m})\leq\frac{1}{2n}$.
However, taking the $2^{m}$ roots of the moduli tells us that
$\varrho_{k}(P)\leq e^{t_{m}}$ and $\varrho_{k+1}(P)\geq e^{-t_{m}}$ where
$t_{m}=\frac{\log(2n)}{2^{m}}$. Given our fixed error tolerance $\tau$, we
choose the smallest $m$ such that $t_{m}<\tau$, guaranteeing a sufficient
separation of $\varrho_{k}(P)<1<\varrho_{k+1}(P)$ up to multiplication by
$e^{\tau}$.

Implementing this algorithm in double precision requires us to truncate
coefficients of our polynomials. The following result3 from Schatzle says that
sufficiently small perturbations of the coefficients lead to small changes in
the estimated $k$th modulus.

Theorem 2 Let $P$ and $\hat{P}$ are complex polynomials with degree $n>0$ such
that $|\hat{P}-P|<\varepsilon|\hat{P}|$ and let $\mu>1$ be an upper bound of
$\varrho_{k}(P)$. If $\varepsilon<\mu^{-n}2^{-4n}$ then

$|\varrho_{k}(\hat{P})-\varrho_{k}(P)|\leq\frac{2\mu(1+\mu)\varepsilon^{1/n}}{1-4(1+\mu)\varepsilon^{1/n}}$

From there, we get the following corollary:

Corollary 2 Let $P$ and $\hat{P}$ be complex polynomials with degree $n>0$,
let $\tau>0$ and $k$ an integer such that $0\leq k\leq n$. If

$|\hat{P}-P|<\varepsilon|\hat{P}|,\qquad\varepsilon<2^{-4n}\tau^{n}e^{-3n\tau/2}$

then

$\varrho_{k}(\hat{P})\geq e^{-3\tau/4}\implies\varrho_{k}(P)\geq e^{-\tau}$

and

$\varrho_{k+1}(\hat{P})\leq e^{3\tau/4}\implies\varrho_{k}(P)\leq e^{\tau}$

so that running the algorithm on the truncated approximation for a radius of
$R$ will find the largest modulus $\varrho_{k}<1$ up to multiplication by
$e^{-\tau}$.

Proof. To apply the theorem, set $\varepsilon\leq\mu^{-n}2^{-4n}$ for
$\mu=\tau^{-1}e^{-3\tau/2}$, with $\tau<1$. Then

$|\varrho_{k}(\hat{P})-\varrho_{k}(P)|\leq\frac{2\mu(1+\mu)\varepsilon^{1/n}}{1-4(1+\mu)\varepsilon^{1/n}}\leq$
$\frac{2\mu(1+\mu)\mu^{-1}2^{-4}}{1-4(1+\mu)\mu^{-1}2^{-4}}=$
$\frac{1+\mu}{8-2(1+\mu^{-1})}=$ $\frac{1+\mu}{6-2\mu^{-1}}=$
$\frac{1+\tau^{-1}e^{3\tau/2}}{6-2\tau e^{-3\tau/2}}$

On the other hand, assuming the corollary doesn’t hold, we see that

$|\varrho_{k}(\hat{P})-\varrho_{k}(P)|>e^{-3\tau/4}-e^{-\tau}$

Taking the difference of the two expressions

$\frac{d}{d\tau}\frac{1+\tau^{-1}e^{3\tau/2}+2e^{\tau/2}-2e^{3\tau/4}}{6-2\tau
e^{-3\tau/2}}e^{-3\tau/4}-e^{-\tau}$

tells us that $e^{-3\tau/4}-e^{-\tau}$ is below our bound for $\tau<1$, so by
nonnegativity of norms, $\varrho_{k}(P)\geq e^{-\tau}$.

The corresponding implication

$\varrho_{k+1}(\hat{P})\leq e^{3\tau/4}\implies\varrho_{k}(P)\leq e^{\tau}$

holds by Corollary 1.

We use this to create the subroutine NRD.

Algorithm 1 NRD: Calculates number of roots in a disk

Input: $(P,R,\tau)$ where $P$ is a polynomial of degree $n\geq 2$, $R>0$ is
the radius of the disk, and $\tau>0$ is an error parameter
Output: An integer $k$ with $0\leq k\leq n$ such that
$\varrho_{k}(P)e^{-\tau}<R<\varrho_{k+1}(P)e^{\tau}$

1:Set $\tau_{0}=\tau$ and
$\varepsilon_{0}=2^{-4n}\tau_{0}^{n}e^{-3n\tau_{0}/2}$ and calculate the
coefficients of the scaled polynomial $P_{0}=P(Rz)$ with a relative precision
of $\varepsilon_{0}/2(n+1)$.

2:Round $P_{0}$ to a polynomial $\hat{P}_{0}$ such that
$|P_{0}-\hat{P}_{0}|<\frac{\varepsilon_{0}}{2}|\hat{P}_{0}|$.

3:Set counter $m=1$.

4:while $\frac{3}{4}\tau_{m-1}<\log(2n)$ do

5: Increment counter $m=m+1$.

6: Set $P_{m}=Graeffe(\hat{P}_{m-1})$.

7: Set $\tau_{m}=\frac{3}{4}\tau_{m-1}$ and
$\varepsilon_{0}m=2^{-4n}\tau_{m}^{n}e^{-3n\tau_{m}/2}$

8: Round $P_{m}$ to a polynomial $\hat{P}_{m}$ such that
$|P_{m}-\hat{P_{m}}|<\varepsilon_{m}|\hat{P}_{m}$

9:end while

10:With $\hat{P_{m}}(x)=b_{0}+b_{1}x+...+b_{n}x^{n}$, return the index
$k=\arg\max_{0\leq i\leq n}|b_{i}|$

Let’s consider the correctness of the algorithm. After the last step, we have
the inequalities $\varrho_{k}(\hat{P}_{m})<\exp(3\tau_{m}/4)$ and
$\varrho_{k+1}(\hat{P}_{m})>exp(-3\tau_{m}/4)$. By Corollary 2, we see that
$\varrho_{k}(P_{m})<e^{\tau_{m}}$ and $\varrho_{k+1}(P_{m})>e^{-\tau_{m}}$. As
$\varrho_{j}(P_{i})=\varrho_{j}(\hat{P}_{i-1})^{2}$, we induct backwards using
the definition of $\tau_{m}$ to see that
$\varrho_{k}(\hat{P}_{j})<e^{3\tau_{j}/4}$ and
$\varrho_{k+1}(\hat{P}_{j})>e^{-3\tau_{j}/4}$ for $0\leq j\leq m-1$. In
particular, for our choice of $k$, $\varrho_{k}(\hat{P_{0}})<e^{3\tau/4}$ and
$\varrho_{k+1}(\hat{P_{0}})>e^{-3\tau/4}$ so by Corollary 2, we conclude that

$\varrho_{k}(P(Rx))e^{-\tau}<1<\varrho_{k+1}(P(Rx))e^{\tau}$

as desired.

### Finding the kth modulus

Next, we will use Graeffe’s method to estimate the $k$th modulus
$\varrho_{k}(P)$ to within a multiplicative error of $e^{\tau}$ for some error
tolerance $\tau>0$. Let $P(x)=a_{0}+a_{1}x+...+a_{n}x^{n}$. Then we will scale
$P$ by some factor $\rho>0$ with $\tilde{P}(x)=P(\rho
x)=\tilde{a}_{0}+\tilde{a}_{1}x+...+\tilde{a}_{n}x^{n}$ such that there exist
two integers $\ell,h$ with

$\ell<k\leq
h,\qquad|\tilde{a}_{\ell}|=|\tilde{a}_{h}|,|\tilde{a}_{j}|\leq|\tilde{a}_{\ell}|\quad\text{for
}j=0,1,...,n$

Then setting $c=q=1$ allows $\ell,h$ to satisfy Theorem 1, with

$\frac{1}{2n}<\varrho_{\ell}(\tilde{P})\leq\varrho_{k}(\tilde{P})\leq\varrho_{h}(\tilde{P})<2n$

We then apply Graeffe as before to tighten the bounds.

We can choose $\ell$ by considering the largest $x$-coordinate of the corners
less than $k$ in the lower convex envelope of the points
$M_{j}=(j,\log(|a_{j}|))$. Similarly, we can choose $h$ as the smallest
$x$-coordinate of the corners greater than $k$ in the same envelope. We would
like $|\tilde{a}_{\ell}|=|\tilde{a}_{h}|$, which requires
$|a_{\ell}|\rho^{\ell}=|a_{h}|\rho^{h}$ or

$\rho=(\frac{|a_{\ell}|}{|a_{h}|})^{1/(h-\ell)}$

In addition, by convexity, as shown in the diagram below, the line passing
through $\log(a_{\ell})$ and $\log(a_{h})$ upper bounds $\log(a_{j})$ for
$0\leq j\leq n$. We switch to $\log(|a_{j}|)$ because
$\log(\rho^{x})=x\log(\rho)$ is monotone in $x$. As a result, the linear upper
bound remains an upper bound when composed with a monotone function
$x\rightarrow x\log(\rho)$ for $0\leq x\leq n$. Therefore, all points
$\tilde{M}_{j}=(j,\log(|\tilde{a}_{j}|))$ are below the line connecting
$\tilde{M}_{\ell},\tilde{M}_{h}$. By our argument above, Theorem 1 is
satisfied.

$\ell$$k$$h$$n$$2$$4$$6$$j$$\log(|a_{j}|)$

Figure 2: Convex envelope of log(coefficients)

For convenience, we set $\rho$ to be a power of $2$, with

$\rho=2^{\beta},\qquad\beta=\left\lfloor\frac{1}{h-\ell}\log_{2}\frac{|a_{\ell}|}{|a_{h}|}+1/2\right\rfloor$

We then see that for $m\geq 0$,

$|\tilde{a}_{h-m}|\leq 2^{m/2}|\tilde{a}_{h}|,\qquad|\tilde{a}_{\ell+m}|\leq
2^{m/2}|\tilde{a}_{\ell}|$

so the conditions of Theorem 1 are valid for $c=1,q=\sqrt{2}$. Due to rounding
error, it is better to work with an upper bound of $c=1,q=3/2$, which gives us
looser initial bounds of

$\frac{1}{3n}<\varrho_{\ell}(\tilde{P})\leq\varrho_{k}(\tilde{P})\leq\varrho_{h}(\tilde{P})<3n$

We then use the following result1 from Gourdon to guarantee the accuracy of
our estimates given a truncated initial polynomial.

Theorem 3 Let $P,\hat{P}$ be two polynomials of degree $n>0$, let $\tau>0$ and
$k$ an integer such that $0\leq k\leq n$. If

$|\hat{P}-P|<\varepsilon|\hat{P}|,\qquad\varepsilon<2^{-n-1}(3n)^{-n}\tau^{n}e^{-3n\tau/2}\qquad\frac{1}{3n}<\varrho_{k}(P)<3n$

Then

$\varrho_{k}(\hat{P})e^{-\tau}\leq\varrho_{k}(P)\leq\varrho_{k}(\hat{P})e^{\tau}$

Algorithm 2 MOD: Evaluates the moduli of the roots.

Input: $(R,k,\tau)$ where $P(x)=a_{0}+a_{1}x+...+a_{n}x^{n}$ is a polynomial
of degree $n\geq 2$, $k$ is an integer with $1\leq k\leq n$, and $\tau>0$ is
an error parameter.
Output: A floating point number $R>0$, such that
$Re^{-\tau}\leq\varrho_{k}(P)\leq Re^{\tau}$.

1:If $a_{0}=a_{1}=...=a_{n-1}=0$, return $R=0$.

2:Otherwise, set $P_{0}=P$ and $\tau_{0}=\tau/8$.

3:Calculate $\tilde{P}_{0}(x)=P_{0}(\rho_{0}x)$ for
$\rho_{0}:=2^{\beta},\quad\beta:=\left\lfloor\frac{1}{h-\ell}\log_{2}\frac{|a_{\ell}|}{|a_{h}|}+1/2\right\rfloor$
where $h,\ell$ are defined as in the convex envelope for $P_{0}$.

4:Set $M$ to be the smallest natural number such that $2^{-M}\log(3n)<\tau/2$

5:for $m=1,2,...,M$ do

6: Round the polynomial $\tilde{P}_{m-1}$ to a polynomial $\hat{P}_{m-1}$ such
that
$|\tilde{P}_{m-1}-\hat{P}_{m-1}|<\varepsilon_{m-1},\quad\varepsilon_{m-1}:=2^{-(n+1)}(3n)^{-n}\tau_{m-1}^{n}e^{-3n\tau_{m-1}/2}$

7: Calculate $P_{m}=Graeffe(\hat{P}_{m-1})$

8: Calculate $\tilde{P}_{m}(x)=P_{m}(\rho_{m}x)$ where $\rho_{m}$ is defined
as in the convex envelope for $P_{m}$.

9: Set $\tau_{m}=\frac{3}{2}\tau_{m-1}$

10:end for

11:Return $R=\rho_{0}\rho_{1}^{2^{-1}}\rho_{2}^{2^{-2}}...\rho_{M}^{2^{-M}}$

Let us evaluate the correctness of our algorithm. For each step $m$, by our
change of coordinates we get
$\varrho_{k}(P_{m})=\rho_{m}\varrho_{k}(\tilde{P}_{m})$. After applying the
perturbation result of Theorem 3, we see that

$\varrho_{m}\varrho_{k}(\hat{P}_{m})e^{-\tau}<\varrho_{k}(P_{m})<\varrho_{m}\varrho_{k}(\hat{P}_{m})e^{\tau}$

Using the relation $\varrho_{k}(P_{m+1})=\varrho_{k}(\hat{P}_{m})^{2}$ gives
us

$\rho_{0}...\rho_{M}^{2^{-M}}\varrho_{k}(\tilde{P}_{M})^{2^{-M}}e^{-\tau^{\prime}}<\varrho_{k}(P_{0})<\rho_{0}...\rho_{M}^{2^{-M}}\varrho_{k}(\tilde{P}_{M})^{2^{-M}}e^{\tau^{\prime}},\quad\tau^{\prime}=\tau_{0}+\frac{\tau_{1}}{2}+...+\frac{\tau_{M-1}}{2^{M-1}}$

From Theorem 3, we see that as $\frac{1}{3n}<\varrho_{k}(\tilde{P}_{M})<3n$,
then

$e^{-\tau/2}<\varrho_{k}(\tilde{P}_{M})^{2^{-M}}<e^{\tau/2}$

Moreover,

$0\leq\tau^{\prime}=\frac{\tau}{8}\left(1+\frac{3}{4}+\left(\frac{3}{4}\right)^{2}+...\right)\leq\frac{\tau}{2}$

so we conclude that

$Re^{-\tau}<\varrho_{k}(P_{0})<Re^{\tau}$

for $R=\rho_{0}...\rho_{M}^{2^{-M}}$.

### Upper bounding the moduli of the roots

For the sketch of our splitting circle algorithm, we need to calculate
$\varrho_{n}(P)$ often. It turns out that this can be accomplished more
efficiently than $\varrho_{k}(P)$ for a general $k$. Again, given an error
parameter $\tau>0$, we would like to find an $R>0$ equal to $\varrho_{n}(P)$
up to multiplicative error, such that

$Re^{-\tau}\leq\varrho_{n}(P)\leq Re^{\tau}$

The idea of the algorithm is the same as in the previous case. However, we
would like to use a different scalar $\rho$. Given an initial
$P=a_{0}+a_{1}x+...+a_{n}x^{n}$, we define $\tilde{P}(x)=P(\rho x)$ such that

$\frac{|\tilde{a}_{j}|}{|\tilde{a}_{0}|}\leq 2^{j}{n\choose j}\quad\text{for
}j=1,...,n\quad\text{and there exists }h,1\leq h\leq n\text{ such that
}\frac{|\tilde{a}_{j}|}{|\tilde{a}_{0}|}\geq{n\choose h}$

and $\rho$ can be determined by

$\rho=2^{\beta},\qquad\beta=\max_{1\leq j\leq
n}\left\lfloor\frac{1}{j}\log\frac{|a_{\ell}|}{|a_{0}|{n\choose
j}}\right\rfloor$

By Vieta’s equations, we get the inequality

$|\tilde{a}_{h}|\leq|\tilde{a}_{0}|{n\choose h}\varrho_{n}(\tilde{P})^{h}$

so $\varrho_{n}(\tilde{P})\geq 1$ when we account for the second assumption
for our choice of $\rho$. Setting $c=1,q=2n$, we see that $2^{j}{n\choose
j}\leq(2n)^{j}$, so Theorem 1’s conditions are satisfied and give an upper
bound of $1\leq\varrho_{n}(\tilde{P})\leq 4n$. We continue with a perturbation
result.

Theorem 4 Let $P,\hat{P}$ be polynomials of degree $n>0$ such that

$|P-\hat{P}|<|\beta|e^{-n\tau}\tau^{n}$

where $\beta$ is the leading coefficient of $P$. Then if $\varrho_{n}(P)\geq
1$, we have that

$\varrho_{n}(\hat{P})e^{-\tau}<\varrho_{n}(P)<\varrho_{n}(\hat{P})e^{\tau}$

Proof. First, let $\varrho_{n}(\hat{P})>\varrho_{n}(P)$ and let $z$ designate
a root of $\hat{P}$ with modulus $\varrho_{n}(P)>1$. Then we have that

$|P(z)|=|P(z)-\hat{P}(z)|\leq|z|^{n}|P-\hat{P}|\leq\varrho_{n}(\hat{P})^{n}|\beta|e^{-n\tau}\tau^{n}$

Moreover, we have

$|P(z)|\geq|\beta|\prod_{i=1}^{n}(|z|-\varrho_{i}(P))\geq|\beta|(\varrho_{n}(\hat{P})-\varrho_{n}(P))^{n}$

which gives us the bound

$\varrho_{n}(\hat{P})-\varrho_{n}(P)\leq\varrho_{n}(\hat{P})e^{-\tau}\tau$

Applying the inequality $1+x\leq e^{x}$ gives us the desired bounds. An
analogous argument holds for $\varrho(\hat{P})<\varrho_{n}(P)$ by switching
the roles of the two terms.

We now describe the algorithm to find the largest modulus of the roots. The
code and justification follow largely from the previous example.

Algorithm 3 MODMAX: Evaluates the largest modulus of the roots.

Input: $(R,\tau)$ where $P(x)=a_{0}+a_{1}x+...+a_{n}x^{n}$ is a polynomial of
degree $n\geq 2$ and $\tau>0$ is an error parameter.
Output: A floating point number $R>0$, such that
$Re^{-\tau}\leq\varrho_{n}(P)\leq Re^{\tau}$.

1:If $a_{0}=a_{1}=...=a_{n-1}=0$, return $R=0$.

2:Otherwise, set $P_{0}=P$ and $\tau_{0}=\tau/8$.

3:Calculate $\tilde{P}_{0}(x)=P_{0}(\rho_{0}x)$ for
$\rho=2^{\beta},\qquad\beta=\max_{1\leq j\leq
n}\left\lfloor\frac{1}{j}\log\frac{|a_{\ell}|}{|a_{0}|{n\choose
j}}\right\rfloor$

4:Set $M$ to be the smallest natural number such that $2^{-M}\log(4n)<\tau/2$

5:for $m=1,2,...,M$ do

6: Round the polynomial $\tilde{P}_{m-1}$ to a polynomial $\hat{P}_{m-1}$ such
that
$|\tilde{P}_{m-1}-\hat{P}_{m-1}|<\varepsilon_{m-1},\quad\varepsilon_{m-1}:=|\beta+{m-1}\tau_{m-1}^{n}e^{-n\tau_{m-1}}$

7: where $\beta_{m-1}$ is the leading coefficient of $\tilde{P}_{m-1}$.

8: Calculate $P_{m}=Graeffe(\hat{P}_{m-1})$

9: Calculate $\tilde{P}_{m}(x)=P_{m}(\rho_{m}x)$ where $\rho_{m}$ is defined
in terms of $P_{m}$ for $P_{0}$.

10: Set $\tau_{m}=\frac{3}{2}\tau_{m-1}$

11:end for

12:Return $R=\rho_{0}\rho_{1}^{2^{-1}}\rho_{2}^{2^{-2}}...\rho_{M}^{2^{-M}}$.

Applying the algorithm to the reciprocal $P^{*}(x)=x^{n}P(\frac{1}{x})$ allows
us to also find the smallest modulus of the roots.

Algorithm 4 MODMIN: Evaluates the smallest modulus of the roots.

Input: $(R,\tau)$ where $P(x)=a_{0}+a_{1}x+...+a_{n}x^{n}$ is a polynomial of
degree $n\geq 2$ and $\tau>0$ is an error parameter.
Output: A floating point number $R>0$, such that
$Re^{-\tau}\leq\varrho_{1}(P)\leq Re^{\tau}$.

1:If $a_{0}=0$, return $R=0$.

2:Otherwise, define the reciprocal polynomial
$Q(x)=a_{n}+a_{n-1}(x)+...+a_{0}x^{n}$ of $P$ and return 1/MODMAX($Q,\tau$).

We now have all the tools we need to find the splitting circle of a
polynomial.

### Algorithms

We can use the previous algorithms to to construct a new algorithm to find a
splitting circle. We will break the search into two parts: finding the center
and finding the radius of the splitting circle.

First, assume the center has already been determined. Finding the center
should also give us an annulus $r<|x|<R$ along with two indices $i,j$ such
that

$1\leq i<j\leq n-1,\qquad\varrho_{i}(P)<r<R<\varrho_{j+1}(P)$

We now have to find a radius $\rho$ in the interval $(r,R)$ from which we can
factorize $P$ into two nontrivial factors. The general idea was discussed
above, and we will deal with practical implementation here.

Algorithm 5 RAD: Returns the radius of the splitting circle.

Input: $(P,r,R,i,j)$ where $P$ is a polynomial of degree $n\geq 2$ that
satisfies the centering conditions above
Output: $(\rho,k,\delta)$, where $\rho>0$ is such that the circle $|z|=\rho$
contains $k$ roots of $P$ with $i\leq k\leq j$ and $\delta>0$ such that no
root of $P$ is within the annulus $\rho e^{-\delta}<|z|<\rho e^{\delta}$

1:We calculate the error margin $\Delta:=\log(R/r)$

2:If $i=j$, then the annulus $r<|z|<R$ doesn’t contain any roots of $P$ and we
can do the following steps:

3: We approximate the moduli $\varrho_{i}(P)$ and $\varrho_{i+1}(P)$ by
$m=MOD(P,i,\Delta/8)$ and $M=MOD(P,i+i,\Delta/8)$

4: We return $(\rho,\delta)$ where $\rho=\sqrt{mM}$ and
$\delta=\frac{1}{2}\log(M/m)-\frac{\Delta}{8}$.

5:Otherwise, we set $\rho=\sqrt{rR},\delta=\frac{\Delta}{8(j-i)}$ and
$k=NRD(P,\rho,\delta)$

6: If $k<(i+j)/2$, return $RAD(R,r,\rho e^{-\delta},i,k)$.

7: If $k>(i+j)/2$, return $RAD(R,\rho e^{\delta},R,k,j)$.

8: If $k=(i+j)/2$ and $k<n/2$, return $RAD(R,r,\rho e^{-\delta},i,k)$.

9:Otherwise if $k=(i+j)/2$ and $k\geq n/2$, return $RAD(R,\rho
e^{\delta},R,k,j)$.

We see that $\delta$ is always greater than $c/(j-i+1)$ for some positive
constant $c$. Furthermore, when the condition $i=j$ is reached, then we have
an annulus $r<|z|<R$ that contains no roots of $P$ and so our final radius
does give us a splitting circle.

In later sections, we will come up with an algorithm
$FCS(P,k,\delta,\varepsilon)$ that splits $P$ into two approximate factors
$\hat{F},\hat{G}$ such that $|P-\hat{F}\hat{G}|<\varepsilon|P|$ given a unit
splitting circle. Combining $FCS$ with $RAD$ gives us an algorithm that
factorizes a centered polynomial $P$.

Algorithm 6 HOM: Factorizes from a dilation of the splitting circle

Input: $(P,r,R,i,j,\varepsilon)$ where $P$ is a polynomial of degree $n\geq 2$
that satisfies the centering conditions above, $\varepsilon>0$ is an error
parameter
Output: $(\hat{F},\hat{G})$, which are approximations of factors $F,G$ such
that $|P-\hat{F}\hat{G}|<\varepsilon|P|$

1:Find $(\rho,k,\delta)=RAD(P,r,R,i,j)$.

2:Find an approximation $\hat{Q}$ of $Q(z)=P(\rho z)$ with a relative
precision of $\varepsilon^{\prime}/n$ for
$\varepsilon^{\prime}=\frac{1}{4}\min(\rho^{-n},\rho^{n})\varepsilon$

3:Calculate $(F_{0},G_{0})=FCS(Q,k,\delta,\varepsilon^{\prime})$

4:Rescale the factors $\hat{F}(z)=F_{0}(z/\rho),\hat{G}(z)=G_{0}(z/\rho)$ with
a relative precision of $\varepsilon^{\prime\prime}=2^{-(n+4)}\varepsilon/n$
and return the scaled factors $(\hat{F},\hat{G})$.

We will now prove that the precision specified gives us
$|P-\hat{F}\hat{G}|<\varepsilon|P|$. In the second step, the approximation
$\hat{Q}$ of $Q$ satisfies

$|\hat{Q}-Q|<\varepsilon^{\prime}|Q|$

Moreover, $F_{0},G_{0}$ satisfy

$|\hat{Q}-\hat{F}\hat{G}|<\varepsilon^{\prime}|\hat{Q}|<\varepsilon^{\prime}|\hat{Q}|<2\varepsilon^{\prime}|Q|$

so

$|Q-\hat{F}\hat{G}|<3\varepsilon^{\prime}|Q|$

The dilation scales errors by a factor at most $\max(\rho^{n},\rho^{-n})$ so
scaling tells us that

$|P-FG|<\frac{3}{4}\varepsilon|P|$

where $F(z)=F_{0}(z/\rho),G(z)=G_{0}(z/\rho)$. Next, from our choice of
$\varepsilon^{\prime\prime}$, we have the inequalities
$|F-\hat{F}|<2^{-n-4}\varepsilon|F|$ and $|G-\hat{G}|<2^{-n-4}\varepsilon|G|$.
Multiplying these together tells us that

$|FG-\hat{F}\hat{G}|<2^{-n-4}\varepsilon(|F|\cdot|G|+|F|\cdot|\hat{G}|)\leq
2^{-n-2}\varepsilon(|F|\cdot|G|)$

From Schonhage, we get the inequality $|F|\cdot|G|\leq 2^{n-1}|FG|$, which is
slightly tighter than the classical inequality $|F|\cdot|G|\leq 2^{n}|FG|$.
This leads to the bound

$|FG-\hat{F}\hat{G}|<\frac{\varepsilon}{8}|FG|\leq\frac{\varepsilon}{4}P$

Putting all of our inequalities together,

$|P-\hat{F}\hat{G}|\leq|P-FG|+|FG-\hat{F}\hat{G}|<\varepsilon|P|$

as desired.

We now proceed to find the center of the splitting circle. The center $u$ for
the polynomial $Q(z)=P(z+u)$ has to be chosen so that the modular ratio
$\frac{\varrho_{n}(Q)}{\varrho_{1}(Q)}$ is sufficiently large. In practice, we
need to preprocess to assure the calculations are well-conditioned. This is
realized by our algorithm $CTR0$ which calls on the centering algorithm $CTR$
proper.

Algorithm 7 CTR0: Factorizes a given polynomial of degree $n\geq 2$

Input: $(P,\varepsilon)$ where $P$ is a polynomial of degree $n\geq 2$,
$\varepsilon>0$ is an error parameter
Output: $(\hat{F},\hat{G})$, which are approximations of factors $F,G$ such
that $|P-\hat{F}\hat{G}|<\varepsilon|P|$

1:If the constant coefficient $a_{0}$ of $P$ satisfies
$|a_{0}|<\varepsilon|P|$, return $(z,(P(z)-P(0))/z)$.

2:Otherwise, set $k_{1}=NRD(P,1.9,0.05)$. If $k_{1}=n$, then
$\varrho_{n}(P)\leq 1.9e^{.05}<2$ and we return $CTR(P,\varepsilon)$.

3:Otherwise, consider the reciprocal polynomial $Q=P^{*}z^{n}P(1/z)$ and
calculate $k_{2}=NRD(Q,1.9,0.05)$. If $k_{2}=n$, then calculate
$(F_{0},G_{0})=CTR(Q,\varepsilon)$ and return $(F_{0}^{*},G_{0}^{*})$.

4:Otherwise, set $R=1.9e^{0.05}$, $r=1/R$, and return
$HOM(P,r,R,n-k_{2},k_{1}-1,\varepsilon)$.

The algorithm $CTR0$ allows us to use $CTR$, which requires that the moduli of
the roots are not too high or even directly move to $HOM$ if the constraints
are satisfied, giving us a stability condition.

Algorithm 8 CTR: Finds the center of the splitting circle

Input: $(P,\varepsilon)$ where $P(z)=a_{0}+a_{1}z+...+a_{n}z^{n}$ is a
polynomial of degree $n\geq 2$ with $\varrho_{n}(P)\leq 2$, $\varepsilon>0$ is
an error parameter
Output: $(\hat{F},\hat{G})$, which are approximations of factors $F,G$ such
that $|P-\hat{F}\hat{G}|<\varepsilon|P|$

1:Set $P_{0}(z)=P(z+u)$ where $u=\frac{-a_{1}}{na_{0}}$ is the center of mass
of the roots.

2:Find an approximation $\hat{P}_{0}$ of $P_{0}$ with a relative precision of
$\varepsilon_{0}^{\prime}=\frac{1}{4}\varepsilon(1+|u|)^{-n}$.

3:Set $\varepsilon_{0}=\frac{\varepsilon}{4}(1+|u|)^{-n}\frac{|P|}{|P_{0}|}$

4:If the constant term $c$ of $\hat{P}_{0}$ satisfies
$|c|<\varepsilon_{0}|\hat{P}_{0}|$ then skip to the last step with the
polynomials $\hat{F}_{0}(z)=z,\hat{G}_{0}(z)=(\hat{P}_{0}(z)-c)/z$.

5:Calculate an approximation $r=MODMAX(\hat{P}_{0},0.01)$ of the largest
modulus of the roots of $\hat{P}_{0}$, then check that $\rho=re^{-0.01}<4$.

6:Calculate an approximation $\hat{P}_{1}$ of the dilation
$P_{1}(z)=\rho^{-n}\hat{P}_{0}(\rho z)$ with a relative precision of
$\varepsilon_{1}/n$ for
$\varepsilon_{1}=\frac{1}{4}\varepsilon_{0}\sup(1,\rho)^{-n}$

7:For $j=0,1,2,3$ calculate the polynomial $Q_{j}=\hat{P}_{1}(z+2e^{ij\pi/2})$
and the values $R_{j}=MODMAX(Q_{j},0.01),r_{j}=MODMIN(Q_{j},0.01)$.

8:Let $j_{0}$ be the index that maximizes the ratio $R_{j}/r_{j}$. Set
$P_{2}=Q_{j_{0}}$ and $v=2e^{ij_{0}\pi/2}$.

9:Calculate $(F_{2},G_{2})=HOM(P_{2},r,R,1,n-1,\varepsilon_{2})$ with
$r=r_{j_{n}}e^{0.01}$, $R=R_{j_{n}}e^{-0.01}$ and
$\varepsilon_{2}=\varepsilon_{1}\frac{|\hat{P}_{0}|}{|P_{2}|}3^{-n}$.

10:Calculate $F_{1}(z)=F_{2}(z-v),G_{1}(z)=G_{2}(z-v)$.

11:Find the approximations $\hat{F}_{0},\hat{G}_{0}$ of the scaled polynomials
$F_{0}=\rho^{k}F_{1}(z/p),F_{0}=\rho^{n-k}G_{1}(z/p)$ where $k=deg(F_{1})$.
Work with a relatie precision of $\varepsilon^{\prime}/n$ for
$\varepsilon^{\prime}=2^{-n-4}\varepsilon_{0}$.

12:Calculate the approximations $\hat{F},\hat{G}$ of
$F=\hat{F}_{0}(z-u),G=\hat{G}_{0}(z-u)$ to a relative accuracy of
$\varepsilon^{\prime\prime}/n$ for $\varepsilon=2^{-n-4}\varepsilon$ and
return $(\hat{F},\hat{G})$.

We will now show that the precision of our calculations is sufficient. After
calling HOM in line 9, we have that
$|P_{2}-F_{2}G_{2}|<\varepsilon_{2}|P_{2}|$, and the transformation
$R(z)\rightarrow R(z-v)$ could amplify errors by a factor of at most
$(1+|v|)^{n}=3^{n}$ at most, giving us a bound of

$|\hat{P}_{1}-F_{1}G_{1}|<3^{n}\varepsilon_{2}|P_{2}|=\varepsilon_{1}|\hat{P}_{0}|$

As $|\hat{P}_{1}-P_{1}|<\varepsilon_{1}|\hat{P}_{0}|$, we see that

$|P_{1}-F_{1}G_{1}|<2\varepsilon_{1}|\hat{P}_{0}|$

Meanwhile, the dilation in line 11 gives us

$|\hat{P}_{0}-F_{0}G_{0}|<2\varepsilon_{1}\sup_{1,\rho}^{n}|\hat{P}_{0}|=\frac{\varepsilon_{0}}{2}|\hat{P}_{0}|$

We repeat the argument to show that the precision used at step 11 to find
$\hat{F}_{0},\hat{G}_{0}$ implies

$|\hat{P}_{0}-\hat{F}_{0}\hat{G}_{0}|<\varepsilon_{0}|\hat{P}_{0}|$

Next, as $|\hat{P}_{0}-P_{0}|<\varepsilon_{0}|\hat{P}_{0}|$ we see by the
triangle inequality that

$|P_{0}-\hat{F}_{0}\hat{G}_{0}|<2\varepsilon_{0}|\hat{P}_{0}|$

and therefore

$|P_{0}-FG|<2\varepsilon_{0}(1+|u|)^{n}|\hat{P}_{0}|=\frac{\varepsilon}{2}|P|$

The precision used to calculate $\hat{F},\hat{G}$ ensures that
$|P-\hat{F}\hat{G}|<\varepsilon|P|$, which guarantees the desired output.

Note that for $|u|\leq 2,\rho\leq 4$, we have that

$\varepsilon_{2}=\frac{1}{4}(1+|u|)^{-n}\sup(1,\rho)^{-n}3^{-n}\frac{|P|}{|P_{2}|}\varepsilon\geq\frac{1}{4}36{-n}\frac{|P|}{|P_{2}|}\varepsilon$

In addition, we have that

$|\hat{P}_{0}|\leq 2|P_{0}|\leq 2(1+|u|)^{n}|P|\leq 2\cdot 3^{n}|P|$

Next, all the roots of $P_{1}$ are in the interior of the unit disk, so if
$\alpha$ is the leading coefficient of $P_{1}$, then
$|P_{1}|\leq|\alpha|2^{n}\leq|\hat{P}_{0}|2^{n}$ as $\alpha$ is also the
dominant coefficient of $\hat{P}_{0}$. We then have that
$|\hat{P}_{1}|\leq|\hat{P}_{0}|2^{n+1}$. Using the inequality $|P_{2}|\leq
3^{n}|\hat{P}_{1}|$ we get that

$\frac{|P_{2}|}{|P|}\leq 3^{2n}\cdot 2^{n+2}=4*18^{n}$

Putting all of our inequalities alows us to conclude that
$\varepsilon_{2}\geq\frac{1}{16}648^{-n}\varepsilon$. This lower bound of
$\varepsilon_{2}$ makes sure that it is sufficiently large with respect to
$\varepsilon$ that the calculations are well-conditioned.

## Using a splitting circle to find the factors

### Overview

Assume that we have a splitting circle $C$ that contains roughly half of the
roots of our polynomial $p(x)$. Then by a change of coordinates, we can assume
that our splitting circle is the unit circle. We would like to find the factor
$F(x)$ that consists of the product of the roots in the circle.

We first find an initial approximation $F_{0}$ of $F$ using the residue
theorem. By polynomial division, we get the approximation $G_{0}$ of the
remainder. We then use a form of Newton’s method to get better approximations
$F_{1},G_{1}$ of $F,G$. We repeat the process until $F_{1},G_{1}$ are linear
factors.

### Finding $F_{0}$

By the residue theorem, we note that along the boundary of our unit circle
$C$:

$\frac{1}{2i\pi}\oint_{C}\frac{p^{\prime}(z)}{p(z)}z^{m}dz=\sum_{|z|<1,p(z)=0}z^{m}=W_{m}$

Using Newton’s identities, we can turn the sums of powers from $m=1,...,n$ to
the coefficients $\phi_{i}$ of $F$ using the recursive formulas:

$\phi_{1}=-W_{1}$
$\phi_{m}=-\frac{1}{m}(W_{1}\phi_{m-1}+...+W_{m-1}\phi_{1}+W_{m})$

For a large enough $N$, we can approximate the contour integral $W_{m}$ with
the discrete sum

$W_{m}=\frac{1}{N}\sum_{j=0}^{N-1}\frac{p^{\prime}(\omega^{j})}{p(\omega^{j})}\omega^{(m+1)j}$

where $\omega$ is the $N$th root of unity. Applying Newton’s identities gives
us our estimates $\phi_{i}$ and our initial approximation $F_{0}$ of the
factor $F$, which should be within $O(e^{-\delta N})$ of the true value.

### Newton-Schonhage iteration

For our final approximation, we would like to find correction factors $f,g$
with degrees $deg(f)<deg(F_{0}),deg(g)<deg(G_{0})$ such that
$F_{1}:=F_{0}+f,G_{1}:=G_{0}+g$ are sufficiently good approximations of the
true factors $F,G$. Specifically, we would like to minimize the L1 norm

$\min_{f,g}||P-(F_{0}+f)(G_{0}+g)||_{1}$

To the first order, we have that

$P-(F_{0}+f)(G_{0}+g)\approx P-F_{0}G_{0}-fG_{0}-gF_{0}$

so we would like to find $f,g$ such that

$P-F_{0}G_{0}=fG_{0}+gF_{0}$

We can ignore $g$ in our calculations by taking the modular representation

$fG_{0}\equiv P\mod F_{0}$

Instead of solving directly for $f$, we will calculate an ”auxiliary
polynomial” $H$ with degree $deg(H)<deg(F_{0})$ such that

$HG_{0}\equiv 1\mod F_{0}$

in which case we can simply find $f$ by

$f\equiv f\cdot HG_{0}\equiv HP\mod F_{0}$

From there, we get the updated estimates $F_{1}=F_{0}+f$ and $G_{1}=P/F_{1}$.

It turns out that $H$ can be represented as another complex integral, with the
following result.

Theorem: Given our unit splitting circle $C$ with $P=FG$ such that $F$ only
has roots inside of $C$ and $G$ does not,

$HG_{0}\equiv 1\mod F_{0}$

is uniquely defined by

$H=\frac{1}{2i\pi}\oint_{C}\frac{1}{F_{0}G_{0}}\frac{F_{0}(z)-F_{0}(t)}{z-t}dt$

Proof. Let

$I(z):=\frac{1}{2i\pi}\oint_{C}\frac{H(t)}{F_{0}(t)}\frac{F_{0}(z)-F_{0}(t)}{z-t}dt$

We can decompose our modular equation into the form

$\frac{1}{F_{0}G_{0}}=R+\frac{H}{F_{0}}+\frac{B}{G_{0}}$

where $deg(A)<deg(F_{0}),deg(B)<deg(G_{0})$ and we have that $HG_{0}\equiv
1\mod F_{0}$.

We note that the map
$x\rightarrow(R(x)+\frac{B(x)}{G_{0}(x)})\frac{F_{0}(z)-F_{0}(x)}{z-x}$ is
analytic in the unit disk, so we can use Cauchy’s integral formula to get

$I(z)=\frac{1}{2i\pi}\oint_{|t|=1}\frac{1}{F_{0}(t)G_{0}(t)}\frac{F_{0}(z)-F_{0}(t)}{z-t}dt=$
$F_{0}(z)\frac{1}{2i\pi}\oint_{|t|=1}\frac{H(t)}{F_{0}(t)}\frac{1}{z-t}dt-\frac{1}{2i\pi}\oint_{|t|=1}\frac{H(t)}{z-t}dt$

We note that for $|z|>1$, the second term vanishes, giving us

$\frac{I(z)}{F_{0}(z)}=\frac{1}{2i\pi}\oint_{|t|=1}\frac{H(t)}{F_{0}(t)}\frac{1}{z-t}dt$

Changing the circle from $|t|=1$ to $|t|=R$ and sending $R$ to infinity gives
us a unique residue of $H(z)/F_{0}(z)$ and the integral is 0 for $R=\infty$,
as $H$ has degree less than $F$. We conclude that $H(z)=I(z)$ for all $z$ with
modulus greater than 1, so by the identity theorem, our claim is valid.

Assuming that the polynomial
$F_{0}(z)=\phi_{0}z^{k}+\phi_{1}z^{k-1}+...+\phi_{k}$, we see by linearity of
integrals that $H$ would be equal to

$H(z)=\sum_{\ell=0}^{k-1}\left(\sum_{m=\ell+1}^{k}\phi_{k-m}v_{m-\ell}\right)z^{\ell}$

with

$v_{m}=\frac{1}{2i\pi}\oint_{C}\frac{t^{m-1}}{F_{0}G_{0}(t)}dt$

Again, we can approximate this integral around the unit circle with discrete
sums of values at roots of unity. We define

$U_{m}=\frac{1}{N}\sum_{j=0}^{N-1}\frac{\omega^{mj}}{P(\omega^{j})}$

One can show that $U_{m}$ is sufficiently close to $v_{m}$ given a large
enough $N$ and that therefore we can make an initial guess $H_{0}$ of $H$ by

$H_{0}(z)=\sum_{\ell=0}^{k-1}\left(\sum_{m=\ell+1}^{k}\phi_{k-m}U_{m-\ell}\right)z^{\ell}$

with $|H-H_{0}|=O(N^{n-1}e^{-\delta N})$. We then apply Newton’s method to get
an improved guess $H_{1}$ of $H$ via the following iteration:

Given $H_{m}$, we find $D_{m}$ by taking the modulus

$H_{m}G_{0}\equiv 1-D_{m}\mod F_{0}$

We then find $H_{m+1}$ by taking the modulus

$H_{m+1}G_{0}\equiv H_{m}(1+D_{m})\mod F_{0}$

We note that if $H_{0}$ is a sufficiently good guess, then the sequence
$\\{D_{m}\\}$ will converge quadratically to 0, as

$D_{m+1}\equiv 1-H_{m+1}G_{0}\equiv 1-(1-D_{m})(1+D_{m})\equiv D_{m}^{2}\mod
F_{0}$

We iterate until reaching a desired error tolerance of $H$.

### Discrete Fourier Transforms

For ease of implementation, it turns out that the discrete integrals we use
over the unit circle exactly correspond to Discrete Fourier transforms, which
are defined as follows. Let $P$ be an $L$-dimensional vector of coefficients
for our polynomial $p$, then

$DFT(P)_{m}=\sum_{k=0}^{L-1}P_{k}\omega^{mk}$

where $\omega$ is the $L$th root of unity. Then taking a discretization of the
integral over the unit circle is just $DFT(P)/L$. In addition, we can use the
Fast Fourier Transform when $L$ is a power of 2 for an $O(n\log n)$
calculation of polynomial multiplication. We note that for $N=KL$ as described
above:

$W_{m}=\sum_{u=0}^{K-1}Y_{m,u}\omega^{(m+1)u}$
$U_{m}=\sum_{u=0}^{K-1}X_{m,u}\omega^{mu}$

with

$Y_{m,u}=\sum_{v=0}^{L-1}\frac{P^{\prime}(\omega^{u+vK})}{P(\omega^{u+vK})}(\omega^{K})^{(m+1)v}$
$X_{m,u}=\sum_{v=0}^{L-1}\frac{1}{P(\omega^{u+vK})}(\omega^{K})^{mv}$

where we can represent the double sums $W_{m},U_{m}$ as sums of compositions
of discrete Fourier transforms. Note that setting $K=1$ just gives us a direct
composition of DFTs. For tuning purposes, it can be helpful to have $K>1$.

We are now ready to reveal the algorithms for estimating the factors
$F(x),G(x)$ of $p(x)$ given a unit splitting circle and the number of roots
inside the circle.

### Algorithms

We first want to take discrete Fourier transforms to estimate our integrals
$W_{m},U_{m}$:

Algorithm 9 DFT: Find the Fourier transforms $W_{m},U_{m}$

Input: $(N,k,P)$ where $P$ is a polynomial having $k<n$ roots in the unit
disk, $N$ is an integer of the form $N=KL$ with $L$ is a power of $2$ such
that $n<L\leq 2n$.
Output: $(W,U)$ where $W_{m},U_{m}$ are defined for $1\leq m\leq k$

1:Calculate $K,L$ where $N=KL$ and $L$ is a power of $2$ such that $n<L\leq
2n$.

2:Set $W_{m}=0,U_{m}=0$ for $1\leq m\leq k$.

3:Set $\omega=e^{2i\pi/N}$ and $\omega_{0}=\omega^{K}$.

4:for $u=0,1,\ldots,K-1$ do

5: Calculate $\alpha=FFT(P,L)$, where
$\alpha_{v}=\sum_{j=0}^{n}(p_{j}\omega^{uj})\omega_{0}^{vj}$ for $v\in 0:L-1$.

6: Calculate $\beta=FFT(P^{\prime},L)$, where
$\beta_{v}=\sum_{j=0}^{n-1}((j+1)p_{j+1}\omega^{uj})\omega_{0}^{vj}$ for $v\in
0:L-1$.

7: Calculate $\gamma_{v}=1/\alpha_{v}$ for $v\in 0:L-1$.

8: Calculate $x=FFT(\gamma,L)$, where
$x_{m}=\sum_{v=0}^{L-1}\gamma_{v}\omega_{0}^{mv}$ for $m\in 0:L-1$.

9: Set $U_{m}=U_{m}+x_{m}\omega^{mu}$ for $m\in 0:k-1$.

10: Calculate $\delta_{v}=\beta_{v}/\alpha_{v}$ for $v\in 0:L-1$.

11: Calculate $y=FFT(\delta,L)$, where
$y_{m}=\sum_{v=0}^{L-1}\delta_{v}\omega_{0}^{mv}$ for $m\in 0:L-1$.

12: Set $W_{m}=W_{m}+y_{m+1}\omega^{(m+1)u}$ for $m\in 0:k-1$.

13:end for

14:Set $W_{m}=W_{m}/N$ and $U_{m}=U_{m}/N$.

Now that we have the appropriate integrals, we would like to plug them in to
get initial guesses for our factor $F$ and auxiliary polynomial $H$.

Algorithm 10 RES: Initial approximation of $F,H$ given the respective contour
integrals

Input: $(N,k,P,\delta,\varepsilon)$ where $P$ is a polynomial having $k<n$
roots in the unit disk with no roots in the annulus
$(e^{-\delta},e^{\delta})$, $N$ is an integer of the form $N=KL$ with $L$ is a
power of $2$ such that $n<L\leq 2n$, $\varepsilon$ is a preset error
tolerance.
Output: $(F_{0},H_{0})$ which are approximations of $F,H$ respectively.

1:Calculate the contour integrals $(W,U)$ approximately using DFT($P,k,N)$.

2:Approximate the factors $\phi_{m}$ in
$F_{0}=z^{k}+\phi_{1}z^{k-1}+...\phi_{k}$ using the recursive formulas
$\phi_{1}=-W_{i}$
$\phi_{m}=-\frac{1}{m}(W_{1}\phi_{m-1}+...+W_{m-1}\phi_{1}+W_{m})$

3:Set the polynomial
$H_{0}(z)=\sum_{\ell=0}^{k-1}\left(\sum_{m=\ell+1}^{k}\phi_{k-m}U_{m-\ell}\right)z^{\ell}$
and keep the coefficients as a $k$-dimensional vector.

We now improve $H_{0}$ for use in Newton’s method:

Algorithm 11 AUX: Improves approximation of auxiliary polynomial $H$

Input: $(F_{0},G_{0},H_{0},\varepsilon)$ where $F_{0},G_{0},H_{0}$ are
approximations of $F,G,H$, $\varepsilon$ is a preset error tolerance.
Output: $(H_{1})$, which is an approximation of $H$ such that
$H_{1}G_{0}\equiv 1-D\mod F_{0}$ with $|D|<\varepsilon$

1:Set $H_{1}:=H_{0}$.

2:for $k=0,1,2,\ldots$ do

3: Find $D$ by the modular relation $H_{1}G_{0}\equiv 1-D\mod F_{0}$.

4: if $||D||_{1}<\varepsilon$ then

5: Return $H_{1}$.

6: else if $||D||_{1}>1$ then

7: Return 0.

8: end if

9: Set $H_{1}\equiv H_{1}(1+D)\mod F_{0}$

10:end for

We now plug in our improved estimate of $H$ to the Newton-Schonhage method:

Algorithm 12 NS: Improves approximation of factor $F$

Input: $(P,F_{0},H_{0},\varepsilon)$ where $P$ is a polynomial, $F_{0},H_{0}$
are approximations of $F,H$, $\varepsilon$ is a preset error tolerance.
Output: $(F_{1},G_{1})$, which are approximations of $F,G$ such that
$|P-F_{1}G_{1}|<\varepsilon|P|$.

1:Set $F_{1}:=F_{0}$

2:Set $G_{1}:=P/F_{1}$ via Euclidean division.

3:Calculate $\varepsilon_{0}=|P-F_{1}G_{1}|/|P|$.

4:for k = 0,1,2…. do

5: if $\varepsilon_{0}<\varepsilon$ then

6: Return $(F_{1},G_{1})$.

7: else if $\varepsilon_{0}>1$ then

8: Return (0,0).

9: end if

10: Set $H:=$ AUX$(F_{1},G_{1},H_{1},\varepsilon_{0})$.

11: Calculate $f:=H_{1}P\mod F_{1}$.

12: Set $F_{1}:=F_{1}+f$, $G_{1}=P/F_{1}$.

13:end for

Note that the above algorithm could fail to converge if the approximation
$F_{0}$ of $F$ is too far off. We combine all of our steps for the final
algorithm:

Algorithm 13 FCS: Factorization given a splitting circle

Input: $(P,k,\delta,\varepsilon)$ where $P$ is a polynomial of degree $n\geq
2$ possessing $k<n$ roots in the unit disk and none in the annulus
$(e^{-\delta},e^{\delta})$, $\varepsilon>0$ being the preset error tolerance
Output: $(F_{1},G_{1})$, which are approximations of $F,G$ such that
$|P-F_{1}G_{1}|<\varepsilon|P|$.

1:Set $L$ to be the power of 2 such that $n<L\leq 2n$

2:Set $K=\max(1/2\delta,2)$ and set $N=KL$.

3:for k = 0,1,2…. do

4: Find the initial guesses $(F_{0},G_{0})=RES(P,k,\delta,N)$

5: Find the updated factors $F_{1},G_{1}=NS(P,F_{0},H_{0},\varepsilon)$

6: if $(F_{1},G_{1})=(0,0)$ then

7: Set $N:=2N$.

8: else

9: Return $(F_{1},G_{1})$.

10: end if

11:end for

## Complete Factorization

Now that we are capable of finding a splitting circle and using it to factor
our degree $n$ polynomial $P$ into two nontrivial factors $F,G$, we can use
our routine $CTR0$ repeatedly until we’ve reduced the problem to linear
factors $L_{1},L_{2},...L_{n}$, with

$|P-L_{1}L_{2}...L_{n}|<\varepsilon|P|$

This would occur if each intermediary call to $CTR0$ left us with $k$ factors
$P_{1},...,P_{k}$ such that

$|P-P_{1}P_{2}...P_{k}|<\frac{k}{n}\varepsilon|P|$

While $k<n$, we can factorize one of the polynomials $P_{i}$ with degree $>1$.
WLOG, let this be $P_{1}$, with a factorization of

$|P-FG|<\varepsilon_{k}|P_{1}|$

which would give

$|P-FGP_{2}...P_{k}|<\frac{k}{n}\varepsilon|P|+\varepsilon_{k}|P_{1}|\cdot|P_{2}...P_{k}|$

This would be upper bounded by $\frac{k+1}{n}\cdot\varepsilon|P|$ if
$\varepsilon_{k}<\frac{\varepsilon}{n}\frac{|P|}{n|P_{1}|\cdot|P_{2}...P_{k}|}$.
From Schonhage, we have the inequality

$|Q|\cdot|R|<2^{\deg(QR)-1}|QR|$

which implies

$|P_{1}|\cdot|P_{2}...P_{k}|<2^{n-1}|P_{1}...P_{k}|<2^{n-1}(1+\frac{k}{n}\varepsilon)|P|<2^{n}|P|$

so our constraint is modified if we choose

$\varepsilon_{k}=2^{-n}\frac{\varepsilon}{n}$

In summary, if each factorization step of the algorithm $CTR0$ uses a
precision of $2^{-n}\varepsilon/n$, then we have our approximate linear
factorization. That gives us the following algorithm:

Algorithm 14 FACT: Approximate complete factorization

Input: $(P,\varepsilon)$ where $P$ is a polynomial of degree $n\geq 1$,
$0<\varepsilon<1$ being the preset error tolerance
Output: $(L_{1},...,L_{n})$, which are linear factors such that such that
$|P-L_{1}L_{2}...L_{n}|<\varepsilon|P|$.

1:If $n=1$, return $P$.

2:Otherwise, we calculate $(F,G)=CTR0(P,\frac{2^{-n}\varepsilon}{n})$.

3:Let $k=deg(F)$. Then we calculate $(L_{1},...,L_{k})=FACT(F,\varepsilon)$
and $(L_{k+1},...,L_{n})=FACT(G,\varepsilon)$

4:Return $(L_{1},...,L_{n})$.

## References

1\. Gourdon, Xavier. Combinatoire, Algorithmique et Géométrique des Polynômes.
PhD Thesis, École Polytechnique, Paris, 1996.

2\. Schonhage, A. Equation solving in terms of computational complexity. In
Proceedings of the International Congress of Mathematicians (1987).

3\. Schatzle, R. Diploma Thesis Universitat Bonn, 1996.