众所周知,生成函数是一个十分强大的东西,许多与多项式相关的算法也就应运而生了,在这里选取几种较为简单的算法做一个介绍.
p.s. 这篇文章在去年NOI前已经完成了一半,现在笔者将其补充完整后发出,同时也为了纪念那一段美好的时光。
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}\rceil})\)
\(G(x)-A(x)\equiv0(mod\ x^{\lceil \frac{n}{2}\rceil})\)
注意到这个式子可以平方一下\(G(x)^2-2A(x)G(x)+A(x)^2\equiv0(mod\ x^n)\)
两边乘上\(F(x),G(x)-2A(x)+G(x)A(x)^2\equiv0(mod\ x^n)\)
最后得到\(G(x)\equiv A(x)(2-G(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}\rceil})\).
考虑将\(F(G(x))\)在\(G_1(x)\)处Taylor展开,有:
注意到在\(\mod n\)的意义下,若\(k\geq 2\), 那么\((G(x)-G_1(x))^k\equiv 0(\mod x^n)\).
于是有
移项之后得到:
已知\(n\)次多项式\(F(x)\),求多项式\(G(x)\)满足\(G(x)^2\equiv 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\)(一个常数)
复合函数求导:记\(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)}\),多项式求逆之后再积分即可。
已知\(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;
原文:https://www.cnblogs.com/encodetalker/p/12663969.html