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})；状态转移方程：
其中 (ext(i,j)) 表示第 (i) 天城市 (j) 通过美食节额外获得的愉悦值。
暴力 DP 肯定是不行的。这个 (T) 这么大，一脸矩阵快速幂的样子。
考虑用矩阵的形式将状态转移方程写出来。我们重定义矩阵乘法，将 ( imes) 改为 (+)，将 (+) 改为 (max)，学过 DDP 的都知道这样满足原矩阵乘法的一切性质（我没学过 DDP 都知道）。需要注意的是，这种乘法意义下的单位元是
具体广义乘法这套理论等到我学 DDP 再研究吧。于是：
其中 (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;
}