NTT&FFT(快速？变换,以及扩展)
FFT&NTT(以及扩展)
预备知识：用于NTT
NTT/FFT其实本质相同，用途是快速求解 多项式乘积
前言
FT: 傅里叶变换:
这是一个工程上的概念，可以简述为：一个周期性的信号波段可以用 若干个正弦曲线 的带权和表示
DFT: 离散傅里叶变换，这是傅里叶变换在离散情况下的变种
FFT: 快速傅里叶变换
NTT: 快速数论变换
谈及核心思想
1.单位根:
构造\(\omega_n\)为\(n\)阶单位根(不知道\(\omega_n\)的值域)，满足性质\(\omega_n^n=\omega_n^0=1\)
对于\(2|n\),\(\omega _n^{\frac{n}{2}}=-1\)
显然\(\omega_n\)满足一个非常简单的性质：折半引理 \(\begin{aligned} \forall 2|i\and 2|n , \omega_n^i=\omega_{\frac{n}{2}}^{\frac{i}{2}}\end{aligned}\)
\(\omega_n\)实际上是一个在幂次上呈现\(n\)元循环的数值
2.多项式与点值式的转化
一个\(n\)阶多项式最普通的表示就是\(F(x)=\sum_{i=0}^{n-1} a_ix^i\)
然而，多项式也可以用\(n\)个互不相关的点表示，即\((x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\)
两者可以互相转化
对于同\(x_i\)的点值，两个多项式卷积时，其\(y_i\)可以直接对应相乘
FFT/NTT的核心过程是
多项式\(\longrightarrow\) 点值式\(\longrightarrow\)点值式对应相乘\(\longrightarrow\)多项式
而用单位根来构造快速的多项式与点值式的转化
3.分治思想
用于降低多项式与点值式转换的复杂度
FFT的单位根
\((x,y)\)指复数\(i=\sqrt{-1},(x,y)=x+yi\)
基本运算\((x,y)+(a,b)=(x+a,y+b),(x,y)\cdot (a,b)=(ax-by,ay+bx)\)
FFT的单位根是:\(\omega_n\)=\((cos(\frac{2\pi}{n}),sin(\frac{2\pi}{n}))\)
而\(\omega _n^i=(cos(\frac{2\pi}{n}\cdot i),sin(\frac{2\pi}{n}\cdot i))\) (展开发现就是三角函数求和公式)
显然满足单位根的性质
(实际上可以发现，这个说是点值其实就是信号序列的三角函数表示)
NTT
相信您已经了解了原根的一些性质，\(\text{NTT}\)的单位根常用原根构造
\(\text{NTT}\)的单位根实际有较大的局限性，对于质数\(P\)只能构造出\(n|P-1,\omega_n=g^{\frac{P-1}{n}}\)
计算在模意义下就能满足单位根的性质
通常我们\(P\)取\(998244353\)，\(2^{23}|(P-1)\)，它的一个原根是3
实际上，为了满足下面分治需要，构造的模数通常满足\(P-1=s\cdot 2^t\)的\(t\)较大，这类模数我们常称作\(\text{NTT}\)模数
以上部分均为基础知识，相对来说应该不会太难，下面是主要难点
多项式转点值式
接下来我们考虑如何将多项式转化为点值式
对于点值式，我们构造的点横坐标为\(x_i=\omega_n^i\)
具体目标是对于函数\(F(x)\)，求出在\(x_0,x_1,\cdots ,x_{n-1}\)上的函数值
即求出\(F(x_i)=a_0\omega_n^0+a_1\omega_n^{i}+a_2\omega_n^{2i}+\cdots\)
接下来就是核心的分治思想，注意，这里的分治是子问题严格等大的
对于当前问题，分成两部分子问题求解(实际是可以分成多部分的，但是这个是特殊情况暂时不予讨论)，即求解
令\(m=\frac{n}{2}\)
\(\begin{aligned} 2|i,G(x_i)=a_0\omega_{m}^0+a_2\omega_{m}^{\frac{i}{2}}+a_4\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\end{aligned}\)
\(\begin{aligned} 2|i,H(x_i)=a_1\omega_{m}^0+a_3\omega_{m}^{\frac{i}{2}}+a_5\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\end{aligned}\)
更简洁的描述为
\(i<m,G(x_i')=a_0\omega_{m}^0+a_2\omega_{m}^{i}+a_4\omega_{m}^{2i}+\cdots\)
\(i<m,H(x_i')=a_1\omega_{m}^0+a_3\omega_{m}^{i}+a_5\omega_{m}^{2i}+\cdots\)
由于\(G(x'_i),H(x'_i)\)计算的是\([0,m-1]\)项，而求\(F(x_i)\)时用到的是\(0,2,4,\cdots\)项，实际需要访问\(G(x^2_i),H(x^2_i)\)
和\(F(x_i)\)的式子比较，我们得到合并的式子为
\(F(x_i)=G(x^2_i)+x_i H(x^2_i)\)
带入折半引理，实际等价于
\(F(x_i)=G(x'_i)+x_i H(x'_i)\)
注意\(x_i=x'_{i\mod m}\)
为了保证复杂度，尽量使得每次分治的子问题都分为两部分，这样的复杂度为\(O(n\log n)\)
附:实际上，分为\(d\)个子问题时，每次合并的复杂度为\(O(n\cdot d)\)，因此复杂度为
保证每次分治为两个严格等大的子问题，可以从一开始就把\(n\)扩充为\(2\)的幂次
int N=1;
while(N<=n+m) N<<=1;
附:\(d\)个子问题时，设子问题答案为\(G_j(x_i)\)，则合并的式子为
\(\begin{aligned} F(x_i)=\sum_{j=0}^{d-1}x_i^jG_j(x_i^d)=\sum_{i=0}^{d-1}x_i^jG_j(x'_{i\mod \frac{n}{d}})\end{aligned}\)
点值式转多项式
一个简单的性质：单位根反演 \(\sum_{j=0}^{n-1}\omega_n^{ij}= \left\{\begin{aligned} \frac{\omega_n^{in}-1}{\omega_n^i-1}=0 && i\ne 0\\ n && i=0\end{aligned} \right.\)
设点值式对应\(y_i\)的序列为\(b_i\)
则\(n\cdot a_i=\sum_{j=0}^{n-1}\omega_n^{-ij} b_j\)
证明
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j=\sum_{j=0}^{n-1} \omega_n^{-ij}(\sum_{k=0}^{n-1}a_k\omega_n^{jk})\end{aligned}\)
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j= \sum_{k=0}^{n-1}a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)} \end{aligned}\)
由上面的式子，发现只有\(k-i=0\)时右边的求和式有值，所以上式成立
因此点值式转多项式直接把系数改为\(\omega_n^{-i}\)即可
Tips:
1.由于单位根的循环特性，溢出会直接溢出到本来的式子里
因此，如果乘法过后的多项式产生了超过\(>n\)的项\(x^i\)，会溢出到\(x^{i\mod n}\)
2.点值式并不是不满足除法，只是除法得到的多项式并不一定是一个\(n\)元以内的多项式，除了恰好整除的情况，得到的通常是一个无穷级数的式子，如\(\begin{aligned} \frac{1}{1-x}=\frac{1-x^{\infty}}{1-x}=\sum_{i=0}^{\infty}x^i\end{aligned}\)
真正要求除法，通常是求前\(n\)项的结果，即需要用到多项式乘法逆
代码实现与优化
然后我们得到一份优美的代码(FFT)
(Complex是C++库自带的复数,M_PI是C++自带\(\pi\)常量)
void FFT(int n,Complex *a,int f) {
if(n==1) return;
Complex tmp[N];
int m=n/2;
rep(i,0,m-1) tmp[i]=a[i<<1],tmp[i+m]=a[i<<1|1]; // 按照奇偶分类
memcpy(a,tmp,sizeof(Complex) * n);
FFT(m,a,f),FFT(m,a+m,f); // 分两半，算g(x),h(x)
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0); // w=x^1,e=x^i
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m]; // f(x_i)=g(x_i)+e*h(x_i)
e=e*w;
}
rep(i,m,n-1) {
tmp[i]=a[i-m]+e*a[i];
e=e*w;
}
memcpy(a,tmp,sizeof(T)*n);
}
由于\((\omega_n)^{\frac{n}{2}}=-1\)，所以还可以简化为
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0);
rep(i,0,m-1) {
tmp[i]=a[i]+e*a[i+m];
tmp[i+m]=a[i]-e*a[i+m];
e=e*w;
}
由于用了double，最后输出要取整
蝴蝶优化
我们加一点优化，取代递归的分治过程
可以看到，分治时我们按照\(i \mod 2\)分成两组，然后继续分
这个过程中，实际上我们就是将\(i\)的二进制位前后翻转
所以我们可以暴力处理出\(i\)分治底层的位置
rep(i,0,n-1) {
int x=i,s=0;
for(int j=1;(j<<c)<=n;++j) {
s=(s<<1)|(x&1);
x>>=1;
} // s就是最终位置
}
当然也是有\(O(n)\)处理方法的
int N=1,c=-1;
while(N<=n+m) N<<=1,c++;
rep(i,1,N-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
(建议自己模拟一下)
有了这个翻转数组，我们可以直接从分治底层开始解决整个问题，每次合并操作完全相同
每次分治问题的大小，依次合并每一个子问题区间即可
为了在一个数组上完成操作，还需要注意合并顺序
代码解释\(i\):分治子问题大小为\(2i\)，\(l\):合并区间的左端点为\(l\)，右端点为\(l+2i\)，\(j\)枚举合并位置
void FFT(int n,Complex *a){
for(int i=1;i<n;i<<=1) {
Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n));
for(int l=0;l<n;l+=i*2) {
Complex e(1,0);
for(int j=l;j<l+i;++j,e=e*w) {
Complex t=a[j+m]*e; // a'[j]=a[j]+e*a[j+m]
// a'[j+i]=a[j]-e*a[j+m]
a[j+m]=a[j]-t;
a[j]=a[j]+t;
}
}
}
}
事实上我们还有更快的写法，就是将\(\omega_n^i\)预处理出来
(注意这个\(\text{FFT}\)的预处理很考验double精度，不能每次都直接累乘上去，隔几个就要重新调用依次三角函数)
当然如果自己写复数会更快
关于点值式转多项式的优化
由于每次求得点值是\(\omega_n^{-i}=\omega_n^{n-i}\)
所以可以直接用 多项式转点值式的函数， 最后把\([1,n-1]\)这一段翻转，每个数除掉\(n\)即可
对于加减运算取模的优化
三目运算
a+=b,a=a>=P?a-P:a;
a-=b,a=a<0?a+P:a;
逻辑运算优化(原理是逻辑预算会在第一个确定表达式值的位置停下)
a+=b,((a>=P)&&(a-=P));
a-=b,((a<0)&&(a+=P));
关于系数预处理优化(以NTT为例)
带入上面已经提到的优化，无预处理系数的\(\text{NTT}\)大概是这样的
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int l=0;l<n;l+=i*2) {
int e=1;
for(int j=l;j<l+i;++j,e=1ll*e*w%P) {
int t=1ll*a[j+i]*e%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
一种简单的预处理是，每次对于每个分治大小，预处理依次系数
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=1;j<i;++j) e[j]=1ll*e[j-1]*w%P;
//for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
//这个版本是沿用上一次预处理的结果，实际(只有)用这种预处理方法可以极大程度上加强FFT的精度
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
另一种是在一开始就把所有的系数用一个数组存下来，具体过程可以描述为
对于每个分治长度\(n\)，我们只需要访问\(\omega_n^{0},\omega_n^{1},\cdots,\omega_n^{\frac{n}{2}-1}\)
那么对于分治长度\(n\)，我们在\(w\)数组的第\(\frac{n}{2}\) ~ \(n-1\)项依次存储这些值
优化:我们只需要对于最大的分治长度处理，剩下的部分发现可以直接用折半引理访问得到
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
const int N=1<<21;
int a[N],w[N];
void Init(){
w[N>>1]=1;
int t=qpow(3,(P-1)/N);
rep(i,(N>>1)+1,N-1) w[i]=1ll*w[i-1]*t%P;
drep(i,(N>>1)-1,1) w[i]=w[i<<1];
}
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<n;i<<=1) {
int *e=w+i;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
三份代码在duck.ac上的评测结果表明，不预处理系数将近慢一倍
单组数据来看，预处理系数会慢一点
多组来看，预处理系数会快
实际差距不大，都可以使用
但是在某些层面来说，下面这份板子才是最好的(适用NTT,FFT且精度较高)，不需要预处理
ll qpow(ll x,ll k=P-2){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int a[N],e[N];
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
e[0]=1;
for(int i=1;i<n;i<<=1) {
int w=qpow(3,(P-1)/i/2);
for(int j=i-2;j>=0;j-=2) e[j+1]=1ll*w*(e[j]=e[j>>1])%P;
for(int l=0;l<n;l+=i*2) {
for(int j=l;j<l+i;++j) {
int t=1ll*a[j+i]*e[j-l]%P;
a[j+i]=a[j]-t,((a[j+i]<0)&&(a[j+i]+=P));
a[j]+=t,((a[j]>=P)&&(a[j]-=P));
}
}
}
if(f==-1) {
reverse(a+1,a+n);
int Inv=qpow(n);
rep(i,0,n-1) a[i]=1ll*a[i]*Inv%P;
}
}
拓展
1.分治+NTT
常用于处理多个计数背包的快速合并 (实际无权值01背包也是可以的)
我们可以用NTT\(n\log n\)合并两个大小为\(n\)的背包
分治时，每次合并两个分治子问题，总共的时间就是\(\sum size\log n\)
每个背包的\(size\)会被计算\(\log n\)次，所以总共复杂度是\(n \log ^2 n\)
2.CDQ+NTT
对于形如\(dp_i=\sum_{j=0}^{i-1}dp_jg_{i-j}\)的\(dp\)转移(就是dp转移与差值有关)
由于求\(dp_i\)时，需要保证\(dp_0,dp_1,\cdots,dp_{i-1}\)才能卷积，这个限制，我们可以用CDQ分治解决
对于当前分治区间\([L,R]\)
依次考虑\([L,mid]\)内部转移，\([L,mid]\)向\([mid+1,R]\)的转移(用FFT/NTT解决)，\([mid+1,R]\)内部转移
算法流程
void Solve(l,r){
if(l==r) return;
mid=(l+r)>>1;
Solve(l,mid);
(l,mid)->(mid+1,r);
Solve(mid+1,r);
}
3.MTT(任意模数NTT)
4.\(n\)元点值式
练习建议：
1.高精度乘法
6.\(CDQ\)分治+\(FFT\)：HDU-5730 题解
7.\(CDQ\)+NTT/降次前缀和优化\(dp\)：HDU-5332 题解
9.图上\(dp\)：
联通图个数：BZOJ-3456 题解
森林数量和带限制森林数量：HDU - 5279 题解
10.点分治+FFT：CodeChef-PRIMEDST 题解
更多应用和优化参见毛啸2016论文
(如：两次FFT做卷积，4次FFT做MTT。。。)
