Kapitel 2: Algebraische Algorithmen (Effiziente Algorithmen, WS 2019) Gerhard Woeginger

Kapitel 2: Algebraische Algorithmen

(Effiziente Algorithmen, WS 2019) Gerhard Woeginger

WS 2019, RWTH

Algebraische Algorithmen

Multiplikation von ganzen Zahlen Multiplikation von Matrizen Berechnung von inversen Matrizen Multiplikation von Polynomen

Zum Aufw¨ armen:

Multiplikation von ganzen Zahlen

Multiplikation von ganzen Zahlen: Schulmethode

In der Grundschule haben wir gelernt,

wie man zwei positive ganze Zahlenx undy mit einander multipliziert:

1 1 9 5 8 3 3 ∗ 2 8 5 1

2 3 9 1 6 6 6

9 5 6 6 6 6 4

5 9 7 9 1 6 5

1 1 9 5 8 3 3

3 4 0 9 3 1 9 8 8 3

Diese Schulmethode verwendetΘ(n²)Operationen.

(Wir betrachten hier die Bit-Komplexit¨at von Algorithmen.) Frage

Geht das irgendwie besser? Schneller?

Multiplikation von ganzen Zahlen: Eine erste Idee

Wir teilen die Ziffernfolgen der beidenn-stelligen Zahlen

x=10^n/2x1+x0 undy=10^n/2y1+y0in zwei gleich lange Teile auf.

x: x1 x0

y: y1 y0

Divide and Conquer

Wir berechnen rekursiv die vier Produktex₀y₀,x₀y₁, x₁y₀, x₁y₁. Wir geben10ⁿx1y1+ (x0y1+x1y0)10^n/2+x0y0aus.

Ergo: T(n) =4T(n/2) + Θ(n)

Multiplikation von ganzen Zahlen: Die zweite Idee

Idee von Anatoly Alexeevitch Karatsuba (1962) Statt den vier Produktenx0y0,x0y1,x1y0, x1y1,

berechnen wir nur die drei Produkte x0y0, x1y1, und (x0+x1)(y0+y1) Damit k¨onnen wir den gemischten Term in der Form

x0y1+x1y0 = (x0+x1)(y0+y1)−x0y0−x1y1

schreiben, und das gew¨unschte Produkt wie folgt berechnen:

xy = 10ⁿx1y1+ ((x0+x1)(y0+y1)−x0y0−x1y1)10^n/2+x0y0

Multiplikation von ganzen Zahlen: Resultat

Ergo:T(n) =3T(n/2) + Θ(n)

Satz (Anatoly Karatsuba, 1962)

Das Produkt von zwein-stelligen Zahlen

kann mitΘ(n^log23)vielen Bit-Operationen berechnet werden.

Anmerkung: log₂3≈1.585

Anmerkungen

Satz (David Harvey & Joris van der Hoeven, M¨arz 2019) Das Produkt von zwein-stelligen Zahlen

kann mitΘ(nlogn)vielen Bit-Operationen berechnet werden.

Spektrum der Wissenschaft (April 2019) https://www.spektrum.de/news/

die-schnellste-art-zu-multiplizieren/1638472

Matrix-Multiplikation

Matrix-Multiplikation: Schulmethode

Problem: Matrix-Multiplikation

Eingabe:Zwein×nMatrizenAund B Gesucht:Das ProduktC :=AB

Trivialer Algorithmus: Berechne EintragCij alsPn

k=1aikbkj

Satz

Das Produkt von zwein×nMatrizen

kann in kubischer ZeitO(n³)berechnet werden.

Anmerkung: Wir z¨ahlen die Anzahl der Multiplikationen.

Frage

Geht das irgendwie besser? Schneller?

Multiplikation von 2 × 2 Matrizen (1)

Mit acht Multiplikationen und vier Additionen:

a b c d

! w x y z

!

= aw+by ax+bz cw+dy cx+dz

!

Multiplikation von 2 × 2 Matrizen (2)

Mit sieben Multiplikationen und achtzehn Additionen:

m1 = (b−d)(y+z) m2 = (a+d)(w +z) m3 = (a−c)(w+x) m₄ = (a+b)z m₅ = a(x−z) m6 = d(y−w) m7 = (c+d)w

aw+by ax+bz cw+dy cx+dz

!

= m1+m2−m4+m6 m4+m5

m +m m −m +m −m

!

Multiplikation von 2 × 2 Matrizen (3)

Mit sechs Multiplikationen und 8.128 Additionen:

Nein, das geht nicht.

Bei sieben Multiplikationen ist bereits Schluss.

