FFT NTT小结
一个问题
给定\(a_i,b_i\)，定义\(c_n=\sum_{i=0}^n a_ib_{n-i}\),求\(c\)序列。
想法
暴力乘需要\(O(nm)\),有没有办法优化呢？
能不能把多项式转化一下，然后快速的计算再转回来呢？
答案是能的。
多项式有两种表示法：点值表示法和系数表示法。
系数表示法是用n+1个数\(a_i\)来表示，也就是\(F(x)=\sum_{i=0}^na_ix^i\)
点值表示法是用n+1个点\((x_i,y_i)\)来表示，也就是\(x=x_i\)时\(F(x)=y_i\)
表示出来以后,可以直接\(O(n)\)计算：如果x相同的两个点分别是\((x,f(x)),(x,g(x))\)，那么乘积显然是
\((x,f(x)g(x))\)
当然由于乘积的次数是\(n+m\)次的，至少要用\(n+m+1\)个点。
问题就在于如何在系数表示法和点值表示法之间转换。
暴力用\(x=0..n+m\)带入计算然后之后插值计算原函数吗？先不说插值的复杂度（我也不知道没写过），点值计算的复杂度就已经有\(O((n+m)^2)\) ，好垃圾啊
单位复数根
设\(\omega^n=1\)，则\(\omega\)有n个解，在复平面单位圆上，排成了一个正n边形，其中一个顶点是1。
定义\(\omega_n\)为最小的那个解，也就是从1逆时针方向第一个解，比如\(\omega_4=i\)
那么\(\omega_n^i\)其中\(i\)从0到n-1刚好构成了全部的解。(\(\omega_n^0=1\))
这个东西满足一些性质：
- \(\omega_n^i=\omega_n^{i+n}\)
- 如果\(n\)为偶数,（n=2k)，则\(\omega_n^i+\omega_n^{k+i}=0\)
- 如果\(n\)为偶数,（n=2k)，则\(\omega_n^{2i}=\omega_k^i\)
- \(\frac{1}{n}\sum_{i=0}^{n-1}\omega_n^{vi}=[n|v]\)
这些性质很重要
则对于任意\(a,b,c\)，显然
\[ \begin{align} [n|a+b-c]&=\sum_{i=0}^{n-1}\omega_n^{(a+b-c)i}\\ &=\frac{1}{n}\sum_{i=0}^{n-1}\omega_n^{ai}\omega^{bi}\omega^{-ci} \end{align} \]
成立
DFT
令\(n\)是大于c长度的最小的形如\(2^k\)的数。
则 \[ \begin{align} c_m&=\sum_{i}\sum_{j}[i+j=m]a_ib_j\\ &=\sum_{i}\sum_{j}[n|i+j-m]a_ib_j\\ &=\frac{1}{n}\sum_{i}\sum_{j}\sum_{v=0}^{n-1}a_i\omega_n^{iv}b_j\omega_n^{jv}\omega^{-mv}\\ &=\frac{1}{n}\sum_{v=0}^{n-1}\omega^{-mv}\sum_{i}a_i\omega_n^{iv}\sum_{j}b_j\omega_n^{jv}\\ \end{align} \] 令 \[ A(x)=\sum_{i=0}^{n-1}a_ix^i\\ B(x)=\sum_{i=0}^{n-1}b_ix^i\\ c_m=\frac{1}{n}\sum_{i=0}^{n-1}C_i\omega_n^{-im}\\ C(x)=\sum_{i=0}^{n-1}C_ix^i \]
则 \[ C_v=A(\omega_n^v)B(\omega_n^v)\\ c_m=\frac{1}{n}\sum_{v=0}^{n-1}\omega^{-mv}A(\omega_{n}^v)B(\omega_{n}^v)\\ c_m=\frac{1}{n}C(\omega^{-m}) \] 所以我们就只需要从\(a_{0}\dots a_{n-1}\)计算\(A(\omega_n^{0})\dots A(\omega_n^{n-1})\)，以及从\(C_i\)计算\(c_{i}\)，前者是从系数表示法转换为点值表示(DFT)，后者点值表示转换为系数表示(IDFT)。
FFT
一个引理：设\(F(x)=a_0+a_1x+a_2x^2+...\)，\(F_0(x)=a_0+a_2x+a_4x^2+a_6x^3+...\),\(F_1(x)=a_1+a_3x+a_5x^2+...\)，那么\(F(x)=F_0(x^2)+xF_1(x^2)\)
所以
\[ F(\omega_n^m)=F_0(\omega_{\frac{n}{2}}^m)+\omega_n^mF_1(\omega_{\frac{n}{2}}^m)\\ F(\omega_n^{m+\frac{n}{2}})=F_0(\omega_{\frac{n}{2}}^m)+\omega_n^{m+\frac{n}{2}}F_1(\omega_{\frac{n}{2}}^m)=F_0(\omega_{\frac{n}{2}}^m)-\omega_n^{m}F_1(\omega_{\frac{n}{2}}^m) \]
所以我们只需要先给\(F_0\)和\(F_1\)DFT，然后\(O(n)\)合并。
显然\(F_0\)和\(F_1\)都是\(\frac{n}{2}\)次的，那么复杂度为\(T(n)=2T(\frac{n}{2})+O(n)\)根据主定理复杂度为\(T(n\log n)\)
那IDFT怎么计算呢？
根据前面的式子，IDFT只需要把\(\omega_n^m\)改成\(\omega_n^{-m}\)再DFT一遍最后每一项\(/n\)就好了。
或者还有一种写法：直接再DFT一遍，然后把\(c_1\dots c_{n-1}\)反转，然后把每一项\(/n\)
至于单位复数根怎么计算呢？
根据欧拉定理，\(\omega_n^m=\cos(\frac{2\pi m}{n})+i\sin(\frac{2\pi m}{n})\)
当然也可以得到了\(\omega_n^1\)后\(\omega_n^m=\omega_n^{m-1}\omega_n^1\)递推。
但是直接递归有点慢....因为要反复的把元素放进两个数组再合并回来，有没有办法优化呢?
发现递归的过程形如这样（系数表示法）设n=8： \[ (a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7)\\ (a_0,a_2,a_4,a_6)(a_1,a_3,a_5,a_7)\\ (a_0,a_4)(a_2,a_6)(a_1,a_5)(a_3,a_7)\\ (a_0)(a_4)(a_2)(a_6)(a_1)(a_5)(a_3)(a_7) \] 如果我们先按最下面的顺序排好，然后向上一层一层的合并，是不是就节省了很多时间呢？
发现(0,4,2,6,1,5,3,7)和(0,1,2,3,4,5,6,7)二进制恰好相反，这是不是巧合？
仔细思考：
把这n个元素分成两组的时候，是按照最低位为0还是1分的，为0则会分到前面一半（最高位为0），为1则会分到后面一半。
也就是最低位会对应到换之后的最高位，然后次低位也是一样会对应到次高位....
所以就是每个元素编号的二进制位翻转。
这个称作位逆序置换。
如何计算位逆序置换呢？
设\(rev(x)\)表示x的位逆序置换，那么\(rev(x)=rev(\lfloor\frac{x}{2}\rfloor)+(x\bmod 2)\times\frac{n}{2}\)
这个不难证明 。
代码：
typedef complex<double> fs;
void FFT(fs *a,int flag){
for(int i=0;i<s;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int i=1;i<s;i<<=1){
fs wa(cos(pi/i),flag*sin(pi/i));
for(int j=0;j<s;j+=i+i){
fs w(1,0),u,v;
for(int k=0;k<i;++k){
u=a[j+k],v=w*a[j+k+i];
a[j+k]=u+v;
a[j+k+i]=u-v;
w=w*wa;
}
}
}
if(flag<0){
for(int i=0;i<s;++i)a[i]/=s;
}
}FFT(a,1)代表对a进行DFT，FFT(a,-1)代表IDFT
当然还有一种常数优化：手写复数类型而不是用complex。
NTT
还有一些题目，要求膜意义下卷积，也就是卷出来的系数要膜dkw一个数。
一般膜的数可以表示为\(2^ab+1\)的形式，其中\(2^a\)要不小于n。（否则就需要一些奇怪的操作）
思考\(\omega\)用在FFT中的性质 这些在膜意义下有什么东西替代吗？
答案是原根。
由于费马小定理\(g^{p-1}=1\)，我们可以用\(g^\frac{p-1}{n}\)代替FFT的\(\omega_n^1\)，发现也同样满足这些性质。
然后对应的IDFT 用的就是\(g^{-\frac{p-1}{n}}\)
代码：
const ll p=998244353;
const ll g=3;
ll fpm(ll a,ll b){
ll c=1;
while(b){
if(b&1)c=c*a%p;
a=a*a%p;
b>>=1;
}
return c;
}
ll invs;
void NTT(ll *a,ll flag){
for(ll i=0;i<s;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(ll i=1;i<s;i<<=1){
ll w=fpm(g,((flag+(p-1))*(p-1)/i/2));
for(ll j=0;j<s;j+=i+i){
ll aw=1,u,v;
for(ll k=0;k<i;++k){
u=a[j+k],v=a[j+k+i]*aw%p;
a[j+k]=(u+v)%p;
a[j+k+i]=(u-v+p)%p;
aw=aw*w%p;
}
}
}
if(flag<0)for(ll i=0;i<s;++i)a[i]=a[i]*invs%p;
}
例题
UOJ34 多项式乘法
一个模板题。
FFT和NTT都可以
代码：
ZJOI 2014 力
\[ E_i=\sum_{j<i}q_j\frac{1}{(i-j)^2}-\sum_{i<j}q_j\frac{1}{(j-i)^2} \]
设 \[ F_i=\sum_{j<i}q_j\frac{1}{(i-j)^2}\\ G_i=\sum_{i<j}q_j\frac{1}{(j-i)^2}\\ E_i=F_i-G_i \]
\(F_i\)一看就能直接卷积，后面一半感觉也像卷积形式，可是\(j>i\)怎么解决呢？
把\(q_j,G_i\)都翻转就能卷积了。
fs q[1<<19],qr[1<<19],k[1<<19],f[1<<19],F[1<<19];
int rev[1<<19];
void FFT(fs *a,int base,int flag){
for(int i=0;i<base;++i)if(i<rev[i]){
swap(a[i],a[rev[i]]);
// printf("Swap%d,%d\n",i,rev[i]);
}
fs ww,w,u,v;
for(int S=2,s=1;s<base;S<<=1,s<<=1){
ww=fs(cos(pi/s),flag*sin(pi/s));
for(int l=0;l<base;l+=S){
w=fs(1,0);
for(int j=0;j<s;++j){
u=a[l+j];
v=w*a[l+j+s];
a[l+j]=u+v;
a[l+j+s]=u-v;
w=w*ww;
}
}
}
if(flag<0)for(int i=0;i<base;++i)a[i].x/=base,a[i].y/=base;
}
int main(){
#ifdef cnyali_lk
freopen("3527.in","r",stdin);
freopen("3527.out","w",stdout);
#endif
int n;
scanf("%d",&n);
for(int i=1;i<=n;++i){
scanf("%lf",&q[i].x);
qr[n+1-i]=q[i];//翻转
k[i].x=1./i/i;
}
int m=n<<1,base=1,l=-1;
while(base<=m){base<<=1;++l;}
for(int i=0;i<base;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
FFT(q,base,1);
FFT(qr,base,1);
FFT(k,base,1);
for(int i=0;i<base;++i)f[i]=q[i]*k[i],F[i]=qr[i]*k[i];
FFT(f,base,-1);
FFT(F,base,-1);
for(int i=1;i<=n;++i)
printf("%lf\n",f[i].x-F[n+1-i].x);//翻转
return 0;
}