[学习笔记]多项式多点求值与多点插值 - 多项式理论 - 学习笔记

下文及代码中所有提及某个函数是以n为界或者是n次的，意思是其最高次项次数$n-1$。
关于求逆，由牛顿迭代：$$\begin{gathered} F=G^{-1} \\ {F}_{n+1}=2F_n-F_n^{2}G \end{gathered}$$
关于多项式取模，A(x)以n为界，B(x)以m为界，需要求C(x)和D(x)，使得C(x)的界是n-m+1，D(x)的界是m-1：
$$\begin{gathered} A(x)=B(x)C(x)+D(x) \\ {x}^{n-1}A\left(\frac{1}{x}\right)=\left(x^{m-1}B\left(\frac{1}{x}\right)\right)\left(x^{n-m}C\left(\frac{1}{x}\right)\right)+x^{n-1}D\left(\frac{1}{x}\right) \end{gathered}$$
注意到$x^{n-1}D\left(\frac{1}{x}\right)\mathrm{mod}{x}^{n-m+1}=0$，因此（设$A^{T_n}(x)$表示将$A(x)$以$n$为界翻转）：
$$\begin{gathered} A^{T_n}(x)=B^{T_m}(x)C^{T_{n-m+1}}(x)(\mathrm{mod}{x}^{n-m+1}) \\ {C}^{T_{n-m+1}}(x)=A^{T_n}(x){\left(B^{T_m}(x)\right)}^{-1} \end{gathered}$$
因此可以求出$C$，再通过一次乘法求出D即可。

多项式多点求值：
现在有一个以$n$为界的多项式$A(x)$，以及$m$个位置$x_1,\dots ,x_m$，求$A(x_1),\dots ,A(x_m)$。
考虑分治，令$L(x)={\prod }_{i=L}^{mid}(x-x_i),R(x)={\prod }_{i=mid+1}^R(x-x_i)$，那么$A(x_i)=\left(A\mathrm{mod}L\right)(x_i),\forall i\in [L,mid]$，右边同理，递归处理即可。$L,R$可以分治NTT的时候保存下来。

多项式多点插值：
现在告诉你$n$个位置的点值$(x_1,y_1),\dots ,(x_n,y_n)$，求一个以$n$为界的多项式$A(x)$，使得$\forall i\in [1,n],A(x_i)=y_i$。
考虑拉格朗日插值：$$A(x)=\sum _{i=1}^n{y}_i\prod _{j=1,|j-i|>0}^n\frac{x-x_j}{x_i-x_j}$$
首先考虑$v_i={\prod }_{j=1,|j-i|>0}^n(x_i-x_j)$怎么求。
注意到令$F(x)={\prod }_{i=1}^n(x-x_i)$，那么$F^{\prime }(x)={\sum }_{i=1}^n{\prod }_{j=1,|j-i|>0}^n(x-x_j)$，因此$v_i=F^{\prime }(x_i)$，分治NTT出$F$求导再多点求值即可。
接下来考虑求$\mathrm{solve}(L,R)={\sum }_{i=L}^R\frac{y_i}{v_i}{\prod }_{j=L,|j-i|>0}^R(x-x_j)$。
令$L(x)={\prod }_{i=L}^{mid}(x-x_i),R(x)={\prod }_{i=mid+1}^R(x-x_i)$
那么$\mathrm{solve}(L,R)=\mathrm{solve}(L,mid)R(x)+L(x)\mathrm{solve}(mid+1,R)$。
几个卡常小技巧是：预处理单位根，多项式求逆的时候手动展开。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define lint long long
#define p 998244353
#define ull unsigned lint
#define db long double
#define pb push_back
#define mp make_pair
#define fir first
#define sec second
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
typedef pair<int,int> pii;
typedef set<int>::iterator sit;
namespace INPUT_SPACE{
	const int BS=(1<<24)+5;char S[BS],*H,*T;
	char gc() { if(H==T) T=(H=S)+fread(S,1,BS,stdin);return (H==T)?EOF:*H++; }
	inline int inn() {
		int x,ch;while((ch=gc())<'0'||ch>'9');x=ch^'0';
		while((ch=gc())>='0'&&ch<='9') x=(x<<1)+(x<<3)+(ch^'0');return x;
	}
}using INPUT_SPACE::inn;
namespace OUTPUT_SPACE{
	char s[100000*15],t[20];int n,m; inline int PC(char c) { return s[++n]=c; }
	inline void print(int x) { if(!x) s[++n]='0';for(m=0;x;x/=10) t[++m]=char(x%10+'0');for(;m;m--) s[++n]=t[m]; }
	inline void Flush() { fwrite(s+1,sizeof(char),n,stdout),n=0; }
}using OUTPUT_SPACE::print;using OUTPUT_SPACE::PC;using OUTPUT_SPACE::Flush;
inline int pr(int *a,int n) { cerr<<"n = "<<n<<", ";rep(i,0,n-1) cerr<<a[i]sp;cerr ln;return 0; }
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%p) (k&1)?ans=(lint)ans*x%p:0;return ans; }
namespace NTT_space{
	const int N=270000;