Shmuel Winograd hat 1971 gezeigt, dass man zur Multiplikation von 2×2Matrizen mindestens sieben Multiplikationen braucht (auch wenn man noch so viele Additionen und Subtraktionen zur Verf¨ugung hat).

Schnelle Matrix-Multiplikation nach Strassen (1)

Idee von Volker Strassen (1969)

Konstruiere einen Divide-and-Conquer Algorithmus f¨ur die Multiplikation vonn×nMatrizen.

Verwende die sieben Multiplikationen und achtzehn Additionen (aus dem Multiplikationsschema f¨ur2×2Matrizen) im Conquer Schritt.

Schnelle Matrix-Multiplikation nach Strassen (2)

Im Divide Schritt unterteilen wir die beiden Matrizenn×nMatrizen Aund B jeweils in viern/2×n/2Matrizen:

A= A₁₁ A₁₂ A21 A22

!

B = B₁₁ B₁₂ B21 B22

!

Im Conquer Schritt berechnen wir durch zehn Matrix-Additionen und durch sieben rekursive Aufrufe die sieben Matrix-Produkte

M₁, . . . ,M₇f¨ur diesen/2×n/2MatrizenA₁₁,A₂₁, . . . ,B₂₂. Wir f¨ugen die Ergebnisse durch acht Matrix-Additionen zur Ergebnismatrix zusammen:

M1+M2−M4+M6 M4+M5

M₆+M₇ M₂−M₃+M₅−M₇

!

Schnelle Matrix-Multiplikation nach Strassen (3)

Ergo:T(n) =7T(n/2) + Θ(n²)

Satz (Strassen, 1969)

Das Produkt von zwein×nMatrizen

kann mitO(n^log27)Integer-Multiplikationen berechnet werden.

Anmerkung: log₂7≈2.807

Anmerkungen

Der konstante Faktor (in der O-Notation versteckt) von Strassen’s Algorithmus ist gr¨osser als der konstante Faktor in derO(n³) Schulmethode.

In der Praxis: Schulmethode f¨ur kleinen, Strassen f¨ur grossen.

Der Cross-over Punkt liegt meistens um n=20.

In der Praxis: Wenn Matrizen d¨unn besetzt sind (= mit sehr vielen Nullen), dann gibt es schnellere Spezialalgorithmen aus der Numerik Don Coppersmith und Shmuel Winograd (1990) haben den

Exponenten von Strassen’s2.807auf2.376verbessert

Virginia Vassilevska Williams (2011) hat den Exponenten weiter auf 2.373verbessert

Als untere Schranke f¨ur die Multiplikation von zwein×nMatrizen kennen wir nur die SchrankeΩ(n²).

Verifikation von

Matrix-Multiplikationen

Verifikation von Matrix-Multiplikationen

Problem: Matrix-Multiplikation Verifikation Eingabe:Drein×nMatrizenA, B,C Frage:GiltAB =C?

Trivialer L¨osungsansatz:

Berechne das ProduktAB

Vergleiche dien² Eintr¨age mit den Eintr¨agen inC Frage

Geht das irgendwie besser? Schneller? Einfacher?

Ein probabilistischer Ansatz

Idee von Rusins Martins Freivalds (1977) W¨ahle zuf¨allig einen Vektorx ∈ {0,1}ⁿ. Wenn ABx =Cx, dann return “AB=C”.

WennABx 6=Cx, dann return “AB 6=C”.

Anmerkungen:

Wenn Freivalds Ungleichheit behauptet, dann gilt wirklich AB 6=C. Wenn Freivalds Gleichheit behauptet, dann gilt aber nicht

notwendigerweiseAB =C.

Beispiel: Falls der zuf¨allig gew¨ahlte Vektorx=0ist, dann behauptet Freivalds Gleichheit v¨ollig unabh¨angig vonA, B,C

Fehleranalyse (1)

Satz

WennD einen×nMatrix mitD 6=0ist undx ein zuf¨allig gew¨ahlter Vektor in {0,1}ⁿ,

dann ist die Wahrscheinlichkeit vonDx =0h¨ochstens1/2.

Beweisskizze:

W¨ahle Indizesk und`, sodassd<sub>k`</sub>6=0in MatrixD gilt.

Setzey :=Dx und betrachte die k-te Komponenteyk iny. Dann giltyk = dk1x1+dk2x2+· · ·+dknxn = dk`x`+R Zufallsexperiment: Wir setzen der Reihe nach die Komponenten xi

miti 6=`zuf¨allig und unabh¨angig von einander auf0/1, und erhalten so die entsprechenden RestsummeR.

Dann ist mindestens einer der beiden Werte R+d_{k`</sub> mit x<sub>`} = +1 und R mitx_`=0ungleich0.

