Trennkreisverfahren

Das Trennkreisverfahren (englisch splitting circle method) ist eine Methode zum numerischen Faktorisieren von Polynomen in einer Variablen mit komplexen Koeffizienten. Dieses Verfahren wurde 1982 von Arnold Schönhage in dem Artikel The fundamental theorem of algebra in terms of computational complexity (bisher nur im Netz veröffentlicht) vorgeschlagen und 1996 von Xavier Gourdon im Computeralgebrasystem PARI/GP und nachfolgend Magma implementiert. Seit der Mitte der 1990er Jahre wurden u. a. von V. Pan und G. Malajovich Verbesserungen des Algorithmus vorgeschlagen, die jedoch bisher nirgends implementiert wurden.

Durch fortgesetztes Zerlegen eines Polynoms in zwei nichttriviale Faktoren kann letztendlich eine Faktorisierung in Linearfaktoren erreicht werden. Dies ist gleichbedeutend zum Auffinden aller komplexen Nullstellen des Polynoms einschließlich der Angabe ihrer Vielfachheit.

Beim numerischen Rechnen mit einer fixierten endlichen Genauigkeit (s. Gleitkommazahl und Festkommazahl) ist es nicht möglich, zwischen einer mehrfach auftretenden Nullstelle und einer gleichmächtigen Gruppe nahe beieinander liegender Nullstellen zu unterscheiden. In diesem Fall ist das Ergebnis des Verfahrens eine Faktorisierung in

• Linearfaktoren für ausreichend isolierte Nullstelle und
• Faktoren höheren Grades für Gruppen von Nullstellen, die in der gewählten Genauigkeit nicht unterscheidbar sind.

Nach dem verallgemeinerten Satz von Vieta sind die Koeffizienten eines normieren Polynoms

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

bis auf ein Vorzeichen die elementarsymmetrischen Polynome in den Nullstellen ${\displaystyle z_{1},\dots ,z_{n}}$ des Polynoms. Es soll eine Zerlegung von ${\displaystyle p(x)}$ in ein Produkt zweier Faktoren ${\displaystyle p(x)=f(x)g(x)}$ gefunden werden, wobei ${\displaystyle f(x)}$ die ersten ${\displaystyle k}$ Nullstellen von ${\displaystyle p(x)}$ als Nullstellen habe, d. h.

${\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})}$ und

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

Die Koeffizienten von ${\displaystyle f(x)}$ und ${\displaystyle g(x)}$ sind zu bestimmen, ohne als Zwischenschritt die Nullstellen bestimmen zu müssen. Dies ist mittels des Residuenkalküls und einer geeigneten Zerlegung der komplexen Ebene möglich. Eine Art der Zerlegung ist die in das Innere und Äußere eines beliebigen Kreises, der dann Trennkreis genannt wird.

Sei ${\displaystyle K}$ eine beschränkte Teilmenge der komplexen Zahlenebene ${\displaystyle \mathbb {C} }$ mit (stückweise) glatter Randkurve ${\displaystyle C}$. Dann gilt nach dem Residuensatz für jede in ${\displaystyle K}$ holomorphe Funktion ${\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 )}$.

Liegen die Nullstellen ${\displaystyle z_{1},\dots ,z_{k}}$ von ${\displaystyle p(x)}$ im Inneren von ${\displaystyle K}$ und alle anderen Nullstellen außerhalb von ${\displaystyle K}$, liegt insbesondere keine Nullstelle auf der Randkurve ${\displaystyle C}$, so gilt also

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

Die in den Koeffizienten von ${\displaystyle f(x)}$ auftretenden Koeffizienten sind Summen in gemischten Produkten der Nullstellen. Die eben angegebene Residuenformel ist daher nicht direkt anwendbar. Es ist aber möglich, die elementarsymmetrischen Polynome durch "ungemischte" Ausdrücke in den Nullstellen darzustellen. Jedes symmetrische Polynom in ${\displaystyle k}$ komplexen Zahlen ${\displaystyle z_{1},\dots ,z_{k}}$ kann durch einen polynomialen Ausdruck in den elementarsymmetrischen Polynomen dieser ${\displaystyle k}$ Zahlen dargestellt werden. Dies gilt insbesondere für die Potenzsummen ${\displaystyle s_{m}=z_{1}^{m}+z_{2}^{m}+\dots +z_{k}^{m}}$. Umgekehrt können die elementarsymmetrischen Polynome und damit die Koeffizienten des Polynoms ${\displaystyle f(x)=(x-z_{1})(x-z_{2})\cdot (x-z_{k})}$ aus den ersten ${\displaystyle k}$ Potenzsummen gewonnen werden, die Umrechnungsformeln dafür sind die Newtonidentitäten. Die Potenzsummen selbst können mittels der Residuenformel zu ${\displaystyle h(z)=z^{m},m=1,\ldots ,k}$ durch ein Konturintegral gewonnen werden.

