Separation circle method

The splitting circle method (engl. Splitting circle method ) is a method for numerically factoring polynomials in a variable complex coefficient. This method was proposed in 1982 by Arnold Schönhage in the article The fundamental theorem of algebra in terms of computational complexity (previously only published on the Internet) and implemented in 1996 by Xavier Gourdon in the computer algebra system PARI / GP and subsequently Magma . Since the mid-1990s, u. a. V. Pan and G. Malajovich suggested improvements to the algorithm, which, however, have not yet been implemented anywhere.

By continually breaking down a polynomial into two nontrivial factors, a factorization into linear factors can ultimately be achieved. This is equivalent to finding all complex zeros of the polynomial including the specification of their multiplicity .

In numerical calculations with a fixed finite precision (see floating point number and fixed point number ) it is not possible to differentiate between a multiple occurring zero and an equal group of closely spaced zeros. In this case the result of the procedure is a factorization in

• Linear factors for a sufficiently isolated zero and
• Factors of a higher degree for groups of zeros which cannot be distinguished in the selected precision.

Factoring with the help of the residual calculus

According to Vieta's generalized theorem , the coefficients are a normalized polynomial

${\displaystyle p(x)=x^{n}+p_{n-1}x^{n-1}+\dots +p_{1}x+p_{0}}$

the elementary symmetrical polynomials in the zeros of the polynomial except for a sign . We want to find a decomposition of into a product of two factors , where the first zeros of have as zeros, i. H. ${\displaystyle z_{1},\dots ,z_{n}}$ ${\displaystyle p(x)}$ ${\displaystyle p(x)=f(x)g(x)}$ ${\displaystyle f(x)}$ ${\displaystyle k}$ ${\displaystyle p(x)}$

${\displaystyle p(x)=(x-z_{1})(x-z_{2})\cdots (x-z_{n})}$,

${\displaystyle f(x)=(x-z_{1})(x-z_{2})\cdots (x-z_{k})}$ and

${\displaystyle g(x)=(x-z_{k+1})\cdots (x-z_{n})}$.

The coefficients of and can be determined without having to determine the zeros as an intermediate step. This is possible using the residual calculus and a suitable decomposition of the complex plane. One type of decomposition is the inside and outside of any circle, which is then called the separation circle . ${\displaystyle f(x)}$ ${\displaystyle g(x)}$

Let be a bounded subset of the complex number plane with (piecewise) smooth boundary curve . Then, according to the residual theorem, holds for every in holomorphic function ${\displaystyle K}$ ${\displaystyle \mathbb {C} }$ ${\displaystyle C}$ ${\displaystyle K}$ ${\displaystyle h}$

${\displaystyle \oint _{C}h(z){\frac {p'(z)}{p(z)}}\,dz=2\pi i\sum _{\zeta \in K\,:\,p(\zeta )=0}\mathrm {Res} \left({\textstyle {\frac {h(z)p'(z)}{p(z)}}},\zeta \right)=2\pi i\sum _{\zeta \in K\,:\,p(\zeta )=0}h(\zeta )}$.

If the zeros of lie in the interior of and all other zeros outside of , in particular there is no zero point on the edge curve , then the following applies ${\displaystyle z_{1},\dots ,z_{k}}$ ${\displaystyle p(x)}$ ${\displaystyle K}$ ${\displaystyle K}$ ${\displaystyle C}$