Daher gilt y_k 6=0mit Wahrscheinlichkeit mindestens1/2.

Fehleranalyse (2)

Satz

Wenn der Freivalds Algorithmus behauptet,

dassAB 6=C, dann stimmt das auf jeden Fall;

dassAB =C, dann stimmt das mit Wahrscheinlichkeit≥1/2.

Beweis:

Wende Satz von vorhergehender Seite auf MatrixD :=C−AB an

Wiederholt man den Freivalds Algorithmus100-mal (unabh¨angig) so sinkt die Fehlerwahrscheinlichkeit von1/2auf2⁻¹⁰⁰≈10⁻³⁰

Laufzeit

Satz

Der Freivalds Algorithmus kann inO(n²)Zeit implementiert werden.

Beweis: Das ProduktABx berechnet man alsA(Bx).

Ubung ¨

Ein viel einfacherer Ansatz (“Ingenieursmethode”) W¨ahle zuf¨allig Zeilenindexi und Spaltenindexj Berechne Eintrag[AB]ij im ProduktAB Wenn [AB]ij =Cij, dann return “AB =C”.

Wenn[AB]ij 6=Cij, dann return “AB 6=C”.

Lineare LaufzeitO(n)

Frage: Wie gross ist die Fehlerwahrscheinlichkeit?

(Extremfall: Nur ein einziger Eintrag inC ist falsch)

Frage: Wie oft muss man diesen einfachen Algorithmus wiederholen, damit die Fehlerwahrscheinlichkeit unter 1/2sinkt?

Berechnung von inversen Matrizen

Inverse Matrizen

Definition

Dieinverse MatrixA⁻¹ zu einer quadratischenn×nMatrixA erf¨ullt die Gleichungen A⁻¹A=AA⁻¹=I.

Eine MatrixAmitdet(A) =0heisstsingul¨ar und besitzt keineinverse Matrix.

Eine MatrixAmitdet(A)6=0heisstregul¨arund besitzt eine eindeutigeinverse Matrix.

F¨ur jede regul¨are MatrixAgilt: det(A)·det(A⁻¹) =1

Das Gauss-Jordan Eliminationsverfahren (Jiu Zhang Suanshu, 100 v.Chr.) berechnet die inverse Matrix inO(n³)Zeit.

Frage

Inversion → Multiplikation

F¨ur zwein×nMatrizenAundB betrachten wir die3n×3nMatrix M =







In A 0 0 In B 0 0 In





 mit M⁻¹=







In −A AB 0 In −B 0 0 In





.

Satz

Wenn die inverse Matrix einern×nMatrix inO(n^α)Zeit berechnet werden kann (mit fixemα≥2),

dann kann auch das Produkt zweier beliebiger n×nMatrizen inO(n^α)Zeit berechnet werden.

Multiplikation → Inversion

Unser n¨achstes Ziel ist es nun,

aus einemO(n^α)Zeit Algorithmus f¨ur Matrix-Multiplikation einenO(n^α)Zeit Algorithmus f¨ur Matrix-Inversion zu bauen.

Satz

Wenn das Produkt zweier beliebigern×nMatrizen inO(n^α)Zeit berechnet werden kann (mit fixemα≥2),

dann kann auch die inverse Matrix von einern×nMatrix inO(n^α)Zeit berechnet werden.

Gegeben sei also eine regul¨aren×nMatrixA.

Wir wollen die inversen×nMatrixA⁻¹f¨urAfinden.

Multiplikation → Inversion: Erste Vereinfachung

Erste Vereinfachung

Wir k¨onnen/d¨urfen/werden annehmen,

dass die Dimension n=2^q der Eingabematrix eine Zweierpotenzist.

A 0 0 I_k

!−1

= A⁻¹ 0 0 I_k

!

Multiplikation → Inversion: Zweite Vereinfachung

Zweite Vereinfachung

Wir k¨onnen/d¨urfen/werden annehmen,

dass die EingabematrixAsymmetrischundpositiv-definitist.

1. Berechne Hilfsmatrix B=A^TA.

Diese HilfsmatrixB ist symmetrisch und positiv-definit.

2. BerechneB⁻¹= (A^TA)⁻¹.

3. Dann giltA⁻¹= (A^TA)⁻¹A^T =B⁻¹A^T.

Schritt 1 kostet eine Matrix-Transposition und eine Matrix-Multiplikation.

Schritt 2 ist die Inversion einer symmetrischen + positiv-definiten Matrix.

Schritt 3 kostet eine Matrix-Multiplikation.

Multiplikation → Inversion: Algorithmus (1)

