FFT & NTT 及其简单优化
FFT
FFT 是一种高效实现 DFT 和 IDFT 的方式，可以在 \(O(n \log n)\) 的时间内求多项式的乘法。
多项式的点值表示
不同于用每项的系数来表示一个多项式，我们知道对于给定的 \(n+1\) 个点值，可以确定唯一的 \(n\) 次多项式。这种用点值表示多项式的方法叫点值表示法。
如果知道 \(F(x)\) 和 \(G(x)\) 的点值表示，求出 \(F(x)\times G(x)\) 的点值表示是 \(O(n+m)\) 的。复杂度瓶颈变成了如何快速转换多项式的两种表示形式。
把系数转成点值的算法叫 离散傅里叶变换(DFT)，相应地，把点值转成系数的算法叫 逆离散傅里叶变换(IDFT)。
DFT
前置知识：单位根
OIWiki-复数-单位根
已经说过，要选取 \(n\) 个点并求出它们对应 \(F(x)\) 的值。先观察 \(F(x)\) 的性质：
把奇偶项分开，得
似乎看起来对每个 \(x\) 求 \(G\) 和 \(H\)，要求的还是 \(n\) 个值。但这时单位根的性质就能用上了：
令 \(x\) 分别为 \(\omega_n^1,\omega_n^2,...,\omega_n^n\)，那么有：
根据 \(\omega_n^{k}=-\omega_n^{k+n/2}\)，且 \(G,H\) 均为偶函数，可以得到：
发现这是好的，因为只用求前一半的 \(G\) 和 \(H\) 的值就可以了，把它们两个分别递归处理即可。
但是 FFT 只能对次数为 \(2^m\) 的多项式这么做，所以如果多项式不是恰好 \(2^m\) 项，需要在后面用一些系数为 \(0\) 的项把它补满。
IDFT
想完成多项式乘法，还要把点值表示重新转换回系数表示，这个过程叫 IDFT。
构造方法实在过于神秘，以至于菜猫猫并不能理解这个思路的由来，所以这里直接给出做法：
设对 \(F(x)\) 进行 DFT 后得到点值表示 \(y_0,y_1,...y_{n-1}\)，构造 \(G(x)=y_0+y_1x+...+y_{n-1}x^{n-1}\)。
将 \(x=\omega_n^0,\omega_n^{-1}...\omega_n^{-(n-1)}\) 代入，得：
那么当且仅当 \(j=k\)，\(G=a_k\times n\)。因此对 \(G(x)\) 代入 \(x=\omega_n^0,\omega_n^{-1}...\omega_n^{-(n-1)}\) 求 DFT，结果除以 \(n\) 即为多项式的各项系数。
实现上可以设定一个参数为 \(1\) 或 \(-1\) 传入 DFT 函数，减小码量。
蝶形运算
一个很诈骗的东西。不要被网上这种图吓到了：
通俗地讲，其实就是，我们通过 \(G(\omega_{n/2}^k)\) 和 \(H(\omega_{n/2}^k)\) 求 \(F(\omega_{n}^k)\)、\(F(\omega_{n}^{k+n/2})\) 时，本来需要额外开一个数组；
但你发现 \(G,H\) 存储的位置和你要把 \(F\) 求完存储的位置是同一个，且求其他位置 \(F\) 的值时也不再用到这两个数据，一边求一边覆盖是不影响正确性的。
讲完了，跟蝶形确实没什么关系对吧。
实现-递归版
众所周知 c++ 自带的 complex 类常数极大，建议手写。
至此，我们可以写出朴素的递归版本 FFT 了：
递归 FFT
void FFT(int limit,Complex a[],int op)
{
if(limit==1) return;
Complex a1[(limit>>1)+5],a2[(limit>>1)+5];
for(int i=0;i<=limit;i+=2) a1[i>>1]=a[i],a2[i>>1]=a[i+1];
FFT(limit>>1,a1,op),FFT(limit>>1,a2,op);
Complex Wn={cos(2.0*pai/limit),sin(2.0*pai/limit)*op},w={1,0};
for(int i=0;i<(limit>>1);i++,w=w*Wn)
{
a[i]=a1[i]+w*a2[i];
a[i+(limit>>1)]=a1[i]-w*a2[i];
}
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++) a[i].x=read();
for(int i=0;i<=m;i++) b[i].x=read();
int w=1;while(w<=m+n) w<<=1;
FFT(w,a,1),FFT(w,b,1);
for(int i=0;i<=w;i++) a[i]=a[i]*b[i];
FFT(w,a,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/w+0.5));
return 0;
}
位逆序置换
递归的常数太大了，接下来讲讲怎么不递归。
考虑每次把系数按奇偶分开的过程，即把下标以二进制的最后一位为关键字排序。也就是说，把二进制倒过来看，如果 \(i\) 按位翻转后为 \(j\)，那么下标为 \(i\) 的数操作到最后就会跑到位置 \(j\) 上。
这启发我们使用迭代实现，预处理出每个系数最后所在的位置，直接从底层倒推即可。
const int N=4e6+5;
const double pi=acos(-1);
int n,m,limit=1;
struct Complex {double x,y;}a[N],b[N];
il Complex operator +(Complex a,Complex b) {return {a.x+b.x,a.y+b.y};}
il Complex operator -(Complex a,Complex b) {return {a.x-b.x,a.y-b.y};}
il Complex operator *(Complex a,Complex b) {return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
int to[N];
void FFT(Complex *a,int tp)
{
for(int i=0;i<limit;i++) if(i<to[i]) swap(a[i],a[to[i]]);
for(int len=1;len<limit;len<<=1)
{
Complex Wn={cos(pi/len),sin(pi/len)*tp};
for(int i=0;i<limit;i+=(len<<1))
{
Complex w={1,0};
for(int j=0;j<len;j++,w=w*Wn)
{
Complex x=a[i+j],y=w*a[i+len+j];
a[i+j]=x+y,a[i+len+j]=x-y;
}
}
}
}
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++) a[i].x=read();
for(int i=0;i<=m;i++) b[i].x=read();
int k=0; while(limit<=n+m) limit<<=1,k++;
for(int i=0;i<limit;i++) to[i]=(to[i>>1]>>1)|((i&1)<<(k-1));
FFT(a,1),FFT(b,1);
for(int i=0;i<limit;i++) a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/limit+0.5));
return 0;
}
三次变两次优化
上文中，进行两个多项式相乘时，我们需要调用三次 FFT 函数：对 \(F\) 和 \(G\) 分别 DFT，对 \(F\times G\) IDFT。这么做是不够优的。
这里就可以用到三次变两次优化的技巧，使常数减小三分之一。
具体地，我们在计算前直接把 \(G\) 放在 \(F\) 的虚部上，乘法操作变为计算 \(F(x)\) 平方。
答案就是 \(F^2(x)\) 的虚部除以 \(2\)。
正确性证明是容易的：\((a+b\text{i})^2=(a^2-b^2)+2ab\text{i}\)。
在主函数上稍作改动即可：
int main()
{
n=read(),m=read();
for(int i=0;i<=n;i++) a[i].x=read();
for(int i=0;i<=m;i++) a[i].y=read();//here
int k=0; while(limit<=n+m) limit<<=1,k++;
for(int i=0;i<limit;i++) to[i]=(to[i>>1]>>1)|((i&1)<<(k-1));
FFT(a,1);
for(int i=0;i<limit;i++) a[i]=a[i]*a[i];//and here
FFT(a,-1);
for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].y/2/limit+0.5));
return 0;
}
NTT
前置知识-原根
阶：使得 \(a^n\equiv1\pmod{m}\) 的最小正整数 \(n\) 称为 \(a\) 模 \(m\) 的阶。
原根：若 \(\gcd(g,m)=1\)，且 \(\delta_m(g)=\varphi(m)\)，则 \(g\) 是 \(m\) 的原根。
存在定理 一个数 \(m\) 存在原根当且仅当 \(m=2,4,p^{\alpha},2p^{\alpha}\)，其中 \(p\) 为奇素数，\(\alpha \in \mathbb{N}^{*}\)。
判定定理 设 \(m\ge 3,\gcd(a,m)=1\)，则 \(a\) 是模 \(m\) 的原根的充要条件是，对于 \(\varphi (m)\) 的每个质因数 \(p\)，都有 \(a^{\frac{\varphi(m)}{p}}\not\equiv 1\pmod m\)。
原根个数 若一个数有原根，则它原根的个数为 \(\varphi(\varphi(m))\)。
以上结论的证明：link
NTT 用来解决模意义下的多项式乘法问题，其运算均为整数，在常数和精度方面均优于 FFT。（划掉，实测还是 FFT 跑得快。
原根的性质与单位根类似，因此我们把 FFT 中所有 \(\omega_n\) 替换成 \(g^{\frac{p-1}{n}}\)，乘 \(-1\) 换成乘逆元即可。
NTT 的常见模数是 \(998244353\)，它的原根是 \(3\)。
为什么我的 NTT 跑得比 FFT 慢啊 /dk
#define int long long
const int N=4e6+5,mod=998244353;
il int qpow(int n,int k=mod-2)
{
int res=1;
for(;k;n=n*n%mod,k>>=1) if(k&1) res=res*n%mod;
return res;
}
int n,m,a[N],b[N],inv=qpow(3),limit=1,to[N];
il void NTT(int *a,int tp)
{
for(int i=0;i<limit;i++) if(i<to[i]) swap(a[i],a[to[i]]);
for(int len=1;len<limit;len<<=1)
{
int Wn=qpow(tp>0?3:inv,(mod-1)/(len<<1));
for(int i=0;i<limit;i+=(len<<1))
for(int j=0,w=1;j<len;j++,w=w*Wn%mod)
{
int x=a[i+j],y=w*a[i+len+j]%mod;
a[i+j]=(x+y)%mod,a[i+len+j]=(x-y+mod)%mod;
}
}
}
signed main()
{
n=read(),m=read();
for(int i=0;i<=n;i++) a[i]=read();
for(int i=0;i<=m;i++) b[i]=read();
int k=0; while(limit<=m+n) limit<<=1,k++;
for(int i=0;i<limit;i++) to[i]=(to[i>>1]>>1)|((i&1)<<(k-1));
NTT(a,1),NTT(b,1);
for(int i=0;i<limit;i++) a[i]=a[i]*b[i]%mod;
NTT(a,-1);
for(int i=0;i<=n+m;i++) printf("%lld ",a[i]*qpow(limit)%mod);
return 0;
}
本文来自博客园，作者：樱雪喵，转载请注明原文链接：https://www.cnblogs.com/ying-xue/p/17676005.html
