Fast Fourier transform (FFT)

Fast Fourier transform (FFT)

Hang a blog
Main idea: Divide and conquer
Main purpose: optimize the time complexity of two polynomial multiplication ( (O (n^2)->o (NLOGN)) \)

Front-facing

For a polynomial\[a (x) =\sum_{i=0}^{n-1}a_ix^i\]Think of it as a function that has\ (y_k=a (x_k) \), the point on the image\ ((x_k,y_k) \)
theorem: Here you\ (n\)A different point, there must be a unique\ (n\)Sub-polynomial\ (A (x) \)
The complexity of multiplying two points is\ (O (1) \)Of
Now consider specifying the polynomial\ (a,b\)The horizontal axis of each n point, the Ordinate (evaluation) is obtained,
and multiply it up to get new\ (n\)Point, calculate\ (a,b\)Product of\ (c\)(interpolated value)
First, regardless of interpolation, evaluation or\ (O (n^2) \)Of
So we have to choose some ghost animals of the Axis speed evaluation process,

That is\ (n\)Secondary unit root\ (w\)Meet\ (w^n=1\)And finds that it happens to be in a plural field.\ (n\)And all on one unit circle
\[w_n^k=e^{2πki/n}\]
which\ (W_n^1=e^{2πi/n}=cos (2π/n) +i\) \ (sin (2π/n) \)is the main\ (n\)Secondary unit root (abbreviated as\ (w_n\)),
All other\ (n\)The second unit root is its entire power.
three large lemma
To eliminate the lemma:\[w_{dn}^{dk}=w_n^k\]
Prove:\[w_{dk}^{dn}=e^{2dkπi/dn}=w_n^k\]
Binary lemma:\[(W_N^{K+N/2}) ^2= (w_n^k) ^2\]
Prove:\[(W_N^{K+N/2}) ^2=w_n^{2k+n}=w_n^nw_n^{2k}= (w_n^k) ^2\]
which\ (w_n^n=cos2π+i\) \ (sin2π=1\)

Summation lemma:
For any integer \ (n≥1\) and non-negative integer ( k\)that cannot be divisible by \ (n\) , there are:
\[\sum_{j=0}^{n-1} (W_n^{k}) ^j=0\]
Proof:\[s= (w_n^k) ^0+ (w_n^k) ^1+...+ (w_n^k) ^{n-1}\]
\[w_n^ks= (w_n^k) ^1+ (w_n^k) ^2+...+ (w_n^k) ^n\]
\[s= (w_n^k) ^n-(w_n^k) ^0/w_n^k-1=0\]

Energetic

