Trasformate notevoli ed esponenziazione di polinomi con la FFT

Due esercitazioni sulla FFT. Nella prima si sfrutta la struttura del vettore per calcolare la DFT con meno operazioni della FFT generale (Fast Fourier Transform - algoritmo ricorsivo e trasformata inversala FFT calcola DFT_n (n potenza di due) con divide et impera: si separano i coefficienti di indice pari e dispari, p(x) = p0(x^2) + x p1(x^2), e per il lemma di dimezzamento bastano due DFT di ordine n/2; y_i = y0_i + w^i y1_i e y_{i+n/2} = y0_i - w^i y1_i. La ricorrenza T(n) = 2T(n/2) + 2n con T(1) = 0 ha soluzione esatta 2 n log2 n. L'inversa ha matrice F^{-1}_{ij} = w^{-ij}/n: si ottiene con la FFT sulle potenze negative (poi si divide per n fuori dalla ricorsione) oppure con DFT_n(y) rovesciato. LIN_CONV moltiplica due polinomi con n coefficienti in 12 n log2(2n) + 4n operazioni.Fast Fourier Transform - algoritmo ricorsivo e trasformata inversa →); nella seconda si usa la rappresentazione per valori per calcolare i coefficienti della potenza $p(x)^k$ di un polinomio. Il metodo è lo stesso del capitolo: si modifica la proprietà di sottostruttura della FFT con le informazioni aggiuntive sull'istanza, e si ricava la nuova ricorrenza. Una situazione analoga (struttura che abbassa il costo) è nelle matrici circolanti (Matrici circolanti e supercircolanti - algoritmi divide and conquerle matrici supercircolanti (quadranti alternati ricorsivi) hanno solo n valori distinti: somma in n operazioni, prodotto in 4n log2 n + n con il trucco (A1+A2)(B1+B2) e (A1-A2)(B1-B2) invece delle 4 chiamate (Theta(n^2)); gli autovalori sono la trasformata di Hadamard della prima riga (FHT, n log2 n) e il prodotto passa dalla base diagonalizzante; le matrici circolanti sono polinomi nella matrice di rotazione R, il loro prodotto è la convoluzione ciclica delle prime righe e si diagonalizzano con la matrice di Fourier: prodotto in O(n log n) con la FFT.Matrici circolanti e supercircolanti - algoritmi divide and conquer →).

Trasformate note a priori

La componente $i$ di $\mathbf y=F_n\mathbf a$ è il prodotto scalare della riga $i$ di $F_n$ con $\mathbf a$. La prima riga di $F_n$ è fatta di soli $1$.

Vettore costante. Se $\mathbf a=(a,a,\dots,a)\in\mathbb C^n$ allora $\mathrm{DFT}_n(\mathbf a)=(n\,a,\,0,\dots,0)$.

Dimostrazione. $y_0=\sum_ja=na$. Per $0<i\le n-1$, $y_i=a\sum_{j=0}^{n-1}(\omega_n^i)^j=0$ per il lemma di somma, perché $i$ non è un multiplo di $n$. $\square$

Vettore con la sola prima componente. Se $\mathbf a=(a,0,\dots,0)$ allora $\mathrm{DFT}_n(\mathbf a)=(a,a,\dots,a)$.

Dimostrazione. $\mathbf a=a\,\mathbf e_0$, con $\mathbf e_0$ la prima colonna di $I_n$, e $F_n\mathbf e_0$ è la prima colonna di $F_n$, tutta di $1$. In termini di polinomio: $p(x)=a$ vale $a$ in ogni punto. $\square$

Più in generale la DFT di $\mathbf e_j$ è la colonna $j$ di $F_n$, cioè $(\omega_n^{ij})_i$.

Esempio ($n=4$). $\mathrm{DFT}_4(7,7,7,7)=(28,0,0,0)$ e $\mathrm{DFT}_4(7,0,0,0)=(7,7,7,7)$ (controllati con numpy).

Vettori $(n,k)$-sparsi

Definizione. Siano $n$ e $k$ potenze di due, $1\le k\le n$. Un vettore $\mathbf x\in\mathbb C^n$ è $(n,k)$-sparso se $x_i=0$ per ogni indice $i$ con $i\bmod k\ne0$.

Le componenti possono essere diverse da zero solo agli indici multipli di $k$ (sono $n/k$, separate da $k-1$ zeri; non si esclude che anche quelle siano nulle). Esempio: un vettore $(8,4)$-sparso ha la forma $(x_0,0,0,0,x_4,0,0,0)$. Casi estremi: ogni vettore è $(n,1)$-sparso; un vettore $(n,n)$-sparso ha al più la sola componente $x_0$ (la sua DFT è nota).

Obiettivo. Un algoritmo D&C che calcoli $\mathrm{DFT}_n(\mathbf x)$ per $\mathbf x$ $(n,k)$-sparso con $T(n,k)=O\!\left(\frac nk\log\frac nk\right)$ operazioni: per $k=1$ si ritrova $O(n\log n)$, per $k=n$ nessuna operazione.

Sottostruttura rinforzata. Siano $\mathbf x^{[0]}$ e $\mathbf x^{[1]}$ i vettori degli indici pari e dispari di $\mathbf x$, come nella FFT. Se $1<k<n$ (cioè $k\ge2$ pari):

• $\mathbf x^{[0]}$ è $(n/2,k/2)$-sparso: $x^{[0]}_j=x_{2j}$ è nullo se $2j\bmod k\ne0$, cioè se $j\bmod(k/2)\ne0$;
• $\mathbf x^{[1]}$ è il vettore nullo: $x^{[1]}_j=x_{2j+1}$ ha indice dispari, che non può essere un multiplo della potenza di due $k\ge2$.

Quindi $\mathbf y^{[1]}=\mathrm{DFT}_{n/2}(\mathbf 0)=\mathbf 0$ e le formule di ricombinazione della FFT diventano $y_i=y^{[0]}_i$ e $y_{i+n/2}=y^{[0]}_i$: la DFT è la concatenazione di $\mathbf y^{[0]}$ con se stessa, senza nessuna operazione.

SPARSE_FFT(x, k)
n <- x.len
if k = 1 then return FFT(x)
if k = n then return (x[0], x[0], ..., x[0])        { n copie }
x0 <- (x[0], x[2], ..., x[n-2])
y0 <- SPARSE_FFT(x0, k/2)
for i <- 0 to n/2 - 1
    y[i] <- y0[i];  y[i + n/2] <- y0[i]
return y

Correttezza. Segue dalla correttezza della FFT e dalla sottostruttura rinforzata, per induzione su $k$.

Costo. Nel conquer ci sono solo assegnamenti, che non si contano come operazioni aritmetiche. La ricorrenza a due parametri è $$T(n,k)=\begin{cases}2n\log_2n&k=1\ (\text{FFT}),\\0&k=n,\\T\!\left(\frac n2,\frac k2\right)&1<k<n.\end{cases}$$ Srotolando (Albero della ricorrenza e formula generale per le ricorrenzeL'albero della ricorrenza ha la forma dell'albero delle chiamate dell'istanza peggiore e T(n) è la somma dei lavori dei nodi. Per la famiglia T(n) = s(n)T(f(n)) + w(n) con base T0 per n <= n0, contando i nodi per livello con le iterate di f e il livello f*(n,n0) si ottiene la formula chiusa T(n) = somma_l [prodotto_j s(f^(j)(n))] w(f^(l)(n)) + [prodotto s] T0, verificata per induzione e applicata a esempi (merge sort, Hanoi, ricorrenze con radice).Albero della ricorrenza e formula generale per le ricorrenze →) per $1<k<n$: $T(n,k)=T(\frac n2,\frac k2)=T(\frac n4,\frac k4)=\dots=T(\frac nk,1)=2\,\frac nk\log_2\frac nk=\Theta\!\left(\frac nk\log\frac nk\right)$. L'ordine di grandezza vale anche per $k=1$ e per $k=n$ ($\log1=0$). Se si contasse un'unità per ogni assegnamento, il lavoro sarebbe lineare e si otterrebbe $O(\frac nk\log\frac nk+n)$. Il conteggio $2\frac nk\log_2\frac nk$ è stato verificato con un'implementazione per tutte le coppie $n=2^0,\dots,2^7$, $k=2^0,\dots,n$.

Conseguenza. La DFT di un vettore $(n,k)$-sparso è formata da $k$ copie della $\mathrm{DFT}_{n/k}$ del vettore $(x_0,x_k,x_{2k},\dots,x_{n-k})$ che contiene le sue $n/k$ componenti non banali. Lo si vede anche direttamente: $y_i=\sum_{m=0}^{n/k-1}x_{mk}\,\omega_n^{imk}=\sum_mx_{mk}\,\omega_{n/k}^{im}$ (cancellazione), che è periodica in $i$ con periodo $n/k$.

Esempio ($n=8$, $k=4$). $\mathbf x=(2,0,0,0,1,0,0,0)$. Le componenti non nulle $(2,1)$ hanno $\mathrm{DFT}_2(2,1)=(2+1,\,2-1)=(3,1)$, ripetuta $4$ volte: $\mathrm{DFT}_8(\mathbf x)=(3,1,3,1,3,1,3,1)$.

Esempio d'esame. $\mathbf x=(5,3,3,3,4,3,3,3)$, $\mathbf y=\mathrm{DFT}_8(\mathbf x)$. Per linearità ($F_n$ è lineare) si scompone $\mathbf x=(3,3,3,3,3,3,3,3)+(2,0,0,0,1,0,0,0)$. Il primo addendo è costante, con DFT $(24,0,\dots,0)$; il secondo è $(8,4)$-sparso con DFT $(3,1,3,1,3,1,3,1)$. Somma: $\mathbf y=(27,1,3,1,3,1,3,1)$ (verificato con la DFT calcolata dalla definizione).

La strategia, ripetibile in esercizi simili, è: guardare la struttura del vettore, scomporlo (anche con la linearità) in pezzi con trasformata nota, oppure riscrivere la sottostruttura della FFT usando la struttura, e studiare la ricorrenza.

Esponenziazione di polinomi

Dato $p(x)=\sum_{i<n}a_ix^i\equiv\mathbf a\in\mathbb C^n$ e un intero $k\ge1$, si vogliono i coefficienti di $q(x)=p(x)^k$, cioè $\mathbf a^{(k)}=\mathbf a\star\mathbf a\star\dots\star\mathbf a$ ($k$ volte). Il limite di grado di $q$ è $(n-1)k+1$.

Primo approccio: quadrature successive con la convoluzione

La convoluzione è associativa, quindi si applica lo schema di Esponenziazione veloce per quadrature successivepotenza x^n con un operatore associativo con elemento neutro, costo = numero di applicazioni dell'operatore; STUPID_POWER costa n-1, la sottostruttura x^n = (x^floor(n/2))^2 (per n dispari per x) dà POWER con T(n) = floor(log2 n) + (numero di bit a 1 di n) - 1, Theta(log n) (caso peggiore 2 floor(log2 n), n = 2^t - 1); versioni iterative POWER2 e POWER3 dalla rappresentazione binaria di n; esempio con la matrice di Fibonacci; il modello a costo unitario vale solo se l'operatore non fa crescere gli operandi.Esponenziazione veloce per quadrature successive → con $\star$ al posto del prodotto. Siano $d=\lfloor\log_2k\rfloor+1$ e $k_{d-1}\dots k_0$ le cifre binarie di $k$, $k=\sum_ik_i2^i$. Allora $\mathbf a^{(k)}=\mathbf a^{(k_02^0)}\star\mathbf a^{(k_12^1)}\star\dots\star\mathbf a^{(k_{d-1}2^{d-1})}$.

EXPONENTIATE(a, k)
d <- floor(log2 k) + 1
h <- k                         { h = (k_{d-1} ... k_0) in base 2 }
y <- a                         { all'iterazione i: y = a^(2^i) }
z <- (1)                       { elemento neutro della convoluzione; accumula il prodotto }
for i <- 0 to d - 1
    if h è dispari then z <- LIN_CONV(z, y)       { cifra k_i = 1 }
    y <- LIN_CONV(y, y)                            { y = a^(2^(i+1)) }
    h <- floor(h / 2)
return z

LIN_CONV va nella versione per vettori di lunghezza qualsiasi (padding alla potenza di due). Analisi. All'inizio dell'iterazione $i$, y ha lunghezza al più $2^in$ (la lunghezza di $\mathbf a^{(2^i)}$ è $2^i(n-1)+1$) e z ha lunghezza minore di $2^in$ (è la convoluzione di fattori di lunghezze $\le 2^jn$, $j<i$, e $\sum_{j<i}2^j<2^i$). Quindi l'iterazione $i$ costa $T_i=O(2^in\log(2^in))=O(2^in\log n+2^in\,i)$, e poiché $2^d\le2k$: $$T(n,k)=\sum_{i=0}^{d-1}T_i=O\bigl(n2^d\log n+n2^dd\bigr)=O\bigl(nk\log n+nk\log k\bigr)=O\bigl(nk\log(nk)\bigr).$$

Secondo approccio: potenza dei valori

Se si ha una rappresentazione per valori di $p$ su una base estesa $\mathbf x^E$ con $h\ge(n-1)k+1$ punti, i valori di $q=p^k$ sono le potenze dei valori di $p$: $q(x_i)=\bigl(p(x_i)\bigr)^k$, cioè $\mathbf y_q=\bigl((y_{p,i})^k\bigr)_i$. Basta quindi valutare, elevare e interpolare.

EXPONENTIATE2(a, k)
n <- a.len
h <- 2^ceil(log2((n-1)*k + 1))              { prima potenza di due >= (n-1)k+1 }
yp <- FFT((a | 0_{h-n}))
for i <- 0 to h - 1
    yq[i] <- POWER3(yp[i], k)                 { potenza scalare per quadrature, O(log k) }
return INV_FFT(yq)

La correttezza segue dalle considerazioni sopra: $q$ ha limite di grado $\le h$, quindi $\mathbf y_q$ lo determina e INV_FFT restituisce i suoi coefficienti (seguiti da zeri di padding). Costo: FFT e INV_FFT su vettori di taglia $h<2((n-1)k+1)\le2nk$, quindi $\Theta(nk\log(nk))$; più $h$ potenze scalari ciascuna di $\Theta(\log k)$ prodotti, in tutto $\Theta(nk\log k)$. Il totale è $\Theta(nk\log(nk)+nk\log k)=\Theta(nk\log(nk))$. Stesso ordine del primo approccio, ma l'algoritmo è più semplice e con costanti migliori: una trasformata diretta, una inversa e $h$ potenze, invece di $O(\log k)$ convoluzioni di vettori sempre più lunghi.

Esempio ($p(x)=1+x$, $k=4$). $n=2$, il limite di grado di $q$ è $(2-1)\cdot4+1=5$, quindi $h=2^{\lceil\log_25\rceil}=8$. $\mathbf a\,|\,0=(1,1,0,0,0,0,0,0)$ e $y_{p,i}=1+\omega_8^i$. Le potenze quarte:

(per esempio $(1+\omega_8)^2=1+2\omega_8+i=(1+\sqrt2)(1+i)$ e quindi $(1+\omega_8)^4=(3+2\sqrt2)\cdot2i=(6+4\sqrt2)\,i$.) L'antitrasformata di questo vettore è $(1,4,6,4,1,0,0,0)$, cioè $(1+x)^4=1+4x+6x^2+4x^3+x^4$, con i tre zeri finali di padding. Verificato con numpy insieme a $20$ casi casuali per ciascun approccio.

Domande d'esame

• Calcolare $\mathrm{DFT}_8(5,3,3,3,4,3,3,3)$ senza eseguire la FFT. Traccia: $(3,\dots,3)+(2,0,0,0,1,0,0,0)$, DFT $(24,0,\dots,0)+(3,1,3,1,3,1,3,1)=(27,1,3,1,3,1,3,1)$.
• Quante operazioni esegue SPARSE_FFT su un vettore $(n,k)$-sparso con $1<k<n$ e perché? Traccia: nessuna nel conquer, $T(n,k)=T(n/2,k/2)=\dots=T(n/k,1)$, quindi $2\frac nk\log_2\frac nk$.
• Per elevare alla sesta potenza un polinomio con 5 coefficienti usando la FFT, su quanti punti si valuta? Traccia: $(n-1)k+1=4\cdot6+1=25$, $h=2^{\lceil\log_225\rceil}=32$.

Esercizi svolti: Esercizio - polinomio con tutte le radici dell'unità tranne 1 e -1, valutazione e DFT dei coefficienti (compito del 26 giugno 2025, es. 1); Esercizio - domande di prima parte su polinomi, DFT e FFT (simulazioni d'esame); Esercizio - matrice di Fourier per una matrice generica a sinistra e a destra e per aI + bU (compito del 27 gennaio 2025, es. 2).