差别

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

--- 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/23 17:38]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:多项式_4 [2020/08/25 17:43] (当前版本)
jxm2001
@@ 行 9: / 行 9: @@
 ==== 性质 ====
-对序列 $A,B$ 做长度为 $n$ 的 $\text{FFT}$ 等价于求序列 $A,B$ 的循环卷积。
+对序列 $A,B$ 做长度为 $n$ 的 $\text{FFT}$ 等价于求序列 $A,B$ 的长度为 $n$ 的循环卷积。
 考虑单位根反演证明
@@ 行 25: / 行 25: @@
 事实上普通的卷积计算相当于长度为 $2^n$ 的循环卷积计算，只是循环卷积长度大于 $C$ 序列的长度，所以循环卷积结果即为普通卷积结果。
+值得一提的是，循环卷积的求逆、快速幂等操作直接对 $\text{DFT}$ 的点值进行相应运算即可。
+另外注意到 $f(w^{-k})=f(w^{n-k})$，于是循环卷积的 $\text{IDFT}$ 过程可以通过将序列 $[1,n-1]$ 部分翻转再 $\text{DFT}$ 的方式实现。
+==== Cooley–Tukey FFT algorithm ====
+=== 算法实现 ===
+普通 $\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}$ 过程，有
@@ 行 38: / 行 166: @@
 $$
