BZOJ 3583: 杰杰的女性朋友

答案就是 \(\sum \limits_{i=0}^d(OI^T)^i[u][v]\)
可以在 \(O\) 中 \(v \to 0\) 和 \(0 \to 0\) 和 \(I\) 中的 \(0 \to 0\) 都设为 \(1\)，\(d + 1\) 次幂的 \(u \to 0\) 就是前 \(d\) 个矩阵 \(u \to v\) 的方案数了
这样做复杂度是 \(O(mn^3\log d)\) 的，没法过
但是注意到 \((OI^T)^{d+1}= O(I^TO)^d I^T\)
这样复杂度就是 \(O(mk^3\log d)\) 的，问题解决。
矩阵乘法那里会做很多次取模，可以先存为long long，大于 \(2^60\) 再取模，速度快了一倍多。。

#include <bits/stdc++.h>
typedef long long ll;

const int N = 1007;
const int MOD = 1000000007;

inline void M(int &x) {
	if (x >= MOD) x -= MOD;
	if (x < 0) x += MOD;
}

struct Mat {
	int r, c;
	std::vector< std::vector<int> > mat;
	inline void init(int _r, int _c, int x = 0) {
		r = _r; c = _c;
		mat.resize(r + 1);
		for (int i = 0; i <= r; i++) {
			mat[i].resize(c + 1);
			for (int j = 0; j <= c; j++) {
				if (i == j) mat[i][j] = x;
				else mat[i][j] = 0;
			}
		}
	}
	inline Mat operator * (const Mat &p) const {
		static Mat res;
		res.init(r, p.c);
		for (int i = 0; i <= r; i++)
			for (int j = 0; j <= p.c; j++) {
				ll x = 0;
				for (int k = 0; k <= c; k++) {
					x += 1LL * mat[i][k] * p.mat[k][j];
					if (x >= (1LL << 60))
						x %= MOD;
				}
				res.mat[i][j] = x % MOD;
			}
		return res;
	}
} O, I;

int n, k;

inline Mat qp(Mat a, int b) {
	static Mat res;
	res.init(a.r, a.c, 1);
	while (b) {
		if (b & 1) res = res * a;
		a = a * a;
		b >>= 1;
	}
	return res;
}

int main() {
	scanf("%d%d", &n, &k);
	O.init(n, k); I.init(k, n);
	for (int i = 1, x; i <= n; i++) {
		for (int j = 1; j <= k; j++) {
			scanf("%d", &x);
			M(O.mat[i][j] = x);
		}
		for (int j = 1; j <= k; j++) {
			scanf("%d", &x);
			M(I.mat[j][i] = x);
		}
	}
	int m;
	scanf("%d", &m);
	while (m--) {
		int u, v, d;
		scanf("%d%d%d", &u, &v, &d);
		O.mat[v][0] = O.mat[0][0] = I.mat[0][0] = 1;
		static Mat res;
		res = (I * O);
		res = O * qp(res, d);
		int ans = 1LL * res.mat[u][0] * I.mat[0][0] % MOD;
		printf("%d\n", ans);
		O.mat[v][0] = O.mat[0][0] = I.mat[0][0] = 0;
	}
	return 0;
}