3.7 Brza Furijeova transformacija

Brza Furijeova transformacija (ili FFT, što je skraćenica od engl. Fast Fourier transform) od svog otkrića sredinom šezdesetih godina dvadesetog veka spada u najznačajnije algoritme 20. veka.¹ Iako su značajne primene brze Furijeove transformacije nastale nakon rada Kulija i Tukija 1965. godine, na polju obrade signala (u cilju detektovanja nuklearnih proba obradom seizmičkih merenja), ispostavilo se da su ova dva autora nezavisno iznova “izmislili” algoritam koji je osmislio Gaus još davne 1805. godine. Brza Furijeova transformacija ima brojne primene u inženjerstvu, u obradi digitalnih signala poput obrade zvuka i obrade digitalnih slika i u matematici. Furijeovom transformacijom neprekidnog signala vrši se njegovo prevođenje iz vremenskog u frekvencijski domen, tj. signal se razlaže na zbir sinusoidalnih elementarnih talasa (na primer, muzički akord se može razložiti na pojedinačne tonove koja ga sačinjavaju). Ako je signal predstavljen vektorom dobijenim merenjem njegovog intenziteta u diskretnim vremenskim intervalima, primenjuje se diskretna Furijeova transformacija. Na slici 1 prikazana su tri signala, posebno u vremenskom, a posebno u frekvencijskom domenu. Vidi se da su prva dva pravilne sinusoide (imaju samo jednu dominantnu frekvenciju), dok je treći dobijen sabiranjem dve pravilne sinusoide (ima dve dominantne frekvencije).

Brza Furijeova transformacija je zapravo samo efikasan algoritam za izračunavanje diskretne Furijeove transformacije. Mi ćemo se u ovom udžbeniku ograničiti na jednu njenu primenu, a to je množenje polinoma.

Formulacija problema koji treba rešiti je precizna samo na prvi pogled, jer nije preciziran način na koji su polinomi predstavljeni. Obično se polinom \(P(x)=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\ldots+a_1x+a_0,\) stepena \(n-1\), predstavlja nizom svojih \(n\) koeficijenata \(a_0, a_1, a_2, \ldots a_{n-1}\) uz \(1,x, x^2, \ldots, x^{n-1}\), međutim, to nije jedina mogućnost. Alternativno, polinom stepena \(n-1\) moguće je predstaviti njegovim vrednostima u \(n\) različitih tačaka: te vrednosti jednoznačno određuju polinom. Na primer, polinom stepena \(1\), odnosno linearna funkcija, jedinstveno je određena vrednostima u dve različite tačke, a polinom stepena \(2\), odnosno kvadratna funkcija, vrednostima u tri različite tačke.

Množenje dva polinoma koja su zadata nizom svojih koeficijenata može se izvršiti osnovnim algoritmom složenosti \(O(n^2)\) koji se sastoji u množenju svakog monoma prvog polinoma svakim monomom drugog, ili korišćenjem Karacubinog algoritma zasnovanog na pametnoj primeni dekompozicije koji je složenosti \(O(n^{\log_2 3})\). Cilj ovog poglavlja je da primenom algoritma FFT izvedemo algoritam složenosti \(O(n\log n)\), ali za to je potrebno da detaljnije razmotrimo drugu reprezentaciju.

Predstavljanje polinoma vrednostima na skupu tačaka je interesantno zbog jednostavnosti množenja. Naime, proizvod dva polinoma stepena \(n-1\) je polinom stepena \(2n-2\), pa je određen svojim vrednostima u \(2n-1\) tačaka. Ako pretpostavimo da su vrednosti polinoma-činilaca stepena \(n-1\) date na istom skupu od \(2n-1\) različitih tačaka, onda se proizvod ova dva polinoma može izračunati pomoću \(2n-1\), odnosno \(O(n)\) običnih množenja. Naime, vrednost polinoma \(PQ\) u bilo kojoj tački \(x\) jednaka je proizvodu vrednosti polinoma \(P\) u tački \(x\) i vrednosti polinoma \(Q\) u tački \(x\), tj. važi \((P\cdot Q)(x) = P(x)\cdot Q(x)\).

Nažalost, predstavljanje polinoma njegovim vrednostima na skupu tačaka nije pogodno za neke primene. Primer je izračunavanje vrednosti polinoma u proizvoljnoj tački: pri reprezentaciji polinoma vrednostima na skupu tačaka, ovo je mnogo teže u odnosu na to kada je polinom zadat nizom svojih koeficijenata (potrebno je rešiti problem interpolacije). Sa druge strane, vrednost polinoma \(P(x)\) stepena \(n-1\) zadatog nizom svojih koeficijenata u proizvoljnoj tački \(x\) može se pomoću Hornerove šeme izračunati korišćenjem \(n\) množenja, tj. algoritmom složenosti \(O(n)\):

\[P(x)=a_0+x(a_1+x(a_2+\ldots+x\cdot a_{n-1}\ldots))\]

Dakle, svaka od reprezentacija je pogodna za jednu, a nije pogodna za drugu operaciju, kako je ilustrovano u tabeli 1.

Ako bismo mogli da efikasno prevodimo polinome iz jedne u drugu reprezentaciju, dobili bismo mogućnost da obe najznačajnije operacije (množenje i izračunavanje vrednosti) izvršavamo efikasno. Primenom brze Furijeove transformacije, tj. algoritma FFT, složenost prevođenja iz jedne u drugu reprezentaciju polinoma postaje \(O(n\log n)\).

Bez primene algoritma FFT prelaz od predstavljanja polinoma koeficijentima na predstavljanje polinoma vrednostima u tačkama (koji ćemo zvati tabeliranje), rešava se uzastopnim izračunavanjem vrednosti polinoma u svakoj od tih tačaka. Izračunavanje vrednosti polinoma \(P(x)\) stepena \(n-1\) u \(n\) proizvoljnih tačaka (što je potrebno za jednoznačno predstavljanje polinoma) izvodljivo je primenom Hornerove sheme pomoću \(O(n^2)\) množenja. Prelaz od predstavljanja polinoma vrednostima na skupu tačaka na predstavljanje koeficijentima se zove interpolacija. Vrednosti koeficijenta polinoma na osnovu vrednosti \((x_j,y_j), 0\le j\leq n\) polinoma u \(n+1\) različitih tačaka mogu se odrediti korišćenjem Lagranžove interpolacione formule (koja se obično izučava u sklopu numeričke matematike):

