Zaawansowane algorytmy i struktury danych/Wykład 4

Abstrakt

Naturalna metoda dodawania dwóch wielomianów wymaga czasu $\mathrm{\Theta }(n)$, natomiast prosty algorytm mnożenia dwóch wielomianów stopnia $n$ wymaga czasu $\mathrm{\Theta }(n^2)$. W wykładzie tym pokażemy, jak z wykorzystaniem szybkiej transformaty Fouriera (STF), wykonać wszystkie podstawowe operacje na wielomianach w czasie większym niż $\mathrm{\Theta }(n)$ o czynnik polilogarytmiczny. Pokażemy jak dla wielomianów stopnia $n$:

• mnożyć je w czasie $O(n\log n)$,
• obliczać wielomian interpolacyjny w czasie $O(n{\log }^{2}n)$,
• obliczać wartość wielomianu w $n$ punktach w czasie $O(n{\log }^{2}n)$,
• dzielić wielomiany w czasie $O(n{\log }^{3}n)$.

Mnożenie wielomianów w punktach

Niech $A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$ i $B(x)=\sum _{i=0}^{n-1}{b}_i{x}^i$ będzą wielomianami stopnia $n$ nad ciałem $F$. Wielomiany te możemy jednoznaczne reprezentować poprzez ich wartości w $n$ punktach. Następujące twierdzenie zostało sformułowane w ramach wykładu z Metod Numerycznych.

Twierdzenie [Twierdzenie o interpolacji wielomianów]

Dla dowolnego zbioru $n$ par $X=\{(x_0,y_0),(x_1,y_1),\dots ,(x_{n-1},y_{n-1})\}$ takiego, że wszystkie wartości $x_i$ są parami różne, istnieje jedyny wielomian $C(x)$ stopnia $n$ taki, że $C(x_i)=y_i$ dla $i=0,1,\dots ,n-1.$

Niech $X$ będzie ustalonym zbiorem parami różnych punktów $x_0,\dots ,x_{2n-1}\in F$. Dla tego zbioru punktów możemy wyznaczyć zbiory wartości wielomianów:

