Zaawansowane algorytmy i struktury danych/Wykład 4

From Studia Informatyczne

Spis treści

Abstrakt

Naturalna metoda dodawania dwóch wielomianów wymaga czasu \Theta(n), natomiast prosty algorytm mnożenia dwóch wielomianów stopnia n wymaga czasu \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ż \Theta(n) o czynnik polilogarytmiczny. W ramach wykładu pokażemy, jak dla wielomianów stopnia n:

  • mnożyć je w czasie O(n \log n),
  • dzielić wielomiany w czasie O(n \log n).

Natomiast jako ćwiczenie pozostanie nam pokazanie, jak wykorzystać te algorytmy do:

  • obliczania wielomianu interpolacyjnego w czasie O(n \log^2 n),
  • obliczania wartości wielomianu w n punktach w czasie O(n \log^2 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ędą 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 (zobacz).

Twierdzenie 1 [Twierdzenie o interpolacji wielomianowej]

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


Niech X będzie ustalonym zbiorem parami różnych punktów x_0, \ldots, 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), \ldots, (x_{2n-1}, A(x_{2n-1}))\},
X_B = \{(x_0, B(x_0)), (x_1, B(x_1), \ldots, (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\cdot B} = \{(x_0, A(x_0)B(x_0)), (x_1, A(x_1)B(x_1), \ldots, (x_{2n-1}, A(x_{2n-1})B(x_{2n-1})) \},

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




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ż \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, że chcemy obliczyć wartości wielomianu A(x) = \sum_{k=0}^{n-1} a_k x^k oraz n jest parzyste. Jeżeli n jest nieparzyste, to dodajemy na początek A(x) jednomian 0x^{n+1} co nie zmienia wyniku działania algorytmu. Punkty X_n = \{x_0,x_1,\ldots,x_{n-1}\} zdefiniujemy w następujący sposób:

x_k = e^{\frac{2\pi ki}{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_2x + a_4x^2 + \ldots + a_{n-2}x^{\frac{n}{2}-1},
A^{[1]}(x) = a_1 + a_3x + a_5x^2 + \ldots + 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 x_0, x_1,\ldots, x_{n-1} sprowadza się do:

  • ewaluacji wielomianów A^{[0]}(x) i A^{[1]}(x) w punktach
X' =\{x_0^2, x_1^2, \ldots, x_{n-1}^2\}.
  • a następnie obliczenie wartości A(x) zgodnie ze wzorem (1).

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

x_k^2 = \left(e^{\frac{2\pi ki}{n}} \right)^2 = e^{\frac{2\pi k i}{n/2}}.

Możemy teraz zauważyć, że zachodzi x_k^2 = x_{k+\frac{n}{2}}^2, a więc X' = X_{\frac{n}{2}}. Udało nam się więc zredukować problem rozmiaru n - obliczenia wartości wielomianu A(x) stopnia n w n 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ę rekurencyjnie, otrzymując następujący algorytm.

Algorytm Szybkiej Transformaty Fouriera


 STF(a = (a_0,\ldots,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, \ldots, a_{n-2})
 a^{[1]} = (a_1, a_3, \ldots, 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 i k}{n}} = x_k.. Czyli:

\begin{array}{r@{}c@{}l} 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}}\right)^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), \end{array}

oraz

\begin{array}{r@{}c@{}l} 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}}). \end{array}

Gdzie w ostatniej równości skorzystaliśmy ze wzoru (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) =2 T(\frac{n}{2}) + \Theta(n) = \Theta(n \log n).

Poniższa animacja przedstawia działanie algorytmu SFT wywołanego dla wielomianu A(x) = x^3 + 2x^2 + 5x + 3. Wartości wielomianu A(x) wyznaczone są z wartości wielomianów A^{[0]}(x) i A^{[1]}(x), a wartości tych wielomianów z wartości wielomianów A^{[00]}(x), A^{[01]}(x), A^{[10]}(x) i A^{[11]}(x).



Odwrotna transformata Fouriera

Aby zakończyć konstrukcję algorytmu dla szybkiego mnożenia wielomianów, pozostaje nam pokazanie, jak 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), \ldots, A(x_{n-1}))^T = V_n (a_0, a_1, \ldots, a_{n-1})^T, gdzie V_n = V(x_0,\ldots,x_{n-1}) jest macierzą Vandermonde'a zawierającą potęgi x_j:

\left( \begin{matrix} A(x_0)\\ A(x_1)\\ A(x_2)\\ \vdots\\ A(x_{n-1}) \end{matrix} \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} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 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{matrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n-1} \end{matrix} \right).

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

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

Korzystając z definicji zbioru X_n, otrzymujemy