\[P(x)=\sum_{j=0}^{n-1}a_jx^j=\sum_{j=0}^{n-1}{y_j\cdot\frac{\prod_{k\neq j} (x-x_k)}{\prod_{k\neq j} (x_j-x_k)}}\]

Na osnovu ovog izraza koeficijenti polinoma mogu se odrediti algoritmom složenosti \(O(n^2)\).

Postavlja se pitanje kako je moguće ubrzati korake tabeliranja i interpolacije? Ključna ideja da se ne koristi proizvoljnih \(n\) tačaka: mi imamo slobodu da po želji izaberemo pogodan skup od \(n\) različitih tačaka. Brza Furijeova transformacija koristi specijalan skup tačaka, tako da se i tabeliranje i interpolacija mogu efikasnije izvršavati.

3.7.1 Direktna brza Furijeova transformacija

Osnovna ideja na kojoj počiva brza Furijeova transformacija je izložena u narednom primeru.

Brza Furijeova transformacija je algoritam tipa podeli-pa-vladaj. Pretpostavlja se da je dimenzija ulaza stepen dvojke. Ako niz koeficijenata nije dužine \(n=2^k\), on se uvek može dopuniti nulama. Razmotrimo za početak slučaj \(n=8=2^3\) tj.
pretpostavimo da je potrebno izračunati vrednost polinoma