Wir verwenden Divide & Conquer und unterteilen dien×nMatrixA in vier ⁿ₂×ⁿ₂ Untermatrizen:

A = B C^T C D

!

Dann gilt:

A⁻¹ = B⁻¹+B⁻¹C^TS⁻¹CB⁻¹ −B⁻¹C^TS⁻¹

−S⁻¹CB⁻¹ S⁻¹

!

wobei die MatrixS dasSchur Komplement vonB inC ist:

S = D−CB⁻¹C^T

Die MatrizenB undS sind symmetrisch und positiv-definit.

Multiplikation → Inversion: Algorithmus (2)

1. Bestimme die Untermatrizen B, C, C^T,D.

2. BerechneB⁻¹ rekursiv.

3. Berechne das ProduktW =CB⁻¹und transponiereW^T =B⁻¹C^T. 4. Berechne das ProduktX =WC^T und die Differenz S=D−X.

(Comment: Es giltX =CB⁻¹C^T undS=D−CB⁻¹C^T) 5. BerechneS⁻¹rekursiv.

6. Berechne das Produkt Y =S⁻¹W und transponiereY^T =W^TS⁻¹. (Comment: Y =S⁻¹CB⁻¹undY^T =B⁻¹C^TS⁻¹)

7. Berechne das ProduktZ =W^TY. (Comment: Z=B⁻¹C^TS⁻¹CB⁻¹)

Multiplikation → Inversion: Algorithmus (3)

Man pr¨uft leicht nach, dass die gesuchte Inverse wie folgt aussieht:

A⁻¹ = B⁻¹+B⁻¹C^TS⁻¹CB⁻¹ −B⁻¹C^TS⁻¹

−S⁻¹CB⁻¹ S⁻¹

!

= B⁻¹+Z −Y^T

−Y S⁻¹

!

In anderen Worten:

Aus den Hilfsmatrizen, die in den sieben Schritten berechnet werden, l¨asst sich leicht die inverse MatrixA⁻¹ zusammensetzen.

Multiplikation → Inversion: Algorithmus (4)

In den sieben Schritten

wird zweimal eine ⁿ₂×ⁿ₂ Matrix invertiert;

werden vier Produkte von ⁿ₂×ⁿ₂ Matrizen berechnet;

werden verschiedene ⁿ₂×ⁿ₂ Matrizen transponiert, zu einander addiert, von einander subtrahiert.

F¨ur die Zeitkomplexit¨atT(n)des Algorithmus gilt daher T(n) ≤ 2T(n/2) +4·O(n^α) +O(n²)

≤ 2T(n/2) +O(n^α)

≤ O(n^α)

Inverse Matrizen: Zusammenfassung

Satz

Multiplikation vonn×nMatrizen und Inversion vonn×n Matrizen sind zwei exakt gleich schwere Probleme.

Jeder Algorithmus f¨ur das eine Problem ¨ubersetzt sich in einen Algorithmus f¨ur das andere Problem mit der selben asymptotischen Worst Case Zeitkomplexit¨at.

Multiplikation von Polynomen

Polynom-Multiplikation: Schulmethode

In der Mittelschule haben wir gelernt,

wie man zwei PolynomeA(x)undB(x)mit einander multipliziert:

(x³+5x²−2x+1)·(3x³−x²+x+2)

= 3x⁶ +15x⁵ −6x⁴ +3x³

−x⁵ −5x⁴ +2x³ −x² x⁴ +5x³ −2x² +x

+2x³ +10x² −4x +2

= 3x⁶ +14x⁵ −10x⁴ +12x³ +7x² −3x +2 Die Schulmethode verwendetΘ(n²)Operationen,

um zwei Polynome vom Gradnmit einander zu multiplizieren.

Polynom-Multiplikation: Illustration

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

a

b

c = a b + a b + a b + a b + a b + a b + a b

Polynom-Multiplikation: Genaue Problemstellung

A(x) = a_n−1xⁿ⁻¹+a_n−2xⁿ⁻²+· · ·+a₂x²+a₁x+a₀ B(x) = b_n−1xⁿ⁻¹+b_n−2xⁿ⁻²+· · ·+b₂x²+b₁x+b₀ C(x) = A(x)B(x)

= c_2n−2x²ⁿ⁻²+c_2n−3x²ⁿ⁻³+· · ·+c₂x²+c₁x+c₀

Problem: Polynom-Multiplikation

Eingabe:Ganze Zahlena0, . . . ,a_n−1 undb0, . . . ,b_n−1 Gesucht:Zahlenc₀, . . . ,c_2n−2 mitc_k =Pk

