[ads2.git] / 7-fft / 7-fft.tex

\def\scharfs{\char"19}
\input lecnotes.tex
\prednaska{9}{Fourierova transformace}{}

\h{Násobení polynomù}

V~této kapitole se budeme zabývat násobením polynomù. Tento na první pohled
triviální algebraický problém má pøekvapivì efektivní øe¹ení, které nás
dovede a¾ k~Fourierovì transformaci.

\s{Znaèení:} {\I Polynomy} jsou výrazy typu
$$
P(x)=\sum_{i=0}^{n-1} p_i\cdot x^i,
$$
kde $x$ je~promìnná a $p_0$ a¾ $p_{n-1}$ jsou èísla, kterým øíkáme {\I koeficienty} polynomu.
Obecnì budeme znaèit polynomy velkými písmeny a jejich koeficienty pøíslu¹nými
malými písmeny s~indexy.

V~algoritmech polynomy obvykle reprezentujeme pomocí {\I vektoru koeficientù}
$(p_0,\ldots,p_{n-1})$; oproti zvyklostem lineární algebry budeme slo¾ky vektorù v~celé
této kapitole indexovat od~0. Poètu koeficientù~$n$ budeme øíkat {\I velikost
polynomu} $\vert P\vert$. Èasovou slo¾itost algoritmu budeme vyjadøovat vzhledem k~velikostem
polynomù na jeho vstupu.

Pokud pøidáme nový koeficient $p_n=0$, hodnota polynomu se pro ¾ádné~$x$
nezmìní. Stejnì tak je-li nejvy¹¹í koeficient~$p_{n-1}$ nulový, mù¾eme ho
vynechat. Takto mù¾eme ka¾dý polynom zmen¹it na {\I normální tvar,}
v~nìm¾ má buïto nenulový nejvy¹¹í koeficient, nebo nemá vùbec ¾ádné
koeficienty -- to je takzvaný {\I nulový polynom,} který pro ka¾dé~$x$ roven nule.
Nejvy¹¹í mocninì s~nenulovým koeficientem se øíká {\I stupeò polynomu} $\deg R$,
nulovému polynomu pøiøazujeme stupeò~$-1$.

\s{Násobení polynomù:}
Polynomy násobíme jako výrazy:
$$
P(x)\cdot Q(x) = \Biggl( \sum_{i=0}^{n-1} p_i\cdot x^i \Biggr) \cdot \Biggl( \sum_{j=0}^{m-1} q_j\cdot x^j \Biggr).
$$
Po~roznásobení mù¾eme tento souèin zapsat jako polynom~$R(x)$, jeho¾
koeficient u~$x^k$ je roven $r_k = p_0q_k + p_1q_{k-1} + \ldots + p_kq_0$.

Snadno nahlédneme, ¾e polynom~$R$ má stupeò $\deg P + \deg Q$ a velikost
$\vert P\vert + \vert Q\vert - 1$.

Algoritmus, který poèítá souèin dvou polynomù velikosti~$n$ pøímo podle definice,
proto spotøebuje èas $\Theta(n)$ na výpoèet ka¾dého koeficientu, celkem
tedy $\Theta(n^2)$. Pokusíme se nalézt efektivnìj¹í zpùsob.

\s{Rovnost polynomù:}
Odboème na chvíli a uva¾ujme, kdy dva polynomy pova¾ujeme za stejné.
Na to se dá nahlí¾et více zpùsoby. Buïto se na polynomy mù¾eme dívat
jako na výrazy a porovnávat jejich symbolické zápisy. Pak jsou si dva
polynomy rovny právì tehdy, mají-li po~normalizaci stejné vektory koeficientù.
Tomu se øíká {\I identická rovnost} polynomù a obvykle se znaèí $P\equiv Q$.

Druhá mo¾nost je porovnávat polynomy jako reálné funkce. Polynomy $P$ a~$Q$
si tedy budou rovny ($P=Q$) právì tehdy, je-li $P(x)=Q(x)$ pro v¹echna $x\in{\bb R}$.
Identicky rovné polynomy si jsou rovny i jako funkce, ale musí to platit i naopak?
Následující vìta uká¾e, ¾e ano a ¾e dokonce staèí rovnost pro koneèný poèet~$x$.

\s{Vìta:} Buïte $P$ a~$Q$ polynomy stupnì nejvý¹e~$d$. Pokud platí $P(x_i) = Q(x_i)$
pro navzájem rùzná èísla $x_0,\ldots,x_d$, pak jsou $P$ a~$Q$ identické.

\proof
Pøipomeòme nejprve následující standardní lemma o~koøenech polynomù:

{\narrower\s{Lemma:}
Polynom~$R$ stupnì~$t\ge 0$ má nejvý¹e $t$ koøenù (èísel~$\alpha$, pro nì¾ je $P(\alpha)=0$).

\proof
Z~algebry víme, ¾e je-li èíslo $\alpha$ koøenem polynomu~$R(x)$, mù¾eme $R(x)$ beze zbytku vydìlit
výrazem $x-\alpha$. To znamená, ¾e $R(x)\equiv (x-\alpha)\cdot R'(x)$ pro nìjaký polynom $R'(x)$
stupnì $t-1$. Dosazením ovìøíme, ¾e koøeny polynomu~$R'$ jsou pøesnì tyté¾ jako koøeny polynomu~$R$,
s~mo¾nou výjimkou koøene~$\alpha$.

Budeme-li tento postup opakovat $t$-krát, buïto nám v~prùbìhu dojdou koøeny
(a~pak lemma jistì platí), nebo dostaneme rovnost $R(x) \equiv (x-\alpha_1)\cdot\ldots\cdot(x-\alpha_t)\cdot R''(x)$,
kde $R''$ je polynom nulového stupnì. Takové polynomy ov¹em nemohou mít ¾ádný
koøen, a~tím pádem nemù¾e mít ¾ádné dal¹í koøeny ani~$R$.
\qed

}

\smallskip

\>Abychom dokázali vìtu, staèí uvá¾it polynom $R(x)\equiv P(x)-Q(x)$.
Tento polynom má stupeò nejvý¹e~$d$, ov¹em ka¾dé z~èísel $x_0,\ldots,x_d$
je jeho koøenem. Podle lemmatu musí tedy být identicky nulový, a~proto
$P\equiv Q$.
\qed

\s{Zpìt k~násobení:}
Vìta, ji¾ jsme právì dokázali, vlastnì øíká, ¾e polynomy mù¾eme reprezentovat
nejen vektorem koeficientù, ale také {\I vektorem funkèních hodnot} v~nìjakých smluvených
bodech -- tomuto vektoru budeme øíkat {\I graf polynomu.} Pokud zvolíme dostateènì
mnoho bodù, je polynom svým grafem jednoznaènì urèen.

V~této reprezentaci je pøitom násobení polynomù triviální: Souèin polynomù
$P$ a~$Q$ má v~bodì~$x$ hodnotu $P(x)\cdot Q(x)$. Staèí tedy grafy vynásobit
po~slo¾kách, co¾ zvládneme v~lineárním èase. Jen je potøeba dát pozor na to,
¾e souèin má vy¹¹í stupeò ne¾ jednotliví èinitelé, tak¾e si potøebujeme poøídit
dostateèný poèet bodù.

\s{Algoritmus pro násobení polynomù:}
\algo
\:Jsou dány polynomy $P$ a~$Q$ velikosti~$n$, urèené svými koeficienty.
  Bez újmy na obecnosti pøedpokládejme, ¾e horních $n/2$ koeficientù je
  u~obou polynomù nulových, tak¾e souèin $R=P\cdot Q$ bude také polynom
  velikosti~$n$.
\:Zvolíme navzájem rùzná èísla $x_0,\ldots,x_{n-1}$.
\:Spoèítáme grafy polynomù~$P$ a~$Q$, tedy vektory $(P(x_0),\ldots,P(x_{n-1}))$
  a $(Q(x_0),\ldots,Q(x_t))$.
\:Z~toho vypoèteme graf souèinu~$R$ vynásobením po~slo¾kách: $R(x_i)=P(x_i)\cdot Q(x_i)$.
\:Nalezneme koeficienty polynomu~$R$ tak, aby odpovídaly grafu.
\endalgo

\>Krok~4 trvá $\Theta(n)$, tak¾e rychlost celého algoritmu stojí a padá s~efektivitou
pøevodù mezi koeficientovou a hodnotovou reprezentací polynomù. To obecnì neumíme
v~lep¹ím ne¾ kvadratickém èase, ale zde máme mo¾nost volby bodù $x_0,\ldots,x_{n-1}$,
tak¾e si je zvolíme tak ¹ikovnì, aby pøevod ¹el provést rychle.

\h{Pokus o~vyhodnocení polynomu metodou Rozdìl a panuj}

Uva¾ujme polynom $P$ velikosti~$n$, který chceme vyhodnotit v~$n$~bodech.
Body si zvolíme tak, aby byly {\I spárované,} tedy aby tvoøily dvojice li¹ící
se pouze znaménkem: $\pm x_{0}, \pm x_{1}, \ldots, \pm x_{n/2-1}$.

Polynom $P$ mù¾eme rozlo¾it na~èleny se sudými exponenty a ty s~lichými:
$$P(x) = (p_{0}x^{0} + p_{2}x^{2} + \ldots + p_{n-2}x^{n-2}) + (p_{1}x^{1} +
        p_{3}x^{3} + \ldots + p_{n-1}x^{n-1}).$$
Navíc mù¾eme z~druhé závorky vytknout~$x$:
$$P(x) = (p_{0}x^{0} + p_{2}x^{2} + \ldots + p_{n-2}x^{n-2}) + x\cdot (p_{1}x^{0} +
        p_{3}x^{2} + \ldots + p_{n-1}x^{n-2}).$$
V~obou závorkách se nyní vyskytují pouze sudé mocniny~$x$. Proto mù¾eme
ka¾dou závorku pova¾ovat za~vyhodnocení nìjakého polynomu velikosti~$n/2$
v~bodì~$x^2$, tedy:
$$P(x) = P_s(x^2) + x\cdot P_\ell(x^2),$$
kde:
$$\eqalign{
P_s(t) &= p_0t^0 + p_{2}t^{1} + \ldots + p_{n - 2}t^{n - 2\over 2}, \cr
P_\ell(t) &= p_1t^0 + p_3t^1 + \ldots + p_{n - 1}t^{n - 2\over 2}. \cr
}$$
Navíc pokud podobným zpùsobem dosadíme do~$P$ hodnotu~$-x$, dostaneme:
$$P(-x) = P_s(x^2) - x\cdot P_\ell(x^2).$$

Vyhodnocení polynomu~$P$ v~bodech $\pm x_0,\ldots,\pm x_{n-1}$ tedy
mù¾eme pøevést na vyhodnocení polynomù $P_s$ a~$P_\ell$ polovièní velikosti
v~bodech $x_0^2,\ldots,x_{n-1}^2$.

To naznaèuje algoritmus s~èasovou slo¾itostí $T(n) = 2T(n/2) + \Theta(n)$
a z~Kuchaøkové vìty víme, ¾e tato rekurence má øe¹ení $T(n) = \Theta(n\log n)$.
Jediný problém je, ¾e tento algoritmus nefunguje: druhé mocniny, které pøedáme
rekurzivnímu volání, jsou v¾dy nezáporné, tak¾e u¾ nemohou být správnì spárované.
Ouvej.

Tedy \dots{} alespoò dokud poèítáme s~reálnými èísly. Uká¾eme, ¾e v~oboru
komplexních èísel u¾ mù¾eme zvolit body, které budou správnì spárované
i po~nìkolikerém umocnìní na druhou.

\h{Intermezzo: Pøipomenutí komplexních èísel}
\def\ii{\mkern 0.5mu{\bf i}}
\def\im{{\bf i}}
\def\e{{\rm e}}
\def\\{\hfil\break}

\s{Základní operace}

\itemize\ibull
\:Definice: ${\bb C} = \{a + b\ii \mid a,b \in {\bb R}\}$, $\im^2=1$.

\:Sèítání: $(a+b\ii)\pm(p+q\ii) = (a\pm p) + (b\pm q)\ii$.

\:Násobení: $(a+b\ii)(p+q\ii) = ap + aq\ii + bp\ii + bq\ii^2 = (ap-bq) + (aq+bp)\ii$.
Pro $\alpha\in{\bb R}$ je $\alpha(a+b\ii) = \alpha a + \alpha b\ii$.

\:Komplexní sdru¾ení: $\overline{a+b\ii} = a-b\ii$. \\
$\overline{\overline x} = x$, $\overline{x\pm y} = \overline{x} \pm
\overline{y}$, $\overline{x\cdot y} = \overline x \cdot \overline y$, $x
\cdot \overline x  = (a+b\ii)(a-b\ii) = a^2+b^2 \in {\bb R}$.

\:Absolutní hodnota: $\vert x \vert = \sqrt{x\cdot\overline{x}}$, tak¾e $\vert
a+b\ii \vert = \sqrt{a^2+b^2}$. \\
Také $\vert \alpha x \vert = \vert \alpha\vert \cdot \vert x \vert$.

\:Dìlení: $x/y = (x\cdot \overline{y}) / (y \cdot \overline{y})$.
\endlist

\s{Gau{\scharfs}ova rovina a goniometrický tvar}

\itemize\ibull
\:Komplexním èíslùm pøiøadíme body v~${\bb R}^2$: $a+b\ii \leftrightarrow (a,b)$.

\:$\vert x\vert$ je vzdálenost od~bodu $(0,0)$.

\:$\vert x\vert = 1$ pro èísla le¾ící na~jednotkové kru¾nici ({\I komplexní jednotky}). \\
Pak platí $x=\cos\varphi + \im\sin\varphi$ pro nìjaké $\varphi\in\left[
0,2\pi \right)$.

\:Pro libovolné $x\in{\bb C}$: $x=\vert x \vert \cdot (\cos\varphi(x) +
\im\sin\varphi(x))$. \\
Èíslu $\varphi(x)\in\left[ 0,2\pi \right)$ øíkáme {\I argument}
èísla~$x$, nìkdy znaèíme $\mathop{\rm arg} x$.

\:Navíc $\varphi({\overline{x}}) = -\varphi(x)$.
\endlist

\s{Exponenciální tvar}

\itemize\ibull
\:Eulerova formule: $\e^{\im\varphi} = \cos\varphi + \im\sin\varphi$.

\:Ka¾dé $x\in{\bb C}$ lze tedy zapsat jako $\vert x\vert \cdot \e^{\ii \varphi(x)}$.

\:Násobení: $xy = \left(\vert x\vert\cdot \e^{\ii\varphi(x)}\right) \cdot
        \left(\vert y\vert\cdot \e^{\ii\varphi(y)}\right) = \vert x\vert
        \cdot \vert y\vert \cdot \e^{\ii(\varphi(x) + \varphi(y))}$ \\
(absolutní hodnoty se násobí, argumenty sèítají).

\:Umocòování: $x^\alpha = \left(\vert x\vert\cdot \e^{\ii\varphi(x)}\right)^
        \alpha = {\vert x\vert}^\alpha\cdot \e^{\ii \alpha \varphi(x)}$.
\endlist

\s{Odmocniny z~jednièky}

\smallskip

Odmocòování v~komplexních èíslech obecnì není jednoznaèné: jestli¾e tøeba
budeme hledat ètvrtou odmocninu z~jednièky, toti¾ øe¹it rovnici $x^4=1$, nalezneme
hned ètyøi øe¹ení: 1, $-1$, $\im$ a $-\im$.

Prozkoumejme nyní obecnìji, jak se chovají $n$-té odmocniny z~jednièky, tedy
komplexní koøeny rovnice $x^n=1$:

\itemize\ibull
\:Jeliko¾ $\vert x^n \vert = \vert x \vert ^n$, musí být $\vert x \vert = 1$.
  Proto $x=\e^{\ii\varphi}$ pro nìjaké~$\varphi$.
\:Má platit $1 = x^n = \e^{\ii\varphi n} = \cos{\varphi n} + \im\sin\varphi n$.
  To nastane, kdykoliv $\varphi n = 2k\pi$ pro nìjaké $k\in{\bb Z}$.
\endlist

\>Dostáváme tedy $n$ rùzných $n$-tých odmocnin z~1, toti¾
$\e^{2k\pi/n}$ pro $k=0,\ldots,n-1$.
Nìkteré z~tìchto odmocnin jsou ov¹em speciální:

\s{Definice:} Komplexní èíslo~$x$ je {\I primitivní} $n$-tá odmocnina z~1, pokud $x^n=1$
a ¾ádné z~èísel $x^1,x^2,\ldots,x^{n-1}$ není rovno~1.

\s{Pøíklad:} Ze~ètyø zmínìných ètvrtých odmocnin z~1 jsou $\im$ a $-\im$ primitivní
a druhé dvì nikoliv (ovìøte sami dosazením). Pro obecné~$n>2$ v¾dy existují alespoò
dvì primitivní odmocniny, toti¾ èísla $\omega = \e^{2\pi \ii / k}$ a $\overline\omega = \e^{-2\pi\ii/k}$.
Platí toti¾, ¾e $\omega^j = \e^{2\pi\ii j/k}$, a~to je rovno~1 právì tehdy,
je-li $j$ násobkem~$k$ (jednotlivé mocniny èísla~$\omega$ postupnì obíhají
jednotkovou kru¾nici). Analogicky pro~$\overline\omega$.

\s{Pozorování:} Pro sudé~$n$ a libovolné èíslo~$\omega$, které je primitivní
$n$-tou odmocninou z~jednièky, platí:

\itemize\ibull
\:$\omega^j \neq \omega^k$, kdykoliv $0\le j<k<n$. Staèí se podívat na podíl
  $\omega^k / \omega^j = \omega^{k-j}$. Ten nemù¾e být roven jedné,
  proto¾e $0 < k-j < n$ a $\omega$ je primitivní.
\:Pro sudé~$n$ je $\omega^{n/2} = -1$. Platí toti¾ $(\omega^{n/2})^2 = \omega^n = 1$,
  tak¾e $\omega^{n/2}$ je druhá odmocnina z~1. Takové odmocniny jsou jenom dvì:
  1 a $-1$, ov¹em 1 to být nemù¾e, proto¾e $\omega$ je primitivní.
\endlist

\h{Vyhodnocení polynomù podruhé -- algoritmus FFT}

Uká¾eme, ¾e primitivních odmocnin lze vyu¾ít k~záchranì na¹eho párovacího
algoritmu na vyhodnocování polynomù.

Nejprve polynomy doplníme nulami tak, aby jejich velikost~$n$ byla mocninou dvojky.
Poté zvolíme nìjakou primitivní $n$-tou odmocninu z~jednièky~$\omega$ a budeme
polynom vyhodnocovat v~bodech $\omega^0,\omega^1,\ldots,\omega^{n-1}$. To jsou
navzájem rùzná komplexní èísla, která jsou správnì spárovaná -- hodnoty $\omega^{n/2},
\ldots,\omega^{n-1}$ se od $\omega^0,\ldots,\omega^{n/2-1}$ li¹í pouze znaménkem.
To snadno ovìøíme: pro $0\le j<n/2$ je $\omega^{n/2+j} = \omega^{n/2} \omega^j = -\omega^j$.
Navíc $\omega^2$ je primitivní $(n/2)$-tá odmocnina z~jednièky, tak¾e se rekurzivnì opìt
voláme na problém tého¾ druhu, a~ten je opìt správnì spárovaný.

Ná¹ plán pou¾ít metodu Rozdìl a panuj tedy vy¹el, opravdu máme algoritmus o~slo¾itosti
$\Theta(n\log n)$ pro vyhodnocení polynomu. Je¹tì ho upravíme tak, aby místo s~polynomy pracoval
s~vektory jejich koeficientù èi hodnot. Tomuto algoritmu se øíká FFT, vzápìtí
prozradíme, proè.

\s{Algoritmus FFT:}
\algo
\algin Èíslo~$n=2^k$, primitivní $n$-tá odmocnina z~jednièky~$\omega$
  a vektor $(p_0,\ldots,p_{n-1})$ koeficientù polynomu~$P$.
\:Pokud $n=1$, polo¾íme $y_0\leftarrow p_0$ a skonèíme.
\:Jinak se rekurzivnì zavoláme na sudou a lichou èást koeficientù:
\::$(s_0,\ldots,s_{n/2-1}) \leftarrow {\rm FFT}(n/2,\omega^2,(p_0,p_2,p_4,\ldots,p_{n-2}))$.
\::$(\ell_0,\ldots,\ell_{n/2-1}) \leftarrow {\rm FFT}(n/2,\omega^2,(p_1,p_3,p_5,\ldots,p_{n-1}))$.
\:Z~grafù obou èástí poskládáme graf celého polynomu:
\::Pro $j=0,\ldots,n/2-1$:
\:::$y_j \leftarrow s_j + \omega^j\cdot \ell_j$.
\:::$y_{j+n/2} \leftarrow s_j - \omega^j\cdot \ell_j$.
\algout Graf polynomu~$P$, tedy vektor $(y_0,\ldots,y_{n-1})$, kde $y_j = P(\omega^j)$.
\endalgo

\h{Diskrétní Fourierova transformace}

Vyhodnotit polynom v~mocninách èísla~$\omega$ umíme, ale je¹tì nejsme v~cíli.
Potøebujeme umìt provést dostateènì rychle i opaèný pøevod -- z~hodnot na~koeficienty. K~tomu nám pomù¾e
podívat se na vyhodnocování polynomu trochu abstraktnìji jako na nìjaké zobrazení,
které jednomu vektoru komplexních èísel pøiøadí jiný vektor. Toto zobrazení
matematici v~mnoha rùzných kontextech potkávají u¾ nìkolik staletí a nazývají
ho Fourierovou transformací.

\s{Definice:}
{\I Diskrétní Fourierova transformace (DFT)} je zobrazení ${\cal F}: { {\bb C} ^n} \rightarrow { {\bb C} ^n}$,
které vektoru ${\bf x}$ pøiøadí vektor ${\bf y}$, jeho¾ slo¾ky jsou dány
pøedpisem $$y_j = \sum_{k=0}^{n-1} x_k \cdot \omega^{jk},$$
kde $\omega$ je nìjaká pevnì zvolená primitivní $n$-tá odmocnina z~jedné.

\s{Pozorování:}
Pokud oznaèíme ${\bf p}$ vektor koeficientù nìjakého polynomu~$P$, pak jeho
Fourierova transformace ${\cal F}({\bf p})$ není nic jiného ne¾ graf tohoto polynomu
v~bodech $\omega^0,\ldots,\omega^{n-1}$. To snadno ovìøíme dosazením do definice.

Ná¹ algoritmus tedy poèítá diskrétní Fourierovu transformaci v~èase $\Theta(n\log n)$.
Odtud pramení jeho název FFT -- Fast Fourier Transform.

\def\OO{{\bf\Omega}}

Také si v¹imnìme, ¾e DFT je lineární zobrazení. Mù¾eme ho proto zapsat
jako násobení nìjakou maticí~$\OO$, kde $\OO_{jk} = \omega^{jk}$.
Pro pøevod grafu na~koeficienty tedy potøebujeme najít inverzní zobrazení
urèené inverzní maticí $\OO^{-1}$.

Jeliko¾ $\omega^{-1} = \overline{\omega}$, pojïme zkusit, zda $\overline{\OO}$
není hledanou inverzní maticí.

\s{Lemma:} $\OO \cdot \overline{\OO} = n\cdot{\bf E}$, kde ${\bf E}$ je jednotková matice.

\proof
Dosazením do definice a elementárními úpravami:
$$\eqalign{
(\OO\cdot \overline{\OO})_{jk}
&= \sum_{\ell=0}^{n-1} \OO_{j\ell}\cdot\overline{\OO}_{\ell k}
= \sum_{\ell=0}^{n-1} \omega^{j\ell} \cdot \overline{\omega^{\ell k}}
= \sum_{\ell=0}^{n-1} \omega^{j\ell} \cdot \overline{\omega}^{\ell k}\cr
&= \sum_{\ell=0}^{n-1} \omega^{j\ell} \cdot (\omega^{-1})^{\ell k}
= \sum_{\ell=0}^{n-1} \omega^{j\ell} \cdot \omega^{-\ell k}
= \sum_{\ell=0}^{n-1} \omega^{(j-k)\ell}.\cr
}$$
To je ov¹em geometrická øada. Pokud je $j=k$, jsou v¹echny èleny øady
jednièky, tak¾e se seètou na~$n$. Pro $j\ne k$ pou¾ijeme známý vztah
pro souèet geometrické øady s~kvocientem $q=\omega^{j-k}$:
$$
\sum_{\ell=0}^{n-1} q^\ell
= {q^n - 1 \over q-1}
= {\omega^{(j-k)n} - 1 \over \omega^{j-k} - 1}
= 0.
$$
Poslední rovnost platí díky tomu, ¾e
$\omega^{(j-k)n} = (\omega^n)^{j-k} = 1^{j-k} = 1$, tak¾e
èitatel zlomku je nulový; naopak jmenovatel urèitì nulový není, jeliko¾
$\omega$ je primitivní a $0<\vert j-k\vert <n$.
\qed

\s{Dùsledek:} $\OO^{-1} = (1/n)\cdot\overline{\OO}$.

Matice~$\OO$ tedy je regulární a její inverze se kromì vydìlení~$n$
li¹í pouze komplexním sdru¾ením. Navíc $\overline{\omega} = \omega^{-1}$
je také primitivní $n$-tou odmocninou z~jednièky, tak¾e a¾ na faktor $1/n$
se jedná opìt o~Fourierovu transformaci, kterou mù¾eme spoèítat stejným
algoritmem FFT. Shròme, co jsme zjistili, do následující vìty:

\s{Vìta:} Je-li $n$ mocnina dvojky, lze v~èase $\Theta(n\log n)$ spoèítat
diskrétní Fourierovu transformaci v~${\bb C}^n$ i její inverzi.

Tím jsme také doplnili poslední èást algoritmu na násobení polynomù:

\s{Vìta:} Polynomy velikosti~$n$ nad tìlesem ${\bb C}$ lze násobit v~èase
$\Theta(n\log n)$.

V~obou vìtách pøitom èiníme pøedpoklad, ¾e základní operace s~komplexními
èísly umíme provést v~konstantním èase. Pokud tomu tak není, staèí slo¾itost
vynásobit slo¾itostí jedné operace.

\s{Dal¹í pou¾ití FFT:}
Dodejme je¹tì, ¾e Fourierova transformace se hodí i k~jiným vìcem ne¾
k~násobení polynomù. Své uplatnìní nachází i v~dal¹ích algebraických algoritmech,
tøeba v~násobení velkých èísel (lze se dostat a¾ ke~slo¾itosti $\Theta(n)$).
Mimo to skýtá i lecjaké fyzikální aplikace -- odpovídá toti¾ spektrálnímu
rozkladu signálu na siny a cosiny o~rùzných frekvencích. Na tom jsou zalo¾eny
napøíklad algoritmy pro filtrování zvuku, také pro kompresi zvuku a obrazu (MP3, JPEG)
nebo rozpoznávání øeèi.

\h{Paralelní implementace FFT}

Zkusme si je¹tì prùbìh algoritmu FFT znázornit graficky
Na~levé stranì následujícího obrázku se nachází vstupní vektor $x_0,\ldots,x_{n-1}$
(v~nìjakém poøadí), na~pravé stranì pak výstupní vektor $y_0,\ldots,y_{n-1}$.
Sledujme chod algoritmu pozpátku: Výstup spoèítáme z~výsledkù \uv{polovièních}
transformací vektorù $x_0,x_2,\ldots,x_{n-2}$ a $x_1,x_3,\ldots,x_{n-1}$.
Èerné krou¾ky pøitom odpovídají výpoètu lineární kombinace $a+\omega^kb$,
kde $a,b$ jsou vstupy krou¾ku a $k$~nìjaké pøirozené èíslo závislé na poloze
krou¾ku. Ka¾dá z~polovièních transformací se poèítá analogicky z~výsledkù
transformace velikosti $n/4$ atd. Celkovì výpoèet probíhá v~$\log_2 n$ vrstvách
po~$\Theta(n)$ operacích.

\figure{img.eps}{Pøíklad prùbìhu algoritmu pro~vstup velikosti 8}{3in}

Na obrázek se také mù¾eme dívat jako na schéma hradlové sítì pro výpoèet DFT.
Krou¾ky jsou pøitom \uv{hradla} pracující s~komplexními èísly. V¹echny operace
v~jedné vrstvì jsou na sobì nezávislé, tak¾e je sí» poèítá paralelnì. Sí» proto
pracuje v~èase $\Theta(\log n)$ a prostoru $\Theta(n)$, opìt pøípadnì násobeno
slo¾itostí jedné operace s~komplexními èísly

\s{Cvièení:} Doka¾te, ¾e permutace vektoru $x_0,\ldots,x_{n-1}$ na levé stranì
hradlové sítì odpovídá bitovému zrcadlení, tedy ¾e na pozici~$b$ shora se vyskytuje
prvek $x_d$, kde~$d$ je èíslo~$b$ zapsané ve~dvojkové soustavì pozpátku.

\h{Nerekurzivní FFT}

Obvod z~pøedchozího obrázku také mù¾eme vyhodnocovat po~hladinách zleva doprava,
èím¾ získáme elegantní nerekurzivní algoritmus pro výpoèet FFT v~èase $\Theta(n\log n)$
a prostoru $\Theta(n)$:

\algo
\algin $x_0,\ldots,x_{n-1}$
\:Pro $k=0,\ldots,n-1$ polo¾íme $y_k\= x_{r(k)}$, kde $r$ je funkce bitového zrcadlení.
\:Pøedpoèítáme tabulku hodnot $\omega^0,\omega^1,\ldots,\omega^{n-1}$.
\:$b\= 1$ \cmt{velikost bloku}
\:Dokud $b<n$, opakujeme:
\::Pro $j=0,\ldots,n-1$ s~krokem~$2b$ opakujeme: \cmt{zaèátek bloku}
\:::Pro $k=0,\ldots,b-1$ opakujeme: \cmt{pozice v~bloku}
\::::$\alpha\=\omega^{nk/2b}$
\::::$(y_{j+k},y_{j+k+b}) \= (y_{j+k}+\alpha\cdot y_{j+k+b}, y_{j+k}-\alpha\cdot y_{j+k+b})$.
\::$b\= 2b$
\algout $y_0,\ldots,y_{n-1}$
\endalgo

\h{FFT v~koneèných tìlesech}

Nakonec dodejme, ¾e Fourierovu transformaci lze zavést nejen nad tìlesem
komplexních èísel, ale i v~nìkterých koneèných tìlesech, pokud zaruèíme
existenci primitivní $n$-té odmocniny z~jednièky. Napøíklad v~tìlese ${\bb
Z}_p$ pro prvoèíslo $p=2^k+1$ platí $2^k=-1$, tak¾e $2^{2k}=1$
a $2^0,2^1,\ldots,2^{2k-1}$ jsou navzájem rùzná. Èíslo~2 je tedy primitivní
$2k$-tá odmocnina z~jedné. To se nám ov¹em nehodí pro algoritmus FFT, nebo»
$2k$ bude málokdy mocnina dvojky.

Zachrání nás ov¹em algebraická vìta, která øíká, ¾e multiplikativní grupa\foot{To je
mno¾ina v¹ech nenulových prvkù tìlesa s~operací násobení.}
libovolného koneèného tìlesa je cyklická, tedy ¾e v¹echny nenulové prvky tìlesa lze
zapsat jako mocniny nìjakého èísla~$g$ (generátoru grupy). Napøíklad pro $p=2^{16}+1=65\,537$
je jedním takovým generátorem èíslo~$3$. Jeliko¾ mezi èísly $g^1,g^2,\ldots,g^{p-1}$
se musí ka¾dý nenulový prvek tìlesa vyskytnout právì jednou, je~$g$ primitivní
$2^k$-tou odmocninou z~jednièky, tak¾e mù¾eme poèítat FFT pro libovolný vektor,
jeho¾ velikost je mocnina dvojky men¹í nebo rovná $2^k$.

Bli¾¹í prùzkum
na¹ich úvah o~FFT dokonce odhalí, ¾e není ani potøeba tìleso. Postaèí libovolný
komutativní okruh, ve~kterém existuje pøíslu¹ná primitivní odmocnina
z~jednièky, její multiplikativní inverze (ta ov¹em existuje v¾dy, proto¾e
$\omega^{-1} = \omega^{n-1}$) a multiplikativní inverze èísla~$n$. To nám
poskytuje je¹tì daleko více volnosti ne¾ tìlesa, ale není snadné takové okruhy
hledat.

Výhodou tìchto podob Fourierovy transformace je, ¾e na rozdíl od~té klasické komplexní
nejsou zatí¾eny zaokrouhlovacími chybami (komplexní odmocniny z~jednièky mají obì slo¾ky
iracionální). To se hodí napøíklad ve~zmiòovaných algoritmech na násobení velkých èísel.

\h{Cvièení}

\itemize\ibull
\:O~jakých vlastnostech vektoru vypovídá nultý a $(n/2)$-tý koeficient jeho
  Fourierova obrazu (výsledku Fourierovy transformace)?
\:Spoèítejte Fourierovy obrazy následujících vektorù z~${\bb C}^n$:
    \itemize\ibull
    \:$(x,\ldots,x)$
    \:$(1,-1,1,-1,\ldots,1,-1)$
    \:$(\omega^0,\omega^1,\omega^2,\ldots,\omega^{n-1})$
    \:$(\omega^0,\omega^2,\omega^4,\ldots,\omega^{2n-2})$
    \endlist
\:Roz¹íøením výsledkù z~pøedchozího cvièení najdìte pro ka¾dé~$j$ vektor, jeho¾
  Fourierova transformace má na $j$-tém místì jednièku a v¹ude jinde nuly. Z~toho lze pøímo
  sestrojit inverzní transformaci.
\:Uka¾te, ¾e je-li $\bf x$ reálný vektor z~${\bb R}^n$, je jeho Fourierova transformace ${\bf y}={\cal F}({\bf x})$
  {\I antisymetrická:} ${\bf y}_j = \overline{{\bf y}_{n-j}}$ pro v¹echna~$j$.
\:Podobnì uka¾te, ¾e Fourierova transformace ka¾dého antisymetrického vektoru je reálná.
\:Uva¾ujme reálnou funkci~$f$ definovanou na intervalu $[0,2\pi)$. Pokud její
  hodnoty {\I navzorkujeme} v~$n$ pravidelnì rozmístìných bodech, získáme vektor
  ${\bf f}\in {\bb R}^n$ o~slo¾kách ${\bf f}_j = f(2\pi j/n)$. Jak vypadá Fourierova
  transformace tohoto vektoru pro následující funkce?
    \itemize\ibull
    \:$\e^{\im kx}$ pro $k\in{\bb N}$
    \:$\cos kx$
    \:$\sin kx$
    \endlist
  Nápovìda: sinus a cosinus mù¾ete zapsat jako lineární kombinaci dvou komplexních exponenciál.
\:Pomocí pøedchozího cvièení doka¾te, ¾e libovolnou reálnou funkci na $[0,2\pi)$
  existuje lineární kombinace funkcí $\sin kx$ a $\cos kx$, která pøi vzorkování
  v~$n$ bodech není od zadané funkce rozli¹itelná.

  Pøesnìji øeèeno, pro ka¾dý vektor~${\bf f}\in {\bb R}^n$ existují vektory ${\bf a},{\bf b}\in {\bb R}^n$
  takové, ¾e platí:
  $$
  {\bf f}_j = \sum_{k=0}^{n-1} {\bf a}_k \sin {2jk\pi \over n} + {\bf b}_k \cos {2jk\pi \over n}.
  $$
  Koeficienty $a_k$ a $b_k$ lze pøítom snadno získat z~Fourierova obrazu vektoru~${\bf f}$.
\endlist

\bye