Robert Sedgewick: Algorithmen

[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]


41. Die schnelle Fourier-Transformation



Implementation

Nunmehr verfügen wir über alle Bestandteile für einen Algorithmus vom Typ »Teile und Herrsche« für die Multiplikation zweier Polynome unter Verwendung von nur ungefähr N lg N Operationen. Das allgemeine Schema hat folgendes Aussehen:

Die obige Beschreibung läßt sich unmittelbar in ein Programm übertragen, in dem eine Prozedur benutzt wird, mit der ein Polynom vom Grad N - 1 für die N-ten Einheitswurzeln berechnet werden kann. Leider handelt es sich bei der gesamten Arithmetik in diesem Algorithmus um komplexe Arithmetik, und in Pascal existiert kein Standardtyp complex. Obwohl es möglich ist, einen anwenderdefinierten Typ für die komplexen Zahlen zu verwenden, ist es dann auch notwendig, Prozeduren oder Funktionen für alle Rechenoperationen mit diesen Zahlen zu definieren, und dadurch wird der Algorithmus unnötig kompliziert. Für das folgende Programm wird ein Typ complex vorausgesetzt, für den die offensichtlichen arithmetischen Funktionen definiert sind:

    eval(p,outN,0);
    eval(q,outN,0);
    for i:=0 to outN do r[i]:=p[i]*q[i];
    eval(r,outN,0);
    for i:=1 to N do
      begin t:=r[i]; r[i]:=r[outN+1-i]; r[outN+1-i]:=t end;
    for i:=0 to outN do r[i]:=r[i]/(outN+1);

Für dieses Programm wird angenommen, daß die globale Variable outN auf 2N - 1 gesetzt worden ist, und daß p, q und r Felder mit Indizes von 0 bis 2N - 1 sind, die komplexe Zahlen enthalten. Die beiden zu multiplizierenden Polynome p und q sind vom Grad N - 1, und die anderen Koeffizienten in diesen Feldern werden anfangs auf null gesetzt. Die Prozedur eval ersetzt die Koeffizienten des als ersten Parameter angegebenen Polynoms durch die Werte, die man erhält, wenn das Polynom für die Einheitswurzeln berechnet wird. Der zweite Parameter gibt den Grad des Polynoms an (um eins kleiner als die Anzahl der Koeffizienten und der Einheitswurzeln), der dritte Parameter wird weiter unten beschrieben. Das obige Programm berechnet das Produkt von p und q und speichert das Ergebnis in r.

Nun benötigen wir noch die Implementation von eval. Wie wir bereits festgestellt haben, kann die Implementation rekursiver Programme, in denen Felder auftreten, sehr schwierig sein. Es zeigt sich, daß es für den vorliegenden Algorithmus möglich ist, das für gewöhnlich auftretende Problem der Speicherplatzverwaltung zu umgehen, indem der Speicherplatz auf geschickte Weise wiederholt genutzt wird. Es wäre wünschenswert, über eine rekursive Prozedur zu verfügen, die als Eingabegröße ein N + 1 Koeffizienten enthaltendes zusammenhängendes Feld verwendet und die N + 1 Werte im gleichen Feld zurückgibt. Doch der rekursive Schritt erfordert die Verarbeitung zweier nicht zusammenhängender Felder: der ungeraden und der geraden Koeffizienten. Der Leser kann sich leicht überlegen, daß das perfekte Mischen aus dem vorangegangenen Kapitel genau das ist, was hier benötigt wird. Wir können die ungeraden Koeffizienten in einem zusammenhängenden Teilfeld (der ersten Hälfte) und die geraden Koeffizienten in einem weiteren zusammenhängenden Teilfeld (der zweiten Hälfte) zusammenfassen, indem wir ein »perfektes Entmischen« der Eingabedaten vornehmen, wie dies in Abbildung 41.1 für N = 15 schematisch dargestellt ist.

Abbildung 41.1
Abbildung 41.1 Perfektes Entmischen für die schnelle Fourier-Transformation.

Natürlich werden für die Realisierung die konkreten Werte der komplexen Einheitswurzeln benötigt. Bekanntlich gilt

Formel

diese Werte lassen sich unter Verwendung der üblichen trigonometrischen Funktionen leicht berechnen. Im nachfolgenden Programm wird angenommen, daß das Feld w die (outN + 1)-ten Einheitswurzeln enthält.

Dies führt zu der folgenden Implementation der schnellen Fourier-Transformation:

    procedure eval(var p: poly; N,k: integer);  
      var i,j: integer;
      begin
      if N=1 then
        begin
        t:=p[k]; p1:=p[k+1];
        p[k]:=t+p1; p[k+1]:=t-p1
        end
      else
        begin
        for i:=0 to N div 2 do
          begin
          j:=k+2*i;
          t[i]:=p[j]; t[i+1+N div 2]:=p[j+1]
          end;
        for i:=0 to N do p[k+i]:=t[i];
        eval(p,N div 2,k);
        eval(p,N div 2,k+1+N div 2);
        j:=(outN+1) div (N+1);
        for i:=0 to N div 2 do
          begin
          t:=w[i*j]*p[k+(N div 2)+1+i];
          t[i]:=p[k+i]+t; t[i+(N div 2)+1]:=p[k+i]-t
          end;
        for i:=0 to N do p[k+i]:=t[i]
        end
      end;

Dieses Programm bewirkt die Transformation des Polynoms vom Grad N an Ort und Stelle im Teilfeld p[k..k + N] unter Benutzung des oben beschriebenen rekursiven Verfahrens. (Der Einfachheit halber wird für das Programm vorausgesetzt, daß N + 1 eine Zweierpotenz ist, obwohl es nicht schwer ist, auf diese Forderung zu verzichten.) Falls N = 1 ist, wird die einfache Berechnung für 1 und -1 ausgeführt. Andernfalls bewirkt die Prozedur zuerst ein Mischen, ruft sich dann rekursiv selbst auf, um die beiden Hälften zu transformieren und kombiniert schließlich die Ergebnisse dieser Berechnungen in der oben beschriebenen Weise. Um die benötigten Einheitswurzeln zu erhalten, wählt das Programm aus dem Feld w mit einer durch die Variable i bestimmten Schrittweite Werte aus. Wenn zum Beispiel outN den Wert 15 hat, werden die vierten Einheitswurzeln in w[0], w[4], w[8] und w[12] gefunden. Dadurch entfällt die Notwendigkeit, Einheitswurzeln jedesmal, wenn sie benötigt werden, neu zu berechnen.

Eigenschaft 41.3 Zwei Polynome vom Grad N können mit Hilfe von 2N lg N + O(N) komplexen Multiplikationen multipliziert werden.

Diese Tatsache folgt unmittelbar aus den Eigenschaften 41.1 und 41.2. w.z.b.w.

Wie zu Beginn bereits erwähnt, sind die Anwendungsmöglichkeiten der schnellen Fourier-Transformation weit größer, als hier beschrieben werden kann, und der Algorithmus wurde für viele Gebiete gründlich untersucht und vielfältig genutzt. Trotzdem ist die prinzipielle Arbeitsweise bei komplizierteren Anwendungen die gleiche wie bei dem hier betrachteten Problem der Multiplikation von Polynomen. Die schnelle Fourier-Transformation ist ein klassisches Beispiel für die Entwicklung von Algorithmen auf der Grundlage des »Teile und Herrsche«-Prinzips zur Erzielung von beträchtlichen Einsparungen bezüglich des Rechenaufwands.


[ Inhaltsverzeichnis ] [ vorhergehende Seite ] [ nächste Seite ] [ Stichwort ]