题解 P5158【【模板】多项式快速插值】

\(\text{Link}\)

题意

给出 \(n-1\) 次多项式 \(F(x)\) 的 \(n\) 个点值 \(F(x_i)=y_i\)，求出这个多项式 \(F(x)\)。

\(1\le n\le10^5\)。

思路

\(\text{Update 2024.03.19}\)：更改了代码与部分说明。

以下默认 \(mid=\lfloor\frac {l+r} 2\rfloor\)。

回顾拉格朗日插值的式子：

\[F(x)=\sum_{i=1}^n y_i\prod_{i\ne j}\dfrac{x-x_j}{x_i-x_j} \]

将后面的分母提到前面：

\[F(x)=\sum_{i=1}^n \dfrac{y_i}{\prod_{i\ne j}{x_i-x_j}}\prod_{i\ne j}{(x-x_j)} \]

设 \(\delta(x)=\displaystyle\prod_{i=1}^n(x-x_i)\)，则：

\[\prod_{i\ne j}({x_i-x_j})=\lim_{x\to x_i}\dfrac{\delta(x)}{x-x_i}=\delta'(x_i) \]

\[F(x)=\sum_{i=1}^n \dfrac{y_i}{\delta'(x_i)}\prod_{i\ne j}{(x-x_j)} \]

\(\delta(x)\) 容易分治 \(\text{NTT}\) 求出，然后使用多点求值求出 \(\delta'(x_i)\)，不会多点求值可以看这里。

接下来直接分治 \(\text{NTT}\) 即可。具体地，设 \(G_{l,r}(x)=\prod_{i=l}^r(x-x_i)\)，\(H_{l,r}(x)\) 为 \((x_l,y_l),(x_{l+1}y_{l+1}),\cdot\cdot\cdot (x_r,y_r)\) 插出来的多项式，即 \(\sum_{i=l}^r \frac{y_i}{\delta'(x_i)}\prod_{i\ne j}^{l\le j\le r}{(x-x_j)}\)，不难导出：

\[G_{l,r}(x)=G_{l,mid}(x)\cdot G_{mid+1,r}(x) \]

\[H_{l,r}(x)=H_{l,mid}(x)\cdot G_{mid+1,r}(x)+H_{mid+1,r}(x)\cdot G_{l,mid}(x) \]

每一步的时间复杂度均为 \(O(n\log^2 n)\)，于是总时间复杂度为 \(O(n\log^2 n)\)。

核心代码：

namespace Interpolation{
	#define ls (rt<<1)
	#define rs (rt<<1|1)
	inline Poly MulT(const Poly &a,const Poly &b){
		Poly F=a,G=b;
		int n=a.size(),m=b.size();
		reverse(G.begin(),G.end());
		init(n);
		F.resize(lim),G.resize(lim);
		NTT(F,1),NTT(G,1);
		for(int i=0;i<lim;i++)
			G[i]=1ll*F[i]*G[i]%mod;
		NTT(G,-1);
		for(int i=m-1;i<n;i++)
			F[i-m+1]=G[i];
		F.resize(max(0,n-m+1));
		return F;
	}
	Poly TR[N],T[N];
	Poly Q,QY,ans;
	inline void build(int rt,int l,int r){
		if(l==r){
			TR[rt]=(Poly){1,dec(0,Q[l])};
			T[rt]=TR[rt],reverse(T[rt].begin(),T[rt].end());
			return ;
		}
		int mid=l+r>>1;
		build(ls,l,mid),build(rs,mid+1,r);
		TR[rt]=TR[ls]*TR[rs];
		T[rt]=TR[rt],reverse(T[rt].begin(),T[rt].end());
	}
	inline void solve(int rt,int l,int r,Poly F){
		if(l==r){
			ans[l]=F[0];
			return ;
		}
		int mid=l+r>>1;
		solve(ls,l,mid,MulT(F,TR[rs]));
		solve(rs,mid+1,r,MulT(F,TR[ls]));
	}
	inline Poly solve(int rt,int l,int r){
		if(l==r)
			return (Poly){1ll*QY[l]*qpow(ans[l],mod-2)%mod};
		int mid=l+r>>1;
		Poly L=solve(ls,l,mid),R=solve(rs,mid+1,r);
		return L*T[rs]+R*T[ls];
	}
	inline Poly solve(Poly X,Poly Y){
		Q=X,QY=Y;
		int n=X.size();
		build(1,0,n-1);
		Poly F=Deriv(T[1]);
		ans.resize(n),F.resize(n*2);
		solve(1,0,n-1,MulT(F,Inv(TR[1])));
		return solve(1,0,n-1);
	}
}