i=0a_ib_k−i Anmerkung: Wir nehmen an, dassa_i=b_i=0f¨uri ≥n

Unser Ziel: Polynom-Multiplikation in sub-quadratischer Zeit

Darstellung von Polynomen (1)

Koeffizienten-Darstellung

Das PolynomA(x) =an−1xⁿ⁻¹+· · ·+a1x+a0wird durch die Folgea₀, . . . ,a_n−1 der Koeffizienten spezifiziert.

Punkt-Wert-Darstellung

Das PolynomA(x) =a_n−1xⁿ⁻¹+· · ·+a1x+a0wird durch nPunkte(x0,y0), . . . ,(x_n−1,y_n−1)spezifiziert, wobeiy_i=A(x_i)f¨ur0≤i ≤n−1gilt.

Darstellung von Polynomen (2)

Zur Erinnerung

F¨urnpaarweise verschiedene Zahlen x0, . . . ,x_n−1, und f¨urnPunkte(x0,y0), . . . ,(x_n−1,y_n−1)

existiert genau ein PolynomA(x) =a_n−1xⁿ⁻¹+· · ·+a₁x+a₀ mity_i =A(x_i)f¨ur0≤i ≤n−1.

Beweis f¨ur “es existiert h¨ochstens ein PolynomA(x)”:

Wenn ein Polynom vom Gradn−1mehr alsn−1Nullstellen hat, so ist es identisch gleich 0

Beweis f¨ur “es existiert mindestens ein PolynomA(x)”:

Lagrange Interpolation A(x) =

n−1

X

k=0

yk

Q

j6=k(x−xj) Q

j6=k(xk−xj)

Darstellung von Polynomen (3): ¨ Ubersetzung/Hin

Koeffizienten→Punkt-Wert

Gegeben: Koeffizienten-Darstellunga₀, . . . ,a_n−1eines PolynomsA(x);

npaarweise verschiedene Zahlen x0, . . . ,x_n−1

Gesucht: Punkt-Wert-Darstellung vonA(x)f¨ur St¨utzstellenx0, . . . ,x_n−1

Jeder Wertyi=A(xi)kann inO(n)Zeit berechnet werden:

1 y= a[n-1];

2 for i= n-2 downto 0 do 3 y= x*y + a[i]

4

5 return y;

Darstellung von Polynomen (4): ¨ Ubersetzung/R¨ uck

Punkt-Wert→Koeffizienten

Gegeben: Punkt-Wert-Darstellung(x₀,y₀), . . . ,(x_n−1,y_n−1)f¨urA(x);

Gesucht: Koeffizienten-Darstellunga₀, . . . ,a_n−1vonA(x)

Die Lagrange’sche Interpolationsformel kann inO(n²)Zeit ausgewertet werden

Zusammenfassend:

Koeffizienten-Darstellung und Punkt-Wert-Darstellung k¨onnen in quadratischer ZeitO(n²)in einander ¨ubergef¨uhrt werden

Zur¨ uck zur Polynom-Multiplikation

Unser Hauptziel: Polynom-Multiplikation in sub-quadratischer Zeit Problem: Polynom-Multiplikation (in Koeffizienten-Darstellung) Eingabe:Ganze Zahlena₀, . . . ,a_n−1 undb₀, . . . ,b_n−1

Gesucht:Zahlenc₀, . . . ,c_2n−2 mitc_k =Pk

i=0a_ib_k−i

Nebenproblem: Polynom-Multiplikation in sub-quadratischer Zeit Problem: Polynom-Multiplikation (in Punkt-Wert-Darstellung) Eingabe:Punkt-Wert-Darstellung(x0,y0), . . . ,(x_n−1,y_n−1)f¨urA(x)

Punkt-Wert-Darstellung(x₀,y₀⁰), . . . ,(x_n−1,y_n−1⁰ )f¨urB(x) Gesucht:Punkt-Wert-Darstellung(x0,y₀⁰⁰), . . . ,(x_n−1,y_n−1⁰⁰ )f¨urA(x)B(x) Einfach in linearer ZeitO(n)zu l¨osen

Arbeitsplan

Schritt 1:

Ubersetze die beiden Polynome¨ A(x)undB(x)

von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung Schritt 2:

MultipliziereA(x)undB(x)in Punkt-Wert-Darstellung Schritt 3:

Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung

Intermezzo:

Rechnen mit komplexen Zahlen

Komplexe Zahlen

Eine komplexe Zahlz kann dargestellt werden:

Algebraisch durch Zerlegung in Realteil und Imagin¨arteil:z=a+i b Polar durch Radius r und Winkelφ: z =rcosφ+i rsinφ

