浅谈多项式乘法，FFT与NTT
闲话
老师，数论会和 dp 有关吗？
“会的，而且如果有关的话也会比较简单”
开篇
FFT 和 NTT 都是用来求多项式的，他们可以求所有的 加法卷积，也就是形如 \(h_x=\sum\limits_{i+j=x}a_ib_j\)。那么显然我们的多项式乘法、以及一些 dp 式都具有卷积的特性。
初中老师都教过我们，如果是多个多项式相乘，我们可以拆括号，那么如果两个序列长度为 \(n\)，\(m\)，那么我们的计算量为 \(\mathcal O(nm)\)。
那么到 OI 里，有没有我们更快的算法呢？前面我们都讲过，OI 最大的乐趣，莫过于 \(\mathcal O(n^2)\rightarrow \mathcal O(n\log n)\)。
那么我们如何优化这种算法呢？
点值表示
我们假设我们有一个 \(n-1\) 次多项式，那么我们假设 \(y=\sum\limits_{i=0}^{n-1} a_ix^i\)，那么我们显然有 \(n\) 个解 \(x_0,x_1,\cdots,x_{n-1}\)，并且有对应的 \(y_0,y_1,\cdots,y_n\)，那么我们将 \(<x_0,y_0>,<x_1,y_1>,\cdots,<x_{n-1},y_{n-1}>\) 记作点值。
那么我么假设我们随了几组特别强的 \(x_i\)，那么当然 \(y_i\) 也会很强，我们考虑怎么确定这个多项式呢？
点值唯一定理
假设我们有不同的 \(x_i\)，那么根据这些 \(<x_i,y_i>\) 我们可以逆推原多项式的 \(a_i\)，即根据数值求系数，并且情况是 唯一的。
可以补充一下怎么求解。
我们写成矩阵形式，即：
以上为范德蒙德卷积形式，只要 \(<x_0,x_1,\cdots,x_{n-1}>\) 互不相同，就有唯一解，这就是 DFT 的过程。
我们可以给个简单的证明。我们将上面的矩阵写为行列式。
我们每一次用前一列乘上 \(x_0\) 然后消去，可以得到：
很好！让我们紧接着提取公因式，得到：
以此类推我们就找到了递推式！然后我们就可以发现，只要 \(x_i\) 互不相同，就不会有 \(0\)，然后就有唯一解了。
非常好，所以现在我们只需要找到若干个点值就可以得到这个多项式的系数了，这就叫做 IDFT。
但是这样是 \(\mathcal O(n^2)\) 的复杂度傻大，所以咋办呢？
DFT
单位根
我们发现上面的做法瓶颈是啥呢？是确定完了 \(n\) 个 \(x_i\) 后求数值，这样就是 \(\mathcal O(n)\times O(n)\)，怎么办呢？
我们可以把目光放到 \(x_i\) 上，我们想有没有一些具有性质的数，能够降低复杂度呢？
答案是 单位根。
先来讲讲复数，我们定义 \(i=\sqrt{-1}\)，而我们设 \(z=a+bi\)，这就是一个复数了。复数的加减乘除和实数完全一致。
我们如果把普通平面直角坐标系放到复数系，那么我们就可以得到复平面的一个向量来表示复数。
我们可以表示 \(z=|z|(\cos\theta+i\sin\theta)\)，\(|z|\) 表示模长。
我们可以简写为 \(z=e^{i\theta}\)
那么单位根是个啥东东呢？我们假设 \(\omega_n\) 为 \(n\) 次单位根，那么 \(\omega_n=\cos{\frac{2\pi}{n}}+i\sin{\frac{2\pi}{n}}\)。
然后我们可以发现一些规律。
- \(\omega_n^{\frac{n}{2}}=-\omega_n^0\)
- \(\omega_n=\omega_0\)
- \(\omega_{nk}^{kd}=\omega_{n}^d\)
- 折半定理：\(\forall n,k\in\{k|k=2y,y\in N\},\omega_n^k=\omega_{\frac{n}{2}}^{\frac{k}{2}}\)
- \(\sum\limits_{i=0}^{n-1} \omega_n^i=0\)
- \(\omega_n^k=\omega_n^{k\% n}\)
根据以上性质，我们就可以快速计算单位根的幂次，这样就不用 \(\mathcal O(n)\) 算值了，而是 \(\mathcal O(\log n)\) 次！
具体的，我们设 \(<\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}>\)，\(y=\sum\limits_{i=0}^{n-1}a_ix^i\)，那么我们方便起见令 \(n=2^L-1\)，也就是比它大的第一个 \(2\) 的幂次减一。
我们把 \(y\) 的多项式分为奇偶项（根据次数）。我们得到 \(y_i=(a_0+a_2\omega_n^{2i}+a_4\omega_n^{4i}+\cdots+a_{n-2}\omega_n^{(n-2)i})+\omega_n^i(a_1+a_3\omega_n^{2i}+\cdots+a_{n-1}\omega_n^{(n-2)i})\)。
根据折半定理，\(y_i=(a_0+a_2\omega_{\frac{n}{2}}^{i}+a_4\omega_{\frac{n}{2}}^{2i}+\cdots+a_{n-2}\omega_{\frac{n}{2}}^{\frac{(n-2)i}{2}})+\omega_n^i(a_1+a_3\omega_{\frac{n}{2}}^{i}+\cdots+a_{n-1}\omega_{\frac{n}{2}}^{\frac{(n-2)i}{2}})\)。
我们令 \(Y_1(x)=a_0+a_2x^{i}+a_4x^{2i}+\cdots+a_{n-2}x^{\frac{(n-2)i}{2}}\)，\(Y_2(x)=a_1+a_3x^{i}+\cdots+a_{n-1}x^{\frac{(n-2)i}{2}}\)，那么 \(y_i=Y_1(\omega_{\frac{n}{2}}^{i})+\omega_n^iY_2(\omega_{\frac{n}{2}}^{i})\)。
但是这样我们的 \(i\) 还是 \(n\) 的级别，怎么与 \(\frac{n}{2}\) 建立联系呢？
我们可以发现当 \(i\in [0,\frac{n}{2})\)，我们就用 \(y_i=Y_1(\omega_{\frac{n}{2}}^{i})+\omega_n^iY_2(\omega_{\frac{n}{2}}^{i})\)。
对于 \(i\in[\frac{n}{2},n)\)，我们可以 \(y_i=Y_1(\omega_{\frac{n}{2}}^{i-\frac{n}{2}})+\omega_n^{i-\frac{n}{2}}\omega_n^{\frac{n}{2}}Y_2(\omega_{\frac{n}{2}}^{i-\frac{n}{2}})\)。
大功告成！这样无论 \(i\) 取何值，我们都有 \(n\) 到 \(\frac{n}{2}\) 的映射了！
而我们发现 \(\omega_n^{\frac{n}{2}}=-1\)，所以对于 \(i\in[\frac{n}{2},n)\)，我们可以 \(y_i=Y_1(\omega_{\frac{n}{2}}^{i-\frac{n}{2}})-\omega_n^{i-\frac{n}{2}}Y_2(\omega_{\frac{n}{2}}^{i-\frac{n}{2}})\)。
我们对于 \(i\in [0,\frac{n}{2})\)，可以这么写：
这样我们原来要求的是 \(y_0,\cdots,y_{n-1}\)，现在我们只需要求 \({Y_1}_0,\cdots,{Y_1}_{\frac{n}{2}-1}\) 和 \({Y_2}_0,\cdots,{Y_2}_{\frac{n}{2}-1}\)。
根据主定理，这显然是 \(\mathcal O(n\log n)\)。
那么问题来了，IDFT 怎么做呢？高斯消元的话就会得到比 \(\mathcal O(n^2)\) 还要恐怖的 \(\mathcal O(n^3)\)。
我们可以考虑带入单位根的负数次方，考虑 \(<\omega_n^0,\cdots,\omega_n^{n-1}>\)，取 \(<\frac{y_0}{n},\cdots,\frac{y_{n-1}}{n}>\)，那么：
所以我们只需要给单位根的指数集体乘上 \(-1\) 即可！
代码如下：
void DFT(cp *f,int l,int op)
{
if(!l) return 0;
cp fl[l],fr[l];
for(int i=0;i<l;i++)
{
fl[i]=f[i<<1];
fr[i]=f[i<<1|1];
}
DFT(fl,l>>1,op),DFT(fr,l>>1,op);
cp Wn=(cp){cos(PI/l),sin(PI/l)*op};
cp w=(cp){1,0};
for(int i=0;i<l;i++)
{
cp t=w*fr[i];
f[l]=fl[i]+t;
f[i+l]=fl[i]-t;
w=w*Wn;
}
}
FFT
很可恶的是递归常数傻大，那么怎么避免递归呢？开始我们的主角了：FFT。
非递归首先要解决奇偶分类问题。我们用一种方法：蝴蝶变换。
我们发现，我们如果从下往上合并，那么我们每一个数在递归最底层是什么样子呢？
比如 \(0\) 到 \(7\)，那么我们就可以发现递归最底层的顺序为 \(0,4,2,6,1,5,3,7\)。
然后我们写出所有的二进制 \(000,100,010,110,001,101,011,111\)。
如果把二进制反转一下，可以得到 \(000,001,010,011,100,101,110,111\)。
发现了什么？是不是反转完了就是二进制从小到大！
然后至于这个二进制反转为啥是人能想到的，我也不知道了，反正就是如下反转：
for(int i=0;i<len;i++)//len=2^L
r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
然后就是蝴蝶变换的 FFT 板子了：
void FFT(cp *f,int op)//op表示是FFT还是IFFT
{
for(int i=0;i<len;i++)
if(r[i]>i) swap(a[r[i]],a[i]);//交换
for(int i=1;i<len;i<<=1)
{
cp Wn=(cp){cos(PI/i),sin(PI/i)*op};
for(int j=0;j<len;j+=(i<<1))
{
cp w=(cp){1,0};
for(int k=0;k<i;k++)
{
cp Y1=f[j+k],Y2=w*f[j+i+k];
f[j+k]=Y1+Y2,f[j+k+i]=Y1-Y2;
w=w*Wn;
}
}
}
}
P3803 【模板】多项式乘法（FFT）
这道题是 FFT 板子。
#include<bits/stdc++.h>
using namespace std;
const double PI=acos(-1.0);
const int N=4e6+5;
struct cp
{
double x,y;
cp(double xx=0,double yy=0){x=xx,y=yy;}
friend cp operator +(cp a,cp b) {return cp(a.x+b.x,a.y+b.y);}
friend cp operator -(cp a,cp b) {return cp(a.x-b.x,a.y-b.y);}
friend cp operator *(cp a,cp b) {return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[N],g[N];
int n,m,pos[N],len=1,L;
void FFT(cp *f,int op)
{
for(int i=0;i<len;i++)
if(pos[i]>i) swap(f[pos[i]],f[i]);
for(int i=1;i<len;i<<=1)
{
cp Wn=cp(cos(PI/i),sin(PI/i)*op);
for(int j=0;j<len;j+=(i<<1))
{
cp W=cp(1,0);
for(int k=0;k<i;k++)
{
cp Y1=f[j+k],Y2=W*f[j+k+i];
f[j+k]=Y1+Y2,f[j+k+i]=Y1-Y2;
W=W*Wn;
}
}
}
}
int main()
{
cin.tie(0)->sync_with_stdio(0);
cout.tie(0)->sync_with_stdio(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;
for(;len<=n+m;len<<=1,L++);
for(int i=0;i<len;i++)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1),FFT(g,1);
for(int i=0;i<=len;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<=n+m;i++)
cout<<(int)(f[i].x/len+0.5)<<" ";
return 0;
}
P3338 [ZJOI2014] 力
这道题看上去很复杂，实际上很唐，我们来一步步分析。
\(E_i=\sum\limits_{j=1}^{i-1}\frac{q_j}{(i-j)^2}-\sum\limits_{j=i+1}^n\frac{q_j}{(i-j)^2}\)。
我们先考虑维护前半部分。对于前一半，我们可以看成 \(q_j\) 与 \((i-j)^2\) 的卷积。
而后半部分我们可以将 \(q\) 反转，然后我们就可以进行卷积加速了。
#include<bits/stdc++.h>
using namespace std;
const double PI=acos(-1.0);
const int N=4e5+5;
struct cp
{
double x,y;
cp (double xx=0,double yy=0) {x=xx,y=yy;};
friend cp operator +(const cp &a,const cp &b) {return cp(a.x+b.x,a.y+b.y);}
friend cp operator -(const cp &a,const cp &b) {return cp(a.x-b.x,a.y-b.y);}
friend cp operator *(const cp &a,const cp &b) {return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[N],g[N],h[N];
int n,len=1,L=0,pos[N];
void FFT(cp *f,int op)
{
for(int i=0;i<len;i++)
if(pos[i]>i) swap(f[pos[i]],f[i]);
for(int i=1;i<len;i<<=1)
{
cp Wn(cos(PI/i),sin(PI/i)*op);
for(int j=0;j<len;j+=(i<<1))
{
cp w(1,0);
for(int k=0;k<i;k++)
{
cp Y1=f[j+k],Y2=w*f[j+k+i];
f[j+k]=Y1+Y2,f[j+k+i]=Y1-Y2;
w=w*Wn;
}
}
}
}
int main()
{
cin>>n;
for(int i=1;i<=n;i++) cin>>f[i].x;
for(int i=1;i<=n;i++) g[i].x=1.0/i/i;
for(int i=1;i<=n;i++) h[n-i+1].x=f[i].x;
for(;len<=n*2+2;len<<=1,L++);
for(int i=0;i<len;i++)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1),FFT(g,1),FFT(h,1);
for(int i=0;i<len;i++)
{
f[i]=f[i]*g[i];
h[i]=h[i]*g[i];
}
FFT(f,-1),FFT(h,-1);
for(int i=1;i<=n;i++)
printf("%.4lf\n",f[i].x/len-h[n-i+1].x/len);
return 0;
}
P3723 [AH2017/HNOI2017] 礼物
这道题很简单，主要是看你数学能力如何。
我们假设你对第一个序列整体加上了 \(k\)。
我们可以化简一下。
那么相信大家都是学过二次函数的，显然关于 \(k\) 的是一个二次函数，我们求个最小值即可，\(2\sum\limits_{i=1}^nx_iy_i\) 是一个卷积，显然可以 FFT 优化，然后两个平方和也可以预处理，这道题就做完了。
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=5e5+5;
const double PI=acos(-1.0);
struct cp
{
long double x,y;
cp (double xx=0,double yy=0) {x=xx,y=yy;}
friend cp operator +(cp a,cp b) {return cp(a.x+b.x,a.y+b.y);}
friend cp operator -(cp a,cp b) {return cp(a.x-b.x,a.y-b.y);}
friend cp operator *(cp a,cp b) {return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[N],g[N];
int n,m,pos[N],len=1,ans,L,cnt,res[N<<2];
void FFT(cp *f,int op)
{
for(int i=0;i<len;i++)
if(pos[i]>i) swap(f[i],f[pos[i]]);
for(int i=1;i<len;i<<=1)
{
cp Wn(cos(PI/i),sin(PI/i)*op);
for(int j=0;j<len;j+=(i<<1))
{
cp w(1,0);
for(int k=0;k<i;k++)
{
cp Y1=f[j+k],Y2=w*f[i+j+k];
f[j+k]=Y1+Y2,f[i+j+k]=Y1-Y2;
w=w*Wn;
}
}
}
}
signed main()
{
cin>>n>>m;
for(int i=1;i<=n;i++) cin>>f[i].x;
for(int i=1;i<=n;i++) cin>>g[i].x;
for(int i=n+1;i<=(n<<1);i++) g[i].x=g[i-n].x;
for(int i=1;i<=n;i++)
{
ans+=f[i].x*f[i].x+g[i].x*g[i].x;
cnt+=g[i].x-f[i].x;
}
int c1=floor(cnt*1.0/n),c2=ceil(cnt*1.0/n);
ans+=min(n*c1*c1-2*c1*cnt,n*c2*c2-2*c2*cnt);
reverse(f+1,f+1+n);
for(;len<=n*3;len<<=1,L++);
for(int i=1;i<len;i++)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1),FFT(g,1);
for(int i=0;i<len;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=1;i<=n*3;i++) res[i]=(int)(f[i].x/len+0.5);
int mx=0;
for(int i=0;i<=n;i++) mx=max(mx,res[i+n+1]);
cout<<ans-mx*2;
return 0;
}
P4199 万径人踪灭
实在是不好评价这道题。
首先我们正难则反，我们考虑用总的回文子序列减去回文子串，那么回文子串可以用 manacher 线性解决。
我们来思考如何求出总的回文子序列。我们假设 \(f_i\) 表示以 \(i\) 为中心对称轴的回文串方法。那么我们可以得到：
这个怎么办呢？我们考虑这样一件事，我们设 \(g_k=[s_k='a']\)，\(h_k=[s_k='b']\)，那么我们可以发现：
然后就可以开心的 FFT 了。
那么我是哪里错了呢？
for(int k=0;k<i;k++,w=w*Wn)
{
cp Y1=f[j+k],Y2=w*f[j+k+i];//一开始写为Y2=f[j+k+i]，如何评价调了一晚上（（（（（（（
f[j+k]=Y1+Y2,f[j+k+i]=Y1-Y2;
}
完整代码：
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MOD=1000000007,N=1e5+5;
const double PI=acos(-1.0);
char s[N],c[N<<1];
struct cp
{
double x,y;
cp (double xx=0,double yy=0) {x=xx,y=yy;}
friend cp operator +(const cp &a,const cp &b) {return cp(a.x+b.x,a.y+b.y);}
friend cp operator -(const cp &a,const cp &b) {return cp(a.x-b.x,a.y-b.y);}
friend cp operator *(const cp &a,const cp &b) {return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[N<<2],g[N<<2];
int n,pos[N<<2],LEN,L,ans[N<<2],cnt,b[N<<2];
void init()
{
for(LEN=1;LEN<=n*2+2;LEN<<=1,L++);
for(int i=0;i<LEN;i++)
pos[i]=(pos[i>>1]>>1)|((i&1)<<(L-1));
}
void FFT(cp *f,int op)
{
for(int i=0;i<LEN;i++)
if(pos[i]>i) swap(f[pos[i]],f[i]);
for(int i=1;i<LEN;i<<=1)
{
cp Wn=cp(cos(PI/i),sin(PI/i)*op);
for(int j=0;j<LEN;j+=(i<<1))
{
cp w=cp(1,0);
for(int k=0;k<i;k++,w=w*Wn)
{
cp Y1=f[j+k],Y2=w*f[j+k+i];
f[j+k]=Y1+Y2,f[j+k+i]=Y1-Y2;
}
}
}
}
int qsm(int x,int a)
{
int res=1;
for(;a;a>>=1,x=x*x%MOD)
if(a&1) res=res*x%MOD;
return res;
}
int manacher()
{
c[0]=c[1]='#';
for(int i=1;i<=n;i++)
{
c[i<<1]=s[i];
c[i<<1|1]='#';
}
n=(n<<1)+2,c[n]=')';
for(int i=1,mr=0,mid=0;i<n;i++)
{
if(i<mr) b[i]=min(b[(mid<<1)-i],b[mid]+mid-i);
else b[i]=1;
for(;c[i+b[i]]==c[i-b[i]];b[i]++);
if(i+b[i]>mr) mr=i+b[i],mid=i;
}
int res=0;
for(int i=0;i<n;i++)
res=(res+(b[i]>>1))%MOD;
return res;
}
signed main()
{
cin>>(s+1);n=strlen(s+1);init();
for(int i=1;i<=n;i++)
if(s[i]=='a') f[i].x=g[i].x=1.0;
FFT(f,1),FFT(g,1);
for(int i=0;i<LEN;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<LEN;i++) ans[i]=(ans[i]+(int)(f[i].x/LEN+0.5));
memset(f,0,sizeof(f)),memset(g,0,sizeof(g));
for(int i=1;i<=n;i++)
if(s[i]=='b') f[i].x=g[i].x=1.0;
FFT(f,1),FFT(g,1);
for(int i=0;i<LEN;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=0;i<LEN;i++) ans[i]=(ans[i]+(int)(f[i].x/LEN+0.5));
for(int i=0;i<LEN;i++)
cnt=(cnt+qsm(2,(ans[i]+1)>>1)-1)%MOD;
cnt=(cnt-manacher()+MOD)%MOD;
cout<<cnt;
return 0;
}
NTT
如果你掌握了 FFT，那就掌握 NTT 了。这是为什么呢？
NTT 的本质就是把单位根变为原根。
那么原根和单位根有什么相似性呢？
假设模数 \(MOD\) 的原根为 \(g\)，那么我们假设 \(g_n=g^{\frac{p-1}{n}}\)，那么我们发现对于 \({g_n}^0,\cdots,{g_n}^{n-1}\)，我们 \({g_n}^n=g^{p-1}=g^{p-1}=g^0=1\)，这样我们就有了和单位根差不多的性质。
原根还有如下性质：
- \({g_n}^0=1,{g_n}^{\frac{n}{2}}=-1=-{g_n}^0\)。
- 相消定理：\(\forall n\ge 0,k\ge 0,d\ge 0,{g_{nd}}^{kd}={g_n}^k\)。
- 折半定理：\(\forall n\ge 0,k\ge 0\)，\(n,k\) 为偶数，\({g_{nd}}^{kd}={g_n}^k\)。
然后其他的就跟 FFT 一样了。
由于 NTT 涉及到原根，这里给出一些常见模数以及原根。
- \(998244353=7\times 17\times 2^{23}+1,g=3\)。
- \(167772161=5\times 2^{25}+1,g=3\)。
- \(469762049=7\times 2^{26}+1,g=3\)。
- \(1004535809=479\times 2^{21}+1,g=3\)。
- \(754974721=45\times 2^{24}+1,g=11\)。
void NTT(int *f,int op)//MOD=998244353
{
for(int i=0;i<K;i++)
if(i<r[i])
swap(f[i],f[r[i]]);
for(int i=2,i2=1;i<=K;i2=i,i<<=1)
{
int base=qpow(op?INV:M,(MOD-1)/i);
for(int j=0;j<K;j+=i)
{
int w=1;
for(int k=0;k<i2;k++,w=w*base%MOD)
{
int x=f[j+k],y=w*f[i2+j+k]%MOD;
f[j+k]=(x+y)%MOD,f[i2+j+k]=(x-y+MOD)%MOD;
}
}
}
int X=qpow(K,MOD-2);
if(op) for(int i=0;i<K;i++) f[i]=f[i]*X%MOD;
}
NTT 的优点在于，他能避免精度问题，但是不能处理系数为实数的多项式乘法，所以两个均要掌握。
浙公网安备 33010602011771号