### [题目链接](https://hydro.ac/d/bzoj/p/3583) 居然显示清华集训的题，再是福建省选的题，真是太奇妙了QwQ 设 $f[t][i]$ 表示 $t$ 步到 $i$ 的方案数 $f[t][i]=\sum\limits_{j=1}^{n} (f[t-1][j]*\sum\limits_{l=1}^{k} O_{j,l}*I_{i,l})$ 显然转移可以写成矩阵的形式。 $f[t]=f[t-1]*(O*I)$ $I$ 表示原 $I$ 的转置（我太懒了） 发现 $f[t]=f[0]*(O*I)^t$，$O*I$ 是 $n\times n$ 的矩阵，跑起来很慢，可以转化一下。 $f[t]=f[0]*O*(I*O)^{t-1}*I$ ，这样 $I*O$ 是 $k\times k$ 的，会快很多。 但是我们发现 $f[t]$ 表示的是恰好 $t$ 时刻的，但我们要求的是 $\sum\limits_{i=0}^d f[i]=f[0]*O*(\sum\limits_{i=0}^{d-1}D^{i})*I$，只要求 $\sum\limits_{i=0}^{d-1}D^{i}$ 这个就好了QwQ。 倍增求解就好了QwQ ### CODE ```cpp #include <bits/stdc++.h> using namespace std; int read() { int x=0, f=0; char c=getchar(); while (!isdigit(c)) f|=c=='-', c=getchar(); while (isdigit(c)) x=(x<<3)+(x<<1)+(c^48), c=getchar(); return f ? -x : x; } const int N=1e3+10, p=1e9+7; int n, k, q; struct matrix{ int n, m; vector<vector<int> > a; matrix(int _n=0, int _m=0) { n=_n, m=_m; if (!m) { m=n; a=vector<vector<int> >(n+1, vector<int> (m+1)); for (int i=1; i<=n; i++) a[i][i]=1; } else a=vector<vector<int> >(n+1, vector<int> (m+1)); } vector<int>& operator[](int x) & {return a[x];} const vector<int>& operator[](int x) const& {return a[x];} matrix operator + (const matrix b) const { assert(n==b.n && m==b.m); matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=(a[i][j]+b[i][j])%p; return r; } matrix operator - (const matrix b) const { assert(n==b.n && m==b.m); matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=(a[i][j]-b[i][j]+p)%p; return r; } matrix operator * (const int k) const { matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=1ll*a[i][j]*k%p; return r; } matrix operator * (const matrix b) const { assert(m==b.n); matrix r(n, b.m); for (int i=1; i<=n; i++) for (int k=1; k<=m; k++) for (int j=1; j<=b.m; j++) r[i][j]=(r[i][j]+1ll*a[i][k]*b[k][j])%p; return r; } matrix operator += (const matrix b) {return *this=*this+b;} matrix operator -= (const matrix b) {return *this=*this-b;} matrix operator *= (const int k) {return *this=*this*k;} matrix operator *= (const matrix b) {return *this=*this*b;} bool operator == (const matrix b) const { if (n!=b.n || m!=b.m) return 0; for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) if (a[i][j]!=b[i][j]) return 0; return 1; } bool operator != (const matrix b) const { return !(*this==b); } matrix operator ^ (int k) const { assert(n==m || k<=1); matrix r(m), a=*this; while (k) { if (k&1) r*=a; a*=a; k>>=1; } return r; } matrix operator ^= (int k) {return *this=*this^k;} } f[33], g[33]; matrix calc(int x) { matrix r(k), t(k); for (int i=0; i<32; i++) if (x>>i&1) { r=r+t*g[i]; t*=f[i]; } return r; } int main() { n=read(), k=read(); matrix O(n, k), I(k, n); for (int i=1; i<=n; i++) { for (int j=1; j<=k; j++) O[i][j]=read(); for (int j=1; j<=k; j++) I[j][i]=read(); } matrix QWQ(k); f[0]=g[0]=I*O; for (int i=1; i<32; i++) { f[i]=f[i-1]*f[i-1]; g[i]=g[i-1]*(QWQ+f[i-1]); } q=read(); while (q--) { int x=read(), y=read(), d=read(); if (d==0) printf("%d\n", x==y); else { matrix a(1, n), b(n, 1); a[1][x]=b[y][1]=1; printf("%d\n", (a*O*calc(d-1)*I*b)[1][1]+(x==y)); } } return 0; } ``` Loading... ### [题目链接](https://hydro.ac/d/bzoj/p/3583) 居然显示清华集训的题，再是福建省选的题，真是太奇妙了QwQ 设 $f[t][i]$ 表示 $t$ 步到 $i$ 的方案数 $f[t][i]=\sum\limits_{j=1}^{n} (f[t-1][j]*\sum\limits_{l=1}^{k} O_{j,l}*I_{i,l})$ 显然转移可以写成矩阵的形式。 $f[t]=f[t-1]*(O*I)$ $I$ 表示原 $I$ 的转置（我太懒了） 发现 $f[t]=f[0]*(O*I)^t$，$O*I$ 是 $n\times n$ 的矩阵，跑起来很慢，可以转化一下。 $f[t]=f[0]*O*(I*O)^{t-1}*I$ ，这样 $I*O$ 是 $k\times k$ 的，会快很多。 但是我们发现 $f[t]$ 表示的是恰好 $t$ 时刻的，但我们要求的是 $\sum\limits_{i=0}^d f[i]=f[0]*O*(\sum\limits_{i=0}^{d-1}D^{i})*I$，只要求 $\sum\limits_{i=0}^{d-1}D^{i}$ 这个就好了QwQ。 倍增求解就好了QwQ ### CODE ```cpp #include <bits/stdc++.h> using namespace std; int read() { int x=0, f=0; char c=getchar(); while (!isdigit(c)) f|=c=='-', c=getchar(); while (isdigit(c)) x=(x<<3)+(x<<1)+(c^48), c=getchar(); return f ? -x : x; } const int N=1e3+10, p=1e9+7; int n, k, q; struct matrix{ int n, m; vector<vector<int> > a; matrix(int _n=0, int _m=0) { n=_n, m=_m; if (!m) { m=n; a=vector<vector<int> >(n+1, vector<int> (m+1)); for (int i=1; i<=n; i++) a[i][i]=1; } else a=vector<vector<int> >(n+1, vector<int> (m+1)); } vector<int>& operator[](int x) & {return a[x];} const vector<int>& operator[](int x) const& {return a[x];} matrix operator + (const matrix b) const { assert(n==b.n && m==b.m); matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=(a[i][j]+b[i][j])%p; return r; } matrix operator - (const matrix b) const { assert(n==b.n && m==b.m); matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=(a[i][j]-b[i][j]+p)%p; return r; } matrix operator * (const int k) const { matrix r(n, m); for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) r[i][j]=1ll*a[i][j]*k%p; return r; } matrix operator * (const matrix b) const { assert(m==b.n); matrix r(n, b.m); for (int i=1; i<=n; i++) for (int k=1; k<=m; k++) for (int j=1; j<=b.m; j++) r[i][j]=(r[i][j]+1ll*a[i][k]*b[k][j])%p; return r; } matrix operator += (const matrix b) {return *this=*this+b;} matrix operator -= (const matrix b) {return *this=*this-b;} matrix operator *= (const int k) {return *this=*this*k;} matrix operator *= (const matrix b) {return *this=*this*b;} bool operator == (const matrix b) const { if (n!=b.n || m!=b.m) return 0; for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) if (a[i][j]!=b[i][j]) return 0; return 1; } bool operator != (const matrix b) const { return !(*this==b); } matrix operator ^ (int k) const { assert(n==m || k<=1); matrix r(m), a=*this; while (k) { if (k&1) r*=a; a*=a; k>>=1; } return r; } matrix operator ^= (int k) {return *this=*this^k;} } f[33], g[33]; matrix calc(int x) { matrix r(k), t(k); for (int i=0; i<32; i++) if (x>>i&1) { r=r+t*g[i]; t*=f[i]; } return r; } int main() { n=read(), k=read(); matrix O(n, k), I(k, n); for (int i=1; i<=n; i++) { for (int j=1; j<=k; j++) O[i][j]=read(); for (int j=1; j<=k; j++) I[j][i]=read(); } matrix QWQ(k); f[0]=g[0]=I*O; for (int i=1; i<32; i++) { f[i]=f[i-1]*f[i-1]; g[i]=g[i-1]*(QWQ+f[i-1]); } q=read(); while (q--) { int x=read(), y=read(), d=read(); if (d==0) printf("%d\n", x==y); else { matrix a(1, n), b(n, 1); a[1][x]=b[y][1]=1; printf("%d\n", (a*O*calc(d-1)*I*b)[1][1]+(x==y)); } } return 0; } ``` 最后修改：2025 年 07 月 02 日 © 允许规范转载 赞 如果觉得我的文章对你有用，请随意赞赏