It's written in 2016,and I just carry it from another blog(zybuluo) of mine.
正文
链接：快速傅里叶变换 易懂版
简单总结一下FFT.
要计算 \(C(x)=A(x)*B(x) \),（A，B，C为带x的多项式），暴力复杂度是o(n^2).
点值表示：设\(y_i=A(x_i)\)，\({ (x_0, y_0), (x_1, y_1), … , (x_{n - 1}, y_{n - 1}) }\)为A的一个点值表示.
假设我们有一种\(O(nlogn)\)算法可以求出A和B的点值表示，那么可以\(o(n)\)求出C的点值表示.
$$\left[
\begin{array} {ccccc}
1 & x_0 & x_0^2 & ... & x_0^{n - 1} \\
1 & x_1 & x_1^2 & ... & x_1^{n - 1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & x_{n-1} & x_{n-1}^2 & ... & x_{n-1}^{n - 1} \\
\end{array}
\right]
\left[
\begin{array} {c}
c_0 \\
c_1 \\
\vdots \\
c_{n-1} \\
\end{array}
\right]
=
\left[
\begin{array} {c}
y_0 \\
y_1 \\
\vdots \\
y_{n-1} \\
\end{array}
\right]$$
如图，已知C的点值表示\({y_i}\),只要将等式右边的列向量乘上左边矩阵的逆就可以得到C的系数.可以证明，这个矩阵的逆存在当且仅当x两两不等.
求出C的系数的过程：系数表示\(\rightarrow\)点值表示\(\rightarrow\)乘法\(\rightarrow\)系数表示
单位复数根：
n次单位复数根是满足\(x^n=1的复数\omega\),记作\(\omega_n\).
几个重要的引理：
\((\omega^0_n, \omega^1_n, … , \omega^{n - 1}_n)\)将复平面平均地分成了n个部分，从\(\omega^0_n=1\)开始，每一个\(\omega\)逆时针旋转\(2\pi/n^\circ\)得到下一个\(\omega\)；\(\omega^p_n=\omega^{p\%n}_n\)； \(\omega^{dp}{dk}=\omega^p_k\).
设
$$A^{[0]}(x) = a_0 + a_2x + a_4x^2 + … + a_{n - 2}x^{n/2 - 1} \\
A^{[1]}(x) = a_1 + a_3x + a_5x^2 + … + a_{n - 1}x^{n/2 - 1}$$
则
$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$$
代入\(\omega^k_n\),得
$$A(\omega^k_n) = A^{[0]}(\omega^{2k}_n) + \omega^k_nA^{[1]}(\omega^{2k}_n)$$
为方便，取n为2的整数次幂，由引理知\(\omega^{(n/2+k)2}_n=\omega^{k2}_n(k>=0)\),故平方之后只有n/2个不同的\(\omega\).
可以得到一个分治算法：从第一层\(\omega^{0..n-1}_n\)开始,每次向下递归询问\(\omega^{0..n/2-1}_n\),边界条件为\(A(\omega^0_n)=a_k\);回溯时\(O(n)\)更新,总复杂度为\(O(nlogn)\).
值得一提的是更新过程（偷懒copy神犇blog）:
$$y_k = A(\omega^k_n) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\
\;\;\; = A^{[0]}(\omega^{2k}_n) + \omega^{k}_n A^{[1]}(\omega^{2k}_n) \\
= y^{[0]}_k + \omega^k_n y^{[1]}_k \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$
$$y_{k + n / 2} = A(\omega^{k + n / 2}_n) \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = A^{[0]}(\omega^{2k + n}_n) + \omega^{k + n / 2}_n A^{[1]}(\omega^{2k + n}_n) \\
\;\;\; = A^{[0]}(\omega^{2k}_n) - \omega^k_n A^{[1]}(\omega^{2k}_n) \\
= y^{[0]}_k - \omega^k_n y^{[1]}_k \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$
上述操作被称作蝴蝶操作.
如果用递归实现，每次要重新计算\({a_i}\)，且递归常数较大,这里有一种快捷生成a队列的雷德算法:
void Rader(Complex *y, int n)
{
for(int i = 1, j = n >> 1, k; i < n - 1; ++ i)
{
if(i < j) swap(y[i], y[j]);
k = n >> 1;
while(j >= k)
{
j -= k;
k >>= 1;
}
if(j < k) j+=k;
}
}
其中y为系数序列（实际上我并没有看懂..).
怎样求
$$\left[
\begin{array} {ccccc}
1 & 1 & 1 & ... & 1 \\
1 & \omega_n & \omega^2_n & ... & \omega^{n - 1}_n \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_n^{n - 1} & \omega_n^{2(n - 1)} & ... & \omega_n^{(n - 1)(n - 1)} \\
\end{array}
\right]$$
的逆矩阵? 有一个公式:
$$a_i = \frac{1}{n} \sum_{k = 0}^{n - 1} y_k \; \omega_n^{-ki}$$
与原FFT模型很相似:
$$y_k = \sum^{n - 1}_{i = 0} a_i \omega^{ki}_n$$
在实际代码中求逆过程直接套用FFT过程.
绕了一大圈，终于总复杂度就变成了伟大的\(O(nlogn)\)!!!（掌声）
还不懂看代码和神犇链接qaq(反正我也是copy的).
const double PI=acos(-1.0);
struct complex{
double r,i;//real && imag
void set(double _,double __){r=_,i=__;}
}a[nn],b[nn];
complex operator +(complex x,complex y) {return (complex){x.r+y.r,x.i+y.i};}
complex operator -(complex x,complex y) {return (complex){x.r-y.r,x.i-y.i};}
complex operator *(complex x,complex y) {return (complex){x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r};}
int len,res[nn];
void rader(complex *y)
{
for(int i=1,j=len>>1,k;i<len-1;++i)
{
if(i<j) swap(y[i],y[j]);
k=len>>1;
while(j>=k) j-=k,k>>=1;
if(j<k) j+=k;
}
}
void FFT(complex *y,int op)
{
rader(y);
complex wn,w,u,t;
for(int h=2;h<=len;h<<=1)
{
wn.set(cos(op*2*PI/h),sin(op*2*PI/h));
for(int i=0;i<len;i+=h)
{
w.set(1,0);
for(int j=i;j<i+h/2;++j)
{
u=y[j];
t=w*y[j+h/2];
y[j]=u+t;
y[j+h/2]=u-t;
w=w*wn;
}
}
}
if(op==-1)
for(int i=0;i<len;++i) y[i].r/=len;
}
void work()
{
FFT(a,1); FFT(b,1);
for(int i=0;i<len;++i) a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<len;++i) res[i]=(int)(a[i].r+0.5);
}
另外补充一下NTT(Number Theory Transform):
当题目给出的模数满足\(2^k|P-1\)，其中\(2^k>n\)时就可以用这个式子：
$$g^k\equiv \omega_{m} \pmod{P}$$
其中g是P的原根,k=\({\frac{P-1}{m}}\).
原根的定义是满足\(\gcd(P,g)=1 且 \; g^{\varphi(p)} \equiv 1 \pmod{P}\)且$$g^{j} \ne 1 \;mod(P), for \; j < \varphi(p)$$
g比较小,直接枚举g即可,一般为3或5.
Comments | NOTHING