	int *dwg[N],*dwgi[N];
	struct PRELUDE_DWG {
		PRELUDE_DWG() {
			int n=N-2;
			for(int i=2,t=1;i<=n;i<<=1,t++)
			{
				dwg[t]=new int[i>>1],dwgi[t]=new int[i>>1];
				int *d=dwg[t],*di=dwgi[t];d[0]=di[0]=1;
				int w=fast_pow(3,(p-1)/i),wi=fast_pow(3,p-1-(p-1)/i);
				rep(j,1,(i>>1)-1) d[j]=(lint)d[j-1]*w%p,di[j]=(lint)di[j-1]*wi%p;
			}
		}
	} __PRELUDE_DWG__;

	namespace NTT_SPACE{
		int r[N];
		inline int pre(int m) {
			int n=1,L=0;while(n<m) n<<=1,L++;
			rep(i,1,n-1) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
			return n;
		}
		// a[n] = 0
		inline int NTT(int *a,int n,int s) {
			rep(i,1,n-1) if(i<r[i]) swap(a[i],a[r[i]]);
			for(int i=2,c=1;i<=n;i<<=1,c++) {
				int *d=(s>0?dwg[c]:dwgi[c]);
				for(int j=0,t=i>>1,x,y;j<n;j+=i) rep(k,0,t-1)
					x=a[j+k],y=(lint)d[k]*a[j+k+t]%p,
					a[j+k]=(x+y>=p?x+y-p:x+y),
					a[j+k+t]=(x-y<0?x-y+p:x-y);
			}
			if(s<0) for(int i=0,v=fast_pow(n,p-2);i<n;i++) a[i]=(lint)a[i]*v%p;
			return 0;
		}
	}using NTT_SPACE::NTT;

	inline void clr(int *a,int n) { memset(a,0,sizeof(int)*n); }
	inline void cpy(int *a,int *b,int n) { memcpy(a,b,sizeof(int)*n); }
	inline void rev(int *a,int n) { reverse(a,a+n); }
	inline int *newInt(int n) { int *t=new int[n];memset(t,0,sizeof(int)*n);return t; }

	namespace TMS_space{
		int A[N],B[N];
		inline int trans(int *a,int *b,int m,int n,int s) { return cpy(a,b,m),clr(a+m,n-m),NTT(a,n,s); }
		// a[m1] = b[m2] = 0
		inline void tms(int *a,int m1,int *b,int m2,int *c,int m3=-1) {
			if(m3<0) m3=m1+m2-1;
			int n=NTT_SPACE::pre(m1+m2-1);
			trans(A,a,m1,n,1),trans(B,b,m2,n,1);
			rep(i,0,n-1) A[i]=(lint)A[i]*B[i]%p;
			NTT(A,n,-1),cpy(c,A,m3);
		}
	}using TMS_space::trans;using TMS_space::tms;

	namespace INV_space{
		int f[N],g[N],F[N],G[N];
		// a[m1] = 0
		inline void polyinv(int *a,int m1,int *b,int m2=-1) {
			if(m2<0) m2=m1;int n=1;while(n<m2) n<<=1;m1=min(m1,m2);
			cpy(g,a,m1),clr(g+m1,n-m1),f[0]=fast_pow(g[0],p-2);
			for(int i=2;i<=n;i<<=1) {
				int l=NTT_SPACE::pre(i<<1);trans(F,f,i>>1,l,1),trans(G,g,i,l,1);
				rep(j,0,l-1) F[j]=(F[j]+F[j]-(lint)F[j]*F[j]%p*G[j])%p,(F[j]<0?F[j]+=p:0);
				NTT(F,l,-1),cpy(f,F,i);
			}
			cpy(b,f,m2);
		}
	}using INV_space::polyinv;

	inline void polydif(int *a,int n,int *b) { rep(i,0,n-2) b[i]=(i+1ll)*a[i+1]%p; }
	inline void polyint(int *a,int n,int *b) { rep(i,0,n-2) b[i+1]=(lint)a[i]*fast_pow(i+1,p-2)%p;b[0]=0; }
	namespace POLYLN_space{
		int b[N],c[N];
		inline void polyln(int *a,int n,int *res) { polyinv(a,n,b),polydif(a,n,c),tms(b,n,c,n,b),polyint(b,n,c),cpy(res,c,n); }
	}using POLYLN_space::polyln;

	namespace POLYEXP_space{
		int b[N],c[N],d[N];
		// a[m] = 0
		inline void polyexp(int *a,int m,int *res) {
			int n=1;while(n<(m<<1)) n<<=1;b[0]=1;
			clr(d,n),cpy(d,a,m);
			for(int i=2;i<=n;i<<=1) {
				polyln(b,i>>1,c);
				rep(j,0,(i>>1)-1) c[j]=d[j]-c[j],(c[j]<0?c[j]+=p:0);
				c[0]++,(c[0]>=p?c[0]-=p:0),tms(b,i>>1,c,i>>1,b);
			}
			cpy(res,b,m);
		}
	}using POLYEXP_space::polyexp;