Der theoretische Faktorisierungsalgorithmus lautet also:

1. Finde eine glatt berandete beschränkte Teilmenge ${\displaystyle K}$, die einige, aber nicht alle Nullstellen von ${\displaystyle p(x)}$ enthält.
2. Bestimme die Konturintegrale, welche die Potenzsummen der Nullstellen ergeben. Mit dem konstanten Polynom ${\displaystyle hh(z)=1}$ kann auch die Anzahl der enthaltenen Nullstellen bestimmt werden.
3. Bestimme mittels der Newton-Identitäten die Koeffizienten von ${\displaystyle f(x)}$, mittels Polynomdivision die Koeffizienten von ${\displaystyle g(x)}$.

In der numerischen Anwendung können die Konturintegrale nicht exakt bestimmt werden. Jedoch kann die numerische Integration mit beliebiger Genauigkeit vorgenommen werden, indem eine genügend kleine Schrittweite gewählt wird. Die mittels der Newton-Identitäten bestimmten genäherten Faktoren seien mit ${\displaystyle F_{0}(x)}$ und ${\displaystyle G_{0}(x)}$ bezeichnet. Für eine schnelle Ausführung der numerischen Integration bietet es sich an, sich auf Kreise in der komplexen Ebene zu beschränken, da dann die numerische Integration, d. h. die Bestimmung der Werte der Polynome an den Stützstellen sowie die Bestimmung der Integrale aus den Werten der Quotienten, mit Hilfe der schnellen Fourier-Transformation ausgeführt werden kann.

Bei genügend hoher Genauigkeit der numerischen Integration werden die Koeffizienten des "Fehlerpolynoms" ${\displaystyle \textstyle p(x)-F_{0}(x)\,G_{0}(x)}$ beliebig klein sein. Ist dieser Fehler von der Größenordnung ${\displaystyle \varepsilon >0}$, so hat der Abstand der Nullstellen von p(x) zu den entsprechenden Nullstellen der Faktoren im ungünstigsten Fall die Größenordnung ${\displaystyle \textstyle {\sqrt[{n}]{\varepsilon }}}$. Die numerische Integration muss so ausgeführt werden, dass mit den Nullstellen von p(x) auch die entsprechenden Nullstellen von ${\displaystyle \textstyle F_{0}(x)}$ innerhalb von K und die von ${\displaystyle \textstyle G_{0}(x)}$ außerhalb von K liegen.

Ist die letzte Forderung erfüllt, so kann die Faktorisierung mittels einer Variante des Newtonverfahrens verbessert werden. Es folgt aus der letztgenannten Forderung, dass sowohl f(x) und g(x) als auch ${\displaystyle F_{0}(x)}$ und ${\displaystyle G_{0}(x)}$ teilerfremd sind. Es gibt nach dem erweiterten euklidischen Algorithmus Polynome a(x) und b(x) mit 1=af+bg. Seien ${\displaystyle A_{0}(x),B_{0}(x)}$ Polynome, für welche die Koeffizienten des Fehlerausdrucks ${\displaystyle 1-(A_{0}F_{0}+B_{0}G_{0})}$ ebenfalls die Größenordnung ${\displaystyle \varepsilon >0}$ haben. Dann können verbesserte Polynome

• ${\displaystyle \textstyle F_{1}=F_{0}+\Delta F}$ mit ${\displaystyle \textstyle \Delta F=A_{0}(p-F_{0}G_{0})\;\mod \;F_{0}}$;
• ${\displaystyle \textstyle G_{1}=G_{0}+\Delta G}$ mit ${\displaystyle \textstyle \Delta G=B_{0}(p-F_{0}G_{0})\;\mod \;G_{0}}$;
• ${\displaystyle \textstyle A_{1}=A_{0}+\Delta A}$ mit ${\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}$ mit ${\displaystyle \textstyle \Delta B=B_{0}(1-A_{0}F_{1}-B_{0}G_{1})\;\mod \;G_{0}}$

bestimmt werden, für welche die Koeffizienten der Fehlerausdrücke ${\displaystyle p(x)-F_{1}(x)\,G_{1}(x)}$ und ${\displaystyle \textstyle 1-(A_{1}F_{1}+B_{1}G_{1})}$ die Größenordnung ${\displaystyle \varepsilon ^{2}}$ besitzen.

