Mnozenje polinomov z metodo FFT

Iz MaFiRaWiki

Vsebina

[spremeni]

Hitra Fouriejeva transformacija

Hitra Fourierjeva transformacija (FFT) je učinkovit algoritem za izračun diskretne Fouriejeve transformacije (DFT) in njenega inverza. FFT igra pomembno vlogo pri mnogih aplikacijah, vse od digitalnega procesiranja signalov, reševanja parcialnih diferencialnih enačb, do algoritmov za hitro množenje velikih števil.

Naj bodo x0, ...., xN-1 kompleksna števila. DFT je definirana kot:

X_k = \sum_{n=0}^{N-1} x_n e^{-{2\pi i \over N} nk } \qquad k = 0,\dots,N-1.

Če bi te vsote izračunali direktno, bi za to potrebovali [[O(N2)]] aritmetičnih opreacij. FFT je algoritem, ki izračuna isto vsoto s samo O(N log N) operacijami.

Veliko FFT algoritmov sloni na dejstvu, da je e^{-{2\pi i \over N}} praštevilski koren enote in ga zato lahko uporabimo v analognih transformacijah nad vsakih končnim poljem, kot številske transformacije.

Ker je inverz DFT enak DFT, ampak z različnim predznakom v eksponentu in 1/N faktor, lahko vsak FFT algoritem zlahka uporabimo tudi pri tem.

Inverz DFT dobimo po formuli:

x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{\frac{2\pi i}{N} k n} \quad \quad n = 0,\dots,N-1.
[spremeni]

Množenje polinomov z uporabo FFT

Naprimer, da želimo izračunati produkt polinomov:c(x) = a(x) · b(x).

Pri standardnem izračuno le tega uporabimo linearno konvolucijo. To lahko drugače zapišemo, če vzamemo najprej vektorja koeficientov za a(x) in b(x), tako da vzamemo najprej konstanten člen, nato pa dodoajamo ničle, tako da imata na koncu vektorja a in b dimenziji d > deg(a(x)) + deg(b(x)). Nato:

\mathbf{c} = \mathbf{a} * \mathbf{b}

kjer je c vektor koeficientov c(x) in konvolucijski operator *\, definiran kot

c_n = \sum_{m=0}^{d-1}a_m b_{n-m\ \mathrm{mod}\ d} \qquad\qquad\qquad n=0,1,...,d-1

Ampak konvolucija postane množenje, če uporabimo DFT:

\mathcal{F}(\mathbf{c}) = \mathcal{F}(\mathbf{a})\mathcal{F}(\mathbf{b})

Tukaj vzamemo vektorski produkt po elementih. Tako postanejo koeficienti produkta polinomov c(x) preprosto členi 0, ..., deg(a(x)) + deg(b(x)) vektorja koeficientov.

\mathbf{c} = \mathcal{F}^{-1}(\mathcal{F}(\mathbf{a})\mathcal{F}(\mathbf{b})).
[spremeni]

Metoda za množenje polinomov s FFT

Ideja za množenje je naslednja: izberemo največji w, tako da ne bo prišlo do prekoračitve v nadaljnem procesu, ki je opisan spodaj. Nato razdelimo naši dve števili v skupini z w biti:

a=\sum_{i=0}^m {2^{wi}a_i} and b=\sum_{j=0}^m {2^{wj}b_j}.

Potem lahko rečemo:

ab=\sum_{i=0}^m \sum_{j=0}^m 2^{w(i+j)}a_i b_j = \sum_{k=0}^{2m} 2^{wk} \sum_{i=0}^m a_i b_{k-i} = \sum_{k=0}^{2m} 2^{wk} c_k

če postavimo bj=0 za j > m, k=i+j in {ck} kot konvolucija {ai} in {bj}. Z uporabo konvolucijskega izreka lahko izračunamo ab:

  1. Izračunamo FFT {ai} in {bj},
  2. Pomnožimo rezultata (enega za drugim)
  3. Izračunamo inverz FFT
  4. Dodamo del ck ki je večji kot 2w do ck+1
[spremeni]

Implementacija v Javi

