差别

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

--- 2020-2021:teams:legal_string:jxm2001:多项式_3 [2020/08/11 23:43]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:多项式_3 [2020/08/25 12:43] (当前版本)
jxm2001
@@ 行 1: / 行 1: @@
 ====== 多项式 3 ======
+===== 倍增 FFT =====
+====  算法简介 ====
+主要用于加速可以倍增转移的动态规划，时间复杂度 $O\left(n\log^2 n\right)$。
+==== 算法例题 ====
+[[https://www.luogu.com.cn/problem/CF755G|CF755G]]
+=== 题意 ===
+将 $n$ 个球排成一排。规定一个组至少包含一个球，至多包含两个相邻的球，且一个球至多属于一个组。
+问从这 $n$ 个球中取 $k$ 组有多少方案，答案对 $998244353$ 取模。
+=== 题解 ===
+设 $\text{dp}(i,j)$ 表示从 $i$ 个球里选 $j$ 组的方案数。
+对于第 $n$ 个球，要么不选，要么单独构成一个组，要么与第 $n-1$ 个球构成一个组，于是有状态转移方程
+$$\text{dp}(i,j)=\text{dp}(i-1,j)+\text{dp}(i-1,j-1)+\text{dp}(i-1,j-2)$$
+同时将 $a+b$ 个球划分为两组 $1\sim a$ 和 $a+1\sim a+b$，根据 $a,a+1$ 是否共同为一组可得状态转移方程
+$$\text{dp}(a+b,k)=\sum_{i=0}^k\text{dp}(a,i)\text{dp}(b,k-i)+\sum_{i=0}^{k-1}\text{dp}(a-1,i)\text{dp}(b-1,k-i-1)$$
+设 $F_n(x)=\sum_{i=0}^{\infty} \text{dp}(n,i)x^i$，于是有
+$$F_n(x)=F_{n-1}(x)+xF_{n-1}(x)+x^2F_{n-2}(x)$$
+$$F_{a+b}(x)=F_a(x)F_b(x)+xF_{a-1}(x)F_{b-1}(x)$$
+$$F_{2n-2}(x)=F_{n-1}^2(x)+xF_{n-2}^2(x)$$
+$$F_{2n-1}(x)=F_n(x)F_{n-1}(x)+xF_{n-1}(x)F_{n-2}(x)$$
+$$F_{2n}(x)=F_n^2(x)+xF_{n-1}^2(x)$$
+于是可以 $\left(F_{n-2},F_{n-1},F_n\right)\to \left(F_{n-1},F_n,F_{n+1}\right),\left(F_{2n-2},F_{2n-1},F_{2n}\right)$
+于是用类似快速幂的算法可以 $O(k\log n\log k)$ 解决上述问题。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=1<<16,Mod=998244353,G=3;
+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;
+}
+int rev[MAXN],Wn[30][2];
+void init(){
+	int m=Mod-1,lg2=0;
+	while(m%2==0)m>>=1,lg2++;
+	Wn[lg2][1]=quick_pow(G,m);
+	Wn[lg2][0]=quick_pow(Wn[lg2][1],Mod-2);
+	while(lg2){
+		m<<=1,lg2--;
+		Wn[lg2][0]=1LL*Wn[lg2+1][0]*Wn[lg2+1][0]%Mod;
+		Wn[lg2][1]=1LL*Wn[lg2+1][1]*Wn[lg2+1][1]%Mod;
+	}
+}
+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=0;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;
+				f[k]=(t1+t2)%Mod,f[k+i]=(t1-t2)%Mod;
+				cur=1LL*cur*w%Mod;
+			}
+		}
+	}
+	if(!type){
+		int div=quick_pow(n,Mod-2);
+		_for(i,0,n)f[i]=(1LL*f[i]*div%Mod+Mod)%Mod;
+	}
+}
+int f[3][MAXN],temp[5][MAXN];
+void Add(int n){
+	_rep(i,0,n)
+	f[0][i]=f[1][i],f[1][i]=f[2][i];
+	f[2][0]=1;
+	_rep(i,1,n)
+	f[2][i]=((f[1][i]+f[1][i-1])%Mod+f[0][i-1])%Mod;
+}
+void Mul(int n){
+	int _n=build(n<<1);
+	_for(i,0,3)NTT(f[i],_n,true);
+	_for(i,0,_n){
+		temp[0][i]=1LL*f[1][i]*f[1][i]%Mod,temp[1][i]=1LL*f[0][i]*f[0][i]%Mod;
+		temp[2][i]=1LL*f[2][i]*f[1][i]%Mod,temp[3][i]=1LL*f[1][i]*f[0][i]%Mod;
+		temp[4][i]=1LL*f[2][i]*f[2][i]%Mod;
+	}
+	_for(i,0,5)NTT(temp[i],_n,false);
+	f[0][0]=temp[0][0],f[1][0]=temp[2][0],f[2][0]=temp[4][0];
+	_rep(i,1,n){
+		f[0][i]=(temp[0][i]+temp[1][i-1])%Mod;
+		f[1][i]=(temp[2][i]+temp[3][i-1])%Mod;
+		f[2][i]=(temp[4][i]+temp[0][i-1])%Mod;
+	}
+	_for(i,n+1,_n)f[0][i]=f[1][i]=f[2][i]=0;
+}
+int main()
+{
+	init();
+	int n=read_int(),k=read_int(),pos=30;
+	while(n<(1<<pos))pos--;
+	f[1][0]=f[2][0]=f[2][1]=1;
+	while(pos--){
+		Mul(k);
+		if(n&(1<<pos))Add(k);
+	}
+	_rep(i,1,k)space(f[2][i]);
+	return 0;
+}
+</code>
+</hidden>
 ===== 分治 FFT =====