\[P(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4 + a_5x^5 + a_6x^6 + a_7x^7\]

u nekih 8 različitih tačaka \(x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\).

Pregrupišimo sabirke na sledeći način:

\[\begin{align*} P(x) &= a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4 + a_5x^5 + a_6x^6 + a_7x^7\\ &= (a_0 + a_2x^2 + a_4x^4 + a_6x^6) + x(a_1 + a_3x^2 + a_5x^4 + a_7x^6). \end{align*}\]

Na ovaj način dobijamo dva manja polinoma \(P_0(x) = a_0 + a_2x + a_4x^2 + a_6x^3\) i \(P_1(x) = a_1 + a_3x + a_5x^2 + a_7x^3\). Važi da je \(P(x_i) = P_0(x_i^2) + x_iP_1(x_i^2)\). Potrebno je, dakle, da izračunamo sledećih 8 vrednosti polinoma \(P\):

Da bi se ovo direktno izračunalo, potrebno je izračunati ukupno 16 vrednosti polinoma \(P_0\) i \(P_1\), međutim, pametnim izborom tačaka moguće je taj broj smanjiti na 8. Naime mi imamo slobodu izbora tačaka \(x_k\) (dok god ih biramo tako da su različite) i možemo odabrati tačke za koje, na primer, važi \(x_0^2 = x_4^2\), \(x_1^2 = x_5^2\), \(x_2^2 = x_6^2\), \(x_3^3 = x_7^2\). Ne smemo da uzmemo iste vrednosti tačaka, ali možemo suprotne. Naime, način da se postigne jednakost kvadrata je da se vrednosti izaberu tako da važi \(x_4=-x_0\), \(x_5=-x_1\), \(x_6=-x_2\) i \(x_7=-x_3\). Tada se izračunavanje svodi na:

Potrebno je, dakle, rekurzivno izračunati samo sledećih 8 vrednosti polinoma \(P_0\) i \(P_1\) (umesto polaznih 16), a onda ih iskombinovati pomoću prethodnih formula:

\[ \begin{array}{ll@{\hspace{1cm}}ll} {\color{red}{P_0(x_0^2)}} &= {\color{red}{a_0 + a_2x_0^2 + a_4(x_0^2)^2 + a_6(x_0^2)^3}} & {\color{red!70!black}{P_1(x_0^2)}} &= {\color{red!70!black}{a_1 + a_3x_0^2 + a_5(x_0^2)^2 + a_7(x_0^2)^3}}\\ {\color{blue}{P_0(x_1^2)}} &= {\color{blue}{a_0 + a_2x_1^2 + a_4(x_1^2)^2 + a_6(x_1^2)^3}} & {\color{blue!70!black}{P_1(x_1^2)}} &= {\color{blue!70!black}{a_1 + a_3x_1^2 + a_5(x_1^2)^2 + a_7(x_1^2)^3}}\\ {\color{green}{P_0(x_2^2)}} &= {\color{green}{a_0 + a_2x_2^2 + a_4(x_2^2)^2 + a_6(x_2^2)^3}} & {\color{green!70!black}{P_1(x_2^2)}} &= {\color{green!70!black}{a_1 + a_3x_2^2 + a_5(x_2^2)^2 + a_7(x_2^2)^3}}\\ {\color{yellow}{P_0(x_3^2)}} &= {\color{yellow}{a_0 + a_2x_3^2 + a_4(x_3^2)^2 + a_6(x_3^2)^3}} & {\color{yellow!70!black}{P_1(x_3^2)}} &= {\color{red!70!yellow}{a_1 + a_3x_3^2 + a_5(x_3^2)^2 + a_7(x_3^2)^3}} \end{array} \]

Da li možemo još uštedeti na broju operacija pametnim izborom tačaka \(x_0\), \(x_1\), \(x_2\) i \(x_3\)? Razmotrimo izračunavanje vrednosti polinoma \(P_0\) (vrednosti polinoma \(P_1\) se računaju analogno). Grupišimo ponovo koeficijente na parnim i na neparnim pozicijama. Time dobijamo još manje polinome \(P_{00}(x) = a_0 + a_4x\) i \(P_{01}(x) = a_2 + a_6x\).

\[ \begin{array}{llll} {\color{red}{P_0(x_0^2)}} &= {\color{Salmon}{P_{00}(x_0^4)}} + x_0^2{\color{Salmon!70!black}{P_{01}(x_0^4)}} &= {\color{Salmon}{(a_0 + a_4x_0^4)}} + x_0^2{\color{Salmon!70!black}{(a_2 + a_6x_0^4)}}\\ {\color{green}{P_0(x_1^2)}} &= {\color{ForestGreen}{P_{00}(x_1^4)}} + x_1^2{\color{ForestGreen!70!black}{P_{01}(x_1^4)}} &= {\color{ForestGreen}{(a_0 + a_4x_1^4)}} + x_1^2{\color{ForestGreen!70!black}{(a_2 + a_6x_1^4)}}\\ {\color{blue}{P_0(x_2^2)}} &= {\color{Turquoise}{P_{00}(x_2^4)}} + x_2^2{\color{Turquoise!70!black}{P_{01}(x_2^4)}} &= {\color{Turquoise}{(a_0 + a_4x_2^4)}} + x_2^2{\color{Turquoise!70!black}{(a_2 + a_6x_2^4)}}\\ {\color{yellow}{P_0(x_3^2)}} &= {\color{Apricot}{P_{00}(x_3^4)}} + x_3^2{\color{Apricot!70!black}{P_{01}(x_3^4)}} &= {\color{Apricot}{(a_0 + a_4x_3^4)}} + x_3^2{\color{Apricot!70!black}{(a_2 + a_6x_3^4)}} \end{array} \]

Za izračunavanje \(4\) vrednosti polinoma \(P_0\) potrebno je izračunati ukupno \(8\) vrednosti polinoma \(P_{00}\) i \(P_{01}\) i želimo to da smanjimo. Međutim, sada imamo problem. Da bismo mogli da smanjimo broj izračunavanja na isti način kao u prvom koraku, potrebno je da odaberemo \(x_0\) i \(x_2\) tako da važi \(x_0^4 = x_2^4\). Ne možemo da odaberemo da je \(x_0^2=x_2^2\), jer je tada ili \(x_0 = x_2\) ili \(x_0 = -x_2 = x_6\), što nije dopušteno, jer sve tačke moraju biti različite. Dakle, mora da važi da je \(x_2^2 = -x_0^2\). Međutim to nije moguće ako su \(x_0\) i \(x_2\) realni različiti brojevi, jer su tada i \(x_0^2\) i \(x_2^2\) nenegativni i bar jedan od njih je pozitivan. Da li je onda uopšte moguće izabrati različite vrednosti \(x_0\) i \(x_2\) tako da je \(x_2^2 = -x_0^2\)? Jeste, ako dopustimo da vrednosti \(x_i\) budu kompleksni brojevi. Naime, tražimo kompleksan koeficijent \(k\) takav da je \(x_2 = kx_0\) i \(x_2^4 = x_0^4\) tj. \(x_2^2 = (kx_0)^2 = k^2x_0^2 = -x_0^2\). Pošto je \(k^2=-1\), za \(k\) se može uzeti vrednost imaginarne jedinice \(i\), za koju važi \(i^2 = -1\). Zaista, ako je \(x_2=ix_0\), tada je \(x_2^2 = (ix_0)^2 = i^2 x_0^2 = -x_0^2\) i \(x_2^4 = (x_2^2)^2 = (-x_0^2)^2 = x_0^4\). Brojeve \(x_3\) i \(x_1\) možemo izabrati tako da je \(x_3=ix_1\) iz čega sledi da je \(x_3^2=-x_1^2\) i \(x_3^4 = x_1^4\). Nakon toga, izračunavanje potrebnih vrednosti poinoma \(P_0\) se svodi na izračunavanje samo \(4\) vrednosti polinoma \(P_{00}\) i \(P_{01}\):

\[ \begin{array}{llll} {\color{red}{P_0(x_0^2)}} &= {\color{Salmon}{P_{00}(x_0^4)}} + x_0^2{\color{Salmon!70!black}{P_{01}(x_0^4)}} &= {\color{Salmon}{(a_0 + a_4x_0^4)}} + x_0^2{\color{Salmon!70!black}{(a_2 + a_6x_0^4)}}\\ {\color{green}{P_0(x_1^2)}} &= {\color{ForestGreen}{P_{00}(x_1^4)}} + x_1^2{\color{ForestGreen!70!black}{P_{01}(x_1^4)}} &= {\color{ForestGreen}{(a_0 + a_4x_1^4)}} + x_1^2{\color{ForestGreen!70!black}{(a_2 + a_6x_1^4)}}\\ {\color{blue}{P_0(x_2^2)}} &= {\color{Salmon}{P_{00}(x_0^4)}} - x_0^2{\color{Salmon!70!black}{P_{01}(x_0^4)}} &= {\color{Salmon}{(a_0 + a_4x_0^4)}} + x_2^2{\color{Salmon!70!black}{(a_2 + a_6x_0^4)}}\\ {\color{yellow}{P_0(x_3^2)}} &= {\color{ForestGreen}{P_{00}(x_1^4)}} - x_1^2{\color{ForestGreen!70!black}{P_{01}(x_1^4)}} &= {\color{ForestGreen}{(a_0 + a_4x_1^4)}} + x_3^2{\color{ForestGreen!70!black}{(a_2 + a_6x_1^4)}} \end{array} \]

Prethodni način izbora tačaka \(x_2\) i \(x_3\) donosi sličnu uštedu i prilikom računanja vrednosti polinoma \(P_1\). Ovim se ceo problem svodi na izračunavanje sledećih \(8\) vrednosti polinoma \(P_{00}\), \(P_{01}\), \(P_{10}\) i \(P_{11}\):

\[ \begin{array}{ll@{\hspace{1cm}}ll} {\color{Salmon}{P_{00}(x_0^4)}} &= {\color{Salmon}{a_0 + a_4x_0^4}} & {\color{Salmon!70!black}{P_{01}(x_0^4)}} &= {\color{Salmon!70!black}{a_2 + a_6x_0^4}}\\ {\color{ForestGreen}{P_{00}(x_1^4)}} &= {\color{ForestGreen}{a_0 + a_4x_1^4}} & {\color{ForestGreen!70!black}{P_{01}(x_1^4)}} &= {\color{ForestGreen!70!black}{a_2 + a_6x_1^4}}\\ {\color{Thistle}{P_{10}(x_0^4)}} &= {\color{Thistle}{a_1 + a_5x_0^4}} & {\color{Thistle!70!black}{P_{11}(x_0^4)}} &= {\color{Thistle!70!black}{a_3 + a_7x_0^4}}\\ {\color{SkyBlue}{P_{10}(x_1^4)}} &= {\color{SkyBlue}{a_1 + a_5x_1^4}} & {\color{SkyBlue!70!black}{P_{11}(x_1^4)}} &= {\color{SkyBlue!70!black}{a_3 + a_7x_1^4}} \end{array} \]

Posednja ušteda se dobija time što se vrednosti \(x_0\) i \(x_1\) odaberu tako da je \(x_1^8 = x_0^8\) tj. \(x_1^4 = -x_0^4\). To se može postići ako je \(x_1 = kx_0\), gde se koeficijent \(k\) određuje tako da je \(k^4 = -1\). Pošto je \(e^{i\pi} = -1\), jedan takav broj je broj \(k=e^{i\frac{\pi}{4}}=\cos(\frac{\pi}{4}) + i\sin(\frac{\pi}{4}) = \frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}\). Podsetimo se Ojlerove formule:

\[e^{i\phi} = \cos\phi + i\sin\phi\]

Važi \((e^{i\phi})^n = e^{in\phi}\), što odgovara Moavrovoj formuli:

\[(\cos\phi + i\sin\phi)^n = \cos(n\phi) + i\sin(n\phi)\]

Tada se prethodno izračunavanje svodi na:

Mogli bismo reći da se vrednosti ovih linearnih polinoma izračunavaju kombinovanjem polinoma stepena \(0\) (u pitanju su \(8\) polinoma koji odgovaraju koeficijentima \(a_k\): \(P_{000} = a_0\), \(P_{001} = a_4\), \(P_{010} = a_2\), \(P_{011} = a_6\), \(P_{100} = a_1\), \(P_{101} = a_5\), \(P_{110} = a_3\) i \(P_{111} = a_7\)) i čija se vrednost trivijalno izračunava (bez daljih rekurzivnih poziva). Da bi se sve uklopilo u opštu shemu, možemo reći da se vrednost tih konstantnih \(8\) polinoma izračunava u tački \(x_0^8\).

Ako rezimiramo sve odnose između promenljivih koje smo do sada uspostavili, vidimo da još jedino možemo slobodno da biramo vrednost \(x_0\), dok se sve ostale vrednosti tačaka izračunavaju na osnovu vrednosti \(x_0\). Imamo niz tačaka

\[\begin{align*} x_0 & & x_4 &= -x_0\\ x_1 &= (\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}) x_0 & x_5 & = -x_1 = (-\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2})x_0\\ x_2 &= ix_0 & x_6 &= -x_2 = -ix_0\\ x_3 &= ix_1 = (-\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2})x_0 & x_7 &= -x_3 = (\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2})x_0 \end{align*}\]

Najjednostavnije formule se dobijaju, naravno, ako se odabere \(x_0 = 1\). Tada se može primetiti da su tačke \(x_i\) osmi koreni iz jedinice. Nadalje ćemo ove tačke obeležavati sa \(w_k\) umesto \(x_k\) (da bismo naglasili da su u pitanju kompleksni brojevi). Dakle, izračunavaćemo vrednosti polinoma redom u tačkama \(w_0=1=e^0\), \(w_1=\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}=e^{i\frac{\pi}{4}}\), \(w_2=i=e^{i\frac{\pi}{2}}\), \(w_3=-\frac{\sqrt{2}}{2}+i\frac{\sqrt{2}}{2}=e^{i\frac{3\pi}{4}}\), \(w_4=-1=e^{i\pi}\), \(w_5=-\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2}=e^{i\frac{5\pi}{4}}\), \(w_6=-i=e^{i\frac{3\pi}{2}}\), \(w_7=\frac{\sqrt{2}}{2}-i\frac{\sqrt{2}}{2}=e^{i\frac{7\pi}{4}}\), tj. u tačkama

\[w_k = e^{\frac{k\pi i}{4}} = e^{\frac{2k\pi i}{8}}, \quad 0 \leq k < 8.\]

U opštem slučaju ćemo izračunavati vrednosti polinoma stepena \(n\) (gde se \(n\) bira tako da bude stepen broja \(2\)) i izračunavaćemo vrednosti polinoma u tačkama \(w_k\) definisanim na sledeći način:

\[w_k = e^{\frac{2k\pi i}{n}}, \quad 0 \leq k < n\]

Na slici 3 je prikazano da se svi brojevi \(w_k\) mogu izraziti kao neki stepen broja \(w_1 = e^{\frac{2\pi i}{n}}\), koji ćemo obeležiti sa \(w\). Broj \(w\) je primitivni \(n\)-ti koren iz jedinice (u našem primeru, osmi), što znači da je \(w^n = 1\), dok je \(w^k \neq 1\) za sve \(0 < k < n\). Na primer, imaginarna jedinica \(i\) je primitivni četvrti koren iz jedinice jer je \(i^1 = i \neq 1\), \(i^2 = -1\neq 1\), \(i^3 = -i \neq 1\) i \(i^4 = 1\), dok broj \(-1\) jeste četvrti koren iz jedinice (jer je \((-1)^4=1\)), ali nije primitivan četvrti koren, jer je \((-1)^2 = 1\) (\(-1\) jeste primitivni drugi koren iz jedinice). Može se pokazati i da je recipročna vrednost broja \(w\) tj. broj \(w^{-1} = e^\frac{-2\pi i}{n}\) takođe primitivni \(n\)-ti koren iz jedinice (što će nam biti važno za inverznu Furijeovu transformaciji o kojoj će biti reči kasnije). Može se pokazati da je svaki broj oblika \(e^{ik\frac{2\pi}{n}}\), za uzajamno proste brojeve \(k\) i \(n\) primitivan \(n\)-ti koren iz jedinice. Kada je \(n\) stepen broja \(2\), tada se primitivni \(n\)-ti koreni dobijaju za sve neparne vrednosti \(k\). Svi brojevi \(w_k\) se mogu izraziti kao stepeni broja \(w\) (kao i bilo kog drugog primitivnog \(n\)-tog korena iz jedinice). Naime, očigledno važi da je

\[w_k = e^{\frac{2k\pi i}{n}} = \left(e^{\frac{2\pi i}{n}}\right)^k = w^k.\]

Još jedna važna činjenica u vezi sa primitivnim \(n\)-tim korenima iz jedinice je sledeća. Ako je \(w\) primitivni \(n\)-ti koren iz jedinice, za neki paran broj \(n > 0\), tada je \(w^2\) primitivni \(\frac{n}{2}\)-ti koren iz jedinice (videćemo da se ovo koristi tokom rekurzivnih poziva u brzoj Furijeovoj transformaciji). Zaista, važi:

\[w^2 = \left(e^\frac{2\pi i}{n}\right)^2 = e^\frac{2\cdot 2\pi i}{n} = e^\frac{2\pi i}{\frac{n}{2}}\]

Naglasimo da u izvođenjima poput prethodnih treba biti jako obazriv. Naime, koristili smo nekoliko puta zakon \((a^b)^c = a^{bc}\) koji ne mora biti tačan za proizvoljne kompleksne (pa čak ni realne brojeve). Na primer, \(1 = \sqrt{1} = \sqrt{(-1)^2} = ((-1)^2)^\frac{1}{2} \neq (-1)^{2 \cdot\frac{1}{2}} = (-1)^1 = -1\). Ipak, može se pokazati da su prethodna izvođenja korektna (pre svega jer je spoljni izložilac stepenovanja uvek bio celobrojan).