Exponentiell: z=r·e^iφ

Rechenoperationen mitz =r ·e^iφ undz⁰=r⁰·e^iφ0 Multiplikation: (r·e^iφ) (r⁰·e^iφ0) = rr⁰·e^i(φ+φ0) Potenzierung: zⁿ = rⁿ·e^{i nφ}

Einheitswurzeln (1)

Dien-ten Einheitswurzelnω⁰_n, ω¹_n, . . . , ωⁿ⁻¹_n

sind dien(komplexen) L¨osungen der Gleichung ωⁿ=1.

1 1

i

i

!₈⁰D!⁸₈

!₈¹

!₈²

!₈³

!⁴₈

!₈⁵

!₈⁶

!₈⁷

1 2πi/n k 2kπi/n

Einheitswurzeln (2)

ωn=e^2πi/nist dien-te Haupt-Einheitswurzel

Alle Einheitswurzeln sind Potenzen der Haupt-Einheitswurzel ωn

Es gelten die ¨ublichen Rechenregelnω^k_nω_n^{`</sup> =ω<sup>k+`}_n

Die n-ten Einheitswurzeln bilden eine multiplikative Gruppe, die zur zyklischen GruppeZn (Restklassengruppe modulon) isomorph ist Wenn nZweierpotenz, dann erzeugt jedes Elementω^k_n mit ungeradem Exponentenk die gesamte Gruppe.

Einheitswurzeln (3)

Halbierungslemma

F¨ur eine Zweierpotenz nfallen die Quadrate dern-ten

Einheitswurzeln mit den(n/2)-ten Einheitswurzeln zusammen.

Summenlemma

Jeden-te Einheitswurzelω_n^k 6=1erf¨ullt die Gleichung (ω^k_n)ⁿ⁻¹+ (ω_n^k)ⁿ⁻²+· · ·+ (ω_n^k)²+ (ω_n^k) +1 = 0

Arbeitsplan: Schritt 1

Schritt 1: Zielsetzung

Zur Erinnerung:

Schritt 1:

Ubersetze die beiden Polynome¨ A(x)undB(x)

von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung

Von jetzt an nehmen wir an, dassn=2^q eine Zweierpotenz ist.

Die Punkt-Wert-Darstellung vonA(x)und B(x)werden wir mit den n-ten Einheitswurzelnhx0, . . . ,x_n−1i=hω⁰_n, . . . , ωⁿ⁻¹_n i als St¨utzstellen bestimmen.

Wir beschreiben das Verfahren nur f¨urA(x).

Das PolynomB(x)wird analog behandelt.

Wir verwenden Divide & Conquer.

Divide & Conquer (1)

Divide & Conquer Ansatz:

Aeven(x) = a0+a2x+a4x²+a6x³+· · ·+a_n−2x^n/2−1 Aodd(x) = a1+a3x+a5x²+a7x³+· · ·+a_n−1x^n/2−1

Koeffizienten-Darstellung:

A_even(x) : ha₀, a₂, a₄, a₆, . . . , a_n−2i Aodd(x) : ha1, a3, a5, a7, . . . , a_n−1i

A(x) = Aeven(x²) +x·Aodd(x²)

Divide & Conquer (2)

Wir beobachten, dass die Zahlen(ω_n⁰)²,(ω_n¹)², . . . ,(ω_nⁿ⁻¹)² exakt mit den(n/2)-ten Einheitswurzeln zusammen fallen

Wir bestimmen rekursiv die Punkt-Wert-Darstellung vonAeven(x) f¨ur dien/2St¨utzstellen(ω⁰_n)²,(ω¹_n)², . . . ,(ωⁿ⁻¹_n )²

Wir bestimmen rekursiv die Punkt-Wert-Darstellung vonA_odd(x) f¨ur dien/2St¨utzstellen(ω⁰_n)²,(ω¹_n)², . . . ,(ωⁿ⁻¹_n )²

Wir berechnen aus diesen beiden Punkt-Wert-Darstellungen die Punkt-Wert-Darstellung vonA(x)f¨ur dienSt¨utzstellen ω_n⁰, . . . , ω_nⁿ⁻¹ mit Hilfe vonA(x) = Aeven(x²) +x·Aodd(x²) Ergo: T(n) =2T(n/2) + Θ(n)

Schritt 1: Zusammenfassung

Ergo: T(n) =2T(n/2) + Θ(n) Ergo: T(n)∈O(nlogn)

Satz

Schritt 1 mit denn-ten Einheitswurzelnhx0, . . . ,x_n−1i=hω⁰_n, . . . , ωⁿ⁻¹_n i als St¨utzstellen kann inO(nlogn)Zeit durchgef¨uhrt werden.

