Modulare Polynomarithmetik: Evaluation und Interpolation

Zwei Polynome NiMtJSJmRzYjJSJYRw== und NiMtJSJnRzYjJSJYRw== vom Grad <= 3 sollen miteinander multipliziert werden.

Das Produkt

NiMvLSUiaEc2IyUiWEcqJi0lImZHRiYiIiItJSJnR0YmRis= ist ein Polynom vom Grad <=7.

f := X -> 5*X^3+3*X^2-4*X+3;

g := X -> 2*X^3-5*X^2+7*X-2;

h := unapply(expand(f(X)*g(X)),X);

Um NiMtJSJoRzYjJSJYRw== nach dem Schema von Evaluation und Interpolation zu berechnen, gen\374gen also 8 Interpolationsstellen, die im Prinzip beliebig gew\344hlt werden k\366nnen

Die Interpolationspunkte:

points := [seq(X,X=0..7)];

Die Werte der Polynome NiMtJSJmRzYjJSJYRw== und NiMtJSJnRzYjJSJYRw== an den Interpolationspunkten:

valuesf := map(f,points);

valuesg := map(g,points);

Die Werte des Polynoms NiMtJSJoRzYjJSJYRw== sind die Produkte der entsprechenden Werte von NiMtJSJmRzYjJSJYRw== und NiMtJSJnRzYjJSJYRw==

zip((x,y)-> x*y,valuesf,valuesg);

Durch Interpolation erh\344lt man das gesuchte Polynom:

interp(points,%,X);

Man h\344tte die Interpolationspunkte auch zuf\344llig w\344hlen k\366nnen.

Beispielsweise mit Hilfe eines Zufallsgenerators.

Man muss sich dann nur vergewissern, dass die Punkte paarweise verschieden sind.

die := rand(-99..99);

randpoints := [seq(die(),X=0..7)];

valuesf := map(f,randpoints);

valuesg := map(g,randpoints);

zip((x,y)-> x*y,valuesf,valuesg);

interp(randpoints,%,X);

Bei der Diskreten Fouriertransformation der Ordnung n w\344hlt man als Interpolationspunkte die sogenannten n-ten Einheitswurzeln in der komplexen Ebene. Das sind die n L\366sungen der Gleichung NiMvKSUiWEclIm5HIiIi,

diese liegen auf dem Einheitskreis; es handelt sich um die n (verschiedenen) Potenzen der

sogenannten "primitiven" n-ten Einheitwurzel NiMpSSJlRzYiKioiIiMiIiJJI3BpR0YlRihJIklHRiVGKEkibkdGJSEiIg==.

Im Falle des Beispiels arbeitet man mit $n=8$.

dftpoints := [solve(X^8=1)];

Die Liste der 8 Punkte ist nicht in fortlaufenen Potenzen geordnet. Um das explizit zu machen:

omega8 := (1+I)/sqrt(2);

dftpoints := [seq(evalc(omega8^k),k=0..7)];

Die Polynome NiMtJSJmRzYjJSJYRw== und NiMtJSJnRzYjJSJYRw== und des Produktes NiMtJSJoRzYjJSJYRw== haben an diesen Stellen komplexe Werte.

valuesf := map(evalc,map(f,dftpoints));

valuesg := map(evalc,map(g,dftpoints));

valuesh := map(evalc,map(h,dftpoints));

Die Werte von NiMtJSJoRzYjJSJYRw== erh\344lt man auch durch Bildung der Produkte der entsprechenden Werte von NiMtJSJmRzYjJSJYRw== und NiMtJSJnRzYjJSJYRw==.

zip((x,y)-> evalc(x*y),valuesf,valuesg);

Die Berechnung (der Koeffizienten) von NiMtJSJoRzYjJSJYRw== kann man in diesem Fall statt durch Interpolation viel einfacher durch Anwendung der inversen Fouriertransformation machen. Dazu interpretiert man die Folge der Funktionswerte als Koeffizienten eines Polynoms NiMtJSJIRzYjJSJYRw==. (Jetzt kommt es auf die richtige Reihenfolge an!)

convert(zip((x,y)->x*y,valuesh,[seq(X^k,k=0..7)]),`+`);

H := unapply(%,X);

Die Auswertung geschieht wiederum an den achten Einheitswurzeln, nun aber im umgekehrten Durchlaufungssinn:

inversedftpoints :=[seq(evalc(omega8^(-k)),k=0..7)];

valuesH := map(expand,map(evalc,map(H,inversedftpoints)));

Nach Divisions durch den Skalierungsfaktor 8 erh\344lt man in der Tat die Koeffizienten des Polynoms NiMtJSJoRzYjJSJYRw==.

map(x->1/8*x,%);

Im Beispielfall n=8 kann man exakt rechnen, da NiMlJ29tZWdhOEc= mittels Wurzeln exakt ausgedr\374ckt werden kann. Das ist im allgemeinen nicht mehr der Fall und man muss die Rechnungen (sofern man nicht in endliche K\366rper und Ringe ausweicht) mit Fliesskomma-Arithmetik, also mit Rundungsfehlern behaftet, durchf\374hren.

Zur Illustration folgt das gleiche Beispiel mit der in Maple eingebauten FFT-Funktion.

Komplexe FFT in Maple

Liste der Koeffizienten des Polynoms NiMtJSJmRzYjJSJYRw==:

fcoeffs := array([seq(coeff(f(X),X,k),k=0..7)]);

Die Koeffizienten m\374ssen in Real- und Imagin\344rteile zerlegt werden:

reell := fcoeffs;

imag := array([0$8]);

Durchf\374hrung der DFT der Ordnung NiMvIiIpKiQpIiIjIiIkIiIi mittels FFT:

FFT(3,reell,imag);

Man erh\344lt im Frequenzbereich die Liste der Werte auch nach Real- und Imagin\344rteil getrennt:

Freell := copy(reell): print(reell);

Fimag := copy(imag): print(imag);

Mittels inverser Fouriertransformation sollten sich die Koeffizienten von NiMtJSJmRzYjJSJYRw== wiederherstellen lassen.

Das gelingt nur bis auf Rundungsfehler!

iFFT(3,reell,imag);

op(reell);

op(imag);

FFT der Koeffizientenliste von NiMtJSJnRzYjJSJYRw==:

gcoeffs := array([seq(coeff(g(X),X,k),k=0..7)]);

reell := gcoeffs;

imag := array([0$8]);

FFT(3,reell,imag);

Greell := reell: print(reell);

Gimag := imag: print(imag);

Punktweise Multiplikation im Frequenzbereich:

Fcoeffs := zip((x,y) -> x+I*y,Freell,Fimag);

Gcoeffs := zip((x,y) -> x+I*y,Greell,Gimag);

Hcoeffs := zip((x,y)-> evalc(x*y),Fcoeffs,Gcoeffs);

Vergleich mit der Fouriertransformierten (in floats) NiMtJSJIRzYjJSJYRw== der Produktpolynoms NiMtJSJoRzYjJSJYRw==:

evalf(H(X));

Hreal := map(x->Re(x),Hcoeffs);

Himag := map(x->Im(x),Hcoeffs);

R\374cktransformation mittels der Werte von NiMtJSJIRzYjJSJYRw==:

iFFT(3,Hreal,Himag);

op(Hreal);

op(Himag);

Vergleich mit dem Produktpolynom NiMtJSJoRzYjJSJYRw==:

h(X);