洛谷 P6772

NOI 签到题我都不会哦~

洛谷题目页面传送门

题意见洛谷。（以下用 (K) 表示 (k) 以避免与迭代变量冲突）

首先有个很显然的 DP：(dp_{i,j}) 表示第 (i) 天小 W 在城市 (j) 的情况下获得的最大愉悦值之和。边界：(dp_{0,i}=egin{cases}c_1&i=1\-infty&i eq1end{cases})；目标：(dp_{T,1})；状态转移方程：

[dp_{i,j}=max_{(k,j,len)in E}{dp_{i-len,k}+c_j+ext(i,j)} ]

其中 (ext(i,j)) 表示第 (i) 天城市 (j) 通过美食节额外获得的愉悦值。

暴力 DP 肯定是不行的。这个 (T) 这么大，一脸矩阵快速幂的样子。

考虑用矩阵的形式将状态转移方程写出来。我们重定义矩阵乘法，将 ( imes) 改为 (+)，将 (+) 改为 (max)，学过 DDP 的都知道这样满足原矩阵乘法的一切性质（我没学过 DDP 都知道）。需要注意的是，这种乘法意义下的单位元是

[I=egin{bmatrix}0&-infty&-infty &cdots&-infty\-infty&0&-infty&cdots&-infty\-infty&-infty&0&cdots&-infty\vdots&vdots&vdots&ddots&-infty\-infty&-infty&-infty&-infty&0end{bmatrix} ]

具体广义乘法这套理论等到我学 DDP 再研究吧。于是：

[egin{bmatrix}dp_{i,1}&dp_{i,2}&cdots&dp_{i,n}&cdots&dp_{i-4,n}end{bmatrix}Delta_i=egin{bmatrix}dp_{i+1,1}&dp_{i+1,2}&cdots&dp_{i+1,n}&cdots&dp_{i-3,n}end{bmatrix} ]

其中 (Delta) 是啥就不详细说了，大概就是前 (n) 列是转移方程里的 (c_j+ext(i,j))（如果转移不到就是 (-infty)），后 (4n) 列是一些 (0) 和 (-infty) 表示移位。

注意到大多数时候 (ext(i,j)=0)，只存在 (mathrm O(K)) 个 (i) 使得不等于零。而所有 (ext) 为零的时候的 (Delta) 是相等的，于是考虑按有美食节的时间分段，每段快速幂，有美食节的时间暴力修改矩阵单独转移。

这样复杂度为 (mathrm O!left(KN^3log T ight))，其中 (N=5n)。显然过不去，瓶颈在于快速幂。

如何优化呢？注意到这里是同底数幂，想到进制光速幂，即设一个 (x) 进制，预处理出 (Delta^0,Delta^{x^y},cdots,Delta^{x^y(x-1)})，其中 (yin[0,lceillog_xT ceil])，然后每次算幂的时候将指数 (x) 进制拆分乘即可。时间复杂度 (mathrm O!left(N^3xlog_xT+KN^3log_xT ight))。本来就卡的很紧，似乎并没啥卵用。

但是！这样一来，算幂的时候每次乘法（共 (mathrm O(log_xT)) 次）可以利用矩阵乘法的结合律，将矩阵乘到向量上去，时间复杂度 (mathrm O(N^2))！然鹅直接使用快速幂的话，一次快速幂是一个整体，无法在其过程中乘以向量；光速幂的优点是，它将乘的过程拆开了。

总时间复杂度 (mathrm O!left(N^3xlog_xT+KN^2log_xT ight))。xjb 取个 (x=2) 即可得到 (mathrm O!left(N^3log T+KN^2log T ight)) 的复杂度。由于矩阵乘法常数非常小，可以通过。

这里「xjb 取个 (x=2)」，换句话说就是光速幂的优化的地方不在「光速」，而在用矩阵乘法的结合律降低单次乘法的复杂度，这大概是个矩阵乘法比较惯用的套路吧。

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define pb push_back
#define mp make_pair
#define X first
#define Y second
const int inf=0x3f3f3f3f3f3f3f3f;
const int N=50,T=200,LOG_T=32; 
int n,m,s,t;
int c[N+1];
vector<pair<int,int> > rnei[N+1];
pair<int,pair<int,int> > spc[T+2];
struct matrix{
	int a[5*N][5*N];
	int *operator[](int x){return a[x];}
	matrix(){
		for(int i=0;i<5*N;i++)for(int j=0;j<5*N;j++)a[i][j]=i==j?0:-inf;
	}
	friend matrix operator*(matrix x,matrix y){
		matrix res;
		for(int i=0;i<5*n;i++)for(int j=0;j<5*n;j++)res[i][j]=-inf;
		for(int i=0;i<5*n;i++)for(int j=0;j<5*n;j++)for(int k=0;k<5*n;k++)
			res[i][j]=max(res[i][j],x[i][k]+y[k][j]);
		return res;
	}
	friend vector<int> operator*(vector<int> x,matrix y){
		vector<int> res(5*n,-inf);
		for(int j=0;j<5*n;j++)for(int k=0;k<5*n;k++)
			res[j]=max(res[j],x[k]+y[k][j]);
		return res;
	}
};
matrix pw[LOG_T];
vector<int> dp;
signed main(){
	cin>>n>>m>>s>>t;
	for(int i=1;i<=n;i++)cin>>c[i];
	while(m--){
		int x,y,z;
		cin>>x>>y>>z;
		rnei[y].pb(mp(x,z));
	}
	for(int i=1;i<=t;i++)cin>>spc[i].X>>spc[i].Y.X>>spc[i].Y.Y;
	for(int i=0;i<5*n;i++)for(int j=0;j<5*n;j++)pw[0][i][j]=-inf;
	for(int i=1;i<=n;i++)for(int j=0;j<rnei[i].size();j++)
		pw[0][(rnei[i][j].Y-1)*n+rnei[i][j].X-1][i-1]=c[i];
	for(int i=n;i<5*n;i++)pw[0][i-n][i]=0;
	for(int i=1;i<LOG_T;i++)pw[i]=pw[i-1]*pw[i-1];
	sort(spc+1,spc+t+1);
	spc[++t]=mp(inf,mp(0,0));
	dp.resize(5*n,-inf);dp[0]=c[1];
	for(int i=1;i<=t;i++){
		if(spc[i].X>s){
			int e=s-spc[i-1].X;
			for(int j=0;j<LOG_T;j++)if(e&1<<j)dp=dp*pw[j];
			break;
		}
		int e=spc[i].X-1-spc[i-1].X;
		for(int j=0;j<LOG_T;j++)if(e&1<<j)dp=dp*pw[j];
		matrix delta=pw[0];
		for(int j=0;j<5*n;j++)delta[j][spc[i].Y.X-1]+=spc[i].Y.Y;
		dp=dp*delta;
	}
	if(dp[0]<0)puts("-1");
	else cout<<dp[0];
	return 0;
}