多项式exp & 牛顿迭代

牛顿迭代解决的是这样一个问题：已知 \(g(f(x))\equiv 0\pmod {x^n}\) 与 \(g(x)\)，求 模 \(x^n\) 意义下的 \(f(x)\)
这个问题可以用倍增的方式解决。首先假设你知道了 \(g(f(x))=0\) 的常数项（一般都能很方便的知道）。
然后，我们假设 \(f_0(x)\) 是模 \(x^{\lceil\frac{n}{2}\rceil}\) 意义下的答案，这个是已经求出来了的答案。
将 \(f(x)\) 在 \(f_0(x)\) 的位置上进行泰勒展开（附赠通用展开公式），得到 \(\sum_{i=0}^{\infty} \frac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i\equiv 0\pmod {x^n}\)，其中 \(f^{(i)}\) 指的是对函数 \(f(x)\) 进行 \(i\) 次求导。
考虑到 \(f(x)-f_0(x)\) 在 \(\lceil\frac{n}{2}\rceil\) 之前的位置上（不包含）一定是没有值的（因为是取模意义下的，所以 \(f(x)\) 和 \(f_0(x)\) 在这些项上的值一定相同），所以对于所有 \(2\le i\)，\((f(x)-f_0(x))^i\equiv 0\pmod{x^n}\)，然后上面那一大坨求和式就只会保留前两项，也就是 \(g(f_0(x))-g'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod{x^n}\)，然后把这个东西解掉就可以得到 \(f(x)=f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n}\)，然后这个东西就可以在倍增里直接做了！

然后考虑这个东西怎么去实现多项式 \(exp\)。
首先我们看到模板题，假设给定的多项式是 \(h(x)\)，那么两边同时取 \(\ln\) 就有 \(g(f(x))=\ln f_0(x)-h(x)\)，求一次导之后就是 \(\frac{1}{f_0(x)}\)，然后我们把这个东西带入到上面的公式里就有 \(f(x)=f_0(x)\times (1-\ln f_0(x)+h(x))\)，然后递归处理即可。

Code Here

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e6+5,mod=998244353,g=3,gi=332748118;
int n,rev[N<<2];
int qpow(int a,int b){
	int ans=1,base=a;
	while(b){
		if(b&1)ans=ans*base%mod;
		base=base*base%mod;
		b>>=1; 
	}
	return ans;
}
#define Poly vector<int>
Poly operator+(Poly x,Poly y){
	Poly z;
	z.resize(max(x.size(),y.size()));
	x.resize(max(x.size(),y.size()));y.resize(max(x.size(),y.size()));
	for(int i=0;i<z.size();++i)z[i]=(x[i]+y[i])%mod;
	return z;
}
Poly operator-(Poly x,Poly y){
	Poly z;
	z.resize(max(x.size(),y.size()));
	x.resize(max(x.size(),y.size()));y.resize(max(x.size(),y.size()));
	for(int i=0;i<z.size();++i)z[i]=(x[i]-y[i]+mod)%mod;
	return z;
}
void NTT(Poly &A,int n,int opt=1){
	A.resize(n+1);
	for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int mid=1;mid<n;mid<<=1){
		int wn=qpow(opt==1?g:gi,(mod-1)/(mid<<1));
		for(int i=0;i<n;i+=(mid<<1)){
			int w=1;
			for(int j=0;j<mid;++j,w=(w*wn)%mod){
				int tmp0=A[i+j],tmp1=w*A[i+j+mid]%mod;
				A[i+j]=(tmp0+tmp1)%mod;
				A[i+j+mid]=(tmp0-tmp1)%mod;
			}
		}
	}
} 
Poly operator*(Poly x,Poly y){
	rev[0]=0;
	int lim=1,L=0;
	while(lim<=x.size()+y.size()+3)lim<<=1ll,L++;
	for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	x.resize(lim+1);y.resize(lim+1);
	NTT(x,lim,1);NTT(y,lim,1);
	for(int i=0;i<lim;++i)x[i]=x[i]*y[i]%mod;
	NTT(x,lim,-1);
	x.resize(lim+1);
	int inv=qpow(lim,mod-2);
	for(int i=0;i<lim;++i)x[i]=(x[i]*inv%mod+mod)%mod;
	return x;
}
Poly _(Poly x,int w){
	for(int i=0;i<x.size();++i)x[i]=(x[i]*w)%mod;
	return x; 
} 
Poly Inv(Poly x,int n){
	x.resize(n+1);
	if(n==1){
//		cout<<"Check: "<<x[0]<<' '<<qpow(x[0],mod-2)<<'\n';
		Poly o;
		o.push_back(qpow(x[0],mod-2));
		return o;
	}
	Poly tmp=Inv(x,n>>1);
	tmp.resize(n+1); 
//	cout<<"Check2: "<<n<<' '<<tmp.size()<<' '<<tmp[0]<<'\n';
	tmp=_(tmp,2)-x*tmp*tmp;
	tmp.resize(n+1);
	return tmp;
}
Poly inv(Poly x){
	int n=x.size();
	int lim=1;
	while(lim<=n)lim<<=1;
	x.resize(lim+1);
	Poly ans=Inv(x,lim);
	ans.resize(lim+1); 
	return ans;
}
Poly Derive(Poly x){
	Poly ans;
	for(int i=1;i<x.size();++i)ans.push_back(i*x[i]%mod);
	int n=x.size();
	int lim=1;
	while(lim<=n)lim<<=1;
	ans.resize(lim+1);
	return ans;
}
Poly Integral(Poly x){
	Poly ans;
	ans.push_back(0);
	for(int i=0;i<x.size();++i)ans.push_back(x[i]*qpow(i+1,mod-2)%mod);
	int n=x.size();
	int lim=1;
	while(lim<=n)lim<<=1;
	ans.resize(lim+1);
	return ans;
}
Poly ln(Poly x){
	Poly ans=Integral(Derive(x)*inv(x));
	return ans;
}
Poly E(Poly x,int n){
	x.resize(n+1);
	if(n==1){
		Poly now;
		now.push_back(1);
		return now;
	}
	Poly tmp=E(x,n+1>>1);
	tmp.resize(n+1);
	Poly tmp2=x-ln(tmp);
	tmp2[0]=(tmp2[0]+1)%mod;
	tmp2.resize(n+1);
	Poly ans=tmp*tmp2;
	ans.resize(n+1);
	return ans;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	Poly a;
	int n,lim=1;
	cin>>n;
	while(lim<=n)lim<<=1;
	a.resize(lim+1);
	for(int i=0;i<n;++i)cin>>a[i];
	a=E(a,n+1);
	for(int i=0;i<n;++i)cout<<a[i]<<' ';
	return 0;
} 

这个东西基本上也是我目前的多项式全家桶
特别注意一下乘法要把两个多项式 resize 成接近的大小，然后 rev 尽可能开大一点，不然很容易 RE。