A moment of reflection

Was haben wir in Schritt 1 eigentlich getan? (1)









































1 1 1 1 1 · · · 1

1 ω ω² ω³ ω⁴ · · · ωⁿ⁻¹

1 ω² ω⁴ ω⁶ ω⁸ · · · ω²⁽ⁿ⁻¹⁾

1 ω³ ω⁶ ω⁹ ω¹² · · · ω³⁽ⁿ⁻¹⁾

1 ω⁴ ω⁸ ω¹² ω¹⁶ · · · ω⁴⁽ⁿ⁻¹⁾

... ... ... ... ... ... ...

... ..

. ..

. ..

. ..

. ..

. ..

. 1 ωⁿ⁻¹ ω²⁽ⁿ⁻¹⁾ ω³⁽ⁿ⁻¹⁾ ω⁴⁽ⁿ⁻¹⁾ · · · ω^{(n−1)(n−1)}















































































 a0

a1

a2

a3

a4

... ... an−1









































=







































 y0

y1

y2

y3

y4

... ... yn−1









































Wir haben den Koeffizientenvektoramit einer MatrixV multipliziert,

und als Resultat den Vektory mit den gew¨unschten y-Koordinaten erhalten.

Was haben wir in Schritt 1 eigentlich getan? (2)

Wir haben den Koeffizientenvektora mit einer MatrixV multipliziert, und als Resultat den Vektory erhalten.

Die MatrixV =V(ω)mit Eintrag ω^rs in Zeiler und Spaltes ist eine so-genannteVandermondeMatrix.

Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelωgenannt.

Analog f¨ur Einheitswurzelnα:=ω^k mit ungeradem Exponentenk: Die Vandermonde MatrixV =V(α)hat den Eintragα^rs in Zeile r und Spaltes stehen.

Der Vektory =VawirdDiskrete Fourier Transformierte (DFT)des Vektorsa bez¨uglich der Einheitswurzelαgenannt.

Was haben wir in Schritt 1 eigentlich getan? (3)

Satz

F¨ur jede Zweierpotenznund

f¨ur jede Einheitswurzelα:=ω^k mit ungeradem Exponentenk, kann die DFT f¨ur einenn-dimensionalen Vektorabez¨uglich der Einheitswurzelα inO(nlogn)Zeit berechnet werden.

Arbeitsplan: Schritt 3

Schritt 3: Zielsetzung

Zur Erinnerung:

Schritt 3:

Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung

Unser Startpunkt ist die Punkt-Wert-Darstellung

eines Polynoms C(x)annSt¨utzstellen(x0,y0), . . . ,(xn−1,yn−1), wobeix_k =ω_n^k f¨ur0≤k ≤n−1

Wir suchen die Koeffizienten-Darstellung c0, . . . ,cn−1von C(x)

Illustration (1)









































1 1 1 1 1 · · · 1

1 ω ω² ω³ ω⁴ · · · ωⁿ⁻¹

1 ω² ω⁴ ω⁶ ω⁸ · · · ω²⁽ⁿ⁻¹⁾

1 ω³ ω⁶ ω⁹ ω¹² · · · ω³⁽ⁿ⁻¹⁾

1 ω⁴ ω⁸ ω¹² ω¹⁶ · · · ω⁴⁽ⁿ⁻¹⁾

... ... ... ... ... ... ...

... ... ... ... ... ... ...

1 ωⁿ⁻¹ ω²⁽ⁿ⁻¹⁾ ω³⁽ⁿ⁻¹⁾ ω⁴⁽ⁿ⁻¹⁾ · · · ω^{(n−1)(n−1)}















































































 c0

c1

c2

c3

c4

... ... cn−1









































=







































 y0

y1

y2

y3

y4

... ... yn−1









































Lineares GleichungssystemVc=y (und wir suchen den Vektorc)

Illustration (2)







































 c0

c1

c2

c3

c4

... ... cn−1









































=









































1 1 1 1 1 · · · 1

1 ω ω² ω³ ω⁴ · · · ωⁿ⁻¹

1 ω² ω⁴ ω⁶ ω⁸ · · · ω²⁽ⁿ⁻¹⁾

1 ω³ ω⁶ ω⁹ ω¹² · · · ω³⁽ⁿ⁻¹⁾

1 ω⁴ ω⁸ ω¹² ω¹⁶ · · · ω⁴⁽ⁿ⁻¹⁾

... ... ... ... ... ... ...

... ... ... ... ... ... ...

