luogu P4726 【模板】多项式指数函数 多项式 exp 牛顿迭代 泰勒展开
LINK:多项式 exp
做多项式的题 简直在嗑药。
前置只是 泰勒展开
这个东西用于 对于一个函数f(x) 我们不好得到 其在x处的取值。
所以另外设一个函数g(x) 来在x点处无限逼近f(x).
具体的 \(f(x) ≈ g(x)=g(0)+\frac{f^1(0)}{1!}x+\frac{f^2(0)}{2!}x^2+...+\frac{f^n(0)}{n!}x^n\)
牛顿迭代:
常用来求一个函数的零点:假设我们已经求得一个近似值x0 那么我们只需要过(x0,f(x0))这个点做函数图像的切线 取切线与x轴的交点作为新的x0.
迭代几次就可以比较精确。
如:现在要求一个函数f(x) 近似值为x0 y=f'(x0)(x-x0)+f(x0);
y=0时 \(x=x0-\frac{f(x0)}{f'(x0)}\)
当然也可以放到多项式上 现在要求一个G(x) 我们想要求出F(G(x))=0的零点G(x).
\(G(x)=G0(x)-\frac{F(G0(x))}{F'(G0(x))}\)
本质上每迭代一次都可以迅速逼近真实值。
如果\(F(G0(x))≡0(\bmod x^{2n})\)那么\(F(G(x))\equiv 0(\bmod{x^{n}})\)
关于这道题的推导过程:
\(B(x)\equiv e^{A(x)}(\bmod x^n)\)
\(InB(x)-A(x)\equiv 0(\bmod x^n)\)
现在设一个函数\(F(G(x))=InG(x)-A(x)\equiv 0\)
其实在求F这个函数的零点.
两边直接求导可得 \((F(G0(x)))'=\frac{G'0(x)}{G0(x)}\)
带入牛顿迭代的式子里。
\(G(x)=\frac{G0(x)(1-InG0(x)+A(x))}{G'0(x)}\)
每次迭代需要 求逆 做多项式In 再来一遍多项式乘法即可。
废什么话 码!
求导 和 积分需要仔细熟悉一下.
const int MAXN=600010,G=3;
int n;
int A[MAXN],B[MAXN],E[MAXN],F[MAXN],C[MAXN],D[MAXN],g[MAXN],rev[MAXN],inv[MAXN],O[MAXN];
inline int mul(int a,int b){return (ll)a*b%mod;}
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline int mus(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int ksm(int b,int p)
{
int cnt=1;
while(p)
{
if(p&1)cnt=mul(cnt,b);
p=p>>1;b=mul(b,b);
}
return cnt;
}
inline void NTT(int *a,int op,int ww)
{
int lim=1;while(lim<ww)lim=lim<<1;
rep(0,lim-1,i)
{
rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
if(i<rev[i])swap(a[i],a[rev[i]]);
}
for(int len=2;len<=lim;len=len<<1)
{
int mid=len>>1;
int wn=ksm(G,op==1?(mod-1)/len:mod-1-(mod-1)/len);
O[0]=1;rep(1,mid-1,i)O[i]=(ll)O[i-1]*wn%mod;
for(int j=0;j<lim;j+=len)
{
for(int i=0;i<mid;++i)
{
int x=a[i+j],y=(ll)a[i+j+mid]*O[i]%mod;
a[i+j]=(x+y)%mod;a[i+j+mid]=(x-y+mod)%mod;
}
}
}
if(op==-1)for(int i=0,inv=ksm(lim,mod-2);i<lim;++i)a[i]=(ll)a[i]*inv%mod;
}
inline void Direv(int *a,int *b,int len)//求导
{
rep(0,len-2,i)b[i]=mul(a[i+1],i+1);b[len-1]=0;
}
inline void Inv(int *a,int *b,int len)
{
if(len==1)return b[0]=ksm(a[0],mod-2),void();
Inv(a,b,len>>1);rep(0,len-1,i)C[i]=a[i],D[i]=b[i];
NTT(C,1,len<<1);NTT(D,1,len<<1);
rep(0,(len<<1)-1,i)D[i]=mul(D[i],mul(D[i],C[i]));
NTT(D,-1,len<<1);
rep(0,len-1,i)b[i]=((ll)2*b[i]-D[i]+mod)%mod;
rep(0,(len<<1)-1,i)C[i]=D[i]=0;
}
inline void Inter(int *a,int *b,int len)
{
rep(1,len-1,i)b[i]=mul(a[i-1],inv[i]);b[0]=0;
}
inline void Ln(int *a,int *b,int len)
{
Inv(a,E,len);Direv(a,F,len);
NTT(E,1,len<<1);NTT(F,1,len<<1);
rep(0,(len<<1)-1,i)E[i]=mul(E[i],F[i]);
NTT(E,-1,len<<1);Inter(E,b,len);
rep(0,(len<<1)-1,i)E[i]=F[i]=0;
}
inline void Exp(int *a,int *b,int len)
{
if(len==1)return b[0]=1,void();
Exp(a,b,len>>1);Ln(b,g,len);
g[0]=(a[0]+1-g[0]+mod)%mod;
rep(1,len-1,i)g[i]=mus(a[i],g[i]);
NTT(g,1,len<<1);NTT(b,1,len<<1);
rep(0,(len<<1)-1,i)b[i]=mul(b[i],g[i]);
NTT(b,-1,len<<1);rep(len,(len<<1)-1,i)g[i]=b[i]=0;
}
int main()
{
//freopen("1.in","r",stdin);
get(n);rep(0,n-1,i)get(A[i]);
int len=1;while(len<n)len=len<<1;
inv[1]=1;rep(2,(len<<1),i)inv[i]=(ll)inv[mod%i]*(mod-mod/i)%mod;
Exp(A,B,len);rep(0,n-1,i)put_(B[i]);return 0;
}