差别

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

--- 2020-2021:teams:legal_string:jxm2001:多项式_3 [2020/08/19 22:22]
jxm2001
+++ 2020-2021:teams:legal_string:jxm2001:多项式_3 [2020/08/25 12:43] (当前版本)
jxm2001
@@ 行 227: / 行 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;
 }
@@ 行 543: / 行 676: @@
 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;
 		}
 	}
@@ 行 566: / 行 700: @@
 		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;
@@ 行 604: / 行 736: @@
 	void Div(const int *f,int _n,const int *g,int _m,int *q,int *r){
 		static int temp[MAXN<<2];
-		_for(i,0,_m)temp[i]=g[_m-1-i];
+		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);
-		_for(i,0,_n)temp[i]=f[_n-1-i];
+		_rep(i,0,_n)temp[i]=f[_n-i];
-		Mul(q,_n-_m+1,temp,_n,_n-_m+1);
+		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]);
-		_for(i,0,_m)r[i]=g[i];_rep(i,0,_n-_m)temp[i]=q[i];
+		int __m=min(_n-_m+1,_m);
-		Mul(r,_m,temp,_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;
 	}

CVBB ACM Team

用户工具

站点工具

差别

页面工具