差别

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

--- 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/24 14:46]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/25 17:43] (当前版本)
jxm2001
@@ 行 32: / 行 32: @@
 ==== Cooley–Tukey FFT algorithm ====
-普通 $\text{DFT}$ 过程是将序列长度根据奇偶幂次分成两段经行递归。现考虑将序列分成 $d$ 段进行递归。
+=== 算法实现 ===
-$$f(w_n)$$
+普通 $\text{DFT}$ 过程是将序列根据奇偶幂次分成两段经行递归，现考虑将序列分成 $d$ 段进行递归，设
+$$f(x)=a_0+a_1x^1+a_2x^2+\cdots a_{n-1}x^{n-1},f_k(x)=a_kx^k+a_{d+k}x^{d+k}+a_{2d+k}x^{2d+k}+\cdots+a_{n-d+k}x^{n-d+k}(0\le k\lt d),m=\frac nd$$
+于是有
+$$f(w_n^{im+j})=\sum_{k=0}^{d-1}w_n^{(im+j)k}f_k(w_n^{(im+j)d})=\sum_{k=0}^{d-1}w_n^{(im+j)k}f_k(w_n^{in}w_n^{jd})=\sum_{k=0}^{d-1}w_n^{(im+j)k}f_k(w_m^j)$$
+计算出每个点值的时间为 $O(d)$，每层共 $O(n)$ 个点值。设 $n=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_k^{\alpha_k}$，于是总时间复杂度 $O\left(n\sum_{i=1}^k p_i\alpha_i\right)$。
+=== 算法例题 ===
+[[https://www.luogu.com.cn/problem/P4191|洛谷p4191]]
+给定长度为 $n$ 的序列 $A,B$，计算 $AB^C$ 的长度为 $n$ 的循环卷积模 $n+1$ 意义下的值。
+数据保证 $n+1$ 为素数，$n\le 5\times 10^5$，$n$ 的最大质因子不超过 $10$。
+== 题解 ==
+$\text{Cooley–Tukey FFT algorithm}$ 板子题，$\text{DFT}$ 后直接对点值快速幂即可，时间复杂度 $O(7n\log n)$
+<hidden 查看代码>
+<code cpp>
+const int MAXN=5e5+5,MAXM=20,MAX_div=7;
+int p;
+int quick_pow(int a,int b){
+	int ans=1;
+	while(b){
+		if(b&1)ans=1LL*ans*a%p;
+		a=1LL*a*a%p;
+		b>>=1;
+	}
+	return ans;
+}
+vector<int> pdiv,frac,Wn[MAXM];
+bool check(int x){
+	_for(i,0,pdiv.size()){
+		if(quick_pow(x,(p-1)/pdiv[i])==1)
+		return false;
+	}
+	return true;
+}
+void get_G(int n){
+	int temp=p-1,g;
+	for(int i=2;i*i<=temp;i++){
+		if(temp%i==0){
+			pdiv.push_back(i);
+			while(temp%i==0)temp/=i;
+		}
+	}
+	if(temp!=1)pdiv.push_back(temp);
+	_for(i,2,p){
+		if(check(i)){
+			g=i;
+			break;
+		}
+	}
+	temp=n;
+	for(int i=2;i*i<=temp;i++){
+		while(temp%i==0){
+			frac.push_back(i);
+			temp/=i;
+		}
+	}
+	if(temp!=1)frac.push_back(temp);
+	int len=n,wn;
+	_for(i,0,frac.size()){
+		Wn[i].resize(len);
+		wn=quick_pow(g,(p-1)/len),Wn[i][0]=1;
+		_for(j,1,len)Wn[i][j]=1LL*Wn[i][j-1]*wn%p;
+		len/=frac[i];
+	}
+}
+int temp[MAXN];
+void DFT(int *f,int n,int dep=0){
+	if(n==1)return;
+	memcpy(temp,f,sizeof(int)*n);
+	int m=n/frac[dep],d=frac[dep],*g[MAX_div];
+	_for(i,0,d)g[i]=f+i*m;
+	for(int i=0,k=0;i<n;i+=d,k++){
+		_for(j,0,d)
+		g[j][k]=temp[i+j];
+	}
+	_for(i,0,d)DFT(g[i],m,dep+1);
+	_for(i,0,d){
+		_for(j,0,m){
+			int pos=i*m+j;
+			temp[pos]=0;
+			_for(k,0,d)
+			temp[pos]=(temp[pos]+1LL*Wn[dep][pos*k%n]*g[k][j])%p;
+		}
+	}
+	memcpy(f,temp,sizeof(int)*n);
+}
+void IDFT(int *f,int n){
+	reverse(f+1,f+n);
+	DFT(f,n);
+	int div=quick_pow(n,p-2);
+	_for(i,0,n)f[i]=1LL*f[i]*div%p;
+}
+int a[MAXN],b[MAXN];
+int main()
+{
+	int n=read_int(),k=read_int();
+	p=n+1;
+	get_G(n);
+	_for(i,0,n)a[i]=read_int();
+	_for(i,0,n)b[i]=read_int();
+	DFT(a,n);DFT(b,n);
+	_for(i,0,n)a[i]=1LL*a[i]*quick_pow(b[i],k)%p;
+	IDFT(a,n);
+	_for(i,0,n)enter(a[i]);
+	return 0;
+}
+</code>
+</hidden>
 ==== Bluestein's Algotithm ====
+=== 算法实现 ===
 考虑 $\text{DFT}$ 过程，有
@@ 行 62: / 行 180: @@
 $\text{ }$
-==== 算法例题 ====
+=== 算法例题 ===
 [[https://www.luogu.com.cn/problem/P5293|洛谷p5293]]
-=== 题意 ===
 给定一个二维 $[L+1]\times n$ 的空间，其中 $(u_1,v_1)\to (u_2,v_2)$ 有 $w_{v_1,v_2}$ 条重边。
@@ 行 76: / 行 192: @@
 数据保证 $p$ 为素数，$10^8\le p\le 2^{30},k\mid p-1,1\le k\lt 65536,1\le n\le 3,L\le 10^9$。
-=== 题解 ===
+== 题解 ==
 假设 $f_{a,b}$ 表示 $m=a$ 且 $y=b$ 的路径数，$g_{a,b}$ 表示将空间的 $X$ 维消去后 $m=a$ 且 $y=b$ 的路径数。
@@ 行 291: / 行 407: @@
 </hidden>
-===== 线性递推 =====
+===== 常系数齐次线性递推 =====
+==== 算法简介 ====
+给定 $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)$。
+==== 算法实现 ====
+[[https://www.luogu.com.cn/problem/P4723|洛谷p4723]]
+考虑求斐波那契数列过程，$f_5=f_4+f_3=2f_3+f_2=3f_2+2f_1=5f_1+3f_0$。
+从多项式角度上考虑该过程，$f_n=f_{n-1}+f_{n-2}$ 的特征多项式为 $x^2-x-1$，并且有
+$$
+\begin{matrix}
+&0x^0&+0x^1&+0x^2&+0x^3&+0x^4&+1x^5\\
+\equiv&0x^0&+0x^1&+0x^2&+1x^3&+1x^4&+0x^5\\
+\equiv&0x^0&+0x^1&+1x^2&+2x^3&+0x^4&+0x^5\\
+\equiv&0x^0&+3x^1&+2x^2&+0x^3&+0x^4&+0x^5\\
+\equiv&3x^0&+5x^1&+0x^2&+0x^3&+0x^4&+0x^5&\pmod {x^2-x-1}
+\end{matrix}
+$$
+于是 $f_5=\left(x^5\bmod {x^2-x-1}\right)\cdot (f_0,f_1)=(3,5)\cdot (f_0,f_1)$。
+以此类推，对 $a_n=f_1a_{n-1}+f_2a_{n-2}+\cdots +f_ka_{n-k}$，其特征多项式为 $x^k-f_1x^{k-1}-f_2x^{k-2}-\cdots -f_k$。
+于是 $a_n= \left(x^n\bmod {x^k-f_1x^{k-1}-f_2x^{k-2}-\cdots -f_k}\right)\cdot (a_0,a_1\cdots a_{k-1})$。
+其中 $\left(x^n\bmod {x^k-f_1x^{k-1}-f_2x^{k-2}-\cdots -f_k}\right)$ 可以通过快速幂与多项式带余除法 $O(k\log k\log n)$ 计算。
+发现每次带余除法计算时 $g(x)$ 都是固定的，可以考虑先计算出 $\frac 1{g^R(x)}$ 减小常数。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=64005,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 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(int *f,int _n,const int *g,int _m,int *gR){
+		static int temp[MAXN<<2],q[MAXN<<2],r[MAXN<<2];
+		_rep(i,0,_n-_m)q[i]=gR[i];
+		_rep(i,0,_n)temp[i]=f[_n-i];
+		Mul(q,_n-_m+1,temp,_n+1,_n-_m+1);
+		for(int i=0,j=_n-_m;i<j;i++,j--)swap(q[i],q[j]);
+		int __m=min(_n-_m+1,_m);
+		_for(i,0,_m)r[i]=g[i];_for(i,0,__m)temp[i]=q[i];
+		Mul(r,_m,temp,__m,_m);
+		_for(i,0,_m)f[i]=(f[i]+Mod-r[i])%Mod;
+	}
+	void Pow_2(int *f,int _n,const int *g,int _m){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		_rep(i,0,_m)temp1[i]=g[_m-i];
+		Inv(temp1,temp2,_m-1);
+		mem(temp1,0);
+		temp1[0]=1;
+		while(_n){
+			int n=build(_m+_m-2);
+			_for(i,_m,n)f[i]=0;
+			NTT(f,n,true);
+			if(_n&1){
+				_for(i,_m,n)temp1[i]=0;
+				NTT(temp1,n,true);
+				_for(i,0,n)temp1[i]=1LL*f[i]*temp1[i]%Mod;
+				NTT(temp1,n,false);
+			}
+			_for(i,0,n)f[i]=1LL*f[i]*f[i]%Mod;
+			NTT(f,n,false);
+			if(_n&1)Div(temp1,_m+_m-2,g,_m,temp2);
+			Div(f,_m+_m-2,g,_m,temp2);
+			_n>>=1;
+		}
+		mem(f,0);
+		_for(i,0,_m)f[i]=temp1[i];
+	}
+}
+int f[MAXN<<2],g[MAXN<<2],a[MAXN];
+int main()
+{
+	Poly::init();
+	int n=read_int(),k=read_int();
+	g[k]=1;
+	for(int i=k-1;i>=0;i--)g[i]=-read_int();
+	_for(i,0,k)a[i]=read_int();
+	f[1]=1;
+	Poly::Pow_2(f,n,g,k);
+	int ans=0;
+	_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

用户工具

站点工具

差别

页面工具