[Notes][多项式]杂记 · 多项式算法—多项式求逆 多项式取模 多项式开根…

多项式

由若干个单项式相加组成的代数式叫做多项式
形如:f(x)=i=0naixi
deg f(x)称为f的度,是f(x)最高次项的次数。

生成函数

形如i=0aixi
生成函数又称母函数,往往和多项式算法联系起来起到优化转移的作用

多项式算法

加减法

多项式加减法比较简单
f(x)=i=0naixig(x)=i=0nbixi

(f±g)(x)=i=0n(ai±bi)xi

代码实现也比较简单


乘法

多项式乘法是所有多项式算法的基础,也是生成函数的卷积运算

f(x)=i=0naixig(x)=i=0nbixi

(f×g)(x)=i=02nj=0iaj×bijxi

朴素的多项式乘法是O(n2)的,但FFT(快速傅里叶变换)运用复数根的特性可以在O(nlogn)的复杂度内实现多项式的系数表达和点值表达的转化,运用点值表达实现O(n)的相乘,总复杂度就是O(nlogn),但是因为常数过大,往往小范围会使用暴力。

inline void FFT(E *a,int r){
  for(int i=0;i<n;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);
  for(int i=1;i<n;i<<=1){
    E wn(cos(M_PI/i),r*sin(M_PI/i));
    for(int j=0;j<n;j+=(i<<1)){
      E w(1,0);
      for(int k=0;k<i;k++,w=w*wn){
    E x=a[j+k],y=w*a[j+k+i];
    a[j+k]=x+y; a[j+k+i]=x-y;
      }
    }
  }
  if(r==-1) for(int i=0;i<n;i++) a[i].real/=n;
}

运用生成函数的计数问题往往会对一个质数取模,如果这个质数可以表示成k×2n+1,其中k为奇数n大于多项式的度数,那么就可以用这个质数的原根代替复数根,实现多项式的系数对一个质数取模,这个算法就是NTT(快速数论变换),这种模数称为NTT模数。

inline void Pre(int n){
  num=n;
  int g=Pow(3,(P-1)/num);
  w[0][0]=w[1][0]=1; for(int i=1;i<num;i++) w[0][i]=1LL*w[0][i-1]*g%P;
  for(int i=1;i<num;i++) w[1][i]=w[0][num-i];
}

inline void NTT(int *a,int n,int r){
  for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
  for(int i=1;i<n;i<<=1)
    for(int j=0;j<n;j+=i<<1)
      for(int k=0;k<i;k++){
        int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
        a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
      }
  if(!r) for(int i=0,inv=Pow(n,P-2);i<n;i++) a[i]=1LL*a[i]*inv%P;
}

不过有些丧心病狂的题的模数不是NTT模数,这个时候需要CRT合并。


多项式求逆

比如BZOJ4555

最后会得到f(x)=f(x)g(x)+1,这是一个非齐次线性递推,因为求的第n项不大,只要直接求出f(x)就行了

很容易可以得到f(x)=11g(x),如果令h(x)=1g(x),那么f(x)=(h(x))1,所以就要求h(x)的逆元
求逆元的话网上讲解也有很多,鏼爷集训队论文里也有,大概是用了倍增的做法,摘一下
A(x)的逆元,即是求B(x),使A(x)B(x)1(modx2k)

如果我们已经求出G(x)A(x)G(x)1(modxn)
那么

A(x)G(x)1(modxn)

A(x)G(x)10(modxn)

(A(x)G(x)1)20(modx2n)

A2(x)G2(x)2A(x)G(x)+10(modx2n)

A(x)(2G(x)A(x)G2(x))1(modx2n)

n=0时,求一下A(x)常数项的乘法逆元就可以了

那么用O(nlogn)的复杂度把B(x)赋成2G(x)A(x)G2(x)进入下一层运算就可以了
总复杂度T(n)=T(n/2)+O(nlogn)=O(nlogn)

void Inv(int *a,int *b,int n){
    if(n==1){
        b[0]=Pow(a[0],P-2); return ;
    }
    Inv(a,b,n>>1);
    for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
    FFT(tmp,n<<1,1); FFT(b,n<<1,1);
    for(int i=0;i<(n<<1);i++)
        tmp[i]=(1LL*b[i]*2+P-1LL*tmp[i]*b[i]%P*b[i]%P)%P;
    FFT(tmp,n<<1,0);
    for(int i=0;i<n;i++) b[i]=tmp[i],b[i+n]=0;
}

应用的话,就是一些解一些现行递推式,或者求伯努利数
预处理伯努利数:http://blog.csdn.net/coldef/article/details/72876762

学习自:http://blog.miskcoo.com/2015/05/polynomial-inverse


多项式开根

就是对A(x),求B(x),使B2(x)A(x)(modxn)
同样是用倍增的方法。
已求得G(x),使G2(x)A(x)(modxn)
那么

B2(x)G2(x)(modxn)

[B2(x)G2(x)]20(modx2n)

B4(x)2B2(x)G2(x)+G4(x)0(modx2n)

B4(x)+2B2(x)G2(x)+G4(x)4B2(x)G2(x)(modx2n)

[B2(x)+G2(x)]2[2B(x)G(x)]2(modx2n)

B2(x)+G2(x)2B(x)G(x)(modx2n)

A(x)+G2(x)2B(x)G(x)(mod22n)