DFT definition
For polynomial \ (a\) coefficient expression \ (\vec{a}= (a_0,a_1,..., a_{n?1})
Definition \[y_k=a (w_n^k) =\sum_{i=0}^{n-1}a_i (w_n^k) ^i\]
Then \ (\vec{y}= (y_0,y_1,..., y_{n?1}) \) is the discrete Fourier transform of \ (\vec{a}\)
Remember as \ (\vec{y}=dft_n (\vec{a}) \)
FFT evaluation
\[a (x) =\sum_{i=0}^{n-1}a_ix^i\] original polynomial
\[a^{[0]} (x) =a_0+a_2x+...+a_{n-2}x^{n/2-1}\] \[a^{[1]} (x) =a_1+a_3x+...+a_{n-1}x^{n/2-1}\] call it sub-polynomial

Can see\ (A (x) =a^{[0]} (x^2) +xa^{[1]} (x^2) \), equivalent to the original polynomial is opened according to the subscript parity
At this point the problem becomes a request\ (x= (w_n^k) ^2\)When\ (A^{[0]} (x) \)And\ (A^{[1]} (x) \)The value\ ((k=0,1,..., n-1) \)
According to the binary lemma, there are\[(w_n^1) ^2=w_n^2=w_{n/2}^1\]\[(W_n^{n-1}) ^2=w_n^{-2} (w_n^n) ^2=w_n^{n-2}=w_{n/2}^{n/2-1}\]Actually asked for\ (x= (w_{n/2}^k) (k=0,1,..., n/2-1) \)
Thus we will issue the\ (size\)Halved, will\ (dft_n\)Split into two\ (dft_{n/2}\)To do it, and keep it in hand.
Until\ (size\)When =1,\ (y_0=a_0w_1^0=a_0\), the complexity of the evaluation becomesO (NLOGN)
Specific practice, set\ (\vec{y^{[0]}}=a^{[0]} (x) \)Of\ (dft\)Results\ (\vec{y^{[1]}}=a^{[1]} (x) \)Of\ (dft\)Results
The following are\ (k= (0,1,..., n/2-1) \)
For the first half part\[y^{[0]}_k=a^{[0]} (W^K_{N/2}) =a^{[0]} (w^{2k}_n) \] \[y^{[1]}_k=a^{[1]} (W^K_{N/2}) =a^{[1]} (w^{2k}_n) \] \[y_k=y^{[0]}_k+w_n^ky^{[1]}_k\]
For the second part\[y_{k+n/2}=y^{[0]}_{k+n/2}+w_n^{k+n/2}y^{[1]}_{k+n/2}\]Because\ (w_n^{k+n/2}=w_2^1w_n^k=e^{πi}w_n^k= (cosπ+i\) \ (sinπ) w_n^k=-w_n^k\)
Yes\[y_{k+n/2}=y^{[0]}_{k+n/2}-w_n^ky^{[1]}_{k+n/2}\]
So we can use variables to store\ (w_n^k\), reducing the amount of computation
As for\ (w_n\)Power, we also choose iterative computation to reduce complexity
FFT interpolation
But the complexity of the interpolation is still wrong, and we'll\ (dft\)In the form of a matrix product\ (\vec{y}=v_n\vec{a}\)

For its inverse operation, select\ (\vec{a}=dft^{-1} (\vec{y}) \)
To find out its inverse matrix\ (v_n^{-1}\), the product of a matrix and its inverse matrix is a\ (n*n\)The Unit matrix
Theorem\[[v_n^{-1}]_{ij}=w_n^{-ij}/n\]
Prove
Consider\ (v_n*v_n^{-1}\)In\ ((i,j) \)The element value, according to the summation lemma
\[[v_nv_n^{-1}]_{ij}=\sum_{k=0}^{n-1} (w_n^{-ki}/n) (w_n^{kj}) =\sum_{k=0}^{n-1}w_n^{k (j-i)}/n\]When\ (i=j\)When\ ([v_nv_n^{-1}]_{ij}=\sum_0^{n-1}w_n^0/n=\sum_0^{n-1}e^0/n=1\)
Otherwise\ ([v_nv_n^{-1}]_{ij}=0\), the proposition can be proven
We found it and\ (dft\)The definition is very similar! The only difference lies in the number of previous\ (1/n\)And\ (w_n\)The extra minus sign on the exponent, that is, we can modify it by slightly\ (fft\)It can be used to calculate\ (dft^{?1}\)The!
So far you can already be in the theory\ (O (NLOGN) \)Under the Halo\ (fft\)The

Optimization

Recursive constants seem a little big.
Then push!!!, from one node to the next, merge two each, until the complete polynomial is counted.
But what are the two of them each time they merge?
Look for a pattern.

\ (tm\) turned the binary inverse of the sequence?!!
\ (init\) the number of each subscript in reverse order, merging adjacent two (see Code for details)

There's no?!!!.

#define IL inline#define ri register int#define de double#include<cstdio> #include <iostream> #include < cmath>using namespace Std;const int n=3e6+5;const de Pi=acos ( -1.0), IL int re () {ri x=0,w=1;register char ch=getchar (    ); while (ch< ' 0 ' | |    Ch> ' 9 ') {if (ch== '-') W=-1;ch=getchar ();}    while (ch>= ' 0 ' &&ch<= ' 9 ') x=x*10+ch-' 0 ', Ch=getchar (); return x*w;}    int n,m,l,s=1,r[n];struct complex{de x, y; Complex (de xx=0,de yy=0) {x=xx,y=yy;}} A[n],b[n];complex operator + (complex A,complex b) {return complex (A.X+B.X,A.Y+B.Y);} Complex operator-(complex A,complex B) {return complex (A.X-B.X,A.Y-B.Y);} Complex operator * (complex A,complex b) {return complex (a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} Il void FFT (complex A[],ri op) {for (RI i=0;i<s;i++) if (I<r[i]) swap (a[i],a[r[i]]);//Group by Parity for (RI I=1;I&L T;s;i<<=1) {//The half length of the enumeration interval (i.e. K/2) complex wn (cos (pi/i), Op*sin (pi/i));//primary n units root for (ri len= (i<<1), J=0;J&L T;s;j+=len) {//len is an interval long coMplex W (1,0);                For (ri k=0;k<i;k++,w=w*wn) {//Group processing complex x=a[j+k],y=w*a[j+i+k per cell];                A[j+k]=x+y;            A[j+i+k]=x-y;    }}}}int Main () {n=re (), M=re ();    For (RI i=0;i<=n;i++) a[i].x=re ();    For (RI i=0;i<=m;i++) b[i].x=re ();    while (s<=n+m) s<<=1,l++; For (RI i=0;i<s;i++) r[i]= (r[i>>1]>>1) |    ((i&1) << (L-1));//The number of FFT (a,1) processed in reverse order; FFT (b,1);//evaluation for (RI i=0;i<=s;i++) a[i]=a[i]*b[i];//multiplies FFT (a,-1);//interpolation for (RI i=0;i<=n+m;i++) printf (" %d ", (int) (a[i].x/s+0.5));//+0.5-off accuracy (remember to/s) return 0;}

Fast Fourier transform (FFT)