boundary C and if the zeros of p(x) are pairwise distinct and not on the boundary C, then from the residue theorem of residual calculus one gets The identity of the left to the right side of this equation also holds for zeros with multiplicities. By using the Newton identities one is able to compute from those sums of powers the factor of p(x) corresponding to the zeros of p(x) inside G. By polynomial division one also obtains the second factor g(x) in p(x) = f(x)g(x). The commonly used regions are circles in the complex plane. Each circle gives raise to a split of the polynomial p(x) in factors f(x) and g(x). Repeating this procedure on the factors using different circles yields finer and finer factorizations. This recursion stops after a finite number of proper splits with all factors being nontrivial powers of linear polynomials. The challenge now consists in the conversion of this analytical procedure into a numerical algorithm with good running time. The integration is approximated by a finite sum of a numerical integration method, making use of the fast Fourier transform for the evaluation of the polynomials p(x) and p'(x). The polynomial f(x) that results will only be an approximate factor. To ensure that its zeros are close to the zeros of p inside G and only to those, one must demand that all zeros of p are far away from the boundary C of the region G. Basic numerical observation (Schönhage 1982) Let be a polynomial of degree n which has k zeros inside the circle of radius 1/2 and the remaining n-k zeros outside the circle of radius 2. With N=O(k) large enough, the approximation of the contour integrals using N points results in an approximation of the factor f with error , where the norm of a polynomial is the sum of the moduli of its coefficients. Since the zeros of a polynomial are continuous in its coefficients, one can make the zeros of as close as wanted to the zeros of f by choosing N large enough. However, one can improve this approximation faster using a Newton method. Division of p with remainder yields an approximation of the remaining factor g. Now , so discarding the last second order term one has to solve using any variant of the extended Euclidean algorithm to obtain the incremented approximations and . This is repeated until