@@ 行 90: / 行 227: @@
 	_for(i,0,n)
 	space(f[i]);
+	return 0;
+}
+</code>
+</hidden>
+==== 算法练习 ====
+[[https://www.luogu.com.cn/problem/CF553E|CF553E]]
+=== 题意 ===
+  - 给定一张 $n$ 个点 $m$ 条边的无重边无自环的有向图，你要从 $1$ 号点到 $n$ 号点去
+  - 如果你在 $t$ 时刻之后才到达 $n$ 号点，你要交 $x$ 元的罚款
+  - 经过每条边需要支付 $w_e$ 费用，且经过该边消耗 $k$ 个单位时间的概率为 $p_{e,k}(1\le k\le t)$
+求最优策略下到达点 $n$ 需要花费的最小费用的期望值，数据范围 $n\le 50,m\le 100,t\le 2\times 10^4$。
+=== 题解 ===
+设 $\text{dp}(i,j)$ 表示走到点 $i$ 且已经花费 $j$ 个单位时间，达到点 $n$ 还需要花费的最小费用的期望值。
+$$\text{dp}(u_e,j)=\min \left(w_e+\sum_{k=1}^tp_{e,k}\text{dp}(v_e,j+k)\right)$$
+边界条件为 $\text{dp}(u_e,j)=x+\text{dis}(u,n)(j\gt t),\text{dp}(n,j)=0(j\le t)$。
+设 $g(e,j)=\sum_{k=1}^tp_{e,k}\text{dp}(v_e,j+k)$，于是有 $\text{dp}(u_e,j)=\min (w_e+g(e,j))$。
+考虑分治 $\text{FFT}$，先计算出 $\text{dp}(u_e,\text{mid}+1\sim \text{rig})$，在计算他们对 $g(e,\text{lef}\sim \text{mid})$ 的贡献，再计算出 $\text{dp}(u_e,\text{lef}\sim \text{mid})$。
+其中，$j\in [\text{lef},\text{mid}]$，而 $j+k$ 需要不遗漏地覆盖 $[\text{mid}+1,\text{rig}]$，于是 $k\in [1,\text{rig}-\text{lef}]$。
+记 $a_i=p_{e,\text{rig}-\text{lef}-i},b_i=\text{dp}(v_e,i+\text{mid}+1)$，于是
+$$g(e,j)\gets \sum p_{e,k}\text{dp}(v_e,j+k)=\sum a_{\text{rig}-\text{lef}-k}b_{k+j-\text{mid}-1}$$
+递归边界为 $\text{lef}=\text{rig}$，此时用 $g(e,j)$ 更新 $\text{dp}(u_e,j)$。由于没有必要通过分治计算出 $\text{dp}(u_e,t+1\sim 2t)$，故可以分治前处理。
+总时间复杂度 $O(mt\log^2 t)$。
+<hidden 查看代码>
+<code cpp>
+const int MAXN=55,MAXM=105,MAXT=2e4+5;
+const double pi=acos(-1.0),Inf=1e9;
+struct complex{
+	double x,y;
+	complex(double x=0.0,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[MAXT<<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;
+}
+complex t1[MAXT<<3],t2[MAXT<<3];
+struct Edge{
+	int u,v;
+	double w;
+}edge[MAXM];
+int edge_cnt;
+double dp[MAXN][MAXT<<1],dp2[MAXM][MAXT<<1],dis[MAXN][MAXN],p[MAXM][MAXT<<1];
+void cal(int lef,int rig){
+	int mid=lef+rig>>1;
+	_for(k,0,edge_cnt){
+		int n1=rig-lef,n2=rig-mid,n=build(n1+n2-2);
+		_for(i,0,n1)t1[i].x=p[k][n1-i],t1[i].y=0.0;_for(i,n1,n)t1[i].x=t1[i].y=0.0;
+		_for(i,0,n2)t2[i].x=dp[edge[k].v][i+mid+1],t2[i].y=0.0;_for(i,n2,n)t2[i].x=t2[i].y=0.0;
+		FFT(t1,n,1);FFT(t2,n,1);
+		_for(i,0,n)t1[i]=t1[i]*t2[i];
+		FFT(t1,n,-1);
+		_rep(i,lef,mid)dp2[k][i]+=t1[i+rig-lef-mid-1].x;
+	}
+}
+void cdq(int lef,int rig){
+	int mid=lef+rig>>1;
+	if(lef==rig){
+		_for(i,0,edge_cnt)
+		dp[edge[i].u][mid]=min(dp[edge[i].u][mid],dp2[i][mid]+edge[i].w);
+		return;
+	}
+	cdq(mid+1,rig);
+	cal(lef,rig);
+	cdq(lef,mid);
+}
+int main()
+{
+	int n=read_int(),m=edge_cnt=read_int(),t=read_int();
+	double x=read_int();
+	_rep(i,1,n)_rep(j,1,n)dis[i][j]=Inf;
+	_rep(i,1,n)dis[i][i]=0.0;
+	_for(i,0,m){
+		edge[i].u=read_int(),edge[i].v=read_int(),edge[i].w=read_int();
+		_rep(j,1,t)p[i][j]=read_int()/1e5;
+		dis[edge[i].u][edge[i].v]=edge[i].w;
+	}
+	_rep(k,1,n)_rep(i,1,n)_rep(j,1,n)
+	dis[i][j]=min(dis[i][j],dis[i][k]+dis[k][j]);
+	_rep(i,0,t)dp[n][i]=0;_rep(i,t+1,t<<1)dp[n][i]=x;
+	_for(i,1,n){
+		_rep(j,0,t)dp[i][j]=Inf;
+		_rep(j,t+1,t<<1)dp[i][j]=x+dis[i][n];
+	}
+	cal(0,t<<1);
+	cdq(0,t);
+	printf("%.10lf",dp[1][0]);
 	return 0;
 }
@@ 行 126: / 行 396: @@
 //递归版
 int temp[MAXN<<2];
-void ployinv(int *f,int *g,int n){
+void polyinv(int *f,int *g,int n){
 	if(n==1)
 	return g[0]=quick_pow(f[0],Mod-2),void();
@@ 行 139: / 行 409: @@
 //递推版
 int temp[MAXN<<2];
-void ployinv(int *f,int *g,int n){
+void polyinv(int *f,int *g,int n){
 	g[0]=quick_pow(f[0],Mod-2);
 	int n1=2,n2=4,pos=2;
@@ 行 197: / 行 467: @@
 }
 int temp[MAXN<<2],inv_f[MAXN<<2];
-void ploysqrt(int *f,int *g,int n){
+void polysqrt(int *f,int *g,int n){
 	f[0]=quick_pow(3,bsgs(3,g[0])/2);
 	int n1=2,n2=4,pos=2,inv2=quick_pow(2,Mod-2);
@@ 行 236: / 行 506: @@
 <code cpp>
 int inv_f[MAXN<<2];
-void ployln(int *f,int n){
+void polyln(int *f,int n){
 	mem(inv_f,0);
 	ployinv(f,inv_f,n);
@@ 行 286: / 行 556: @@
 考虑牛顿迭代法，设 $F(x)\equiv \exp f(x)\pmod {x^n}$，于是有 $g(F(x))\equiv \ln F(x)-f(x)\equiv 0\pmod {x^n}$。
-$$F(x)\equiv F_0(x)-\frac {g(F_0(x))}{g^{\prime}(F_0(x))}\equiv F_0(x)-\frac {\ln F_0(x)-f(x)}{\frac 1{F_0(x)}}\equiv F_0(x)\left(1+f(x)-F_0(x)\right)\pmod {x^n}$$
+$$F(x)\equiv F_0(x)-\frac {g(F_0(x))}{g^{\prime}(F_0(x))}\equiv F_0(x)-\frac {\ln F_0(x)-f(x)}{\frac 1{F_0(x)}}\equiv F_0(x)\left(1+f(x)-\ln F_0(x)\right)\pmod {x^n}$$
 <code cpp>
 int ln_g[MAXN<<2];
-void ployexp(int *f,int *g,int n){
+void polyexp(int *f,int *g,int n){
 	g[0]=1;
 	int n1=2,n2=4,pos=2;
@@ 行 332: / 行 602: @@
 <code cpp>
 int ln_f[MAXN<<2];
-void ploypow(int *f,int n,int k1,int k2){
+void polypow(int *f,int n,int k1,int k2){
 	LL pos=0,posv;
 	while(!f[pos]&&pos<n)pos++;
@@ 行 360: / 行 630: @@
 ==== 算法实现 ====
-构造函数 $f^{R}(x)=x^{\text{deg}(f)}f(\frac 1x)$，易知 $f^{R}(x)$ 与 $f(x)$ 系数恰好颠倒，可以 $O(n)$ 相互转化。
+构造函数 $f^{R}(x)=x^{\text{deg}(f)}f(\frac 1x)$，易知 $f^{R}(x)$ 与 $f(x)$ 系数恰好反转，可以 $O(n)$ 相互转化。
 根据已知，有
@@ 行 378: / 行 648: @@
 <code cpp>
 int temp1[MAXN<<2],temp2[MAXN<<2];
-void ploydiv(int *f,int *g,int *q,int *r,int n,int m){
+void polydiv(int *f,int *g,int *q,int *r,int n,int m){
 	_for(i,0,n)temp1[i]=f[n-1-i];
 	_for(i,0,m)temp2[i]=g[m-1-i];
@@ 行 399: / 行 669: @@
 }
 </code>
+===== 多项式板子汇总 =====
+<hidden 查看代码>
+<code cpp>
+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(const int *f,int _n,const int *g,int _m,int *q,int *r){
+		static int temp[MAXN<<2];
+		if(_n<_m){
+			_rep(i,0,n)r[i]=f[i];
+			return;
+		}
+		_rep(i,0,_m)temp[i]=g[_m-i];
+		Inv(temp,q,_n-_m+1);
+		_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)r[i]=(f[i]+Mod-r[i])%Mod;
+	}
+	void Ln(const int *f,int *g,int _n){
+		static int temp[MAXN<<2];
+		Inv(f,g,_n);
+		_for(i,1,_n)temp[i-1]=1LL*f[i]*i%Mod;
+		temp[_n-1]=0;
+		Mul(g,_n,temp,_n-1,_n);
+		for(int i=_n-1;i;i--)g[i]=1LL*g[i-1]*quick_pow(i,Mod-2)%Mod;
+		g[0]=0;
+	}
+	void Exp(const int *f,int *g,int _n){
+		static int temp[MAXN<<2];
+		if(_n==1)return g[0]=1,void();
+		Exp(f,g,(_n+1)>>1);
+		_for(i,(_n+1)>>1,_n)g[i]=0;
+		Ln(g,temp,_n);
+		temp[0]=(1+f[0]-temp[0])%Mod;
+		_for(i,1,_n)temp[i]=(f[i]-temp[i])%Mod;
+		Mul(g,(_n+1)>>1,temp,_n,_n);
+	}
+	void Pow(const int *f,int *g,int _n,int k1,int k2){
+		static int temp[MAXN<<2];
+		int pos=0,posv;
+		while(!f[pos]&&pos<_n)pos++;
+		if(1LL*pos*k2>=_n){
+			_for(i,0,_n)g[i]=0;
+			return;
+		}
+		posv=quick_pow(f[pos],Mod-2);
+		_for(i,pos,_n)g[i-pos]=1LL*f[i]*posv%Mod;
+		_for(i,_n-pos,_n)g[i]=0;
+		Ln(g,temp,_n);
+		_for(i,0,_n)temp[i]=1LL*temp[i]*k1%Mod;
+		Exp(temp,g,_n);
+		pos=pos*k2,posv=quick_pow(posv,1LL*k2*(Mod-2)%(Mod-1));
+		for(int i=_n-1;i>=pos;i--)g[i]=1LL*g[i-pos]*posv%Mod;
+		_for(i,0,pos)g[i]=0;
+	}
+	void Sqrt(const int *f,int *g,int _n){
+		static int temp1[MAXN<<2],temp2[MAXN<<2];
+		if(_n==1)return g[0]=quick_pow(G,bsgs(3,f[0])/2),void();
+		Sqrt(f,g,(_n+1)>>1);
+		_for(i,(_n+1)>>1,_n)g[i]=0;
+		_for(i,0,_n)temp1[i]=f[i];
+		Inv(g,temp2,_n);
+		Mul(temp1,_n,temp2,_n);
+		int div2=quick_pow(2,Mod-2);
+		_for(i,0,_n)g[i]=1LL*(g[i]+temp1[i])*div2%Mod;
+	}
+}
+</code>
+</hidden>

CVBB ACM Team

用户工具

站点工具

差别

页面工具