众所周知,生成函数是一个十分强大的东西,许多与多项式相关的算法也就应运而生了,在这里选取几种较为简单的算法做一个介绍.
p.s. 这篇文章在去年NOI前已经完成了一半,现在笔者将其补充完整后发出,同时也为了纪念那一段美好的时光。
多项式乘法(FFT/NTT)
https://www.cnblogs.com/encodetalker/p/10212036.html
https://www.cnblogs.com/encodetalker/p/10285657.html
多项式求逆
已知(F(x)),求(G(x))使得(F(x)G(x)equiv 1(mod x^n))
倍增,假设已经求出(A(x)F(x)equiv1(mod x^{lceil frac{n}{2} ceil}))
(G(x)-A(x)equiv0(mod x^{lceil frac{n}{2} ceil}))
注意到这个式子可以平方一下(G(x)^2-2A(x)G(x)+A(x)^2equiv0(mod x^n))
两边乘上(F(x),G(x)-2A(x)+F(x)A(x)^2equiv0(mod x^n))
最后得到(G(x)equiv A(x)(2-F(x)A(x))(mod x^n))
多项式取模
咕咕咕
多项式牛顿迭代
假设已知(n)次多项式(F(x)),求多项式(G(x))使得(F(G(x))equiv 0 (mod x^n)).
同样考虑倍增,假设(F(G_1(x))equiv 0(mod x^{lceil frac{n}{2} ceil})).
考虑将(F(G(x)))在(G_1(x))处Taylor展开,有:
注意到在(mod n)的意义下,若(kgeq 2), 那么((G(x)-G_1(x))^kequiv 0(mod x^n)).
于是有
移项之后得到:
多项式开方
已知(n)次多项式(F(x)),求多项式(G(x))满足(G(x)^2equiv F(x)(mod x^n)).
考虑牛顿迭代,构造(A(G(x))=G(x)^2-F(x).)显然最后需要(A(G(x))equiv 0(mod x^n)).
带入牛顿迭代的结论中可以得到:
多项式求导&&积分
对(F(x)=sum_{i=0}^n a_ix^i)
求导:(F'(x)=sum_{i=1}^na_iix^{i-1})
积分:(int F(x)=sum_{i=0}^n(frac{a_i}{i+1}x^{i+1})+C)(一个常数)
多项式Ln
复合函数求导:记(H(x)=F(G(x))),则(H'(x)=F'(G(x))G'(x))
已知多项式(F(x)),求(G(x)=ln(F(x)))
对两边求导:(G'(x)=frac{F'(x)}{F(x)})
之后再积分回去:(G(x)=int frac{F'(x)}{F(x)}),多项式求逆之后再积分即可。
多项式Exp
已知(n)次多项式(F(x)),求(G(x))满足(G(x)equiv e^{F(x)}(mod x^n)).
两边同时取(ln)得(ln G(x)equiv F(x)(mod x^n)).
令(A(x)=ln G(x)-F(x)), 继续套用牛顿迭代得到:
总代码如下:
namespace Polynomial{
int A[N<<2],r[N<<2],B[N<<2],C[N<<2],D[N<<2];
int calcr(int len)
{
int lim=1,cnt=0;
while (lim<len) {lim<<=1;cnt++;}
rep(i,0,lim-1)
r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
return lim;
}
void ntt(int lim,int *a,int typ)
{
rep(i,0,lim-1)
if (i<r[i]) swap(a[i],a[r[i]]);
for (int mid=1;mid<lim;mid<<=1)
{
int len=(mid<<1),gn=qpow(3,(maxd-1)/len);
if (typ==-1) gn=getinv(gn);
for (int sta=0;sta<lim;sta+=len)
{
int g=1;
for (int j=0;j<mid;j++,g=mul(g,gn))
{
int x=a[sta+j],y=mul(a[sta+j+mid],g);
a[sta+j]=add(x,y);a[sta+j+mid]=dec(x,y);
}
}
}
if (typ==-1)
{
int inv=getinv(lim);
rep(i,0,lim-1) a[i]=mul(a[i],inv);
}
}
void Inv(int *a,int *b,int n)
{
if (n==1)
{
b[0]=getinv(a[0]);
return;
}
Inv(a,b,(n+1)>>1);
int lim=calcr(n<<1);
rep(i,0,n-1) A[i]=a[i];
ntt(lim,A,1);ntt(lim,b,1);
rep(i,0,lim-1) b[i]=mul(b[i],dec(2,mul(A[i],b[i])));
ntt(lim,b,-1);
rep(i,0,lim-1) A[i]=0;
rep(i,n,lim-1) b[i]=0;
}
void Direv(int *a,int *b,int n)
{
rep(i,1,n-1) b[i-1]=mul(a[i],i);
b[n-1]=0;
}
void Inter(int *a,int *b,int n)
{
rep(i,1,n-1) b[i]=mul(a[i-1],getinv(i));
b[0]=0;
}
void Ln(int *a,int *b,int n)
{
rep(i,0,(n<<1)-1) B[i]=C[i]=0;
Direv(a,B,n);Inv(a,C,n);
int lim=calcr(n<<1);
ntt(lim,B,1);ntt(lim,C,1);
rep(i,0,lim-1) B[i]=mul(B[i],C[i]);
ntt(lim,B,-1);Inter(B,b,n);
rep(i,n,lim-1) b[i]=0;
}
void Exp(int *a,int *b,int n)
{
if (n==1) {b[0]=1;return;}
Exp(a,b,(n+1)>>1);
Ln(b,D,n);
rep(i,0,n-1)
{
if (i) D[i]=add(maxd-D[i],a[i]);
else D[i]=add(maxd-D[i],a[i]+1);
}
int lim=calcr(n<<1);
ntt(lim,D,1);ntt(lim,b,1);
rep(i,0,lim-1) b[i]=mul(b[i],D[i]);
ntt(lim,b,-1);
rep(i,n,lim-1) b[i]=0;
}
void Sqrt(int *a,int *b,int n)
{
if (n==1) {b[0]=1;return;}
Sqrt(a,b,(n+1)>>1);
rep(i,0,(n<<1)-1) B[i]=C[i]=0;
Inv(b,B,n);
int lim=calcr(n<<1);
rep(i,0,n-1) C[i]=a[i];
ntt(lim,B,1);ntt(lim,C,1);
rep(i,0,lim-1) B[i]=mul(C[i],B[i]);
ntt(lim,B,-1);
rep(i,0,n-1) b[i]=mul(add(b[i],B[i]),inv2);
}
}
using namespace Polynomial;