FFT 相关
FFT 相关
FFT 数组记得开两倍！
前言
update 2025.09.28 学习了 FFT 自己卷自己。
参考博客
FFT
简介
用于求卷积（\(a,b\) 已知）：
或者多项式乘法（\(A,B\) 均是多项式且已知）：
\(A=\sum_{i=0}^{n} a_i x^i\\ B=\sum_{i=0}^{m} b_i x^i\)
可见 \(C\) 是 \(n+m\) 次多项式。
如果我们把卷积的 \(a_i,b_i\) 看成多项式的系数，卷积就变成求：
求卷积或者多项式乘法的时间复杂度是 \(O(n^2)\) 的。使用 FFT 可以做到 \(O(n\log n)\)。
大体思路
设 \(C\) 的项数为 \(n\)（不是次数），若 \(A,B\) 不足 \(n\) 项就补系数 \(0\)。
显然这个过程很难优化，我们从另一个角度去想。
对于一个多项式，求其在 \(x\) 处的值的时间复杂度是 \(O(n)\) 的，我们把这个操作叫做点值（DFT）。
\(n\) 个点可以唯一确定一个 \(n\) 项多项式（即 \(n-1\) 次多项式），证明不显然但是略过我不会，大概感受一下吧。拉格朗日插值法给出了一种 \(O(n^2)\) 的求出这个多项式（即求出它的系数）的方法原理我不会，这个由点值求系数的过程叫做插值（IDFT）。
因此求多项式乘法，可以变成求任意 \(n\) 个点在两个多项式 \(A,B\) 的点值，然后由 \(C(x)=A(x)B(x)\),\(O(n)\) 求出多项式 \(C\) 在这 \(n\) 个点的点值，然后做一次插值求出 \(C\) 的系数。当然这个也是 \(O(n^2)\) 的，但是聪明的傅里叶给出了一种基于单位根的特殊性质的分治方法求 DFT 和 IDFT，成为快速傅里叶变换（FFT）。
就是先求出 \(A,B\) 的点值（DFT），进行点值相乘得到 \(C\) 的点值，然后插值（IDFT）得到 \(C\) 的系数。
单位根
写作 \(\omega_n^k\)，读作 \(n\) 次单位根的 \(k\) 次方。
\(\omega_n^k\) 是在复数域上的向量，形如在复平面上分成 \(n\) 等分，其中 \(\omega_n^0=\omega_n^n=1\)。按逆时针分别为 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\)。
后面为了方便，我们会假设 \(n=2^k\)。
有几个重要的性质。
- \(\omega_n^n=1\) 正确性显然
- \(\omega_{an}^{ak}=\omega_n^k\) 在平面上想象一下，显然正确
- \(\omega_n^{k-\frac{n}{2}}=-\omega_n^k\) 相当于把向量转 \(180^。\)。
点值
基于这些性质，我们求 \(A,B\) 在 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\) 处的点值。
以求 \(A(x)\) 为例，我们要求所有 \(A(\omega_n^k),0\le k<n\)。
我们把 \(A\) 按奇偶分为两部分：
设
因此有：
代入 \(x=\omega_n^k\)，有：
我们先求出 \(k< \frac{n}{2}\) 的 \(A_0,A_1\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\) 在 \(k< \frac{n}{2}\) 的点值。
然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{n}{2}\) 发现：
因此你发现，这俩十分地相似，因此你求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值之后，可以直接 \(O(n)\) 求出 \(A\) 的 \(n\) 个点值了。然后就这样分治下去，点值时间复杂度为 \(O(n\log n)\)。
因为要一直对 \(n\) 除以 \(2\) 所以令 \(n=2^k\) 意义就在此。（位数不足高位补系数 \(0\)）
插值
求插值的过程是类似的，改几个参数就行了。
把求点值的过程写成矩阵。
由 \(A\) 得到 \(dft(A)_i = A(\omega_{n}^i)\)。有
把 \(p_{i,j}\) 写成矩阵就是：
只要把参数改成逆矩阵，其他不变，就可以由 \(dft(A)\) 求出 \(A\) 了。
根据单位根的性质，矩阵的逆恰好是指数取反：
改进
然后你会发现过不了板子……
因为这个递归的过程常数很大，假设 FFT 的常数本身就大，然后就超时了。
考虑从下往上递推分治的过程。
然后你惊奇地发现最下层的顺序是 \(000,100,010,110,001,101,011,111\)，刚好是 \(0\sim 7\) 的二进制的 reverse。
然后就有了如下板子：
Code
struct fushu {
double x,y;
fushu (double _x=0,double _y=0):x(_x),y(_y){}
}a[N],b[N];
fushu operator + (fushu a,fushu b) { return {a.x+b.x,a.y+b.y}; }
fushu operator - (fushu a,fushu b) { return {a.x-b.x,a.y-b.y}; }
fushu operator * (fushu a,fushu b) { return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }
int len;
const double pi=acos(-1.0);
int re[N];
void FFT(fushu *c,int type) {
rep(i,0,(1<<len)-1) {
if(i<re[i]) swap(c[i],c[re[i]]);
}
for(int k=1;k<(1<<len);k<<=1) {
fushu wn(cos(pi/k),type*sin(pi/k));
for(int r=k<<1,j=0;j<(1<<len);j+=r) {
fushu w(1,0);
for(int i=0;i<k;i++,w=w*wn) {
fushu x=c[j+i],y=w*c[j+k+i];
c[j+i]=x+y;
c[j+k+i]=x-y;
}
}
}
}
int n,m;
int main(){
sf("%d%d",&n,&m);
rep(i,0,n) sf("%lf",&a[i].x);
rep(i,0,m) sf("%lf",&b[i].x);
while((1<<len)<=n+m) len++;
rep(i,0,(1<<len)-1) {
re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
}
FFT(a,1),FFT(b,1);
rep(i,0,(1<<len)-1) a[i]=a[i]*b[i];
FFT(a,-1);
rep(i,0,n+m) pf("%d ",(int)(a[i].x/(1<<len)+0.5));
}
算法缺陷
由于复数域是用浮点数计算的，所以会存在掉精度问题。如果答案是对一个特别的指数取模，如著名的 \(998244353\)，可以使用原根代替单位根计算，在剩余系里计算而不是在复数域计算。详见 NTT。
NTT
阶和原根
欧拉定理：若 \(\gcd(a,n)=1\)，则 \(a^{\varphi(n)} \equiv 1(\bmod n)\)。
阶：设 \(m>1\)，且 \(\gcd(a,m)=1\)，根据欧拉定理一定存在正整数 \(d<m\) 使得 \(a^d \equiv 1(\bmod m)\)。把满足该式子的最小的正整数定义为 \(a\) 对模 \(m\) 的阶（指数），写作 \(\text{ord}_m(a)\)。
阶的一个结论：对于正整数 \(d\)，\(a^d \equiv 1 (\bmod m)\) 的一个充要条件是 \(d \mid \text{ord}_m(a)\)。
原根：设 \(m>1\)，\(\gcd(a,m)=1\)，若 \(\text{ord}_m(a)=\varphi(m)\)，则 \(a\) 为 \(m\) 的原根，记做 \(g_m\)。
常见质数 \(998244353\) 的原根是 \(3\)。
定理 1：\(g_m,g_m^2,\dots g_m^{\varphi(m)-1}\) 两两模 \(m\) 不同余且均与 \(m\) 互质。（其实是两条定理）
定理 2：若 \(m\) 是质数，\(g_m,g_m^2,\dots,g_{m}^{m-2}\) 模 \(m\) 的余数恰好形成 \(1 \sim m-1\) 的排列，而 \(g_m^{m-1}\) 模 \(m\) 的余数则是 \(0\)。
由此定理可以发现原根有很多类似单位根的性质。
定理 3：\(m\) 的原根有 \(\varphi(\varphi(m))\) 个。
简介
NTT 的原理是用原根代替单位根。质数 \(p\) 的原根在模 \(p\) 剩余系意义下具有我们利用的单位根的性质的相同性质。因此如果要求的多项式的系数只需要模 \(p\) 意义下的，可以用原根代替单位根。因为原根是整型，所以可以避免精度问题。
Code
const int N=4e6+7,mod=998244353,G=3,invG=332748118;
ll a[N],b[N];
int len;
int re[N];
ll ksm(ll a,ll b=mod-2) {
ll s=1;
while(b) {
if(b&1) s=s*a%mod;
a=a*a%mod;
b>>=1;
}
return s;
}
void NTT(ll *c,int type) {
rep(i,0,(1<<len)-1) {
if(i<re[i]) swap(c[i],c[re[i]]);
}
for(int k=1;k<(1<<len);k<<=1) {
ll wn=ksm(type==1?G:invG,(mod-1)/(k<<1));
for(int r=k<<1,j=0;j<(1<<len);j+=r) {
ll w=1;
for(int i=0;i<k;i++,w=w*wn%mod) {
ll x=c[j+i],y=w*c[j+k+i]%mod;
c[j+i]=(x+y)%mod;
c[j+k+i]=(x-y+mod)%mod;
}
}
}
}
int n,m;
int main(){
#ifdef LOCAL
freopen("in.txt","r",stdin);
freopen("my.out","w",stdout);
#endif
sf("%d%d",&n,&m);
rep(i,0,n) sf("%lld",&a[i]),a[i]%=mod;
rep(i,0,m) sf("%lld",&b[i]),b[i]%=mod;
while((1<<len)<=n+m) len++;
ll inv=ksm(1<<len);
rep(i,0,(1<<len)-1) {
re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
}
NTT(a,1),NTT(b,1);
rep(i,0,(1<<len)-1) a[i]=a[i]*b[i]%mod;
NTT(a,-1);
rep(i,0,n+m) pf("%lld ",a[i]*inv%mod);
}
分治 FFT
形如 \(f_i=\sum_{j=1}^i f_{i-j}g_j\) 的卷积（\(f_0\) 需要初值）。
一般的卷积我们是知道卷起来的两个多项式的系数的，但是这个卷积显然我们只知道 \(g\) 的系数，却不知道 \(f\) 的系数，那么如何卷呢？可以分治 FFT 解决。
它的思想类于 CDQ 分治。
要求 \(f_{1 \sim n}\)，先求 \(f_{1 \sim mid}\)，然后计算左边对右边的贡献，然后再递归计算右边。
左边对右边的贡献，即算 \(\forall i \in [1,n],\sum_{j=1}^{mid} f_{j}g_{i-j}\)，只需要对 \(f,g\) 做一次 \(O(n\log n)\) 的 FTT。
因为这个卷积不是标准卷积（指 \(j\) 不是枚举到 \(i\)），因此我们可以给 \(j>mid\) 的 \(f_j\) 当做 \(0\) 来做卷积。变成求 \(\sum_{j=1}^i f_j g_{i-j}\)。
然后我们再求 \(f_{mid+1\sim r}\)，因为左边已经对右边贡献过了，所以我们递归右边的时候，计算贡献不需要再带上左边。
整理一下，现在我们要计算 \(f_{[l,r]}\)：
- 先递归算好 \(f_{[l,mid]}\)。
- 算 \(f_{[l,mid]}\) 对 \(f_{[mid+1,r]}\) 的贡献。\(f_i \gets \sum_{j=l}^{mid} f_j g_{i-j}\)。
- 递归计算 \(f_{[mid+1,r]}\)。
\(f_i \gets \sum_{j=l}^{mid} f_j g_{i-j}\)，这是一个标准卷积，FFT 时间复杂度是 \(O(len \log len)\)。
像这样分治下去，直到长度为 \(1\) 时，不用卷积了，直接返回值。分治一共 \(\log\) 层，每一层每个区间都要做 \(O(len\log len)\) 的 FTT，总时间复杂度为 \(O(n\log^2 n)\)。
FFT 自己卷自己
前言
感觉它的思想有点像这题 https://vjudge.net/contest/744791#problem/E。我就说吧，不会的东西先放着，说不定之后那一次恰好就理解了它的思想呢！所以我要学多少次 SAM 才能学会？
介绍
形如
参考分治 FFT 的过程分治。
现在我们要计算 \(f_{[l,r]}\)：
- 先递归算好 \(f_{[l,mid]}\)。
- 算 \(f_{[l,mid]}\) 对 \(f_{[mid+1,r]}\) 的贡献。
- 若 \(l \neq 1\)，可以和分治 FFT 一样做。即 \(f_i \gets \sum_{j=l}^{mid} f_j f_{i-j}\)，\(f_{i-j}\) 一定已知。
- 若 \(l = 1\)，\(f_{[1,mid]} \times f_{[1,mid]}\) 的贡献是可以算的，但是 \(f_{[1,mid]} \times f_{[mid+1,r]}\) 的贡献目前没法算。这部分贡献可以在递归 \([mid+1,r]\) 的时候再算。
- 递归计算 \(f_{[mid+1,r]}\)。
总结整个算法流程，目前分治区间为 \([l,r]\)：
- 若 \(l=1\)：
solve(l,mid)。- 将 \(f_{[1,mid]} \times f_{[1,mid]}\) 贡献至 \(f_{[mid+1,r]}\)。
solve(mid+1,r)。
- 若 \(l \neq 1\)：
solve(l,mid)。- 将 \(f_{[l,mid]} \times f_{[1,r-l]}\) 贡献至 \(f_{[mid+1,r]}\)。
- 将 \(f_{[1,r-l]} \times f_{[l,mid]}\) 贡献至 \(f_{[mid+1,r]}\)。
solve(mid+1,r)。
时间复杂度 \(O(n \log^2 n)\)。
本文来自博客园，作者：wing_heart，转载请注明原文链接：https://www.cnblogs.com/wingheart/p/18412826
浙公网安备 33010602011771号