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

下文及代码中所有提及某个函数是以n为界或者是n次的,意思是其最高次项次数 n − 1 n-1 n1
关于求逆,由牛顿迭代: F = G − 1 F n + 1 = 2 F n − F n 2 G F=G^{-1}\\F_{n+1}=2F_n-F_n^2G F=G1Fn+1=2FnFn2G
关于多项式取模,A(x)以n为界,B(x)以m为界,需要求C(x)和D(x),使得C(x)的界是n-m+1,D(x)的界是m-1:
A ( x ) = B ( x ) C ( x ) + D ( x ) x n − 1 A ( 1 x ) = ( x m − 1 B ( 1 x ) ) ( x n − m C ( 1 x ) ) + x n − 1 D ( 1 x ) A(x)=B(x)C(x)+D(x)\\x^{n-1}A\left(\frac 1x\right)=\left(x^{m-1}B\left(\frac 1x\right)\right)\left(x^{n-m}C\left(\frac1x\right)\right)+x^{n-1}D\left(\frac 1x\right) A(x)=B(x)C(x)+D(x)xn1A(x1)=(xm1B(x1))(xnmC(x1))+xn1D(x1)
注意到 x n − 1 D ( 1 x )   m o d   x n − m + 1 = 0 x^{n-1}D\left(\frac 1x\right)\bmod {x^{n-m+1}}=0 xn1D(x1)modxnm+1=0,因此(设 A T n ( x ) A^{T_n}(x) ATn(x)表示将 A ( x ) A(x) A(x) n n n为界翻转):
A T n ( x ) = B T m ( x ) C T n − m + 1 ( x ) ( m o d x n − m + 1 ) C T n − m + 1 ( x ) = A T n ( x ) ( B T m ( x ) ) − 1 A^{T_{n}}(x)=B^{T_{m}}(x)C^{T_{n-m+1}}(x)\pmod {x^{n-m+1}}\\ C^{T_{n-m+1}}(x)=A^{T_{n}}(x)\left(B^{T_{m}}(x)\right)^{-1} ATn(x)=BTm(x)CTnm+1(x)(modxnm+1)CTnm+1(x)=ATn(x)(BTm(x))1
因此可以求出 C C C,再通过一次乘法求出D即可。

多项式多点求值:
现在有一个以 n n n为界的多项式 A ( x ) A(x) A(x),以及 m m m个位置 x 1 , … , x m x_1,\dots,x_m x1,,xm,求 A ( x 1 ) , … , A ( x m ) A(x_1),\dots,A(x_m) A(x1),,A(xm)
考虑分治,令 L ( x ) = ∏ i = L m i d ( x − x i ) , R ( x ) = ∏ i = m i d + 1 R ( x − x i ) L(x)=\prod_{i=L}^{mid}(x-x_i),R(x)=\prod_{i=mid+1}^R(x-x_i) L(x)=i=Lmid(xxi),R(x)=i=mid+1R(xxi),那么 A ( x i ) = ( A   m o d   L ) ( x i ) , ∀ i ∈ [ L , m i d ] A(x_i)=\left(A\bmod L\right)(x_i),\forall i\in[L,mid] A(xi)=(AmodL)(xi),i[L,mid],右边同理,递归处理即可。 L , R L,R L,R可以分治NTT的时候保存下来。

多项式多点插值:
现在告诉你 n n n个位置的点值 ( x 1 , y 1 ) , … , ( x n , y n ) (x_1,y_1),\dots,(x_n,y_n) (x1,y1),,(xn,yn),求一个以 n n n为界的多项式 A ( x ) A(x) A(x),使得 ∀ i ∈ [ 1 , n ] , A ( x i ) = y i \forall i\in[1,n],A(x_i)=y_i i[1,n],A(xi)=yi
考虑拉格朗日插值: A ( x ) = ∑ i = 1 n y i ∏ j = 1 , ∣ j − i ∣ > 0 n x − x j x i − x j A(x)=\sum_{i=1}^ny_i\prod_{j=1,|j-i|>0}^n\frac{x-x_j}{x_i-x_j} A(x)=i=1nyij=1,ji>0nxixjxxj
首先考虑 v i = ∏ j = 1 , ∣ j − i ∣ > 0 n ( x i − x j ) v_i=\prod_{j=1,|j-i|>0}^n (x_i-x_j) vi=j=1,ji>0n(xixj)怎么求。
注意到令 F ( x ) = ∏ i = 1 n ( x − x i ) F(x)=\prod_{i=1}^n(x-x_i) F(x)=i=1n(xxi),那么 F ′ ( x ) = ∑ i = 1 n ∏ j = 1 , ∣ j − i ∣ > 0 n ( x − x j ) F'(x)=\sum_{i=1}^n\prod_{j=1,|j-i|>0}^n(x-x_j) F(x)=i=1nj=1,ji>0n(xxj),因此 v i = F ′ ( x i ) v_i=F'(x_i) vi=F(xi),分治NTT出 F F F求导再多点求值即可。
接下来考虑求 s o l v e ( L , R ) = ∑ i = L R y i v i ∏ j = L , ∣ j − i ∣ > 0 R ( x − x j ) \mathrm{solve}(L,R)=\sum_{i=L}^R\frac {y_i}{v_i}\prod_{j=L,|j-i|>0}^R(x-x_j) solve(L,R)=i=LRviyij=L,ji>0R(xxj)
L ( x ) = ∏ i = L m i d ( x − x i ) , R ( x ) = ∏ i = m i d + 1 R ( x − x i ) L(x)=\prod_{i=L}^{mid}(x-x_i),R(x)=\prod_{i=mid+1}^R(x-x_i) L(x)=i=Lmid(xxi),R(x)=i=mid+1R(xxi)
那么 s o l v e ( L , R ) = s o l v e ( L , m i d ) R ( x ) + L ( x ) s o l v e ( m i d + 1 , R ) \mathrm{solve}(L,R)=\mathrm{solve}(L,mid)R(x)+L(x)\mathrm{solve}(mid+1,R) solve(L,R)=solve(L,mid)R(x)+L(x)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;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值