（笔记）多项式基础 FFT 快速傅里叶变换 NTT 快速数论变换
多项式
对多项式进行乘法，就是对两个多项式进行加法卷积。其中卷积结果 \(C(k)=\sum_{i=0}^kA(i)B(k-i)\)。
分治乘法
将 \(A(x)\) 左右拆半，不足则末尾（最高位）补上 \(0\)，令 \(n=2^k\)。则
同理，\(B(x)\) 左右拆半，则卷积
这样可以转化为 \(3\) 个规模为 \(n/2\) 的子问题，总时间为 \(T(n)=O(n)+3T(n/2),T(n)=O(n^{\log_2 3})=O(n^{1.585})\)。
点值与插值
点值：对于一个 \(x\) 的 \(F(x)\) 的值。
插值：已知 \(F(x)\) 的若干点值，求其系数序列 \(G(x)\)。
根据定义，\(F(x)=\sum_{i=0}^{n-1}x^iG(i)\)。
FFT 中的插值使用单位根反演实现。
复平面单位根
\(z=\cos\theta+i\sin\theta,\omega_n^k=\cos(\frac{2\pi k}{n})+i\sin(\frac{2\pi k}{n})\)，其点集为单位圆均分成 \(n\) 份的点集，且点集包含 \(1\)。
-
性质 1：\(\omega_n^k=\omega_n^{k+n}\)
-
性质 2：\(\omega_{nb}^{kb}=\omega_n^k\)。
-
性质 3：若 \(n\) 为偶数，\(\omega_n^{k+n/2}=-\omega_n^k\)
FFT
DFT
FFT 多项式乘法外层过程如下：
- 选定 \(n(n=2^{k'})\) 个数 \(\omega_n,\omega_n^{2},...\omega_n^{n-1}\)。
- 对于 \(i\in [0,n-1]\) 求出点值 \(A(\omega_n^i),B(\omega_n^i)\)（DFT）
- 根据 \(C(\omega_n^i)=A(\omega_n^i)B(\omega_n^i)\) 求出 \(C\) 的点值。
- 插值求出 \(C(x)\) 的系数。（IDFT，即 DFT 逆过程）
类似分治乘法，我们采用拆半，但是是奇偶拆半。
假设我们已经知道了一堆点值，代入 \(\omega_n^k,k\in[0,\frac{n}{2})\)。
考虑到 \(FL,FR\) 的计算，我们每次实际分治时只需要计算所有 \(k<\frac{n}{2}\) 的部分，\(\frac{n}{2}\le k<n\) 的部分根据性质 1 \(\omega_{n/2}^{k+n/2}=\omega_{n/2}^k\)，用到的 \(FL,FR\) 是和前面一样的，但是根据性质 3，我们在转移式子中用到的 \(\omega_n^{k+n/2}\) 需要变成 \(-\omega_n^k\)，即：
这是 DFT 的部分。
IDFT
根据上面的 DFT：
根据单位根反演：
所以求出点值以后，我们再跑一边 DFT 即可。
这样我们就可以用子问题求出原问题的答案了，时间复杂度 \(T(n)=T(n/2)+O(n),T(n)=O(n\log n)\)。
一些常数优化
迭代版神秘常数问题
- 如果去掉分治过程直接迭代，先枚举起始点再枚举增量，否则会变得异常缓慢。
蝴蝶变换
-
为了避免大规模数组拷贝，我们使用蝴蝶变换。注意到分治的过程中我们先不断地按照奇偶把数组切分然后分到两边，为了避免这个切分我们考虑一开始就把点放到它最下面那层分治所处的位置，然后合并和之前一样左边作为 \(FL\)，右边作为 \(FR\)。
我们发现一个点 \(i\) 的二进制形式决定了其最终所处位置，把 \(i\) 视为 \(n=2^{k'}\) 的 \(k'\) 位二进制数，那么 \(i\in[0,n-1]\)，从上至下每次奇偶分类相当于看二进制最低位，\(0\) 就分到左边，\(1\) 就分到右边，然后右移一位。不妨把这个二进制数
reverse（翻转）一下，那么每次的选择及其权重和就恰好对应了翻转后的二进制数。于是预处理 \(i\) 的reverse\(tr[i]\)，然后若 \(i<tr[i]\) 交换即可得到最下层位置（避免两次交换）。
三次变两次优化
构造出：
即把 \(F(x)\) 和 \(iG(x)\) 系数相加得到 \(P(x)\) 系数后做一次 DFT，得到的点值平方后做 IDFT，最后得到的系数虚部 \(/2\) 就是 \(F(x)G(x)\)。double 运算容易掉精度，承受能力为跨度上线 \(10^{12}\) 左右，如果相乘系数差太大（如一个 \(\in[10^{-6},10^{-5}]\)，一个 \([10^5,10^6]\)）在有平方项会导致 \(10^{24}\) 的精度跨度存在。解决办法就是乘上一个定值，做完 FFT 再除回去即可。
Code
分治版
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double Pi=acos(-1);
const int N=1.35e6+5;
int n,m,tr[N<<1];
struct CP{
CP(double xx=0,double yy=0){x=xx,y=yy;}
double x,y;
CP operator +(const CP &b)const
{return (CP){x+b.x,y+b.y};}
CP operator -(const CP &b)const
{return (CP){x-b.x,y-b.y};}
CP operator *(const CP &b)const
{return (CP){x*b.x-y*b.y,x*b.y+y*b.x};}
}f[N<<1],g[N<<1];
void FFT(CP *f,int len,bool tf){
if(len==1)return ;
FFT(f,len/2,tf);
FFT(f+len/2,len/2,tf);
CP per(cos(2*Pi/len),sin(2*Pi/len)),now(1,0);
if(!tf)per.y=-per.y;
for(int k=0;k<len/2;k++){
CP tt=now*f[k+len/2];
f[k+len/2]=f[k]-tt;
f[k]=f[k]+tt;
now=now*per;
}
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>n>>m;
for(int i=0;i<=n;i++)cin>>f[i].x;
for(int i=0;i<=m;i++)cin>>g[i].x;
m=m+n;for(n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
for(int i=0;i<n;i++)
if(i<tr[i])
swap(f[i],f[tr[i]]),
swap(g[i],g[tr[i]]);
FFT(f,n,1);FFT(g,n,1);
for(int i=0;i<n;i++)f[i]=f[i]*g[i];
for(int i=0;i<n;i++)
if(i<tr[i])
swap(f[i],f[tr[i]]);
FFT(f,n,0);
for(int i=0;i<=m;i++)cout<<(int)(f[i].x/n+0.49)<<' ';
return 0;
}
迭代版
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const double Pi=acos(-1);
const int N=1.35e6+5;
int n,m,tr[N<<1];
struct CP{
CP(double xx=0,double yy=0){x=xx,y=yy;}
double x,y;
CP operator +(const CP &b)const
{return (CP){x+b.x,y+b.y};}
CP operator -(const CP &b)const
{return (CP){x-b.x,y-b.y};}
CP operator *(const CP &b)const
{return (CP){x*b.x-y*b.y,x*b.y+y*b.x};}
}f[N<<1],g[N<<1];
void FFT(CP *f,bool tf){
if(n==1)return ;
for(int i=0;i<n;i++)
if(i<tr[i])swap(f[i],f[tr[i]]);
for(int p=2;p<=n;p<<=1){
int len=p>>1;
CP per(cos(2*Pi/p),sin(2*Pi/p));
if(!tf)per.y=-per.y;
for(int l=0;l<n;l+=p){
CP now(1,0);
for(int k=l;k<l+len;k++){
CP tt=now*f[k+len];
f[k+len]=f[k]-tt;
f[k]=f[k]+tt;
now=now*per;
}
}
}
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>n>>m;
for(int i=0;i<=n;i++)cin>>f[i].x;
for(int i=0;i<=m;i++)cin>>g[i].x;
m=m+n;for(n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
FFT(f,1);FFT(g,1);
for(int i=0;i<n;i++)f[i]=f[i]*g[i];
FFT(f,0);
for(int i=0;i<=m;i++)cout<<(int)(f[i].x/n+0.49)<<' ';
return 0;
}
NTT
找到单位根 \(\omega_n^1\) 的替代品原根。
同余类与剩余系
根据证明，在复数域 \(\mathbb C\) 中不存在与单位根性质相同的其他数了。考虑转移到同余意义下。我们把 \(a\equiv x\pmod p\) 的所有 \(a\in \mathbb Z\) 称为模 \(p\) 意义下 \(x\) 同余类。根据费马小定理，\(p\in\mathbb P\) 时有 \(a^{p-1}\equiv 1\pmod p\)。更进一步地，欧拉定理指出 \(p>0\) 时有 \(a^{\varphi(p)}\equiv 1\pmod p\)。这引出我们对模意义下的同余类在不停地乘一个数（单位根状物）时循环性质的思考。
完全剩余系
对于 \(a_1,a_2,\dots,a_p\)，对于任意 \(z\in\mathbb Z\) 都能找到唯一一个与之在一个同余类的 \(a_u\)，那么这 \(p\) 个数就是模 \(p\) 意义下的一个完全剩余系。需要注意，\(a_1,a_2,\dots,a_p\) 可以不止是 \([0,p-1]\) 之间的数，
既约剩余系
相对地来说，需要保证 \(a_1,a_2,\dots,a_m\) 都有 \((a_i\bmod p,p)=1\)。
我们发现对于 \(p\in\mathbb P\) 完全剩余系和既约剩余系是可以等同的。
阶
最小的满足 \(n=\xi_p(a),a^n\equiv 1\pmod p\) 的正整数 \(n\)。
我们发现 \(n\) 就是一个最小循环节，利用这个性质我们可以做很多事情。
原根
满足 \(\xi_p(a)=\varphi(p)\) 的任意 \(a\)。
在OI-wiki上有关于原根规模（\(O(n^{0.25})\)）和其存在性判断的详细解析（\(a\) 是 \(p\) 的原根当且仅当 \(a^{\varphi(p)}\equiv 1\pmod p\)，且对于所有 \(p\) 的质因数 \(k\)，都有 \(a^{\frac{\varphi(p)}{k}}\not\equiv 1\pmod p\)）。然而 NTT 的应用我们只需要找到一些特殊的模数及其原根即可。
由于 FFT 中 \(n\) 是 \(2^{k'}\) 的形式，我们可以考虑找到一个大质数 \(p=1+a2^b\) 其最小原根为 \(g\)，那么就有 \(g^{p-1}\equiv 1\pmod p\)，且根据原根的性质 \(g^0,g^1,\dots,g^{p-1}\) 均没有任何一对模意义下在一个同余类的数，这当然就可以构成一个 \(p-1\) 次类单位根状物。尝试把它变成 \(n\) 次单位根，把每 \((p-1)/n\) 个数合并即可，得到 \(g^{\frac{p-1}{n}}\) 是一个 \(n\) 次单位根。验证一下关键性质 \((\omega_{2n}^1)^2=\omega_n^1\)，即 \((g^{\frac{p-1}{2n}})^2\equiv g^{\frac{p-1}{n}}\pmod p\)，指数运算是成立的。这里在分治过程中还要保证 \(p-1\) 足够大，覆盖 \(n\) 的所有 \(2\) 的幂次可能。这里我们选用模数 \(998244353=119\times 2^{23}+1\)，原根为 \(3\)，足够使用。
实现
代码中直接替换即可。同等数据规模下，NTT \(1.44s\)，性能远超依赖 double 运算的 FFT（\(2.18s\)）。
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL MOD=998244353;
const double Pi=acos(-1);
const int N=1.35e6+5;
int n,m,tr[N<<1];
LL qkpow(LL x,LL y){
LL res=1;x%=MOD;
while(y){
if(y&1)res=res*x%MOD;
x=x*x%MOD;
y>>=1;
}
return res;
}
LL f[N<<1],g[N<<1];
LL G=3;
void NTT(LL *f,int len,LL per){
if(len==1)return ;
LL nxt=per*per%MOD;
NTT(f,len/2,nxt);
NTT(f+len/2,len/2,nxt);
LL now=1;
for(int k=0;k<len/2;k++){
int tt=now*f[k+len/2]%MOD;
f[k+len/2]=f[k]-tt;
if(f[k+len/2]<0)f[k+len/2]+=MOD;
f[k]=f[k]+tt;
if(f[k]>=MOD)f[k]-=MOD;
now=now*per%MOD;
}
}
int main(){
ios::sync_with_stdio(0);
cin.tie(0);cout.tie(0);
cin>>n>>m;
for(int i=0;i<=n;i++)cin>>f[i];
for(int i=0;i<=m;i++)cin>>g[i];
m=m+n;for(n=1;n<=m;n<<=1);
for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
for(int i=0;i<n;i++)
if(i<tr[i])
swap(f[i],f[tr[i]]),
swap(g[i],g[tr[i]]);
LL w1=qkpow(G,(MOD-1)/n);
NTT(f,n,w1);NTT(g,n,w1);
for(int i=0;i<n;i++)f[i]=f[i]*g[i]%MOD;
for(int i=0;i<n;i++)
if(i<tr[i])
swap(f[i],f[tr[i]]);
NTT(f,n,qkpow(w1,MOD-2));
LL invn=qkpow(n,MOD-2);
for(int i=0;i<=m;i++)cout<<f[i]*invn%MOD<<' ';
return 0;
}
浙公网安备 33010602011771号