B(x)A(x)+G2(x)2G(x)(mod22n)

n=0的时候就是B(x)的常数项就是1(我也不知道为什么,Manchery说的)
时间复杂度T(n)=T(n/2)++O(nlogn)=O(nlogn)

void Sqrt(int *a,int *b,int n){
    if(n==1) return void(b[0]=1);
    Sqrt(a,b,n>>1);
    memset(invb,0,sizeof(int)*n); Inv(b,invb,n);
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
    for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
    FFT(invb,n<<1,1); FFT(tmp,n<<1,1); 
    for(int i=0,inv2=P+1>>1;i<(n<<1);i++) tmp[i]=1LL*tmp[i]*inv2%P*invb[i]%P;
    FFT(tmp,n<<1,0);
    for(int i=0,inv2=P+1>>1;i<n;i++) b[i]=(1LL*b[i]*inv2+tmp[i])%P;
}

学习自: http://www.cnblogs.com/acha/p/6472798.html


多项式取模

此取模不是指系数取模,而是多项式模多项式


UPD:
暴力的多项式取模比较简单,
设模多项式为 M(x)=i=0naixi,化一下可以得到 xn=i=0n1aianxi,将所有大于 n 的项带入就行


对于A(x)B(x),设degA(x)=n,degB(x)=m(mn)
D(x)R(x),满足A(x)=D(x)B(x)+R(x)
只要求出D(x),通过多项式加减乘就可以得到R(x)
1x带入多项式

A(1x)=D(1x)B(1x)+R(1x)

两边乘xn

xnA(1x)=xnmD(1x)xmB(1x)+xnR(1x)

fR(x)=xdegf(x)f(1x)

AR(x)=DR(x)BR(x)+xndegR(x)RR(x)

发现degR(x)<m

那么

AR(x)DR(x)BR(x)(modxnm+1)

AR(x)BR(x)DR(x)(modxnm+1)

这样多项式求个逆就可以求出D(x)
带回去加加减减就可以得到R(x)

inline void Div(int *a,int n,int *b,int m){
  static int tmp[N],A[N],B[N];
  for(int i=0;i<m;i++) t[i]=b[m-i-1];
  for(int i=0;i<n;i++) A[i]=a[n-i-1];
  int nn=1,d=n-m+1; for(;nn<d<<1;nn<<=1);
  for(int i=n;i<nn;i++) A[i]=0;
  for(int i=d;i<nn;i++) t[i]=0;
  for(int i=0;i<nn;i++) B[i]=0;
  Inv(t,B,nn);
  for(int i=d;i<nn;i++) B[i]=0;
  Mul(A,B,max(n,d));
  for(int i=d;i<=n<<1;i++) A[i]=0;
  for(int i=0;i<d;i++) if(i>d-i-1) swap(A[i],A[d-i-1]);
  for(int i=0;i<m;i++) t[i]=b[i];
  Mul(t,A,max(d,m));
  for(int i=0;i<n;i++) a[i]=(a[i]+P-t[i])%P;
}

学习自:http://blog.miskcoo.com/2015/05/polynomial-division


大致就是这些算法

还有多项式求exp什么的…还不会…先挖个坑

多项式exp

可以参考金策的2015年国家集训队论文

B(x)=eA(x)

两边同求导

B(x)=eA(x)A(x)

B(x)=B(x)A(x)

对于每一项的系数

(i+1)bi+1=j+k=ibjak+1(k+1)

也就是说

bi=1ij=1ijajbij

分治FFT

牛顿迭代不会~~

void exp(int *a,int *b,int l,int r){
  if(l==r){
    add(b[l],1LL*inv[l]*a[l]%P); return ;
  }
  int mid=l+r>>1;
  exp(a,b,l,mid); 
  static int tmpa[N],tmpb[N];
  for(int i=l;i<=mid;i++) tmpa[i-l]=b[i];
  for(int i=1;i<=r-l;i++) tmpb[i]=a[i];
  int m,L=0;
  for(m=1;m<r-l+1;m<<=1,L++); m<<=1;
  for(int i=1;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
  NTT(tmpa,m,1); NTT(tmpb,m,1);
  for(int i=0;i<m;i++) tmpa[i]=1LL*tmpa[i]*tmpb[i]%P;
  NTT(tmpa,m,0);
  for(int i=mid+1;i<=r;i++) add(b[i],1LL*inv[i]*tmpa[i-l]%P);
  for(int i=0;i<m;i++) tmpa[i]=tmpb[i]=0;
  exp(a,b,mid+1,r);
}

应用

根据生成函数的卷积运算优化转移什么的都很常见
还有就是多项式的多点求值,这玩意看http://blog.miskcoo.com/2015/05/polynomial-multipoint-eval-and-interpolation
我实现了一下…然而自带大常数,跑的不如暴力…

NOI2017D1T3,就运用了多项式取模来解决线性递推数列的第n项(n很大),也就是f(x)=f(x)g(x)+1的第n项,

具体看叉姐论文和Picks的博客
http://www.docin.com/p-724323397.html
http://picks.logdown.com/posts/197262-polynomial-division
codechef rng是模板题,可以在上面看大神代码

其他的用处…惯例先挖个坑


本人因能力限制,代码仅供参考,如有疑问,欢迎与本人讨论

阅读更多 登录后自动展开
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页