	namespace DIV_space{
		int A[N],B[N],C[N],D[N];
		// A / B = C , ... D, a[n] = b[m] = 0
		inline void polydiv(int *a,int n,int *b,int m,int *c,int *d) {
			if(n<m) { clr(d,m-1),cpy(d,a,n);return; }
			cpy(A,a,n),rev(A,n),cpy(B,b,m),rev(B,m),polyinv(B,m,B,n-m+1);
			tms(A,n-m+1,B,n-m+1,C,n-m+1),rev(C,n-m+1),tms(C,n-m+1,b,m,D);
			rep(i,0,m-2) d[i]=a[i]-D[i],(d[i]<0?d[i]+=p:0);cpy(c,C,n-m+1);
		}
		inline void polymul(int *a,int n,int *b,int m,int *d) { polydiv(a,n,b,m,C,d); }
	}using DIV_space::polydiv;using DIV_space::polymul;

	int *pd[N],*tmp[N],len[N],res[N];
	namespace VAL_space{
		inline void build(int *x,int k,int L,int R) {
			if(L==R) { pd[k]=newInt(len[k]=2),pd[k][0]=(x[L]?p-x[L]:0),pd[k][1]=1;return; }
			int mid=(L+R)>>1,lc=k<<1,rc=k<<1|1;build(x,lc,L,mid),build(x,rc,mid+1,R);
			pd[k]=newInt(len[k]=len[lc]+len[rc]-1),tms(pd[lc],len[lc],pd[rc],len[rc],pd[k]);
		}
		inline void solve(int *a,int n,int k,int L,int R,int *x,int *y) {
			if(L==R) { y[L]=0;rep(i,0,n-1) y[L]=((lint)y[L]*x[L]+a[n-i-1])%p;return; }
			int mid=(L+R)>>1,lc=k<<1,rc=k<<1|1;
			tmp[lc]=newInt(len[lc]-1),polymul(a,n,pd[lc],len[lc],tmp[lc]),solve(tmp[lc],len[lc]-1,lc,L,mid,x,y);
			tmp[rc]=newInt(len[rc]-1),polymul(a,n,pd[rc],len[rc],tmp[rc]),solve(tmp[rc],len[rc]-1,rc,mid+1,R,x,y);
		}
		inline void del(int k,int L,int R) {
			int mid=(L+R)>>1,lc=k<<1,rc=k<<1|1;
			delete[] pd[k];if(k>1) delete[] tmp[k];
			if(L<R) del(lc,L,mid),del(rc,mid+1,R);
		}
		// x from 1 to m, a[n] = 0
		inline void polyval(int *a,int n,int *x,int m,int *y) {
			build(x,1,1,m),solve(a,n,1,1,m,x,res),del(1,1,m),cpy(y+1,res+1,m);
		}
	}using VAL_space::polyval;

	namespace COEF_space{
		int t[N],v[N],tpl[N],tpr[N],*pl[N];
		inline int solve(int k,int L,int R) {
			if(L==R) return pl[k]=newInt(1),pl[k][0]=v[L],1; pl[k]=newInt(R-L+1);
			int mid=(L+R)>>1,lc=k<<1,rc=k<<1|1,ll=solve(lc,L,mid),rl=solve(rc,mid+1,R);
			tms(pl[lc],ll,pd[rc],len[rc],tpl),tms(pl[rc],rl,pd[lc],len[lc],tpr);
			rep(i,0,R-L) pl[k][i]=tpl[i]+tpr[i],(pl[k][i]>=p?pl[k][i]-=p:0);return R-L+1;
		}
		inline void del(int k,int L,int R) {
			delete[] pl[k];if(L==R) return;
			int mid=(L+R)>>1;del(k<<1,L,mid),del(k<<1|1,mid+1,R);
		}
		// a[n] = 0, (x, y) from 1 to n
		inline void polycoef(int *x,int *y,int n,int *a) {
			VAL_space::build(x,1,1,n),polydif(pd[1],len[1],t);
			VAL_space::solve(t,len[1]-1,1,1,n,x,v);
			rep(i,1,n) v[i]=(lint)fast_pow(v[i],p-2)*y[i]%p;
			solve(1,1,n),cpy(a,pl[1],n),VAL_space::del(1,1,n),del(1,1,n);
		}
	}using COEF_space::polycoef;
}
using NTT_space::NTT;
using NTT_space::tms;
using NTT_space::polyinv;
using NTT_space::polyln;
using NTT_space::polyexp;
using NTT_space::polydiv;
using NTT_space::polyval;
using NTT_space::polycoef;
const int N=100010;int x[N],y[N],a[N];
int main() {
	int n=inn();
	rep(i,1,n) x[i]=inn(),y[i]=inn();
	polycoef(x,y,n,a);
	rep(i,0,n-1) print(a[i]),PC(' ');
	return PC('\n'),Flush(),0;
}