1 ωⁿ⁻¹ ω²⁽ⁿ⁻¹⁾ ω³⁽ⁿ⁻¹⁾ ω⁴⁽ⁿ⁻¹⁾ · · · ω^{(n−1)(n−1)}









































−1





































 y0

y1

y2

y3

y4

... ... yn−1









































Lineares Gleichungssystemc=V⁻¹y (und wir suchen den Vektorc)

Das lineare Gleichungssystem (1)

Es seiω=e^2πi/n dien-te Haupt-Einheitswurzel.

Es seiV = (vj,k)dien×nMatrix mitvj,k =ω^jk. Es seieny₀, . . . ,y_n−1 die vorgegebenen Funktionswerte

des PolynomsC(x)an den St¨utzstellenω⁰, . . . , ωⁿ⁻¹

Unser Ziel: L¨ose das GleichungssystemVc =y nach dem Vektorc auf

Das lineare Gleichungssystem (2)

Satz

Die inverse MatrixV⁻¹= (wj,k)ist gegeben durchwj,k =ω^−jk/n.

Beweis: Wir verifizieren, dassV⁻¹V =I gilt.

Der Eintrag in derr-ten Zeile unds-ten Spalte vonV⁻¹V ist [V⁻¹V]_r,s =

n−1

X

k=0

(ω^−rk/n) (ω^ks)

= 1 n

n−1

X

k=0

(ω^s−r)^k

F¨urr =s istω^s−r =1, und daherPn−1

k=0(ω^s−r)^k =n F¨urr 6=s istω^s−r 6=1, und daherPn−1

k=0(ω^s−r)^k =0

Das lineare Gleichungssystem (3)

Altes Ziel:

L¨ose das GleichungssystemVc =y nach dem Vektorc auf

Neues Ziel:

Berechne den Vektorc=V⁻¹y, wobei[V⁻¹]j,k =ω^−jk/ngilt.

Neues Ziel, besser formuliert:

Berechne die Diskrete Fourier Transformierte (DFT)d des Vektorsy bez¨uglich der Einheitswurzelω⁻¹.

Der gesuchte Vektor ist dannc =d/n.

Ergo: Schritt 3 kann inO(nlogn)Zeit erledigt werden

Polynom-Multiplikation:

Zusammenfassung

Arbeitsplan

Schritt 1:

Ubersetze die beiden Polynome¨ A(x)undB(x)

von Koeffizienten-Darstellung nach Punkt-Wert-Darstellung Laufzeit:O(nlogn)

Schritt 2:

MultipliziereA(x)undB(x)in Punkt-Wert-Darstellung Laufzeit:O(n)

Schritt 3:

Ubersetze das in Schritt 2 berechnete Produkt¨ A(x)B(x) von Punkt-Wert-Darstellung nach Koeffizienten-Darstellung

Hauptresultat

Satz

Das Produkt zweier Polynome vom Gradn−1in

Koeffizienten-Darstellung kann inO(nlogn)Zeit berechnet werden.

Anmerkungen:

Da das ErgebnispolynomC(x)den Grad2n−2hat, m¨ussen wir ganz am Anfang die Koeffizientenvektoren a0, . . . ,a_n−1und

b0, . . . ,bn−1durch Hinzuf¨ugen von Komponenten mit Wert0auf die Dimension 2n−2bringen.

Durch Hinzuf¨ugen von noch mehr Komponenten mit Wert0machen wir die Dimension zu einer Zweierpotenz

Anmerkungen

Die Diskrete Fourier Transform wurde von James William Cooley und John Wilder Tukey 1965 algorithmisch untersucht.

(“An algorithm for the machine calculation of complex Fourier series”, Mathematics of Computation 11, pp 297–301)

Der Numeriker Carl David Tolm´e Runge verwendete bereits 1903 die Reduktion der n-dimensionalen DFT auf zwei(n/2)-dimensionale DFTs

Der rekursive Algorithmus wurde bereits von Carl Friedrich Gauss im Jahre 1805 benutzt, um die Flugbahnen der Asteroiden Pallas und Juno zu interpolieren. (“Theoria interpolationis methodo nova tractata”, verfasst in Neu-Latein)

Die FFT (Fast Fourier Transformation) wurde als einer der Top 10 Algorithms of the 20th Centurygelistet (Algorithms with the greatest influence on the development and practice of science and

Die Liste der Top 10 Algorithmen

Metropolis Algorithm for Monte Carlo Simplex Method for Linear Programming Krylov Subspace Iteration Methods

The Decompositional Approach to Matrix Computations The Fortran Optimizing Compiler

QR Algorithm for Computing Eigenvalues Quicksort Algorithm for Sorting

Fast Fourier Transform Integer Relation Detection Fast Multipole Method