Ovo izračunavanje se shematski prikazuje korišćenjem tzv. “leptir” sheme (slika 4). Jedan “leptir” opisuje izračunavanje prikazano na slici 5.

Izračunavanje iz prethodnog primera je prikazano “leptir” shemom na slici 6.

Na slici 4 se jasno vidi rekurzivna struktura algoritma (svaki rekurzivni poziv je predstavljen jednim obojenim pravougaonikom). Struktura rekurzivnih poziva prikazana je i na slici 7.

Svaki rekurzivni poziv prima polinom \(P_n\) tj. niz njegovih koeficijenata \(a_i\) dužine \(n\) i broj \(W\) koji je neki stepen polaznog primitivnog korena iz jedinice (a zapravo je primitivni \(n\)-ti koren iz jedinice). Rekurzivni poziv izračunava vrednosti polinoma \(P_n\) u tačkama \(W^0, W^1, \ldots W^{n-1}\). Ako je \(n=1\), izlazi se iz rekurzije i rezultat je jedini koeficijent \(a_i\). U suprotnom se vrše dva rekurzivna poziva (za polinome dobijene od koeficijenata na parnim i neparnim pozicijama) i rezultati se objedinjavaju na osnovu pravila, čija korektnost sledi iz naredne leme.

Opis i implementacija algoritma brze Furijeove transformacije dati su u poglavlju 3.7.3.

3.7.2 Inverzna Furijeova transformacija

Brza Furijeova transformacija rešava samo pola problema: na osnovu koeficijenata polinoma stepena \(n\) određuju se njegove vrednosti u \(n\) specijalno odabranih tačaka tj. vrši se tabeliranje polinoma. Ostaje pitanje kako da na osnovu vrednosti polinoma u tim tačkama odredimo njegove koeficijente, tj. kako da efikasno izvršimo interpolaciju. Ispostavlja se da je problem interpolacije vrlo sličan problemu tabeliranja i da ga rešava praktično isti algoritam.

Primetimo prvo da se tabeliranje tj. izračunavanje vrednosti polinoma u tačkama \(1, w, \ldots, w^{n-1}\) može predstaviti i matrično.