Predstavili bomo enostavno implementacijo v Javi. Najprej definiramo konstruktorje, ki jih bomo uporabili pri metodi. Več konsturktorjev ustvarimo zato, da s tem poenostavimo klicanje kasneje v programu.

Metodo FFT konstruiramo po formuli zgoraj: 

  1.  
  2. public class FFT
  3. {
  4.   public static Kompleksno[] FFT(Kompleksno[] x)
  5. 	   {
  6. 	      return FFT(x, velikost(x.length));
  7.    }
  8.   public static Kompleksno[] FFT(Kompleksno[] x, int n)
  9. 	   {
  10. 	      Kompleksno w = new Kompleksno (n);
  11. 	      return FFT(x, n, w);
  12.    }
  13.  
  14.    public static Kompleksno[] FFT(Kompleksno[] x, int n, Kompleksno z)
  15.    {
  16.       Kompleksno[] y = new Kompleksno[n];
  17.  
  18.       if (n == 1) {
  19.          y[0] = x[0];
  20.       }
  21.

Definiramo dve tabeli kompleksnih števil (tukaj poimenovani lihi in sodi), ki 

  1.  
  2.       else {
  3.          Kompleksno[] sodi = new Kompleksno[n/2];
  4.          Kompleksno[] lihi = new Kompleksno[n/2];
  5.          for (int i = 0; i <n/2; ++i) {
  6.             lihi[i] = (2 * i + 1 < x.length) ? x[2 * i + 1] : new Kompleksno (0.0);
  7.             sodi[i] = (2 * i < x.length) ? x[2 * i] : new Kompleksno (0.0);
  8.  
  9.          }
  10.  
  11.            lihi = FFT(lihi, n/2, z.kvadrat());
  12.            sodi = FFT(sodi, n/2, z.kvadrat());
  13.  
  14.  
  15.          Kompleksno z1 = new Kompleksno(1.0);
  16.  
  17.          for (int i = 0; i<n/2; ++i) {
  18.             Kompleksno t = z1.produkt(lihi[i]);
  19.  
  20.             y[i] = sodi[i].vsota(t);
  21.             y[i + n/2] = sodi[i].razlika(t);
  22.             z1 = z1.produkt(z);
  23.          }
  24.       }
  25.  
  26.       return y;
  27.    }
  28.  
  29.    public static int velikost(int n)
  30.    {
  31.       --n;
  32.       int k = 0;
  33.       while (n > 0) {
  34.          n >>= 1;
  35.          ++k;
  36.       }
  37.       return 1 << k;
  38.    }
  39.


Sedaj, ko smo ustvarili metodo za FFT, potrebujemo še njegov inverz. To je enostavno za izračunati, saj je glavna razlika v predzankih 

  1.  
  2.  
  3.    public static Kompleksno[] invFFT(Kompleksno[] y)
  4.    {
  5.       return invFFT(y, velikost(y.length));
  6.    }
  7.  
  8.    public static Kompleksno[] invFFT(Kompleksno[] y, int n)
  9.    {
  10.       Kompleksno z = new Kompleksno (n);
  11.  
  12.       Kompleksno[] x = FFT(y, n, z.inverz());
  13.  
  14.       for (int i = 0; i < n; ++i) x[i].deljeno(n);
  15.       return x;
  16.    }
  17.

Ko imamo enkrat metodo za FFT in njen inverz, je enostavno izračunati produkt dveh polinomov, uporabimo isto formulo kot zgoraj. 

  1.  
  2.  
  3.    public static Kompleksno[] produkt(Kompleksno[] p, Kompleksno[] q)
  4.    {
  5.       int n = velikost(p.length + q.length - 1);
  6.       Kompleksno[] p_1 = FFT(p, n);
  7.       Kompleksno[] q_1 = FFT(q, n);
  8.       Kompleksno[] pq = new Kompleksno[n];
  9.       for (int i = 0; i < n; ++i) pq[i] = p_1[i].produkt(q_1[i]);
  10.  
  11.       return invFFT(pq, n);
  12.    }
  13.       
  14. }
  15.
Vzpostavljeno iz "http://wiki.fmf.uni-lj.si/wiki/Mnozenje_polinomov_z_metodo_FFT"
Osebna orodja