差别

这里会显示出您选择的修订版和当前版本之间的差别。

--- 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/24 23:52]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/25 17:43] (当前版本)
jxm2001
@@ 行 411: / 行 411: @@
 ==== 算法简介 ====
-给定 $a_n=f_1a_{n-1}+f_2a_{n-2}+\cdots +f_ka_{n-k}$ 以及 $a_0,a_1\cdots a_{k-1}$，$O(k\log k\log n)$ 时间求 $a_n$。
+给定 $a_n=f_1a_{n-1}+f_2a_{n-2}+\cdots +f_ka_{n-k}$ 以及 $a_0,a_1\cdots a_{k-1}$，询问 $a_n$ 的值。时间复杂度 $O(k\log k\log n)$。
 ==== 算法实现 ====
@@ 行 456: / 行 456: @@
 namespace Poly{
 	const int G=3;
-	int rev[MAXN<<2],Wn[30][2];
+	int rev[MAXN<<2],Pool[MAXN<<3],*Wn[30];
 	void init(){
-		int m=Mod-1,lg2=0;
+		int lg2=0,*pos=Pool,n,w;
-		while(m%2==0)m>>=1,lg2++;
+		while((1<<lg2)<MAXN*2)lg2++;
-		Wn[lg2][1]=quick_pow(G,m);
+		n=1<<lg2,w=quick_pow(G,(Mod-1)/(1<<lg2));
-		Wn[lg2][0]=quick_pow(Wn[lg2][1],Mod-2);
+		while(~lg2){
-		while(lg2){
+			Wn[lg2]=pos,pos+=n;
-			m<<=1,lg2--;
+			Wn[lg2][0]=1;
-			Wn[lg2][0]=1LL*Wn[lg2+1][0]*Wn[lg2+1][0]%Mod;
+			_for(i,1,n)Wn[lg2][i]=1LL*Wn[lg2][i-1]*w%Mod;
-			Wn[lg2][1]=1LL*Wn[lg2+1][1]*Wn[lg2+1][1]%Mod;
+			w=1LL*w*w%Mod;
+			lg2--;n>>=1;
 		}
 	}
@@ 行 479: / 行 480: @@
 		swap(f[i],f[rev[i]]);
 		int t1,t2;
-		for(int i=1,lg2=0;i<n;i<<=1,lg2++){
+		for(int i=1,lg2=1;i<n;i<<=1,lg2++){
-			int w=Wn[lg2+1][type];
 			for(int j=0;j<n;j+=(i<<1)){
-				int cur=1;
 				_for(k,j,j+i){
-					t1=f[k],t2=1LL*cur*f[k+i]%Mod;
+					t1=f[k],t2=1LL*Wn[lg2][k-j]*f[k+i]%Mod;
 					f[k]=(t1+t2)%Mod,f[k+i]=(t1-t2)%Mod;
-					cur=1LL*cur*w%Mod;
 				}
 			}
 		}
 		if(!type){
+			reverse(f+1,f+n);
 			int div=quick_pow(n,Mod-2);
 			_for(i,0,n)f[i]=(1LL*f[i]*div%Mod+Mod)%Mod;
@@ 行 565: / 行 564: @@
 	_for(i,0,k)ans=(ans+1LL*f[i]*a[i])%Mod;
 	enter((ans+Mod)%Mod);
+	return 0;
+}
+</code>
+</hidden>
+===== 多项式多点求值 =====
+==== 算法简介 ====
+给定一个 $n$ 次多项式 $f(x)$，求 $f(a_i)(1\le i\le m)$。时间复杂度 $O(n\log n+m\log^2 m)$，空间复杂度 $O(m\log m)$。
+==== 算法实现 ====
+[[https://www.luogu.com.cn/problem/P5050|洛谷p5050]]
+设 $f(x)=c_0+c_1x+c_2x^2\cdots c_nx^n$，于是有
+$$
+\begin{equation}\begin{split}
+f(x)-f(x_0)&\equiv c_0-c_0+c_1(x-x_0)+c_2(x^2-x_0^2)\cdots +c_n(x^n-x_0^n)\\
+&\equiv c_1(x-x_0)+c_2(x-x_0)(x+x_0)\cdots +c_n(x-x_0)(x^{n-1}+x^{n-2}x_0+\cdots +x_0^{n-1})\\
+&\equiv 0\pmod {x-x_0}
+\end{split}\end{equation}
+$$
+即 $f(x_0)\equiv f(x)\pmod {x-x_0}$，于是如果能快速计算 $f(x)\bmod {x-x_0}$，就可以得到 $f(x_0)$。
+记 $f_{l,r}(x)\equiv f(x)\pmod {\prod_{i=l}^r (x-a_i)}$，于是有 $f_{l,r}(x)\equiv f_{l,m}(x)\pmod{\prod_{i=l}^m (x-a_i)},f_{l,r}(x)\equiv f_{m+1,r}(x)\pmod{\prod_{i=m+1}^r (x-a_i)}$。
+考虑 $O(m\log^2 m)$ 自底向上建立线段树，每个节点存储 $\prod_{i=l}^r (x-a_i)$ 的多项式。
+然后 $O(m\log^2 m)$ 自顶向下进行带余除法，跑到叶子节点就可以得到每个 $f(a_i)$ 的值。
+注意根节点需要提前进行带余除法，同时该算法常数较大，考虑常数优化。
+不妨令 $r-l$ 小于某个固定值 $L$ 后采用 $O(L^2)$ 暴力秦九韶展开，秦九韶复杂度为 $O(\cfrac mLL^2)=O(mL)$。
+于是总时间复杂变为 $O(n\log n+m\log m(\log m-\log L)+mL)$，经检验 $m=64000$ 时取 $L=640$ 效果较佳。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=64005,MINL=640,Mod=998244353;
+int quick_pow(int a,int b){
+	int ans=1;
+	while(b){
+		if(b&1)
+		ans=1LL*ans*a%Mod;
+		a=1LL*a*a%Mod;
+		b>>=1;
+	}
+	return ans;
+}
+namespace Poly{
+	const int G=3;
+	int rev[MAXN<<2],Pool[MAXN<<3],*Wn[30];
+	void init(){
+		int lg2=0,*pos=Pool,n,w;
+		while((1<<lg2)<MAXN*2)lg2++;
+		n=1<<lg2,w=quick_pow(G,(Mod-1)/(1<<lg2));
+		while(~lg2){
+			Wn[lg2]=pos,pos+=n;
+			Wn[lg2][0]=1;
+			_for(i,1,n)Wn[lg2][i]=1LL*Wn[lg2][i-1]*w%Mod;
+			w=1LL*w*w%Mod;
+			lg2--;n>>=1;
+		}
+	}
+	int build(int k){
+		int n,pos=0;
+		while((1<<pos)<=k)pos++;
+		n=1<<pos;
+		_for(i,0,n)rev[i]=(rev[i>>1]>>1)|((i&1)<<(pos-1));
+		return n;
+	}
+	void NTT(int *f,int n,bool type){
+		_for(i,0,n)if(i<rev[i])
+		swap(f[i],f[rev[i]]);
+		int t1,t2;
+		for(int i=1,lg2=1;i<n;i<<=1,lg2++){
+			for(int j=0;j<n;j+=(i<<1)){
+				_for(k,j,j+i){
+					t1=f[k],t2=1LL*Wn[lg2][k-j]*f[k+i]%Mod;
+					f[k]=(t1+t2)%Mod,f[k+i]=(t1-t2)%Mod;
+				}
+			}
+		}
+		if(!type){
+			reverse(f+1,f+n);
+			int div=quick_pow(n,Mod-2);
+			_for(i,0,n)f[i]=(1LL*f[i]*div%Mod+Mod)%Mod;
+		}
+	}
+	void Mul(int *f,int _n,int *g,int _m,int xmod=0){
+		int n=build(_n+_m-2);
+		_for(i,_n,n)f[i]=0;_for(i,_m,n)g[i]=0;
+		NTT(f,n,true);NTT(g,n,true);
+		_for(i,0,n)f[i]=1LL*f[i]*g[i]%Mod;
+		NTT(f,n,false);
+		if(xmod)_for(i,xmod,n)f[i]=0;
+	}
+	void Mul2(const int *f,int _n,const int *g,int _m,int *h){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		int n=build(_n+_m-2);
+		memcpy(temp1,f,sizeof(int)*_n);memcpy(temp2,g,sizeof(int)*_m);
+		_for(i,_n,n)temp1[i]=0;_for(i,_m,n)temp2[i]=0;
+		NTT(temp1,n,true);NTT(temp2,n,true);
+		_for(i,0,n)temp1[i]=1LL*temp1[i]*temp2[i]%Mod;
+		NTT(temp1,n,false);
+		n=_n+_m-1;
+		_for(i,0,n)h[i]=temp1[i];
+	}
+	void Inv(const int *f,int *g,int _n){
+		static int temp[MAXN<<2];
+		if(_n==1)return g[0]=quick_pow(f[0],Mod-2),void();
+		Inv(f,g,(_n+1)>>1);
+		int n=build((_n-1)<<1);
+		_for(i,(_n+1)>>1,n)g[i]=0;
+		_for(i,0,_n)temp[i]=f[i];_for(i,_n,n)temp[i]=0;
+		NTT(g,n,true);NTT(temp,n,true);
+		_for(i,0,n)g[i]=(2-1LL*temp[i]*g[i]%Mod)*g[i]%Mod;
+		NTT(g,n,false);
+		_for(i,_n,n)g[i]=0;
+	}
+	void Div(const int *f,int _n,const int *g,int _m,int *r){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		if(_n<_m){
+			_rep(i,0,_n)r[i]=f[i];
+			return;
+		}
+		_rep(i,0,_m)temp1[i]=g[_m-i];
+		Inv(temp1,temp2,_n-_m+1);
+		_rep(i,0,_n)temp1[i]=f[_n-i];
+		Mul(temp2,_n-_m+1,temp1,_n+1,_n-_m+1);
+		for(int i=0,j=_n-_m;i<j;i++,j--)swap(temp2[i],temp2[j]);
+		int __m=min(_n-_m+1,_m);
+		_for(i,0,_m)temp1[i]=g[i];
+		Mul(temp1,_m,temp2,__m);
+		_for(i,0,_m)r[i]=(f[i]+Mod-temp1[i])%Mod;
+	}
+}
+int a[MAXN],b[MAXN],c[MAXN],pool[MAXN*40],*pos=pool,*g[MAXN<<2];
+void build(int k,int L,int R){
+	g[k]=pos,pos+=R-L+2;
+	int M=L+R>>1;
+	if(L==R)return g[k][0]=Mod-a[M],g[k][1]=1,void();
+	build(k<<1,L,M);build(k<<1|1,M+1,R);
+	Poly::Mul2(g[k<<1],M-L+2,g[k<<1|1],R-M+1,g[k]);
+}
+void query(int k,int L,int R,int *f){
+	if(R-L<MINL){
+		_rep(i,L,R){
+			b[i]=0;
+			for(int j=R-L;~j;j--)
+			b[i]=(1LL*b[i]*a[i]+f[j])%Mod;
+		}
+		return;
+	}
+	int *temp=pos,M=L+R>>1;pos+=R-L+1;
+	Poly::Div(f,R-L,g[k<<1],M-L+1,temp);
+	query(k<<1,L,M,temp);
+	Poly::Div(f,R-L,g[k<<1|1],R-M,temp);
+	query(k<<1|1,M+1,R,temp);
+}
+int main()
+{
+	Poly::init();
+	int n=read_int(),m=read_int();
+	_rep(i,0,n)c[i]=read_int();
+	_for(i,0,m)a[i]=read_int();
+	build(1,0,m-1);
+	int *f=pos;pos+=m;
+	Poly::Div(c,n,g[1],m,f);
+	query(1,0,m-1,f);
+	_for(i,0,m)enter(b[i]);
+	return 0;
+}
+</code>
+</hidden>
+===== 多项式快速插值 =====
+==== 算法简介 ====
+给定 $n$ 个点 $(x_i,y_i)$，求次数不超过 $n-1$ 次的多项式 $f(x)$ 满足 $f(x_i)=y_i$。时间复杂度 $O(n\log^2 n)$，空间复杂度 $O(m\log m)$。
+==== 算法实现 ====
+[[https://www.luogu.com.cn/problem/P5158|洛谷p5158]]
+根据拉格朗日插值法，有
+$$f(x)=\sum_{i=1}^n y_i\prod_{j=1,j\neq i}^n \frac {x-x_j}{x_i-x_j}=\sum_{i=1}^n \cfrac {y_i}{\prod_{j=1,j\neq i}(x_i-x_j)}\prod_{j=1,j\neq i}^n x-x_j$$
+设 $g(x)=\prod_{i=1}^n x-x_i$，根据洛必达法则，有
+$$\prod_{j=1,j\neq i}^n x_i-x_j=\lim_{x\to x_i}\frac {g(x)}{x-x_i}=g'(x_i)$$
+于是有 $f(x)=\sum_{i=1}^n \cfrac {y_i}{g'(x_i)}\prod_{j=1,j\neq i}^n x-x_j$。
+考虑 $O(n\log^2 n)$ 分治乘法计算出 $g(x)$，再通过 $O(n\log^2 n)$ 多项式多点求值计算出所有 $g'(x_i)$。
+设 $f_{l,r}(x)=\sum_{i=l}^r \cfrac {y_i}{g'(x_i)}\prod_{j=l,j\neq i}^r x-x_j$。考虑分治，有
+$$
+\begin{equation}\begin{split}
+f_{l,r}(x)&=\left(\sum_{i=l}^m \cfrac {y_i}{g'(x_i)}\prod_{j=l,j\neq i}^m x-x_j\right)\left(\prod_{j=m+1}^r x-x_j\right)+\left(\sum_{i=m+1}^r \cfrac {y_i}{g'(x_i)}\prod_{j=m+1,j\neq i}^r x-x_j\right)\left(\prod_{j=l}^m x-x_j\right)\\
+&=f_{l,m}(x)\left(\prod_{j=m+1}^r x-x_j\right)+f_{m+1,r}(x)\left(\prod_{j=l}^m x-x_j\right)
+\end{split}\end{equation}
+$$
+于是可以 $O(n\log^2 n)$ 完成分治，总时间复杂度 $O(n\log^2 n)$。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=1e5+5,MINL=640,Mod=998244353;
+int quick_pow(int a,int b){
+	int ans=1;
+	while(b){
+		if(b&1)
+		ans=1LL*ans*a%Mod;
+		a=1LL*a*a%Mod;
+		b>>=1;
+	}
+	return ans;
+}
+namespace Poly{
+	const int G=3;
+	int rev[MAXN<<2],Pool[MAXN<<3],*Wn[30];
+	void init(){
+		int lg2=0,*pos=Pool,n,w;
+		while((1<<lg2)<MAXN*2)lg2++;
+		n=1<<lg2,w=quick_pow(G,(Mod-1)/(1<<lg2));
+		while(~lg2){
+			Wn[lg2]=pos,pos+=n;
+			Wn[lg2][0]=1;
+			_for(i,1,n)Wn[lg2][i]=1LL*Wn[lg2][i-1]*w%Mod;
+			w=1LL*w*w%Mod;
+			lg2--;n>>=1;
+		}
+	}
+	int build(int k){
+		int n,pos=0;
+		while((1<<pos)<=k)pos++;
+		n=1<<pos;
+		_for(i,0,n)rev[i]=(rev[i>>1]>>1)|((i&1)<<(pos-1));
+		return n;
+	}
+	void NTT(int *f,int n,bool type){
+		_for(i,0,n)if(i<rev[i])
+		swap(f[i],f[rev[i]]);
+		int t1,t2;
+		for(int i=1,lg2=1;i<n;i<<=1,lg2++){
+			for(int j=0;j<n;j+=(i<<1)){
+				_for(k,j,j+i){
+					t1=f[k],t2=1LL*Wn[lg2][k-j]*f[k+i]%Mod;
+					f[k]=(t1+t2)%Mod,f[k+i]=(t1-t2)%Mod;
+				}
+			}
+		}
+		if(!type){
+			reverse(f+1,f+n);
+			int div=quick_pow(n,Mod-2);
+			_for(i,0,n)f[i]=(1LL*f[i]*div%Mod+Mod)%Mod;
+		}
+	}
+	void Mul(int *f,int _n,int *g,int _m,int xmod=0){
+		int n=build(_n+_m-2);
+		_for(i,_n,n)f[i]=0;_for(i,_m,n)g[i]=0;
+		NTT(f,n,true);NTT(g,n,true);
+		_for(i,0,n)f[i]=1LL*f[i]*g[i]%Mod;
+		NTT(f,n,false);
+		if(xmod)_for(i,xmod,n)f[i]=0;
+	}
+	void Mul2(const int *f,int _n,const int *g,int _m,int *h){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		int n=build(_n+_m-2);
+		memcpy(temp1,f,sizeof(int)*_n);memcpy(temp2,g,sizeof(int)*_m);
+		_for(i,_n,n)temp1[i]=0;_for(i,_m,n)temp2[i]=0;
+		NTT(temp1,n,true);NTT(temp2,n,true);
+		_for(i,0,n)temp1[i]=1LL*temp1[i]*temp2[i]%Mod;
+		NTT(temp1,n,false);
+		n=_n+_m-1;
+		_for(i,0,n)h[i]=temp1[i];
+	}
+	void Mul3(int *f,int _n,int *g,int _m){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		int n=build(_n+_m-2);
+		memcpy(temp1,f,sizeof(int)*_n);memcpy(temp2,g,sizeof(int)*_m);
+		_for(i,_n,n)temp1[i]=0;_for(i,_m,n)temp2[i]=0;
+		NTT(temp1,n,true);NTT(temp2,n,true);
+		_for(i,0,n)temp1[i]=1LL*temp1[i]*temp2[i]%Mod;
+		NTT(temp1,n,false);
+		n=_n+_m-1;
+		_for(i,0,n)f[i]=temp1[i];
+	}
+	void Inv(const int *f,int *g,int _n){
+		static int temp[MAXN<<2];
+		if(_n==1)return g[0]=quick_pow(f[0],Mod-2),void();
+		Inv(f,g,(_n+1)>>1);
+		int n=build((_n-1)<<1);
+		_for(i,(_n+1)>>1,n)g[i]=0;
+		_for(i,0,_n)temp[i]=f[i];_for(i,_n,n)temp[i]=0;
+		NTT(g,n,true);NTT(temp,n,true);
+		_for(i,0,n)g[i]=(2-1LL*temp[i]*g[i]%Mod)*g[i]%Mod;
+		NTT(g,n,false);
+		_for(i,_n,n)g[i]=0;
+	}
+	void Div(const int *f,int _n,const int *g,int _m,int *r){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		if(_n<_m){
+			_rep(i,0,_n)r[i]=f[i];
+			return;
+		}
+		_rep(i,0,_m)temp1[i]=g[_m-i];
+		Inv(temp1,temp2,_n-_m+1);
+		_rep(i,0,_n)temp1[i]=f[_n-i];
+		Mul(temp2,_n-_m+1,temp1,_n+1,_n-_m+1);
+		for(int i=0,j=_n-_m;i<j;i++,j--)swap(temp2[i],temp2[j]);
+		int __m=min(_n-_m+1,_m);
+		_for(i,0,_m)temp1[i]=g[i];
+		Mul(temp1,_m,temp2,__m);
+		_for(i,0,_m)r[i]=(f[i]+Mod-temp1[i])%Mod;
+	}
+}
+int x[MAXN],y[MAXN],gy[MAXN],pool[MAXN*50],*pos=pool,*t[MAXN<<2];
+void build(int k,int L,int R){
+	t[k]=pos,pos+=R-L+2;
+	int M=L+R>>1;
+	if(L==R)return t[k][0]=Mod-x[M],t[k][1]=1,void();
+	build(k<<1,L,M);build(k<<1|1,M+1,R);
+	Poly::Mul2(t[k<<1],M-L+2,t[k<<1|1],R-M+1,t[k]);
+}
+void query(int k,int L,int R,int *f){
+	if(R-L<MINL){
+		_rep(i,L,R){
+			gy[i]=0;
+			for(int j=R-L;~j;j--)
+			gy[i]=(1LL*gy[i]*x[i]+f[j])%Mod;
+		}
+		return;
+	}
+	int *temp=pos,M=L+R>>1;pos+=R-L+1;
+	Poly::Div(f,R-L,t[k<<1],M-L+1,temp);
+	query(k<<1,L,M,temp);
+	Poly::Div(f,R-L,t[k<<1|1],R-M,temp);
+	query(k<<1|1,M+1,R,temp);
+}
+void query2(int k,int L,int R,int *f){
+	int M=L+R>>1;
+	if(L==R)return f[0]=1LL*y[M]*quick_pow(gy[M],Mod-2)%Mod,void();
+	int *temp=pos;pos+=R-L+1;
+	query2(k<<1,L,M,f);query2(k<<1|1,M+1,R,temp);
+	Poly::Mul3(f,M-L+1,t[k<<1|1],R-M+1);Poly::Mul3(temp,R-M,t[k<<1],M-L+2);
+	int n=R-L+1;
+	_for(i,0,n)f[i]=(f[i]+temp[i])%Mod;
+}
+int main()
+{
+	Poly::init();
+	int n=read_int();
+	_for(i,0,n)x[i]=read_int(),y[i]=read_int();
+	build(1,0,n-1);
+	int *g=pos;pos+=n;
+	memcpy(g,t[1],sizeof(int)*(n+1));
+	_for(i,0,n)g[i]=1LL*g[i+1]*(i+1)%Mod;
+	query(1,0,n-1,g);
+	int *f=pos;pos+=n;
+	query2(1,0,n-1,f);
+	_for(i,0,n)space(f[i]);
 	return 0;
 }
 </code>
 </hidden>

CVBB ACM Team

用户工具

站点工具

差别

页面工具