\[\begin{equation*} \begin{pmatrix} 1 & 1 & 1 & \ldots & 1 \\ 1 & w &w^2 & \ldots & w^{n-1} \\ 1 & w^2 &w^{2\cdot2}&\ldots &w^{2\cdot(n-1)} \\ \ldots&\ldots &\ldots & \ldots & \ldots \\ 1 &w^{n-1}&w^{(n-1)\cdot2}&\ldots& w^{(n-1)\cdot(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ \ldots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} P(1) \\ P(w) \\ \ldots \\ P(w^{n-1}) \end{pmatrix} \end{equation*}\]

Uvedimo sledeće oznake:

• vektor koeficijenata polinoma \(P\) obeležimo sa \({{a}}=(a_0,a_1,\ldots,a_{n-1})^T\),
• vektor vrednosti polinoma sa \({{v}}=(P(1),P({w}),\ldots,P({w}^{n-1}))^T\),
• Vandermondovu matricu koja sadrži stepene broja \(w\) obeležimo sa \(V(w)\).

Ako su zadati koeficijenti polinoma \({{a}}\), njegove vrednosti \({{v}}\) u \(n\) tačaka \(1,{w},\ldots,{w}^{n-1}\) dobijaju se izračunavanjem proizvoda:

\[{{v}}=V({w}){{a}}.\]

Ovaj proizvod se naivnim algoritmom množenja matrice i vektora izvršava u vremenu \(O(n^2)\), a direktnom brzom Furijeovom transformacijom u vremenu \(O(n \log{n})\).

S druge strane, ako su zadate vrednosti polinoma \({{v}}=(P(1),P(w),\ldots,P({w}^{n-1}))^T= (v_0,v_1,\ldots,v_{n-1})^T,\) a potrebno je izračunati njegove koeficijente \({{a}}\), matrična jednakost postaje sistem linearnih jednačina po nepoznatim koeficijentima \({{a}}\). Rešavanje sistema jednačina svodi se na računanje inverzne matrice \(V(w)^{-1}\), što se uvek može uraditi Gausovom eliminacijom, što je algoritam složenosti \(O(n^3)\). S druge strane, ako bi se matrica \(V(w)^{-1}\) unapred znala, sistem bi se mogao rešiti množenjem jednakosti \({{v}}=V(w) {{a}}\) s leve strane matricom \(V(w)^{-1}\):

\[{{a}}=V(w)^{-1}{{v}}.\]

U opštem slučaju, ovo je algoritam složenosti \(O(n^2)\). Primetimo, međutim, da Vandermondova matrica \(V(\omega)\) ima veoma pravilnu strukturu koja se može iskoristiti i za efikasnije izračunavanje njoj inverzne matrice i za efikasnije množenje vektora njome.

Rešavanje sistema jednačina \({{v}}= V({w}){{a}}\) svodi se, dakle, na izračunavanje proizvoda \({{a}}=\frac1n V({w}^{-1}){{v}}.\) Posao se dalje pojednostavljuje zahvaljujući sledećoj teoremi.

Pošto brza Furijeova transformacija efikasno računa proizvod neke Vandermondove matrice i vektora, proizvod \(\frac{1}{n} V({w}^{-1}){v}\) može se izračunati primenom brze Furijeove transformacije, zamenjujući u tom algoritmu polaznu vrednost \({w}\) sa \({w}^{-1}\), izračunavanjem proizvoda \(V({w}^{-1}){v}\) i zatim deljenjem komponenti dobijenog vektora sa \(n\). Ova transformacija zove se inverzna Furijeova transformacija (engl. inverse Fourier transform).

Napomenimo da se u primenama obrade signala terminologija razlikuje tako što se menja uloga direktne i inverzne transformacije.

Uzimajući sve do sada rečeno u obzir, proizvod polinoma \(A\) i \(B\) stepena \(n-1\) (za \(n=2^k\)) može se izračunati korišćenjem \(O(n\log n)\) operacija (sa kompleksnim brojevima) na sledeći način:

• računamo vrednost polinoma \(A\) i \(B\) u \(2n-1\) tačaka (stepeni \(2n\)-tog primitivnog korena iz jedinice \(w\)) brzom Furijeovom transformacijom,
• u svakoj od tačaka množimo vrednost polinoma \(A\) i \(B\),
• računamo brzu Furijeovu transformaciju dobijenog vektora, pri čemu umesto vrednosti \(w\) koristimo vrednost \(w^{-1}\) i rezultat množimo sa \(\frac{1}{2n}\).

3.7.3 Algoritam brze Furijeove transformacije

Prikažimo sada algoritam brze Furijeove tranformacije, najpre u pseudokodu (algoritam 7), a zatim i u jeziku C++.

\begin{algorithm}
\caption{Brza Furijeova transformacija}

\begin{algorithmic}
\State{\textbf{Input}: $n$ -- broj koeficijenata polinoma $P$}
\State{\textbf{Input}: $a_0, \ldots, a_{n-1}$ -- koeficijenti polinoma $P$}
\State{\textbf{Input}: $w$ - primitivni $n$-ti koren iz jedinice}
\State{\textbf{Output}: vrednosti $P_0, \ldots, P_{n-1}$ polinoma $P$ na skupu tačaka $1, w, \ldots, w^{n-1}$}

\Procedure{FFT}{$n$, $[a_0, a_1, \ldots a_{n-1}]$, $w$}
\If{$n = 1$}
\State{$P_0 = a_0$}
\Else
\State{$P^{parno}$ = \Call{FFT}{$n/2$, $[a_0, a_2, \ldots, a_{n-2}]$, $w^2$}}
\State{$P^{neparno}$ = \Call{FFT}{$n/2$, $[a_1, a_3, \ldots, a_{n-1}]$, $w^2$}}
\For{$k = 0$ to $n/2 - 1$}
\State{$P_k = P^{parno}_k + w^k \cdot P^{neparno}_k$}
\State{$P_{k + n/2} = P^{parno}_k - w^k \cdot P^{neparno}_k$}
\EndFor          
\EndIf
\EndProcedure
\end{algorithmic}
\end{algorithm}

Jasno je da algoritam realizovan prethodnom procedurom zadovoljava jednačinu \(T(n) = 2T(n/2) + O(n)\), pa mu je složenost \(O(n\log{n})\). Nakon transformacije (promene reprezentacije izračunavanjem vrednosti) dva polinoma njihovo množenje je moguće u vremenu \(O(n)\), pa je ukupna složenost množenja polinoma pomoću algoritma FFT jednaka \(O(n\log{n})\).

Inverzna transformacija se može izračunati korišćenjem prethodnog algoritma tako što se u glavnom pozivu funkcije umesto argumenta \(w\) zada argument \(w^{-1}\) i tako što se na kraju svaki element rezultata podeli sa \(n\).

Razmotrimo sada C++ implementaciju. U jeziku C++ kompleksne brojeve imamo na raspolaganju u obliku tipova complex<double> i complex<float> (zapis u dvostrukoj i jednostrukoj tačnosti). Direktan način implementacije je sledeći.

typedef complex<double> Complex;
typedef vector<Complex> ComplexVector;

// brza Furijeova transformacija vektora a duzine n=2^k
// bool parametar inverse odredjuje da li je direktna ili inverzna
ComplexVector fft(const ComplexVector& a, bool inverzna) {
  // broj koeficijenata polinoma
  int n = a.size();

  // ako je stepen polinoma 0, vrednost u svakoj tacki jednaka
  // je jedinom koeficijentu
  if (n == 1)
    return ComplexVector(1, a[0]);

  // izdvajamo koeficijente na parnim i na neparnim pozicijama
  ComplexVector Pparno(n / 2), Pneparno(n / 2);
  for (int i = 0; i < n / 2; i++) {
    Pparno[i] = a[2 * i];
    Pneparno[i] = a[2 * i + 1];
  }

  // rekurzivno izracunavamo Furijeove transformacije tih polinoma
  ComplexVector fftParno   = fft(Pparno, inverzna), 
                fftNeparno = fft(Pneparno, inverzna);
  
  // objedinjujemo rezultat
  ComplexVector rezultat(n);
  for (int k = 0; k < n / 2; k++) {
    // odredjujemo primitivni n-ti koren iz jedinice
    double koeficijent = inverzna ? -1.0 : 1.0;
    Complex w = exp((koeficijent * 2 * k * M_PI / n) * 1i);
    // racunamo vrednost polinoma u toj tacki
    rezultat[k] = fftParno[k] + w * fftNeparno[k];
    rezultat[k + n/2] = fftParno[k] - w * fftNeparno[k];
  }
  // vracamo konacan rezultat
  return rezultat;
}

// funkcija vrsi direktnu Furijeovu transformaciju polinoma ciji su
// koeficijenti odredjeni nizom a duzine 2^k
ComplexVector fft(const ComplexVector& a) {
  return fft(a, false);
}

// funkcija vrsi inverznu Furijeovu transformaciju polinoma ciji su
// koeficijenti odredjeni nizom a duzine 2^k
ComplexVector ifft(const ComplexVector& a) {
  ComplexVector rezultat = fft(a, true);
  // nakon izracunavanja vrednosti, potrebno je jos podeliti 
  // sve koeficijente duzinom vektora
  int n = a.size();
  for (int k = 0; k < n; k++)
    rezultat[k] /= n;
  return rezultat;
}

Prethodnu direktnu implementaciju je moguće poboljšati. Najveći problem je to što se \(n\)-ti koreni iz jedinice računaju zasebno u svakom rekurzivnom pozivu. Efikasnost bi se značajno popravila ako bi se niz korena izračunao samo jednom, pre funkcije. Problem je u tome što se tokom rekurzije \(n\) smanjuje, pa su nam u svakom pozivu potrebni različiti koreni. Međutim, ako je neki broj \(k\)-ti koren iz jedinice, onda je on i \(2k\)-ti koren iz jedinice. Zato, ako znamo niz \(n\)-tih korena iz jedinice (\(e^{\frac{2k\pi i}{n}}\), za \(k\) od \(0\) do \(n-1\)), tada su elementi na parnim pozicijama tog vektora \(n/2\)-ti koreni iz jedinice. Zaista, ako je \(k=2k'\), tada je \(e^{\frac{2k\pi i}{n}} = e^{\frac{2k'\pi i}{\frac{n}{2}}}\) i važi da ako je \(0 \leq k < n\), tada je \(0 \leq k' < n/2\). Dakle, u početku možemo izračunati niz svih \(n\)-tih korena iz jedinice i njih koristiti na početnom nivou rekurzije, na narednom nivou rekurzije ćemo koristiti svaki drugi element tog vektora, na narednom svaki četvrti, zatim svaki osmi i tako dalje. Na primer, ako je \(n=8\), početni niz korena je \(1\), \(\frac{\sqrt{2}}{2} + i\frac{\sqrt{2}}{2}\), \(i\), \(-\frac{\sqrt{2}}{2} + i\frac{\sqrt{2}}{2}\), \(-1\), \(-\frac{\sqrt{2}}{2} - i\frac{\sqrt{2}}{2}\), \(-i\), \(\frac{\sqrt{2}}{2} - i\frac{\sqrt{2}}{2}\), na narednom nivou rekurzije je \(n=4\) i koriste se koreni \(1\), \(i\), \(-1\), \(-i\), na narednom nivou rekurzije je \(n=2\) i koriste se koreni \(1\) i \(-1\), a na poslednjem nivou rekurzije je \(n=1\) i tada se izlazi iz rekurzije bez korišćenja ovih korena.

Još jedan problem prethodne implementacije je to što se tokom rekurzije alociraju i popunjavaju pomoćni vektori, što dovodi do gubitka i vremena i memorije. Furijeovu transformaciju je moguće realizovati i bez korišćenja pomoćne memorije. Na početnom nivou rekurzije koeficijenti ulaznog polinoma su svi dati početnim vektorima koeficijenata. Na narednom se posmatraju elementi na pozicijama \(0\), \(2\), \(4\), …, \(n-2\) i na pozicijama \(1\), \(3\), …, \(n-1\). Na narednom se posmatraju elementi na pozicijama \(0\), \(4\), \(8\), \(n-4\), zatim elementi na pozicijama \(1\), \(5\), …, \(n-3\), zatim elemeni na pozicijama \(2\), \(6\), …, \(n-2\), i na kraju elementi na pozicijama \(3\), \(7\), …, \(n-1\). Slično se nastavlja i na daljim nivoima rekurzije. Dakle, umesto formiranja pomoćnog ulaznog vektora sa pogodno odabranim ulaznim koeficijentima, prosleđivaćemo originalni vektor, poziciju početka \(s\) i pomeraj \(d\) i posmatraćemo njegove elemente na pozicijama \(s + dk\), za \(0 \leq k < n\), gde je \(n = n_0 / d\), a \(n_0\) je dužina početnog vektora. Rezultate rekurzivnih poziva možemo smestiti u dve polovine rezultujućeg niza. Nakon toga rezultate objedinjujemo. Vrednosti na pozicijama \(k\) i \(k+n/2\) rezultujućeg vektora određene su vrednostima na poziciji \(k\) u rezultatu prvog i drugog rekurzivnog, međutim, one se nalaze upravo na pozicijama \(k\) i \(k+n/2\) rezultujućeg vektora (jer smo rezultate rekurzivnih poziva smestili u prvu i drugu polovinu rezultata). Moramo voditi računa da te dve vrednosti moramo istovremeno izračunati i ažurirati.

typedef complex<double> Complex;
typedef vector<Complex> ComplexVector;

// funkcija vrši Furijeovu transformaciju (direkntu ili inverznu) elemenata
//   a[start_a], a[start_a + korak], a[start_a + 2korak], ...
// i rezultat smešta u niz
//   rezultat[start_rezultat], rezultat[start_rezultat + 1], rezultat[start_rezultat + 2], ...
// koristeci primitivne korene iz jedinice smestene u niz
//   w[0], w[korak], w[2*korak], ...
void fft(const ComplexVector& a, int start_a,
         const ComplexVector& w, 
         ComplexVector& rezultat, int start_rezultat,
         int korak) {
  // broj elemenata niza koji se transformise
  int n = a.size() / korak;

  // stepen polinoma je nula, pa mu je vrednost u svakoj tacki jednaka
  // konstantnom koeficijentu
  if (n == 1) {
    rezultat[start_rezultat] = a[start_a];
    return;
  }

  // rekurzivno transformisemo niz koeficijenata na parnim pozicijama
  // smestajuci rezultat u prvu polovinu niza rez
  fft(a, start_a, w, rezultat, start_rezultat, korak*2);
  // rekurzivno transformisemo niz koeficijenata na neparnim pozicijama
  // smestajuci rezultat u drugu polovinu niza rez
  fft(a, start_a + korak, w, rezultat, start_rezultat + n/2, korak*2);

  // objedinjujemo dve polovine u rezultujuci niz
  for (int k = 0; k < n/2; k++) {
    Complex r1 = rezultat[start_rezultat + k];
    Complex r2 = rezultat[start_rezultat + (k + n/2)];
    rezultat[start_rezultat + k] = r1 + w[k*korak] * r2;
    rezultat[start_rezultat + (k + n/2)] = r1 - w[k*korak] * r2;
  }
}

// funkcija vrsi direktnu Furijeovu transformaciju polinoma ciji su
// koeficijenti odredjeni nizom a duzine 2^k
ComplexVector fft(const ComplexVector& a) {
  // duzina niza koeficijenata polinoma
  int n = a.size();
  // izracunavamo primitivne n-te korene iz jedinice
  ComplexVector w(n/2);
  for (int k = 0; k < n/2; k++)
    w[k] = exp((2 * k * M_PI / n) * 1i);
  // vektor u koji se smesta rezultat
  ComplexVector rezultat(n);
  // vrsimo transformaciju
  fft(a, 0, w, rezultat, 0, 1);
  // vracamo dobijeni rezultat
  return rezultat;
}

// funkcija vrsi inverznu Furijeovu transformaciju polinoma ciji su
// koeficijenti odredjeni nizom a duzine 2^k
ComplexVector ifft(const ComplexVector& a) {
  // duzina niza koeficijenata polinoma
  int n = a.size();
  // izracunavamo primitivne n-te korene iz jedinice
  ComplexVector w(n);
  for (int k = 0; k < n; k++)
    w[k] = exp((- 2 * k * M_PI / n) * 1i);
  // vektor u koji se smesta rezultat
  ComplexVector rezultat(n);
  // vrsimo transformaciju
  fft(a, 0, w, rezultat, 0, 1);
  // popravljamo rezultat 
  for (int i = 0; i < n; i++)
    rezultat[i] /= n;
  // vracamo dobijeni rezultat
  return rezultat;
}

vector<double> proizvod(const vector<double>& p1,
                        const vector<double>& p2) {
  // duzina niza koeficijenata
  int n = p1.size();
  // kreiramo nizove kompleksnih koeficijenata dopunjavajuci ih nulama
  // do dvostruke duzine
  int N = 2*n;
  ComplexVector a(N, 0.0);
  copy(begin(p1), end(p1), begin(a));
  ComplexVector b(N, 0.0);
  copy(begin(p2), end(p2), begin(b));

  // vrsimo Furijeove transformacije oba vektora koeficijenata 
  // (slozenost je O(n log(n)))
  ComplexVector va = fft(a), vb = fft(b);
  // mnozimo vrednosti u pojedinacnim tackama (slozenost je O(n))
  ComplexVector vc(N);
  for (int i = 0; i < N; i++)
    vc[i] = va[i] * vb[i];
  // inverznom Furijeovom transformacijom rekonstruisemo koeficijente proizvoda
  // (slozenost je O(n log(n))
  ComplexVector c = ifft(vc);

  // realne delove kompleksnih brojeva smestamo u niz rezultat i vracamo ga
  vector<double> rezultat(N);
  transform(begin(c), end(c), begin(rezultat),
            [](Complex x) { return real(x); });
  return rezultat;
}

Pažljivom analizom je moguće ukloniti rekurziju iz prethodne implementacije (tako se dobija Kuli-Tukijev nerekurzivni algoritam za FFT, zasnovan na ranije prikazanoj leptir-shemi). Ta implementacija neće biti prikazana u ovom udžbeniku.

3.7.4 NTT: Furijeova transformacija u modularnoj aritmetici

U algoritmu FFT prelaz sa realnih brojeva na kompleksne je bio neophodan zato što kompleksni brojevi garantuju to da za svako \(n\) postoji primitivan \(n\)-ti koren iz jedinice, tj. to da svaka jednačina oblika \(w^n = 1\) ima \(n\) različitih rešenja, što kod realnih brojeva nije slučaj već za \(n > 2\). Kao ni realni, ni kompleksni brojevi ne mogu biti potpuno precizno predstavljeni u računaru, već se koristi približan zapis (najčešće zapis u pokretnom zarezu), što može dovesti do računskih grešaka. Umesto kompleksnih brojeva moguće je koristiti i neke druge brojevne strukture. Na primer, moguće je koristiti modularnu aritmetiku i konačna polja.

U modularnoj aritmetici \(n\)-ti koren iz jedinice po modulu \(p\) je broj \(w\) takav da važi \(w^n \equiv 1\ (\mathrm{mod}\ p)\). Ako za svako \(1 \leq m < n\) važi \(w^m \not\equiv 1\ (\mathrm{mod}\ p)\), tada je \(w\) primitivan \(n\)-ti koren iz jedinice. Prirodno se postavlja pitanje da li za dati prost broj \(p\) i dati stepen \(n\) postoji primitivni \(n\)-ti koren iz jedinice.

Na osnovu prethodne leme sledi da u modularnim grupama primitivni \(n\)-ti koreni ne postoje za svako \(n\). Međutim, imamo slobodu izbora modula \(p\) na osnovu poznatog stepena polinoma \(n\). Naime, uvek imamo mogućnost izbora prostog broja \(p\) na osnovu poznate vrednosti \(n\) tako da uslov prethodne teoreme bude ispunjen (tj. da je \(p-1\) deljivo sa \(n\)). Pošto nas zanimaju vrednosti \(n\) koje su stepeni dvojke potrebno je da se pronađe neki prost broj \(p\) takav da je \(p-1\) deljivo sa što više stepena dvojke (onda taj isti broj možemo koristiti za različite vrednosti \(n\) koje su stepeni dvojke). Na primer, jedan takav prost broj je \(p = 2^{23}\cdot 7 \cdot 17 + 1 = 998\,244\,353\) i on ima primitivni \(n\)-ti koren za sve stepene dvojke \(n\) do \(2^{23}\) (što je u mnogim primenama i više nego dovoljno, jer nam omogućava množenje polinoma sa po oko \(8\) miliona koeficijenata koji su u opsegu do skoro milijardu).

Ostaje pitanje kako za datu vrednost \(n\) efektivno pronaći primitivni \(n\)-ti koren iz jedinice po modulu \(p\). Pokažimo jedan efektivni postupak za to. Pretpostavimo da je \(p-1 = kn\). Iz dokaza leme 3.7.4 jasno je da je ključni problem pronaći generator \(a\) grupe \((\mathbb{Z}_p, \times_p)\). tj. element \(a\in Z_P \setminus\{0\}\) takav da je \(a^{p-1} \equiv 1\ (\mathrm{mod}\ p)\) (što uvek važi), ali ni za jedno \(1\leq m < p-1\) ne važi \(a^m \equiv 1\ (\mathrm{mod}\ p)\). Da bi ovaj uslov važio za neko \(m\) potrebno je da \(m\) bude delilac broja \(p-1\). Dovoljno je, dakle, da proverimo da je \(a^m \not\equiv 1\ (\mathrm{mod}\ p)\) za sve delioce \(m\) broja \(p-1\). Ovo se može dalje redukovati time što se ispitaju samo oni delioci koji su oblika \(\frac{p-1}{f_i}\), gde je \(f_i\) neki prost činilac broja \(p-1\). Naime, ovi delioci su maksimalni u mreži deljivosti (za svaki od njih, naredni delilac je upravo \(p-1\)) i ako za neki delilac \(m\) ispred njih važi \(a^m \equiv 1\ (\mathrm{mod}\ p)\), tada će bar neki od njih, tj. bar neki delilac oblika \(\frac{p-1}{f_i}\) biti deljiv sa \(m\) i važiće \(a^{\frac{p-1}{f_i}} \equiv 1\ (\mathrm{mod}\ p)\).

Kada je poznat primitivni \(n\)-ti koren iz jedinice \(w\), izvođenje Furijeove transformacije teče na potpuno isti način, jedino što se umesto kompleksnih brojeva koriste prirodni brojevi i sve operacije se izvršavaju po modulu \(p\). Ovaj oblik Furijeove transformacije se nekada naziva NTT (engl. number theoretic transform). Ona može da ima prednosti u odnosu na klasičan pristup sa kompleksnim brojevima, jer ne postoje problemi sa preciznošću zapisa broja.

Zadaci:

• Poravnavanje niski
• Digitalni brojač