Der Kernpunkt des numerischen Verfahrens besteht im Auffinden geeigneter Trennkreise. Schönhage (1982) schlägt vor, den Betrag der größten Nullstelle zu schätzen und auf einen Kreis doppelten Radius drei gleichverteilte Punkte zu wählen. zusammen mit dem Koordinatenursprung werden diese dann als Mittelpunkt des Trennkreises getestet. Zu jedem dieser Testpunkte werden Schätzungen für die Abstände der Nullstellen des Polynoms bestimmt. Ergibt sich aus diesen Schätzungen ein Kreisring um den Testpunkt ohne enthaltene Nullstellen, so ist dies ein Kandidat für einen Trennkreis. Die relative Breite, d. h. das Verhältnis aus äußerem und inneren Radius, bestimmt die minimal notwendige Genauigkeit bei der numerischen Integration. Man wählt den besten Kandidaten nach den Kriterien der größten relativen Breite des Kreisrings und der gleichmäßigsten Aufteilung der Nullstellen auf das Innere und Äußere des Kreisrings.

Eine verbesserte Konstruktion der Menge der Testpunkte, die eine gleichmäßige Aufteilung der Nullstellen garantiert, wurde in Pan (1996,2002) vorgeschlagen. Eine weitere Variante, Gruppen trennbarer Nullstellen aufzufinden, sind Bisektions-Exklusionsverfahren (Weyl, Yakoubsohn).

Das Produkt ${\displaystyle p(x)p(-x)}$ enthält nur gerade Potenzen von x. Ersetzt man darin x² durch x, so hat das entstehende Polynom ${\displaystyle p_{1}(x)=(-1)^{n}p({\sqrt {x}})p(-{\sqrt {x}})}$ Nullstellen in den Quadraten der Nullstellen von p. Dies ist die Grundlage des Dandelin-Gräffe-Verfahrens zur Nullstellenbestimmung. Hier wird es jedoch nur zur Schätzung der Beträge der Nullstellen verwendet. Gleichzeitig mit den Nullstellen werden auch die relativen Breiten nullstellenfreier Kreisringe quadriert. Wiederholt man dieses Quadrieren oft genug, so werden diese Kreisringe auch in den Schätzungen sichtbar.

Es ist möglich, die durch die Gräffe-Iteration verbreiterten Kreisringe zu benutzen, um die Anfangsfaktorisierung von ${\displaystyle p_{1}(x)}$ mit einer wesentlich geringeren Genauigkeit der numerischen Integration als der für ${\displaystyle p_{0}(x)=p(x)}$ notwendigen durchzuführen. Im Extremfall ist keine numerische Integration erforderlich (Malajovich). Mittels der angegebenen Variante des Newton-Verfahrens wird die Anfangsfaktorisierung von ${\displaystyle p_{1}(x)}$ zu einer Faktorisierung mit kleinem Fehler ${\displaystyle p_{1}(x)-f_{1}(x)g_{1}(x)}$ verbessert. Für die Faktoren von p(x) gilt ${\displaystyle f_{1}(x^{2})\approx (-1)^{k}f(x)f(-x)}$ und ${\displaystyle g_{1}(x^{2})\approx (-1)^{n-k}g(x)g(-x)}$, daher gilt

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

Mit der Methode der Padé-Approximation für die aus der linken Seite entstehende (formale) Potenzreihe kann der gemeinsame Faktor ${\displaystyle f(-x)}$ auf der linken Seite numerisch gekürzt werden und damit (Approximationen für) Zähler und Nenner der rechten Seite bestimmt werden.

Entsprechend muss, wenn die Gräffe-Iteration mehrfach angewandt wird, die Hebung der Faktorisierung mehrfach ausgeführt werden.

• 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 Geometrie des Polynomes. Ecole Polytechnique, Paris 1996 (HTML). 
• Pan, Victor Y. (1996). Optimal and nearly optimal algorithms for approximating polynomial zeros. (PDF; 2,8 MB) Comput. Math. Appl., 31, 97–138.
• Malajovich, Gregorio; Zubelli, Jorge P. (1997): A fast and stable algorithm for splitting polynomials. (PDF; 300 kB) Comput. Math Appl., 33, 1–23.
• Pan, Victor (1998). Algorithm for Approximating Complex Polynomial Zeros, Notizen eines Einführungsvortrags
• Pan, Victor (2002). Univariate Polynomials: Nearly Optimal Algorithms for Numerical Factorization and Root-finding (PDF; 528 kB)

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