V(x_0,\ldots,x_{n-1})_{j,k} = \left(e^{\frac{2\pi i j}{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), \ldots, A(x_{n-1}))^T.

Lemat 2

Macierz W_n 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 i j l}{n}} \frac{1}{n} e^{\frac{-2\pi i l k}{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.

image:End_of_proof.gif

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.

Dzielenie wielomianów

W tej części wykładu skupimy się na problemie dzielenia dwóch wielomianów. Niech A(x) będzie wielomianem stopnia m, a B(x) wielomianem stopnia n. Zakładamy bez straty ogólności, że b_{n-1}\neq 0. W problemie dzielenia wielomianów chcemy obliczyć dwa wielomiany D(x) i R(x) takie, że

A(x) = D(x) B(x)  + R(x),      (2)

oraz stopień wielomianu R(x) jest ostro mniejszy niż n. Wielomian D(x) nazywamy wynikiem dzielenia, a wielomian R(x) to reszta z dzielenia. Pierwszy pomysł jaki się od razu nasuwa, to spróbować policzyć odwrotność wielomianu B(x) i przemnożyć przez tę odwrotność strony tego równania. Niestety wielomiany nie mają odwrotności będących wielomianami. Jednak nie pozostajemy tutaj zupełnie bezradni. Możemy rozszerzyć naszą dziedzinę obliczeń tak aby zagwarantować istnienie pewnych odwrotności.

Obliczenia będziemy wykonywać nad zbiorem szeregów formalnych F[[x]] nad ciałem F, patrz Wykład z matematyki dyskretnej o funkcjach tworzących. Dla części elementów F[[x]] istnieją odwrotności. Elementy te są postaci a + x A(x), gdzie a\neq 0 i A(x) \in F[[x]].


Ćwiczenie

Czy umiesz obliczyć odwrotność dla szeregu 1 + x?

Odwrotność to \sum_{i=0}^{\infty} (-x)^{i}.

Ćwiczenie

Czy umiesz obliczyć odwrotność dla szeregu \sum_{j=0}^{\infty}(j+1)x^j?

Odwrotność to (1-x)^2.


Do wzoru (2) wstawmy x = \frac{1}{z}, otrzymamy wtedy:

A^R(z) = D^R(z) B^R(z)  + z^{m-n-1}R^R(z) = B^R(z)D^R(z) \mod  z^{m-n-1},

gdzie A^R(z) = z^m A(\frac{1}{z}), D^R(z) = z^{m-n} D(\frac{1}{z}), B^R(z) = z^{n} B(\frac{1}{z}) i R^R(z) = z^{n-1} R(\frac{1}{z}), oznaczają wielomiany otrzymane poprzez odwrócenie kolejności współczynników. Z założenia, że b_{n-1}\neq 0 wiemy, że wielomian B^R ma odwrotność nad zbiorem szeregów formalnych. Możemy zapisać więc:

D^R(z) = A^R(x) \left(B^R(z)\right)^{-1} \mod z^{m-n-1}.      (3)

Zauważmy, że w celu wyznaczenia D^R(z) potrzebujemy tylko m-n-1 wyrazów szeregu \left(B^R(z)\right)^{-1}. Wyższe wyrazy i tak znikną w wyniku wykonania mnożenia modulo z^{m-n-1}. Pozostało nam teraz tylko pokazać, jak wyznaczyć odwrotność dla szeregu formalnego. Algorytm ten przedstawiony jest poniżej.


Algorytm obliczania pierwszych m wyrazów odwrotnośći szeregu formalnego


 ODWROTNOŚĆ(A(x) = \sum_{i=0}^{n-1} a_i x^i, m)
 if m=1 then return 1/a_0
 A^{[0]}(x) =ODWROTNOŚĆ(A(x), \lfloor \frac{m}{2} \rfloor)
 return (A^{[0]}(x) -(A(x) A^{[0]}(x) -1)A^{[0]}(x)) \mod x^m

Obliczenie to jest poprawne, ponieważ:

\frac{1}{A(x)} - (A^{[0]}(x) -(A(x) A^{[0]}(x) -1) A^{[0]}(x)) = A(x)\left(\frac{1}{A(x)} - A^{[0]}(x)\right)^2,

a to jest równe wielokrotności x^{m}. Jeżeli wykorzystamy szybkie mnożenie wielomianów do obliczenia (A^{[0]}(x) -(A(x) A^{[0]}(x) -1)A^{[0]}(x)), to złożoność tego algorytmu wynosić będzie O(n\log n). Możemy teraz skonstruować algorytm wykonujący dzielenie wielomianu A(x) przez wielomian B(x) w czasie O(m \log m), gdzie m to stopień wielomianu A(x).

Algorytm dzielenia wielomianów


 PODZIEL(A(x), B(x))
  A^R(z) = z^m A(\frac{1}{z})
  B^R(z) = z^{n} B(\frac{1}{z})
  (B^R(z)^{-1})(z) =ODWROTNOŚĆ(B^R(z), m-n-1)
  D^R(z) = A^R(x) \left(B^R(z)\right)^{-1} \mod z^{m-n-1}
  D(x) = z^{m-n-1} D^R(\frac{1}{x})
  R(x) = A(x) - D(x)B(x)
  return (D(x), R(x))

Poprawność tego algorytmu wynika wprost ze wzoru (3), a jego złożoność wynosi naturalnie O(n \log n).