${\displaystyle \sum _{j=1}^{k}h(z_{j})={\frac {1}{2\pi i}}\oint _{C}h(z){\frac {p'(z)}{p(z)}}\,dz}$.

The coefficients occurring in the coefficients of are sums in mixed products of the zeros. The residual formula just given is therefore not directly applicable. However, it is possible to represent the elementary symmetrical polynomials using "unmixed" expressions in the zeros. Any symmetric polynomial in complex numbers can be represented by a polynomial expression in the elementary symmetric polynomials of those numbers. This is especially true for the power sums . Conversely, the elementary symmetric polynomials and thus the coefficients of the polynomial can be obtained from the first power sums, the conversion formulas for this are the Newton identities . The power sums themselves can be obtained from a contour integral using the residual formula . ${\displaystyle f(x)}$ ${\displaystyle k}$ ${\displaystyle z_{1},\dots ,z_{k}}$ ${\displaystyle k}$ ${\displaystyle s_{m}=z_{1}^{m}+z_{2}^{m}+\dots +z_{k}^{m}}$ ${\displaystyle f(x)=(x-z_{1})(x-z_{2})\cdot (x-z_{k})}$ ${\displaystyle k}$ ${\displaystyle h(z)=z^{m},m=1,\ldots ,k}$

So the theoretical factorization algorithm is:

1. Find a smoothly bounded subset that contains some, but not all, roots of . ${\displaystyle K}$ ${\displaystyle p(x)}$
2. Determine the contour integrals which result in the power sums of the zeros. The constant polynomial can also be used to determine the number of zeros it contains. ${\displaystyle hh(z)=1}$
3. Determine the coefficients of using the Newton identities, and the coefficients of using polynomial division . ${\displaystyle f(x)}$ ${\displaystyle g(x)}$

Approximate factorization and its improvement

In the numerical application, the contour integrals cannot be determined exactly. However, the numerical integration can be carried out with any precision by choosing a sufficiently small step size. The approximated factors determined by means of the Newton identities are denoted by and . For a quick execution of the numerical integration, it is advisable to restrict oneself to circles in the complex plane, since then the numerical integration, i. H. the determination of the values ​​of the polynomials at the support points and the determination of the integrals from the values ​​of the quotients can be carried out with the aid of the fast Fourier transformation . ${\displaystyle F_{0}(x)}$ ${\displaystyle G_{0}(x)}$

If the accuracy of the numerical integration is sufficiently high, the coefficients of the "error polynomial" will be arbitrarily small. If this error is of the order of magnitude , then the distance between the zeros of p (x) and the corresponding zeros of the factors is in the worst case the order of magnitude . The numerical integration must be carried out in such a way that with the zeros of p (x) the corresponding zeros from inside K and those from outside of K are also located. ${\displaystyle \textstyle p(x)-F_{0}(x)\,G_{0}(x)}$ ${\displaystyle \varepsilon >0}$ ${\displaystyle \textstyle {\sqrt[{n}]{\varepsilon }}}$ ${\displaystyle \textstyle F_{0}(x)}$ ${\displaystyle \textstyle G_{0}(x)}$

If the last requirement is met, the factorization can be improved using a variant of Newton's method . It follows from the last-mentioned requirement that both f (x) and g (x) and and are coprime. According to the extended Euclidean algorithm, there are polynomials a (x) and b (x) with 1 = af + bg . Let be polynomials for which the coefficients of the error expression also have the order of magnitude . Then you can use improved polynomials ${\displaystyle F_{0}(x)}$ ${\displaystyle G_{0}(x)}$ ${\displaystyle A_{0}(x),B_{0}(x)}$ ${\displaystyle 1-(A_{0}F_{0}+B_{0}G_{0})}$ ${\displaystyle \varepsilon >0}$

• ${\displaystyle \textstyle F_{1}=F_{0}+\Delta F}$with ; ${\displaystyle \textstyle \Delta F=A_{0}(p-F_{0}G_{0})\;\mod \;F_{0}}$
• ${\displaystyle \textstyle G_{1}=G_{0}+\Delta G}$with ; ${\displaystyle \textstyle \Delta G=B_{0}(p-F_{0}G_{0})\;\mod \;G_{0}}$
• ${\displaystyle \textstyle A_{1}=A_{0}+\Delta A}$ With ${\displaystyle \textstyle \Delta A=A_{0}(1-A_{0}F_{1}-B_{0}G_{1})\;\mod \;F_{0}}$
• ${\displaystyle \textstyle B_{1}=B_{0}+\Delta B}$ With ${\displaystyle \textstyle \Delta B=B_{0}(1-A_{0}F_{1}-B_{0}G_{1})\;\mod \;G_{0}}$

can be determined for which the coefficients of the error terms and have the magnitude . ${\displaystyle p(x)-F_{1}(x)\,G_{1}(x)}$ ${\displaystyle \textstyle 1-(A_{1}F_{1}+B_{1}G_{1})}$ ${\displaystyle \varepsilon ^{2}}$

Finding suitable separation circles

The key point of the numerical procedure is to find suitable separating circles. Schönhage (1982) suggests estimating the absolute value of the largest zero and choosing three equally distributed points on a circle with a double radius. Together with the coordinate origin, these are then tested as the center of the separating circle. Estimates for the distances between the zeros of the polynomial are determined for each of these test points. If these estimates result in a circular ring around the test point without contained zeros, then this is a candidate for a separating circle. The relative latitude, i.e. H. the ratio of the outer and inner radius determines the minimum necessary accuracy for the numerical integration. The best candidate is chosen according to the criteria of the greatest relative width of the annulus and the most even distribution of the zeros on the inside and outside of the annulus.

An improved construction of the set of test points, which guarantees an even distribution of the zeros, was proposed in Pan (1996,2002). Another way of finding groups of separable zeros is the bisection exclusion method (Weyl, Yakoubsohn).

Better separation circles using Gräffe iteration

The product contains only even powers of x . If you replace x² by x , the resulting polynomial has zeros in the squares of the zeros of p . This is the basis of the Dandelin-Gräffe method for determining zeros. Here, however, it is only used to estimate the amounts of the zeros. At the same time as the zeros, the relative widths of zero-free circular rings are squared. If you repeat this squaring often enough, these circular rings also become visible in the estimates. ${\displaystyle p(x)p(-x)}$ ${\displaystyle p_{1}(x)=(-1)^{n}p({\sqrt {x}})p(-{\sqrt {x}})}$

It is possible to use the circular rings broadened by the Gräffe iteration in order to carry out the initial factorization of with a significantly lower accuracy of the numerical integration than that required. In the extreme case, no numerical integration is required (Malajovich). By means of the specified variant of Newton's method, the initial factorization is improved from to a factorization with a small error . For the factors of p (x) we have and , therefore we have ${\displaystyle p_{1}(x)}$ ${\displaystyle p_{0}(x)=p(x)}$ ${\displaystyle p_{1}(x)}$ ${\displaystyle p_{1}(x)-f_{1}(x)g_{1}(x)}$ ${\displaystyle f_{1}(x^{2})\approx (-1)^{k}f(x)f(-x)}$ ${\displaystyle g_{1}(x^{2})\approx (-1)^{n-k}g(x)g(-x)}$

${\displaystyle {\frac {f_{1}(x^{2})}{p(-x)}}\approx (-1)^{k}{\frac {f(x)}{g(-x)}}}$.

With the method of the Padé approximation for the (formal) power series resulting from the left side , the common factor on the left side can be numerically reduced and thus (approximations for) numerators and denominators of the right side can be determined. ${\displaystyle f(-x)}$

Correspondingly, if the Gräffe iteration is used several times, the factorization must be raised several times.

literature

• Schönhage, Arnold (1982): The fundamental theorem of algebra in terms of computational complexity. Preliminary Report, Math. Inst. Univ. Tübingen (1982), 49 pages. (ps.gz; 289 kB)
• Xavier Gourdon: Combinatoire, Algorithmique et Geometry des Polynomials . Ecole Polytechnique, Paris 1996 ( HTML ).
• Pan, Victor Y. (1996). Optimal and nearly optimal algorithms for approximating polynomial zeros. (PDF file; 2.8 MB) Comput. Math. Appl., 31, 97-138.
• Malajovich, Gregorio; Zubelli, Jorge P. (1997): A fast and stable algorithm for splitting polynomials. (PDF file; 300 kB) Comput. Math Appl., 33, 1-23.
• Pan, Victor (1998). Algorithm for Approximating Complex Polynomial Zeros , notes from an introductory lecture
• Pan, Victor (2002). Univariate Polynomials: Nearly Optimal Algorithms for Numerical Factorization and Root-finding (PDF file; 528 kB)

Web links

• Magma documentation. Real and Complex Fields: Element Operations .