不知不觉,就开始学习多项式了,好快啊
重新拾起我之前学的百嘛不是的\(FFT\),新学\(NTT/FWT\)
开始吧学习博客
FFT:快速傅里叶变换
总的来说就是快速求取两个多项式的乘积
把板子粘在这里Luogu3803 多项式乘法
code
#include<bits/stdc++.h>
using namespace std;
#define fo(i,x,y) for(int i=(x);i<=(y);i++)
#define fu(i,x,y) for(int i=(x);i>=(y);i--)
const double pi=acos(-1.0);
const int N=1<<21;
struct coex{
double r,i;
coex(){}
coex(double x,double y){r=x;i=y;}
coex operator + (coex a){return coex(r+a.r,i+a.i);}
coex operator - (coex a){return coex(r-a.r,i-a.i);}
coex operator * (coex a){return coex(r*a.r-i*a.i,r*a.i+i*a.r);}
}a[N],b[N],w[N];
int af[N],la,lb,lim,len;
void fft(coex *a,int lim){
fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]);
for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
coex tmp=w[t*j]*a[i+j+d];
a[i+j+d]=a[i+j]-tmp;
a[i+j]=a[i+j]+tmp;
}
}
signed main(){
scanf("%d%d",&la,&lb);
fo(i,0,la)scanf("%lf",&a[i].r);
fo(i,0,lb)scanf("%lf",&b[i].r);
for(lim=1,len=0;lim<=la+lb;lim<<=1,len++);
fo(i,0,lim-1){
af[i]=(af[i>>1]>>1)|((i&1)<<(len-1));
w[i]=coex(cos(2.0*i*pi/lim),sin(2.0*i*pi/lim));
}
fft(a,lim);fft(b,lim);
fo(i,0,lim-1)a[i]=a[i]*b[i],w[i].i=-w[i].i;
fft(a,lim);
fo(i,0,la+lb)printf("%d ",(int)(a[i].r/lim+0.5));
}
需要注意这么几点:
1、\(lim\)要大于多项式次数,严格大于
2、搞清楚\(d\)和\(t\)的关系,以及各种2倍
3、背板子,背就完事了。
4、最后答案别忘了除以\(lim\),预处理单位根别忘了\(2*i*pi/lim\),别忘了除以\(lim\)
NTT:快速数论变换
code
void ntt(int *a,int lim){
fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]);
for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=g[t*j]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod;
a[i+j]=(a[i+j]+tmp)%mod;
}
}
int bas[M],ans[M];
int a[M],b[M];
void mul(int *x,int *y,int *z){
fo(i,0,lim-1)a[i]=y[i];
fo(i,0,lim-1)b[i]=z[i];
g[0]=1;g[1]=ksm(3,(mod-1)/lim,mod);
fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod;
ntt(a,lim);ntt(b,lim);
g[0]=1;g[1]=ksm(g[1],mod-2,mod);
fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod;
fo(i,0,lim-1)a[i]=a[i]*b[i]%mod;
ntt(a,lim);int inv=ksm(lim,mod-2,mod);
fo(i,0,lim-1)a[i]=a[i]*inv%mod;
fo(i,0,m-2)a[i]=(a[i]+a[i+m-1])%mod;
fo(i,0,m-2)x[i]=a[i];
}
注意:
1、注意模数,以及自己找原根
2、注意最后出来的时候除以\(lim\)
3、别忘了变换原根
多项式求逆
就是求某一个多项式的逆元\(mod\ x^n\)意义下
主要吧,这个板子还挺难背的
递推上来的证明过程
所以结论就是\(B=2B'-AB'^2\)(\(B'\)表示\(mod\ x^{\lceil\frac{n}{2}\rceil}\)意义下的逆元)
code
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define fo(i,x,y) for(int i=(x);i<=(y);i++)
#define fu(i,x,y) for(int i=(x);i>=(y);i--)
int read(){
int s=0,t=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')t=-1;ch=getchar();}
while(isdigit(ch)){s=s*10+ch-'0';ch=getchar();}
return s*t;
}
const int mod=998244353;
const int N=1<<18;
int ksm(int x,int y){
int ret=1;
while(y){
if(y&1)ret=ret*x%mod;
x=x*x%mod;y>>=1;
}return ret;
}
int af[N],g[N],lim,len;
void ntt(int *a,int lim){
fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]);
for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=g[t*j]*a[i+j+d]%mod;
a[i+j+d]=(a[i+j]-tmp+mod)%mod;
a[i+j]=(a[i+j]+tmp)%mod;
}
}
int a[N],b[N],c[N];
void calc(int deg){//deg表示当前mod的是x^deg
if(deg==1){return b[0]=ksm(a[0],mod-2),void();}
calc((deg+1)>>1);
for(lim=1,len=0;lim<deg*2;lim<<=1,len++);
fo(i,0,lim-1)af[i]=(af[i>>1]>>1)|((i&1)<<(len-1));
g[0]=1;g[1]=ksm(3,(mod-1)/lim);
fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod;
fo(i,0,deg-1)c[i]=a[i];//注意这里只赋值前deg-1位,因为这次乘法只能用到这些位
//这里a有deg-1位,而b只有deg+1>>1位,因为b是要平方的,并且这里是求a的deg-1位的逆元
ntt(c,lim);ntt(b,lim);
g[0]=1;g[1]=ksm(g[1],mod-2);
fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod;
fo(i,0,lim-1)b[i]=(2-c[i]*b[i]%mod+mod)*b[i]%mod;
ntt(b,lim);int iv=ksm(lim,mod-2);
fo(i,0,lim-1)b[i]=b[i]*iv%mod;
fo(i,deg,lim-1)b[i]=0;//把后边的清空,防止影响下一次运算
fo(i,0,lim-1)c[i]=0;//清空
}
int n;
signed main(){
n=read();
fo(i,0,n-1)a[i]=read();
calc(n);
fo(i,0,n-1)printf("%lld ",b[i]);
}
FWT:快速沃尔什变换
发现我们所做的卷积都是\(c=\sum_{i=1}^{n}\sum_{j=1}^{i}a_j*b_{i-j}\)
也就是\(FFT\)和\(NTT\)能处理的加法卷积,如果是与?或?异或呢?
就要用到\(FWT\),也就是快速位运算卷积
|:或卷积
这个要是直接卷的话,好像有点困难,但是高维前缀和就是或卷积,所以我们可以用高维前缀和解决或卷积的问题
然而高维前缀和还是没有\(FWT\)来的直接
\(FWT\)也和\(FFT\)一样有变换和逆变换,这个直接看最上面我给出的那个博客吧
这次变换不能叫点值表达式了,应该叫它子集和表达式
相当于是枚举每一个二进制位来求子集和,逆变换拆回去的时候是一样的
code
void ort(int *a,int lim,int tp){//tp=1表示顺变换,tp=2表示逆变换
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j+d]=(a[i+j+d]+a[i+j]*tp)%mod;
}
}
&:与卷积
这应该叫做高维后缀和吧,我们可以用高维后缀和解决与卷积问题
变换就是后缀和,逆变换一样
code
void andt(int *a,int lim,int tp){//tp=1表示顺变换,tp=2表示逆变换
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j]=(a[i+j]+a[i+j+d]*tp)%mod;
}
}
^:异或卷积
这个应该叫高维后缀差分
好像有点难以概括,额,不概括了,直接看码吧
code
void xort(int *a,int lim,int tp){//tp=1表示顺变换,tp=0.5表示逆变换
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=a[i+j+d];
a[i+j+d]=(a[i+j]-tmp+mod)*tp%mod;
a[i+j]=(a[i+j]+tmp)*tp%mod;
}
}
高维前缀和的求法,枚举每一维,做一遍前缀和
在位运算上也就是枚举每一个二进制位,将这一位为\(0\)的加到这一位为\(1\)的上面
这样做高维前缀和虽然朴素,但是有正确性的保证,保证了不重不漏
code
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define fo(i,x,y) for(int i=(x);i<=(y);i++)
#define fu(i,x,y) for(int i=(x);i>=(y);i--)
int read(){
int s=0,t=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')t=-1;ch=getchar();}
while(isdigit(ch)){s=s*10+ch-'0';ch=getchar();}
return s*t;
}
const int mod=998244353;
const int N=1<<18;
int ksm(int x,int y){
int ret=1;
while(y){
if(y&1)ret=ret*x%mod;
x=x*x%mod;y>>=1;
}return ret;
}
int n,lim,a[N],b[N],c[N],aa[N],bb[N];
void init(int lim){fo(i,0,lim-1)a[i]=aa[i],b[i]=bb[i];}
void ort(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j+d]=(a[i+j+d]+a[i+j]*tp)%mod;
}
}
void andt(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j]=(a[i+j]+a[i+j+d]*tp)%mod;
}
}
void xort(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=a[i+j+d];
a[i+j+d]=(a[i+j]-tmp+mod)*tp%mod;
a[i+j]=(a[i+j]+tmp)*tp%mod;
}
}
void mix(int lim){fo(i,0,lim-1)c[i]=a[i]*b[i]%mod;}
void pt(int *a,int lim){fo(i,0,lim-1)printf("%lld ",a[i]);printf("\n");}
signed main(){
n=read();lim=1<<n;
fo(i,0,lim-1)aa[i]=read();
fo(i,0,lim-1)bb[i]=read();
init(lim);ort(a,lim,1);ort(b,lim,1);mix(lim);ort(c,lim,mod-1);pt(c,lim);
init(lim);andt(a,lim,1);andt(b,lim,1);mix(lim);andt(c,lim,mod-1);pt(c,lim);
init(lim);xort(a,lim,1);xort(b,lim,1);mix(lim);xort(c,lim,ksm(2,mod-2));pt(c,lim);
}
技巧
1、生成函数时,有关于对称的问题,如果对称的话,那么两个点的位置相加是对称中心的两倍,意思就是多项式卷积之后对称轴是同一条的都在一起
2、对于生成函数的利用,我们只需要系数,而并不是真正的把\(x\)带进去找值,我只需要系数就行了
比如:\(\sum_{i=0}^{n}\sum_{j=0}^{i}j!*2^{i-j}\),这个东西阶乘和幂可以分别看作是两个多项式的系数
这样我们可以通过卷积快速求出后面的和,也就是卷积后多项式的第\(i\)项系数,符合多项式乘法的规则,然后带入求和