多项式
1.FFT
复数
复数乘法的运算法则是模长相乘，幅角相加。
在单位圆中（模长为1），复数相乘只会幅角转动。
定义单位根 \(w_n\)，表示转 \(\frac{1}{n}\) 圆的角。
满足 \((w_n)^n=1\)。
单位根的 \(0\sim n-1\) 次幂均分整个圆，即满足 \((w_n^i)^n=1\).
它们互不相同。
引入
我们现在要求多项式乘法。
求 \(C=A\times B\).
1.取若干个数 \(1,w_n,w_n^2....,w_n^{n-1}\).
2.对这些数求出 \(A(w_n^i),B(w_n^i)\)，这一步称为 DFT.
3.求出 \(C(w_n^i)=A(w_n^i)\cdot B(w_n^i)\).
4.插值求出 \(C\) 的系数，这一步称为 IDFT.
DFT
给出 \(A\) 函数，求 \(A(w_n^i)\).
将 \(n\) 补成 \(2\) 的整数次幂。
分治：将 \(A\) 拆成奇偶两部分
\(A(x)=A_0(x^2)+xA_1(x^2)\)
我们往下分治，假设我们已知 \(A_0(w_{n/2}^i),A_1(w_{n/2}^i)\) 的值。
我们可以发现，\(A(w_n^i)=A_0(w_n^{2i})+w_n^i A_1(w_n^{2i})\).
\(A(w_n^i)=A_0(w_{n/2}^{i})+w_n^i A_1(w_{n/2}^{i})\)
就变成子问题求解。
IDFT
给出 \(C\) 函数在 \(w_n^i\) 的取值，求 \(C\) 的系数。
只需要把 \(w_n^i\) 的取值看成一个多项式，再求出 \(w_n^{-i}\) 的取值，最后除以 \(n\) 即为答案。
笔者尚且无法证明。
上述过程和 DFT 并无二致，只是一个符号的区别。
这样我们求出了多项式乘法。
卡常
我们设一个函数 \(f=A+Bi\).
DFT，再求出 \((f(w_n^i))^2\).
再 IDFT，求出 \(G=(A+Bi)^2=A^2-B^2+2ABi\)
只需取出 \(G\) 的虚部再除以 \(2\) 就是答案。
这样只需一次 DFT 和 IDFT，比原来少了一次 DFT.
实现
正常的写法，由于数组拷贝，会比较慢。
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e6+10;
const double pi=acos(-1);
int read() {
char ch=0;
while(ch<48||ch>57) ch=getchar();
return ch-'0';
}
struct CP {
double a,b;
CP operator + (const CP x) const {return (CP){a+x.a,b+x.b};}
CP operator - (const CP x) const {return (CP){a-x.a,b-x.b};}
CP operator * (const CP x) const {return (CP){a*x.a-b*x.b,x.b*a+b*x.a};}
};
void dft(vector<CP> &F) {
if(F.size()==1) return ;
int n=F.size();
vector<CP> F0,F1;
F0.resize(n/2); F1.resize(n/2);
for(int i=0; i<n; i+=2) F0[i/2]=F[i];
for(int i=1; i<n; i+=2) F1[i/2]=F[i];
dft(F0); dft(F1);
CP w=(CP){cos(2*pi/n),sin(2*pi/n)},buf=(CP){1,0};
for(int i=0; i<n/2; i++) {
CP tmp=buf*F1[i];
F[i]=F0[i]+tmp;
F[i+n/2]=F0[i]-tmp;
buf=buf*w;
}
return ;
}
void idft(vector<CP> &F) {
int n=F.size();
dft(F);
reverse(&F[1],&F[n]);
for(int i=0; i<n; i++) F[i].a/=n,F[i].b/=n;
}
int n,m,lim;
int main() {
scanf("%d%d",&n,&m);
lim=1; while(lim<n+m+1) lim<<=1;
vector<CP> F; F.resize(lim);
for(int i=0; i<=n; i++) F[i].a=read();
for(int i=0; i<=m; i++) F[i].b=read();
dft(F);
for(int i=0; i<lim; i++) F[i]=F[i]*F[i];
idft(F);
for(int i=0; i<=n+m; i++) {
int ret=(int)(F[i].b*0.5+0.49);
printf("%d ",ret);
}
return 0;
}
迭代写法：
我们发现原来的写法需要很多数组拷贝。
如果我们把一个位置的数奇偶变换的最后位置找出，
那么我们就无须拷贝新的数组了。
我们发现一个位置最后在它的二进制反转的位置。
看代码。
code
#include<bits/stdc++.h>
using namespace std;
const int N=4e6+10;
const double pi=acos(-1);
int read() {
char ch=0;
while(ch<48||ch>57) ch=getchar();
return ch-'0';
}
struct CP {
double a,b;
CP operator + (const CP x) const {return (CP){a+x.a,b+x.b};}
CP operator - (const CP x) const {return (CP){a-x.a,b-x.b};}
CP operator * (const CP x) const {return (CP){a*x.a-b*x.b,x.b*a+b*x.a};}
} F[N];
void dft(CP *F,int n) {
if(n==1) return ;
dft(&F[0],n/2); dft(&F[n/2],n/2);
CP w=(CP){cos(2*pi/n),sin(2*pi/n)},buf=(CP){1,0};
for(int i=0; i<n/2; i++) {
CP tmp=buf*F[i+n/2];
F[i+n/2]=F[i]-tmp;
F[i]=F[i]+tmp;
buf=buf*w;
}
return ;
}
void idft(CP *F,int n) {
dft(F,n);
reverse(&F[1],&F[n]);
for(int i=0; i<n; i++) F[i].a/=n,F[i].b/=n;
}
int n,m,lim;
int tr[N];
int main() {
scanf("%d%d",&n,&m);
lim=1; while(lim<n+m+1) lim<<=1;
for(int i=0; i<=n; i++) F[i].a=read();
for(int i=0; i<=m; i++) F[i].b=read();
for(int i=1; i<lim; i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?lim>>1:0);
for(int i=1; i<lim; i++)
if(i<tr[i]) swap(F[i],F[tr[i]]);
dft(F,lim);
for(int i=0; i<lim; i++) F[i]=F[i]*F[i];
for(int i=1; i<lim; i++)
if(i<tr[i]) swap(F[i],F[tr[i]]);
idft(F,lim);
for(int i=0; i<=n+m; i++) {
int ret=(int)(F[i].b*0.5+0.49);
printf("%d ",ret);
}
return 0;
}
2.NTT
原根
如果我们想要提升精度，那么我们可以用取模计算代替复数运算。
对于取模 \(p\)，我们发现了原根 \(g\)，满足 \(g\) 的 \(1\sim p-1\) 次幂互不相同。
原根的性质与单位根的性质是相似的。
我们可以构造 \(w_n=g^{(p-1)/n}\).
但这要求 \(p-1\) 含有很多 \(2\) 的因子。
常见的原根：\(g_{998244353}=3\).
引入
大体与 FFT 没什么不同。
实现
code
#include<bits/stdc++.h>
using namespace std;
const int N=1e6+10;
const int p=998244353,g=3;
int qpow(int a,int b) {
int res=1;
for(; b; b>>=1) {
if(b&1) res=1ll*res*a%p;
a=1ll*a*a%p;
}
return res;
}
void ntt(vector<int> &F,int op) {
if(F.size()==1) return ;
int n=F.size();
vector<int> F0,F1;
F0.resize(n/2); F1.resize(n/2);
for(int i=0; i<n; i+=2) F0[i/2]=F[i];
for(int i=1; i<n; i+=2) F1[i/2]=F[i];
ntt(F0,op); ntt(F1,op);
int w=qpow(g,(p-1)/n),buf=1;
if(op==-1) w=qpow(w,p-2);
for(int i=0; i<n/2; i++) {
int tmp=1ll*F1[i]*buf%p;
F[i]=(tmp+F0[i])%p;
F[i+n/2]=(-tmp+p+F0[i])%p;
buf=1ll*buf*w%p;
}
return ;
}
int n,m,lim;
int main() {
scanf("%d%d",&n,&m);
lim=1; while(lim<n+m+1) lim<<=1;
vector<int> F,G;
F.resize(lim); G.resize(lim);
for(int i=0; i<lim; i++) F[i]=G[i]=0;
for(int i=0; i<=n; i++) scanf("%d",&F[i]);
for(int i=0; i<=m; i++) scanf("%d",&G[i]);
ntt(F,1); ntt(G,1);
for(int i=0; i<lim; i++) F[i]=1ll*F[i]*G[i]%p;
ntt(F,-1);
int inv=qpow(lim,p-2);
for(int i=0; i<=n+m; i++) printf("%d ",(int)(1ll*F[i]*inv%p));
return 0;
}
3.几个模板
任意模数多项式乘法
如果按照原来的值域计算，那么可能会导致运算过程中值超过了 long long 范围。
可以取几个模数来 NTT，最后 CRT，但是需要用到一些科技，不推荐。
我们考虑将一个数拆成 \(a\times base+b\) 的形式，\(base\) 可以取根号值域，取 \(2^{15}\) 即可。
这样最大的值域是可接受的。
对于 \(A,B\) 两个多项式，
最后结果是 \(A_1B_1base^2+(A_1B_2+A_2B_1)base+A_2B_2\).
其中 \(A=A_1base+A_2\)，\(B=B_1base+B_2\).
这样要用到总共 8 次 DFT 和 IDFT，常数过大。
我们设 \(P_1=A_1+A_2 i,P_2=A_1-A_2 i,Q=B_1+B_2 i\)
考虑求 \(P_1 Q\) 和 \(P_2 Q\).
\(P_1 Q =A_1B_1-A_2B_2+(A_1B_2+A_2B_1)i\).
\(P_2 Q =A_1B_1+A_2B_2+(A_1B_2-A_2B_1)i\).
两式相加，相减，分离实，虚部，得解。
这样只需 5 次。
code
#include<bits/stdc++.h>
using namespace std;
const int N=4e6+10;
const long double pi=acos(-1);
struct CP {
long double a,b;
CP operator + (const CP x) const {return (CP){a+x.a,b+x.b};}
CP operator - (const CP x) const {return (CP){a-x.a,b-x.b};}
CP operator * (const CP x) const {return (CP){a*x.a-b*x.b,x.b*a+b*x.a};}
} F1[N],F2[N],Q[N];
void dft(CP *F,int n) {
if(n==1) return ;
dft(&F[0],n/2); dft(&F[n/2],n/2);
CP w=(CP){cos(2*pi/n),sin(2*pi/n)},buf=(CP){1,0};
for(int i=0; i<n/2; i++) {
CP tmp=buf*F[i+n/2];
F[i+n/2]=F[i]-tmp;
F[i]=F[i]+tmp;
buf=buf*w;
}
return ;
}
void idft(CP *F,int n) {
dft(F,n);
reverse(&F[1],&F[n]);
}
int n,m,p,lim;
int tr[N];
const int base=32000;
int main() {
scanf("%d%d%d",&n,&m,&p);
lim=1; while(lim<n+m+1) lim<<=1;
for(int i=0; i<=n; i++) {
int x; scanf("%d",&x);
F1[i]=(CP){1.0*(x/base),1.0*(x%base)};
F2[i]=(CP){1.0*(x/base),1.0*(-(x%base))};
}
for(int i=0,x; i<=m; i++) {
scanf("%d",&x);
Q[i]=(CP){1.0*(x/base),1.0*(x%base)};
}
for(int i=1; i<lim; i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?lim>>1:0);
for(int i=1; i<lim; i++)
if(i<tr[i]) {
swap(F1[i],F1[tr[i]]);
swap(F2[i],F2[tr[i]]);
swap(Q[i],Q[tr[i]]);
}
dft(F1,lim); dft(F2,lim); dft(Q,lim);
for(int i=0; i<lim; i++) Q[i].a=Q[i].a/lim,Q[i].b=Q[i].b/lim;
for(int i=0; i<lim; i++) {
F1[i]=F1[i]*Q[i];
F2[i]=F2[i]*Q[i];
}
for(int i=1; i<lim; i++)
if(i<tr[i]) {
swap(F1[i],F1[tr[i]]);
swap(F2[i],F2[tr[i]]);
}
idft(F1,lim); idft(F2,lim);
for(int i=0; i<=n+m; i++) {
long long a1b1,a2b1,a1b2,a2b2;
a1b1=((long long)floor((F1[i].a+F2[i].a)*0.5+0.49))%p;
a1b2=((long long)floor((F1[i].b+F2[i].b)*0.5+0.49))%p;
a2b1=((long long)floor(F1[i].b+0.49)-a1b2)%p;
a2b2=((long long)floor(F2[i].a+0.49)-a1b1)%p;
printf("%lld ",((1ll*base*base%p*a1b1%p+1ll*base*(a1b2+a2b1)%p+a2b2)%p+p)%p);
}
return 0;
}
多项式乘法逆
求 \(f^{-1}(x)\times f(x)=1\pmod {x^n}\)
我们假设已经求出了 \(f_*^{-1}(x)\times f(x)=1\pmod {x^{n/2}}\) （\(n/2\) 向上取整）
而 \(f^{-1}(x)\times f(x)=1\pmod {x^{n/2}}\)
\(f_*^{-1}(x)-f^{-1}(x)=0\pmod {x^{n/2}}\)
\(f_*^{-2}(x)+f^{-2}(x)-2f_*^{-1}(x)f^{-1}(x)=0\pmod {x^n}\)
\(f_*^{-2}(x)f(x)+f^{-1}(x)-2f_*^{-1}(x)=0\pmod {x^n}\)
\(f^{-1}(x)=(2-f_*^{-1}(x)\times f(x))f_*^{-1}(x)\)
递归 +NTT 即可。
边界： \(n=1\) 时，为 \(f(0)\) 的逆元。
code
#include<bits/stdc++.h>
using namespace std;
const int N=4e6+10;
const int p=998244353,g=3,gi=332748118;
int qpow(int a,int b) {
int res=1;
for(; b; b>>=1) {
if(b&1) res=1ll*res*a%p;
a=1ll*a*a%p;
}
return res;
}
int n,m,lim,tr[N];
void ntt(int *F,int n,int op) {
if(n==1) return ;
ntt(&F[0],n/2,op); ntt(&F[n/2],n/2,op);
int w=qpow(op==1?g:gi,(p-1)/n),buf=1;
for(int i=0; i<n/2; i++) {
int tmp=1ll*F[i+n/2]*buf%p;
F[i+n/2]=(F[i]-tmp+p)%p;
F[i]=(F[i]+tmp)%p;
buf=1ll*buf*w%p;
}
return ;
}
int c[N];
void inv(int *F,int *G,int n) {
if(n==1) {
G[0]=qpow(F[0],p-2);
return ;
}
inv(F,G,(n+1)/2);
for(int i=0; i<n; i++) c[i]=F[i];
for(int i=n; i<lim; i++) c[i]=0;
lim=1;
while(lim<n*2) lim<<=1;
for(int i=1; i<lim; i++)
tr[i]=(tr[i>>1]>>1)|(i&1?(lim>>1):0);
for(int i=1; i<lim; i++) {
if(i<tr[i]) swap(c[i],c[tr[i]]),swap(G[i],G[tr[i]]);
}
ntt(c,lim,1); ntt(G,lim,1);
for(int i=0; i<lim; i++) G[i]=1ll*G[i]*(-1ll*c[i]*G[i]%p+2ll)%p;
for(int i=1; i<lim; i++)
if(i<tr[i]) swap(G[i],G[tr[i]]);
ntt(G,lim,-1);
int iv=qpow(lim,p-2);
for(int i=0; i<lim; i++) G[i]=(1ll*G[i]*iv%p+p)%p;
for(int i=n; i<lim; i++) G[i]=0;
}
int a[N],ans[N];
int main() {
scanf("%d",&n);
for(int i=0; i<n; i++) scanf("%d",&a[i]),a[i]=(a[i]+p)%p;
inv(a,ans,n);
for(int i=0; i<n; i++) printf("%d ",ans[i]);
return 0;
}
3.一些典题
P3723 [AH2017/HNOI2017] 礼物
考虑求 \(\min{\sum (a_i-b_i+c)^2}\),\(b_i\) 可以循环移位。
考虑将 \(b\) 复制一遍放到 \(i+1 \sim 2n\) 处。
拆式子： \(\sum (a_i^2+b_i^2)+nc^2+\sum (a_i-b_i) c -2\sum a_ib_i\)
\(\sum (a_i^2+b_i^2)\) 是定值，\(nc^2+\sum (a_i-b_i) c\) 是二次函数。
最后求 \(\max{\sum a_i b_{i+k}}\)，发现是“差卷积形式”。
将 \(a\) 翻转为 \(a'\)，\(\max a'_{n-i+1} b_{i+k}\) 是标准卷积。
将 \(a',b\) 转为多项式卷积，最后 \(n+k+1\) 就是偏移 \(k\) 位的答案。
浙公网安备 33010602011771号