@KirinBill
2018-07-07T22:28:11.000000Z
字数 9182
阅读 2383
多项式
目录
这里主要写一下基于FFT的多项式运算,应用方面主要是生成函数,和特征多项式优化线性递推式,见那两篇博客。
感觉还是要先看一下生成函数那个博客,可以不会写代码,但是要知道生成函数是什么东西、怎么列式子,不然会不知道下面要求啥以及是干啥用的...
前置技能:[FFT优化的多项式乘法]("快速傅里叶变换")
...
都说了前置技能...
用到了求逆,最后说。
不会的话去问数学老师...
给定多项式,求,。
循环卷积的快速幂:
inline void mod_conv(int a[],int b[],int c[],int len){
static int tmpa[MAXL],tmpb[MAXL],tmpc[MAXL];
memcpy(tmpa,a,len<<2);
memcpy(tmpb,b,len<<2);
poly_mul(tmpa,tmpb,tmpc,len);
memset(c,0,len<<2);
for(int i=0;i<(len<<1);++i)
c[i]=(c[i]+tmpc[i%len])%P;
}
inline void poly_pow(int x[],int y,int ret[],int len){
static int tmpx[MAXL];
memcpy(tmpx,x,len<<2);
memset(ret,0,len<<2);
ret[0]=1;
for(;y;y>>=1){
if(y&1) mod_conv(ret,tmpx,ret,len);
mod_conv(tmpx,tmpx,tmpx,len);
}
}
已知一个多项式,我们要求的是函数,其中函数是一个初等函数,但是自变量为多项式。这个操作直观意义不是很明显,主要是解决一些生成函数问题用的,所以可以理解为将多项式带入函数的泰勒展开式,譬如就是求:,而并非真要算的“多项式次方”。我们都知道生成函数是形式幂级数,在计算时只需要保留有用的前项。而初等函数运算就是将多项式带入一个泰勒展开式,也只能得到无穷级数,而它主要是为生成函数运算服务的,所以我们也可以把它放在意义下求,这样就方便计算了。
注意:下文有时为了突出多项式是函数自变量,会将多项式写成,所有大写的之类的则代表自变量为多项式的函数。
也就是函数,但这个显然无法快速幂,而且这个由于比较特殊,相当于乘法逆元(是后面基本上所有东西的基础),我们要单独想一种做法...
代码:
inline void poly_inv(int a[],int inva[],int len){
static int tmpa[MAXL];
memset(tmpa,0,len<<1<<2);
memset(inva,0,len<<1<<2);
inva[0]=inv(a[0]);
for(int i=2;i<=len;i<<=1){
memcpy(tmpa,a,i<<2);
NTT(tmpa,i<<1,1),NTT(inva,i<<1,1);
for(int j=0;j<(i<<1);++j)
inva[j]=(2-(long long)tmpa[j]*inva[j]%P+P)*inva[j]%P;
NTT(inva,i<<1,-1),div_len(inva,i<<1);
memset(inva+i,0,i<<2);
}
}
主要是配合后面的来用的。
由于在出无定义,我们不方便泰勒展开,所以用了一种同后面牛顿迭代不同的方法。
设:
求导和积分都是的,求个逆再乘一下就好了,复杂度。
代码:
inline void poly_ln(int a[],int lna[],int len){
static int da[MAXL],inva[MAXL],dlna[MAXL];
memset(da,0,len<<1<<2);
poly_diff(a,da,len);
poly_inv(a,inva,len);
poly_mul(da,inva,dlna,len<<1);
poly_integ(dlna,lna,len);
}
求函数的零点常用的一种方法是牛顿迭代法,其公式可以解释为每次拿在当前的近似解处的切线与的交点作为新的解。这样不断迭代,通常能收敛到零点处。但是这个解释是有些片面的,真正的原理是每次在当前解处将进行泰勒展开:,带入我们要求的零点,同时省略后面的项,可以得到一个方程解出一个。由于上面的省略操作,这个是近似的解。
对于多项式的一些运算,我们可以化为求一个自变量为多项式的函数的零点。考虑将牛顿迭代法推广过来:
代码:
inline void poly_sqrt(int a[],int sqrta[],int len){
static const int INV2=inv(2);
static int tmpa[MAXL],invsqrta[MAXL];
memset(tmpa,0,len<<1<<2);
sqrta[0]=sqrt(a[0]); //通常a[0]=1
for(int i=2;i<=len;i<<=1){
memcpy(tmpa,a,i<<2);
poly_inv(sqrta,invsqrta,i);
NTT(tmpa,i<<1,1),NTT(sqrta,i<<1,1),NTT(invsqrta,i<<1,1);
for(int j=0;j<(i<<1);++j)
sqrta[j]=((long long)sqrta[j]*sqrta[j]%P+tmpa[j])%P*INV2%P*invsqrta[j]%P;
NTT(sqrta,i<<1,-1),div_len(sqrta,i<<1);
memset(sqrta+i,0,i<<2);
}
}
代码:
inline void poly_exp(int a[],int expa[],int len){
static int tmpa[MAXL],lnexpa[MAXL];
memset(tmpa,0,len<<1<<2);
memset(expa,0,len<1<<2);
assert(a[0]==0),expa[0]=1;
for(int i=2;i<=len;i<<=1){
memcpy(tmpa,a,i<<2);
poly_ln(expa,lnexpa,i);
NTT(tmpa,i<<1,1),NTT(expa,i<<1,1),NTT(lnexpa,i<<1,1);
for(int j=0;j<(i<<1);++j)
expa[j]=(long long)expa[j]*(tmpa[j]+1-lnexpa[j]+P)%P;
NTT(expa,i<<1,-1),div_len(expa,i<<1);
memset(expa+i,0,i<<2);
}
}
多项式的幂:
inline void poly_pow(int a[],int k[],int powa[],int len){
static int lna[MAXL],klna[MAXL];
poly_ln(a,lna,len);
poly_mul(lna,k,klna,len);
poly_exp(klna,powa,len);
}
给定次多项式和次多项式(通常),求满足:的多项式和。显然可以类比的认为它俩分别是商和余数。
代码:
inline void poly_div(int a[],int n,int b[],int m,int d[]){
static int ra[MAXL],rb[MAXL],invrb[MAXL];
if(m>n){
d[0]=0;
return;
}
int len=poly_len(n-m);
memset(ra,0,len<<1<<2);
for(int i=0;i<=n-m;++i)
ra[i]=a[n-i];
memset(rb,0,len<<1<<2);
for(int i=0;i<=min(n-m,m);++i)
rb[i]=b[m-i];
poly_inv(rb,invrb,len);
poly_mul(ra,n-m,invrb,n-m,d,n-m);
reverse(d,d+n-m+1);
}
inline void poly_mod(int a[],int n,int b[],int m,int r[]){
static int d[MAXL],tmp[MAXL];
if(m>n){
memcpy(r,a,n+1<<2);
return;
}
poly_div(a,n,b,m,d);
poly_mul(b,m,d,n-m,tmp,m-1);
poly_minus(a,tmp,r,m-1);
}
给一个次多项式和个点,求。假设同阶(实际上复杂度瓶颈是关于的)。
O2
要跑快2s了。代码:
void poly_DC(int d,int l,int r,int x[],int dc[][MAXN<<1]){
if(l==r){
dc[d][l<<1]=P-x[l],dc[d][l<<1|1]=1;
return;
}
int mid=l+r>>1;
poly_DC(d+1,l,mid,x,dc),poly_DC(d+1,mid+1,r,x,dc);
poly_mul(dc[d+1]+(l<<1),mid-l+1,dc[d+1]+(mid+1<<1),r-mid,dc[d]+(l<<1),r-l+1);
}
void poly_eval(int d,int l,int r,int f[],int len,int dc[][MAXN<<1],int y[]){
static int eval[MAXD][MAXN];
poly_mod(f,len,dc[d]+(l<<1),r-l+1,eval[d]+l);
if(l==r){
y[l]=eval[d][l];
return;
}
int mid=l+r>>1;
poly_eval(d+1,l,mid,eval[d]+l,r-l,dc,y);
poly_eval(d+1,mid+1,r,eval[d]+l,r-l,dc,y);
}
inline void poly_eval(int f[],int len,int x[],int n,int y[]){
static int dc[MAXD][MAXN<<1];
poly_DC(0,0,n-1,x,dc);
poly_eval(0,0,n-1,f,len,dc,y);
}
给定个点值,要求插值出一个次多项式来。所谓插值,就是把这个值“插到”多项式里,也就是多项式在这个点的取值是给定值。
代码:
void poly_interp(int d,int l,int r,int val[],int dc[][MAXN<<1],int f[]){
if(l==r){
f[l]=val[l];
return;
}
int mid=l+r>>1;
poly_interp(d+1,l,mid,val,dc,f);
poly_interp(d+1,mid+1,r,val,dc,f);
static int g[MAXN],h[MAXN];
poly_mul(f+l,mid-l,dc[d+1]+(mid+1<<1),r-mid,g,r-l);
poly_mul(f+mid+1,r-mid-1,dc[d+1]+(l<<1),mid-l+1,h,r-l);
poly_add(g,h,f+l,r-l);
}
inline void poly_interp(int x[],int y[],int n,int f[]){
static int dc[MAXD][MAXN<<1],tmp[MAXN],val[MAXN];
poly_DC(0,0,n-1,x,dc);
poly_diff(dc[0],n,tmp);
poly_eval(tmp,n-1,x,n,val);
for(int i=0;i<n;++i)
val[i]=(long long)y[i]*inv(val[i])%P;
poly_interp(0,0,n-1,val,dc,f);
}
ps:上面讲的多项式操作复杂度都做到了理论下界,但是FFT和NTT常数那么大,跟多一两个似的也是没谁了。