$X_A=\{(x_0,A(x_0)),(x_1,A(x_1),\dots ,(x_{2n-1},A(x_{2n-1}))\}$

$X_B=\{(x_0,B(x_0)),(x_1,B(x_1),\dots ,(x_{2n-1},B(x_{2n-1}))\}$

Niech $C$ będzie wynikiem mnożenia wielomianów $A$ i $B$, mamy wtedy

$C(x_i)=A(x_i)\cdot B(x_i)$.

Ponieważ stopień wielomianu $C$ jest nie większy niż $2n$ to z Twierdzenia o interpolacji zbiór wartości

$X_{A\times B}=\{(x_0,A(x_0)B(x_0)),(x_1,A(x_1)B(x_1),\dots ,(2x_{n-1},A(x_{2n-1})B(x_{2n-1}))\}$,

jednoznacznie wyznacza wielomian $A\times B$. Mając zbiory $X_A$ i $X_B$ możemy wyznaczyć zbiór $X_C$ w czasie $O(n)$. Procedura ta jest przedstawiona na następującym rysunku:

<flash>file=Zasd_fft1.swf|width=460|height=350</flash>

Jednak aby ostatecznie otrzymać szybszy algorytm niż algorytm naiwny musimy pokazać jak rozwiązać problem obliczania wartości wielomianu w $n$ punktach w czasie szybszym niż $\mathrm{\Theta }(n^2)$. Podobnie musimy umieć obliczać wielomian interpolacyjny dla danego zbioru punktów.

Szybka transformata Fouriera (STF)

Problem obliczania wartości wielomianu w $n$ punktach i problem jego interpolacji rozwiążemy wykorzystując szybką transformatę Fouriera. W poprzednim rozdziale nie zakładaliśmy nic na temat zbioru punktów $X$. Głównym pomysłem w konstrukcji algorytmu STF będzie właśnie wybór odpowiedniego zbioru punktów X tak, aby jak największa ilość wykonywanych obliczeń powtarzała się.

Założymy chcemy obliczyć wartości wielomianu $A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$ oraz $n$ jest parzyste. Jeżeli $n$ jest nieparzyste to dodajemy na początek $A(x)$ jednomian $0x^{n+1}$ co nie zmienia nam wyniku działania algorytmu. Punkty $X_n=\{x_0,x_1,\dots ,x_{n-1}\}$ zdefiniujemy w następujący sposób:

$x_i=e^{\frac{2\pi i}{n}}$.

Dla wielomianu $A(x)$ definiujemy dwa nowe wielomiany $A^{[0]}(x)$ i $A^{[1]}(x)$ poprzez wybranie do nich współczynników $A(x)$ o numerach odpowiednio parzystych i nieparzystych:

$A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+\dots +a_{n-2}{x}^{\frac{n}{2}-1}$,

$A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+\dots +a_{n-1}{x}^{\frac{n}{2}-1}$.

Wielomiany $A^{[0]}(x)$ oraz $A^{[1]}(x)$ są stopnia co najwyżej $\frac{n}{2}$. Co więcej zachodzi:

$A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$      (1)

Widzimy teraz, że problem ewaluacji wielomianu $A(x)$ w punktach ${\omega }_n^0,{\omega }_n^1,\dots ,{\omega }_n^{n-1}$ sprowadza się do:

• ewaluacji wielomianów $A^{[0]}(x)$ i $A^{[1]}(x)$

w punktach

$X^{\prime }=\{x_0^2,x_1^2,\dots ,x_{n-1}^2\}$.

• a następnie obliczenie wartości $A(x)$wyniku zgodnie ze wzorem (1).

Zauważmy, że z definicji punktów $x_i$ mamy:

$x_i^2={\left(e^{\frac{2\pi i}{n}}\right)}^2=e^{\frac{2\pi i}{n/2}}$.

Możemy teraz zauważyć, że zachodzi $x_i^2=x_{i+\frac{n}{2}}^2$, a więc $X^{\prime }=X_{\frac{n}{2}}$. Udało nam się więc zredukować problem rozmiaru $n$ - obliczenia wartości wielomianu $A(x)$ stopnia $n$ w $n$do punktach, do dwóch problemów rozmiaru $\frac{n}{2}$ - obliczenia wartości wielomianów $A^{[0]}(x)$ i $A^{[1]}(x)$ stopnia $\frac{n}{2}$ w $\frac{n}{2}$ punktach. Możemy teraz zastosować tą technikę rekurencyjne otrzymując następujący algorytm.

Algorytm Algorytm Szybkiej Transformaty Fouriera

---

 STF($a=(a_0,\dots ,a_{n-1})$)
 if $n$ nieparzyste then
   dodaj wyraz $a_n$ do $a$
   zwiększ $n$
 if $n=1$ then return a
 ${\omega }_n=e^{\frac{2\pi i}{n}}$
 $\omega =1$
 $a^{[0]}=(a_0,a_2,\dots ,a_{n-2})$
 $a^{[1]}=(a_1,a_3,\dots ,a_{n-1})$
 $y^{[0]}=SFT(a^{[0]})$
 $y^{[1]}=SFT(a^{[1]})$
 for k=0 to $\frac{n}{2}-1$ do
   $y_k=y_k^{[0]}+\omega y_k^{[1]}$
   $y_{k+\frac{n}{2}}=y_k^{[0]}-\omega y_k^{[1]}$
   \omega = \omega \omega_n
 return y

Algorytm ten najpierw oblicza SFT wielomianów $A^{[0]}(x)$ i $A^{[1]}(x)$ a następnie łączy te wyniki w celu wyliczenia SFT dla wielomianu $A(x)$. Przeanalizujmy teraz wykonanie pętli. Zauważmy najpierw, że w $k$'tym kroku pętli mamy $\omega ={\omega }_n^k=e^{\frac{2\pi ik}{n}}=x_k.$. Czyli:

Parser nie mógł rozpoznać (nieznana funkcja „\fight”): {\displaystyle y_k = y_k^{[0]} + x_k y_k^{[1]} = A^{[0]}(e^{\frac{2\pi i k}{n/2}}) + x_k A^{[1]}(e^{\frac{2\pi i k}{n/2}}) = A^{[0]}\left(\left(e^{\frac{2\pi i k}{n}}\fight)^2\right) + x_k A^{[1]}\left(\left(e^{\frac{2\pi i k}{n}}\right)^2\right) = A^{[0]}(x_k^2) + x_k A^{[1]}(x_k^2) = A(x_k), }

oraz

Parser nie mógł rozpoznać (błąd składni): {\displaystyle y_{k+\frac{n}{2}} = y_k^{[0]} - x_k y_k^{[1]} = = A^{[0]}(e^{\frac{2\pi i k}{n/2}}) - e^{\frac{2\pi i k}{n}} A^{[1]}(e^{\frac{2\pi i k}{n/2}}) = A^{[0]}(e^{\frac{2\pi i (k + n/2)}{n/2}}) + e^{\frac{2\pi i k + n/2}{n}} A^{[1]}(e^{\frac{2\pi i (k + n/2)}{n/2}}) = A^{[0]}(x_{k + n/2}^2) + x_{k + n/2}^2 A^{[1]}(x_{k + n/2}^2) = A(x_{k+\frac{n}{2}). }

Gdzie w ostatniej równości skorzystaliśmy ze wzoru wzor_1|(1). Widzimy zatem, że algorytm poprawnie oblicza wartość STF dla wielomianu $A(x)$. Równanie rekurencyjne na czas działania procedury STF wygląda następująco:

$T(n)=2T(\frac{n}{2})+\mathrm{\Theta }(n)=\mathrm{\Theta }(n\log n)$.

Odwrotna transformata Fouriera

Aby zakończyć konstrukcję algorytmu dla szybkiego mnożenia wielomianów pozostaje nam pokazanie jak wykonać obliczyć wielomian interpolujący dla zbioru punktów $X_n$. Obliczenie wykonane w czasie szybkiej transformaty Fouriera możemy przedstawić w postaci macierzowej jako mnożenie macierzy przez wektor $(A(x_0),A(x_1),\dots ,A(x_{n-1}){)}^T=V_n(a_0,a_1,\dots ,a_{n-1}{)}^T$, gdzie $V_n=V(x_0,\dots ,x_{n-1})$ jest macierzą Vandermonde'a zawierającą potęgi $x_j$

$\left(\begin{array}{c}A(x_0)\\ A(x_1)\\ A(x_2)\\ ⋮\\ A(x_{n-1})\end{array}\right)=\left[\begin{array}{cccccc}{x}_0^0& x_0^1& x_0^2& x_0^3& \cdots & x_0^{n-1}\\ x_1^0& x_1^1& x_1^2& x_1^3& \cdots & x_1^{n-1}\\ x_2^0& x_2^1& x_2^2& x_2^3& \cdots & x_2^{n-1}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{n-1}^0& x_{n-1}^1& x_{n-1}^2& x_{n-1}^3& \cdots & x_{n-1}^{n-1}\\ \end{array}\right]\left(\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right)$

Element macierzy $V(x_0,\dots ,x_{n-1})$ dany jest jako

$(V_n{)}_{j,k}=V(x_0,\dots ,x_{n-1}{)}_{j,k}=x_j^k$.

Korzystając z definicji zbioru $X_n$ otrzymujemy

$V(x_0,\dots ,x_{n-1}{)}_{j,k}={\left(e^{\frac{2\pi ij}{n}}\right)}^k=e^{\frac{2\pi ijk}{n}}.$

W celu wykonania operacji odwrotnej do SFT, czyli obliczenia wielomianu interpolacyjnego, musimy wykonać mnożenie $V_n^{-1}(A(x_0),A(x_1),\dots ,A(x_{n-1}){)}^T$.

{{lemat||| Niech macierz $W_n$ będzie zdefiniowana jako

${\left(W_n\right)}_{j,k}=\frac{1}{n}{e}^{\frac{-2\pi ijk}{n}},$

jest macierzą odwrotną do macierzy $V_n$. } Dowód

Pokażemy, że $V_n{W}_n=I$. Rozważmy pozycję $(j,k)$ macierzy $V_n{W}_n$:

${\left(V_n{W}_n\right)}_{j,k}=\sum _{l=0}^{n-1}{\left(V_n\right)}_{j,l}{\left(W_n\right)}_{l,k}=\sum _{l=0}^{n-1}{e}^{\frac{2\pi ijl}{n}}\frac{1}{n}{e}^{\frac{-2\pi ilk}{n}}=\sum _{l=0}^{n-1}\frac{1}{n}{e}^{\frac{2\pi l(j-k)}{n}}=$

Jeżeli $j=k$ to $e^{\frac{2\pi k(j-k)}{n}}=1$ i suma ta jest równa $1$. W przeciwnym przypadku możemy skorzystać ze wzoru na sumę szeregu geometrycznego:

$=\frac{1}{n}\frac{1-e^{\frac{2\pi n(j-k)}{n}}}{1-e^{\frac{2\pi (j-k)}{n}}}=\frac{1}{n}\frac{1-1^{(j-k)}}{1-e^{\frac{2\pi (j-k)}{n}}}=0.$

Czyli rzeczywiście $V_n{W}_n=I$. \qed

Porównując postać macierzy $V_n$ oraz macierzy $W_n$ widzimy, że w celu obliczenia transformaty odwrotnej możemy użyć Algorytmu Szybkiej Transformaty Fouriera, musimy tylko zamienić linijkę ${\omega }_m=e^{\frac{2\pi i}{n}}$ na ${\omega }_m=e^{-\frac{2\pi i}{n}}$ i podzielić otrzymany wynik przez $n$.