-其中 ${n\choose 2}$ 表示组合数。易知上式可以通过普通卷积求解，$\text{IDFT}$ 过程类似。
+其中 ${n\choose 2}$ 表示组合数，易知上式可以通过普通卷积求解。
-值得一提的是，循环卷积的求逆、快速幂等操作直接对 $\text{DFT}$ 的点值进行相应运算即可。
+考虑 $\text{IDFT}$ 过程，同样有
-==== 算法例题 ====
+$$
+\begin{equation}\begin{split}
+a_k&=\frac 1n\sum_{i=0}^{n-1}y_iw_n^{-ki}\\
+&=\frac 1n\sum_{i=0}^{n-1}y_iw_n^{-{k+i\choose 2}+{k\choose 2}+{i\choose 2}}\\
+&=\frac 1nw_n^{{k\choose 2}}\sum_{i=0}^{n-1}\left(y_iw_n^{{i\choose 2}}\right)w_n^{-{k+i\choose 2}}
+\end{split}\end{equation}
+$$
+$\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}$ 条重边。
+假设起点为 $(0,x)$，终点为 $(\ast,y)$($\ast$ 为任意值)，路径长度 $m$ 定义为路径的边数。
+对每个 $0\le t\lt k$，求满足所有 $m\equiv t\pmod k$ 且横坐标单增的路径数模 $p$ 意义下的值。
+数据保证 $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$ 的路径数。
+于是有状态转移方程
+$$g_{a,b}=\sum_{i=1}^n g_{a-1,i}w_{i,b}\text{ },\text{ }f_{a,b}={L\choose a}g_{a,b}$$
+设矩阵
+$$W=\begin{bmatrix}w_{1,1} & \cdots & w_{1,n} \\
+\vdots &\ddots & \vdots \\
+w_{n,1} &\cdots & w_{n,n}\end{bmatrix}$$
+设 $G_i=(g_{i,1}\cdots g_{i,n})$，于是有 $G_i=G_0W^i$。
+考虑单位根反演，设 $w_k\equiv g^{\frac {p-1}k}\pmod p,g$ 为 $p$ 的原根，有
+$$
+\begin{equation}\begin{split}
+\text{ans}_t&=\sum_{i=0}^Lf_{i,y}[i\bmod k=t]\\
+&=\frac 1k\sum_{i=0}^Lf_{i,y}\sum_{j=0}^{k-1}w_k^{(i-t)j}\\
+&=\frac 1k\sum_{j=0}^{k-1}w_k^{-tj}\sum_{i=0}^L f_{i,y}w_k^{ij}
+\end{split}\end{equation}
+$$
+根据二项式定理，有
+$$\sum_{i=0}^L w_k^{ij}(f_{i,1}\cdots f_{i,n})=\sum_{i=0}^L w_k^{ij}{L\choose i}(g_{i,1}\cdots g_{i,n})=G_0\sum_{i=0}^L {L\choose i}\left(w_k^jW\right)^i=G_0\left(w_k^jW+I\right)^L$$
+于是根据矩阵快速幂可以 $O(kn^3\log L)$ 计算出所有 $h_j=\sum_{i=0}^L f_{i,y}w_k^{ij}(0\le j\lt k)$。
+于是有
+$$\text{ans}_t=\frac 1k\sum_{i=0}^{k-1}w_k^{-ti}h_i$$
+发现上式就是 $\text{Bluestein's Algotithm}$ 的 $\text{IDFT}$ 过程，直接求解时间复杂度 $O(k\log k)$。
+总时间复杂度 $O(kn^3\log L)$，主要用于计算矩阵快速幂。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=1<<16|5;
+const long double pi=acos(-1.0);
+int p,sz,gw[MAXN];
+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;
+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 k){
+	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;
+		}
+	}
+	int w=quick_pow(g,(p-1)/k);
+	gw[0]=1;
+	_rep(i,1,k)gw[i]=1LL*gw[i-1]*w%p;
+}
+struct Matrix{
+	int ele[3][3];
+	Matrix(int x=0){
+		mem(ele,0);
+		_for(i,0,sz)ele[i][i]=x;
+	}
+	Matrix operator + (const Matrix b){
+		Matrix c;
+		_for(i,0,sz)
+		_for(j,0,sz)
+		c.ele[i][j]=(ele[i][j]+b.ele[i][j])%p;
+		return c;
+	}
+	Matrix operator * (const int b){
+		Matrix c;
+		_for(i,0,sz)
+		_for(j,0,sz)
+		c.ele[i][j]=1LL*ele[i][j]*b%p;
+		return c;
+	}
+	void operator *= (const Matrix b){
+		Matrix c=*this;
+		_for(i,0,sz)
+		_for(j,0,sz){
+			ele[i][j]=0;
+			_for(k,0,sz)
+			ele[i][j]=(ele[i][j]+1LL*c.ele[i][k]*b.ele[k][j])%p;
+		}
+	}
+};
+Matrix quick_pow(const Matrix &a,int k){
+	Matrix Ans(1),Base;
+	Base=a;
+	while(k){
+		if(k&1)Ans*=Base;
+		k>>=1;
+		Base*=Base;
+	}
+	return Ans;
+}
+struct complex{
+	long double x,y;
+	complex(long double x=0.0,long double y=0.0):x(x),y(y){}
+	complex operator + (const complex &b){
+		return complex(x+b.x,y+b.y);
+	}
+	complex operator - (const complex &b){
+		return complex(x-b.x,y-b.y);
+	}
+	complex operator * (const complex &b){
+		return complex(x*b.x-y*b.y,x*b.y+y*b.x);
+	}
+};
+int rev[MAXN<<2];
+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 FFT(complex *f,int n,int type){
+	_for(i,0,n)if(i<rev[i])
+	swap(f[i],f[rev[i]]);
+	complex t1,t2;
+	for(int i=1;i<n;i<<=1){
+		complex w(cos(pi/i),type*sin(pi/i));
+		for(int j=0;j<n;j+=(i<<1)){
+			complex cur(1.0,0.0);
+			_for(k,j,j+i){
+				t1=f[k],t2=cur*f[k+i];
+				f[k]=t1+t2,f[k+i]=t1-t2;
+				cur=cur*w;
+			}
+		}
+	}
+	if(type==-1)_for(i,0,n)
+	f[i].x/=n,f[i].y/=n;
+}
+void FFT2(complex *f1,complex *f2,int n){
+	FFT(f1,n,1);
+	f2[0].x=f1[0].x,f2[0].y=-f1[0].y;
+	_for(i,1,n)
+	f2[i].x=f1[n-i].x,f2[i].y=-f1[n-i].y;
+	complex t1,t2;
+	_for(i,0,n){
+		t1=f1[i],t2=f2[i];
+		f1[i]=complex((t1.x+t2.x)*0.5,(t1.y+t2.y)*0.5);
+		f2[i]=complex((t1.y-t2.y)*0.5,(t2.x-t1.x)*0.5);
+	}
+}
+void MTT(int *f,int n1,int *g,int n2,int *ans,int mod){
+	static complex f1[MAXN<<2],f2[MAXN<<2],g1[MAXN<<2],g2[MAXN<<2],temp[2][MAXN<<2];
+	int n=build(n1+n2),m=sqrt(mod)+1;
+	_rep(i,0,n1)f1[i].x=f[i]/m,f1[i].y=f[i]%m;
+	_rep(i,0,n2)g1[i].x=g[i]/m,g1[i].y=g[i]%m;
+	FFT2(f1,f2,n);FFT2(g1,g2,n);
+	complex I(0.0,1.0);
+	_for(i,0,n){
+		temp[0][i]=f1[i]*g1[i]+I*f2[i]*g1[i];
+		temp[1][i]=f1[i]*g2[i]+I*f2[i]*g2[i];
+	}
+	FFT(temp[0],n,-1);FFT(temp[1],n,-1);
+	LL a,b,c;
+	_rep(i,0,n1+n2){
+		a=temp[0][i].x+0.5,b=temp[0][i].y+temp[1][i].x+0.5,c=temp[1][i].y+0.5;
+		ans[i]=((a%mod*m%mod*m%mod+b%mod*m%mod+c%mod)%mod+mod)%mod;
+	}
+}
+Matrix W;
+int a[MAXN<<2],b[MAXN<<2],c[MAXN<<2];
+int main()
+{
+	int k,L,x,y;
+	sz=read_int(),k=read_int(),L=read_int(),x=read_int()-1,y=read_int()-1,p=read_int();
+	get_G(k);
+	_for(i,0,sz)
+	_for(j,0,sz)
+	W.ele[i][j]=read_int();
+	_for(i,0,k){
+		Matrix temp=quick_pow(W*gw[i]+Matrix(1),L);
+		a[i]=1LL*temp.ele[x][y]*gw[1LL*i*(i-1)/2%k]%p;
+	}
+	_for(i,0,2*k-1)b[i]=gw[k-1LL*i*(i-1)/2%k]%p;
+	reverse(a,a+k);
+	MTT(a,k-1,b,2*k-2,c,p);
+	int div=quick_pow(k,p-2);
+	_for(i,0,k)
+	enter(1LL*div*gw[1LL*i*(i-1)/2%k]%p*c[k-1+i]%p);
+	return 0;
+}
+</code>
+</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

用户工具

站点工具

差别

页面工具