http://210.33.19.103/problem/2174
很显然是矩阵快速幂的题，设有in和ou矩阵，设in矩阵的转置为in'
显然可以直接暴力求出任意两点间走一步路径条数，然后求其d次幂，但是这样子复杂度不对
注意到设任意两点间走一步路径条数的矩阵为A，那么A=ou*in'，A^d=(ou*in')^d=ou*(in'*ou)^(d-1)*in'（当然d==1时直接特判掉）
in'*ou就是一个K*K的矩阵了，复杂度很对的样子
好像还有点不对。。。题意要求的是前缀和，问题也不大，
设P=in'*ou
题目要求的是ou*P^0*in'+ou*P^1*in'+..+ou*P^(d-1)*in'=ou*(P^0+P^1+..+P^(d-1))*in'
这个P^0+P^1+..+P^(d-1)，搞一个矩阵套矩阵可以算（的确是可以用的23333）
搞两个矩阵X,Y，X=(I,I)，Y是一个2*2的矩阵，第一行为(P,P)，第二行为(0,I)（I表示K阶单位矩阵，0表示K阶零矩阵）
那么求Z=X*Y^(d-1)，Z的第一行第二列的值即为P^0+P^1+..+P^(d-1)
（以上与求值的前缀和的方法是一样的）
明明是思路很清晰的题啊，然而我又是一天被续掉了。。。
鬼畜卡常题啊。。。
此题最好，最正确的方法是不用任何结构体，手写所有矩阵相乘
绝对不能用vector来存不同大小的矩阵!不然会T飞
绝对不要尝试在这类地方用模板类来实现嵌套矩阵！搞不出来的
嵌套矩阵很难写，常数又大，后来我干脆把这个嵌套的矩阵拆开（2*2的矩阵,其中每个元素是K*K的矩阵，就拆成(2*K)*(2*K)的矩阵），目前可以确认是对的
O3卡过代码：
1 #pragma GCC optimize(3) 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 #include<vector> 6 #include<cassert> 7 using namespace std; 8 #define fi first 9 #define se second 10 #define mp make_pair 11 #define pb push_back 12 typedef long long ll; 13 typedef unsigned long long ull; 14 typedef pair<int,int> pii; 15 //#define assert(x) 16 const int md=1000000007; 17 struct M_NN 18 { 19 int d[1010][1010];int x,y; 20 void resize(int a,int b) 21 { 22 x=a;y=b; 23 memset(d,0,sizeof(d)); 24 //d.clear(); 25 //d.resize(x+1); 26 //for(int i=1;i<=x;i++) d[i].resize(y+1); 27 } 28 }; 29 struct M_NK 30 { 31 int d[1010][42];int x,y; 32 void resize(int a,int b) 33 { 34 x=a;y=b; 35 memset(d,0,sizeof(d)); 36 //d.clear(); 37 //d.resize(x+1); 38 //for(int i=1;i<=x;i++) d[i].resize(y+1); 39 } 40 }; 41 struct M_KN 42 { 43 int d[42][1010];int x,y; 44 void resize(int a,int b) 45 { 46 x=a;y=b; 47 memset(d,0,sizeof(d)); 48 //d.clear(); 49 //d.resize(x+1); 50 //for(int i=1;i<=x;i++) d[i].resize(y+1); 51 } 52 }; 53 struct M_KK 54 { 55 int d[42][42];int x,y; 56 void resize(int a,int b) 57 { 58 x=a;y=b; 59 memset(d,0,sizeof(d)); 60 //d.clear(); 61 //d.resize(x+1); 62 //for(int i=1;i<=x;i++) d[i].resize(y+1); 63 } 64 }; 65 /* 66 void init_0(M &p) 67 { 68 p.resize(p.x,p.y); 69 } 70 */ 71 M_NK operator*(const M_NK &a,const M_KK &b) 72 { 73 assert(a.y==b.x); 74 M_NK c;c.resize(a.x,b.y); 75 int i,j,k; 76 for(i=1;i<=a.x;i++) 77 for(j=1;j<=b.x;j++) 78 for(k=1;k<=b.y;k++) 79 c.d[i][k]=(c.d[i][k]+ll(a.d[i][j])*b.d[j][k]%md)%md; 80 return c; 81 } 82 83 M_NN operator*(const M_NK &a,const M_KN &b) 84 { 85 assert(a.y==b.x); 86 M_NN c;c.resize(a.x,b.y); 87 int i,j,k; 88 for(i=1;i<=a.x;i++) 89 for(j=1;j<=b.x;j++) 90 for(k=1;k<=b.y;k++) 91 c.d[i][k]=(c.d[i][k]+ll(a.d[i][j])*b.d[j][k]%md)%md; 92 return c; 93 } 94 M_KK operator*(const M_KN &a,const M_NK &b) 95 { 96 assert(a.y==b.x); 97 M_KK c;c.resize(a.x,b.y); 98 int i,j,k; 99 for(i=1;i<=a.x;i++) 100 for(j=1;j<=b.x;j++) 101 for(k=1;k<=b.y;k++) 102 c.d[i][k]=(c.d[i][k]+ll(a.d[i][j])*b.d[j][k]%md)%md; 103 return c; 104 } 105 M_KK operator*(const M_KK &a,const M_KK &b) 106 { 107 assert(a.y==b.x); 108 M_KK c;c.resize(a.x,b.y); 109 int i,j,k; 110 for(i=1;i<=a.x;i++) 111 for(j=1;j<=b.x;j++) 112 for(k=1;k<=b.y;k++) 113 c.d[i][k]=(c.d[i][k]+ll(a.d[i][j])*b.d[j][k]%md)%md; 114 return c; 115 } 116 M_KK o,tt2,t1,t2; 117 M_NK ou1;M_KN in1; 118 M_KK p; 119 int in[1010][22],ou[1010][22]; 120 M_KK poww(const M_KK &a,int b) 121 { 122 M_KK ans=o,base=a; 123 assert(a.x==a.y); 124 for(;b;base=base*base,b>>=1) 125 if(b&1) 126 ans=ans*base; 127 return ans; 128 } 129 130 int n,K,m; 131 int solve(int u,int v,int d) 132 { 133 if(d==0) return u==v; 134 M_KK tt=t1*poww(t2,d-1); 135 int i,j; 136 for(i=1;i<=K;i++) 137 { 138 for(j=1;j<=K;j++) 139 { 140 tt2.d[i][j]=tt.d[i][j+K]; 141 } 142 } 143 M_NK t2=ou1*tt2;int an=u==v; 144 for(i=1;i<=t2.y;i++) 145 an=(an+ll(t2.d[u][i])*in1.d[i][v]%md)%md; 146 return an; 147 } 148 int main() 149 { 150 //freopen("/tmp/3583/6.in","r",stdin); 151 //freopen("/tmp/3583/6.ans","w",stdout); 152 int i,j,u,v,d; 153 scanf("%d%d",&n,&K); 154 for(i=1;i<=n;i++) 155 { 156 for(j=1;j<=K;j++) 157 { 158 scanf("%d",&ou[i][j]); 159 } 160 for(j=1;j<=K;j++) 161 { 162 scanf("%d",&in[i][j]); 163 } 164 } 165 ou1.resize(n,K); 166 for(i=1;i<=n;i++) 167 { 168 for(j=1;j<=K;j++) 169 { 170 ou1.d[i][j]=ou[i][j]; 171 } 172 } 173 in1.resize(K,n); 174 for(i=1;i<=K;i++) 175 { 176 for(j=1;j<=n;j++) 177 { 178 in1.d[i][j]=in[j][i]; 179 } 180 } 181 p=in1*ou1; 182 t1.resize(K,2*K); 183 for(i=1;i<=K;i++) t1.d[i][i]=t1.d[i][i+K]=1; 184 t2.resize(2*K,2*K); 185 for(i=1;i<=K;i++) 186 { 187 for(j=1;j<=K;j++) 188 { 189 t2.d[i][j]=t2.d[i][j+K]=p.d[i][j]; 190 } 191 } 192 for(i=1;i<=K;i++) t2.d[i+K][i+K]=1; 193 tt2.resize(K,K); 194 o.resize(2*K,2*K); 195 for(i=1;i<=2*K;i++) o.d[i][i]=1; 196 scanf("%d",&m); 197 while(m--) 198 { 199 scanf("%d%d%d",&u,&v,&d); 200 printf("%d ",solve(u,v,d)); 201 //return 0; 202 } 203 return 0; 204 }
一份有分，可以跑的代码：
1 #pragma GCC optimize(3) 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 #include<vector> 6 //#include<cassert> 7 using namespace std; 8 #define fi first 9 #define se second 10 #define mp make_pair 11 #define pb push_back 12 typedef long long ll; 13 typedef unsigned long long ull; 14 typedef pair<int,int> pii; 15 #define vector myvec 16 template<typename T> 17 struct myvec 18 { 19 T *p;int sz; 20 myvec():p(0),sz(0){} 21 ~myvec(){delete[] p;} 22 myvec &operator=(const myvec &b) 23 { 24 delete[] p; 25 sz=b.sz; 26 p=new T[sz]; 27 for(T *i1=p,*i2=b.p;i1!=p+sz;++i1,++i2) *i1=*i2; 28 return *this; 29 } 30 myvec(const myvec &b):p(0),sz(0){*this=b;} 31 void resize(int d){sz=d;delete[] p;p=new T[d]();} 32 T &operator[](int d){return p[d];} 33 const T &operator[](int d)const{return p[d];} 34 }; 35 const ll md=1000000007; 36 struct M1 37 { 38 vector<vector<ll> > d;int x,y; 39 M1(); 40 void resize(int xx,int yy) 41 { 42 x=xx;y=yy; 43 d.resize(x+1); 44 for(int i=1;i<=x;i++) d[i].resize(y+1); 45 } 46 }; 47 struct M2 48 { 49 M1 d[3][3];int x,y; 50 M2(); 51 void resize_sub(int xx,int yy) 52 { 53 int i,j; 54 for(i=1;i<=x;i++) 55 for(j=1;j<=y;j++) 56 d[i][j].resize(xx,yy); 57 } 58 }; 59 void init_0(M1 &p) 60 { 61 //p.d.clear(); 62 p.resize(p.x,p.y); 63 } 64 void init_0(M2 &p) 65 { 66 int i,j; 67 for(i=1;i<=2;i++) 68 for(j=1;j<=2;j++) 69 init_0(p.d[i][j]); 70 } 71 void init_1(M1 &p) 72 { 73 init_0(p);//assert(p.x==p.y); 74 for(int i=1;i<=p.x;i++) p.d[i][i]=1; 75 } 76 void init_1(M2 &p) 77 { 78 init_0(p);//assert(p.x==p.y); 79 for(int i=1;i<=p.x;i++) init_1(p.d[i][i]); 80 } 81 M1::M1(){x=y=0;init_0(*this);} 82 M2::M2(){x=y=0;init_0(*this);} 83 M1 operator+(const M1 &a,const M1 &b) 84 { 85 //assert(a.x==b.x&&a.y==b.y); 86 M1 c;c.resize(a.x,a.y); 87 int i,j; 88 for(i=1;i<=a.x;i++) 89 for(j=1;j<=a.y;j++) 90 c.d[i][j]=(a.d[i][j]+b.d[i][j])%md; 91 return c; 92 } 93 M1 operator*(const M1 &a,const M1 &b) 94 { 95 //assert(a.y==b.x); 96 M1 c;c.resize(a.x,b.y); 97 int i,j,k; 98 for(i=1;i<=a.x;i++) 99 for(j=1;j<=b.x;j++) 100 for(k=1;k<=b.y;k++) 101 c.d[i][k]=(c.d[i][k]+a.d[i][j]*b.d[j][k])%md; 102 return c; 103 } 104 M2 operator*(const M2 &a,const M2 &b) 105 { 106 //assert(a.y==b.x); 107 M2 c;c.x=a.x;c.y=b.y;c.resize_sub(a.d[1][1].x,a.d[1][1].y); 108 int i,j,k; 109 for(i=1;i<=a.x;i++) 110 for(j=1;j<=b.x;j++) 111 for(k=1;k<=b.y;k++) 112 c.d[i][k]=(c.d[i][k]+a.d[i][j]*b.d[j][k]); 113 return c; 114 } 115 M1 ou1,in1,p; 116 int ou[1010][22],in[1010][22]; 117 M2 poww(const M2 &a,ll b) 118 { 119 //printf("tt%d %d ",a.x,a.y); 120 M2 ans,base=a;/*assert(a.x==a.y);*/ans.x=ans.y=a.x; 121 ans.resize_sub(a.d[1][1].x,a.d[1][1].y);init_1(ans); 122 for(;b;base=base*base,b>>=1) 123 if(b&1) 124 ans=ans*base; 125 return ans; 126 } 127 int n,K,m; 128 129 /* 130 void out(const M1 &a) 131 { 132 puts("start"); 133 int i,j; 134 for(i=1;i<=a.x;i++) 135 { 136 for(j=1;j<=a.y;j++) 137 { 138 printf("%lld ",a.d[i][j]); 139 } 140 puts(""); 141 } 142 puts("end"); 143 fflush(stdout); 144 } 145 146 147 void out(const M2 &p) 148 { 149 puts("st"); 150 printf("xy%d %d ",p.x,p.y); 151 for(int i=1;i<=p.x;i++) 152 for(int j=1;j<=p.y;j++) 153 out(p.d[i][j]),puts(""); 154 puts("ed"); 155 } 156 */ 157 M2 t1,t2; 158 int solve(int u,int v,int d) 159 { 160 if(d==0) return u==v; 161 return ((u==v)+(ou1*((t1*poww(t2,d-1)).d[1][2])*in1).d[u][v])%md; 162 } 163 int main() 164 { 165 int i,j,u,v,d; 166 scanf("%d%d",&n,&K); 167 for(i=1;i<=n;i++) 168 { 169 for(j=1;j<=K;j++) 170 { 171 scanf("%d",&ou[i][j]); 172 } 173 for(j=1;j<=K;j++) 174 { 175 scanf("%d",&in[i][j]); 176 } 177 } 178 ou1.resize(n,K); 179 for(i=1;i<=n;i++) 180 { 181 for(j=1;j<=K;j++) 182 { 183 ou1.d[i][j]=ou[i][j]; 184 } 185 } 186 in1.resize(K,n); 187 for(i=1;i<=K;i++) 188 { 189 for(j=1;j<=n;j++) 190 { 191 in1.d[i][j]=in[j][i]; 192 } 193 } 194 p=in1*ou1; 195 t1.x=1;t1.y=2;t1.resize_sub(K,K); 196 //t1.d[1][1]=M1(1,K,K); 197 init_1(t1.d[1][1]);init_1(t1.d[1][2]); 198 t2.x=t2.y=2;t2.resize_sub(K,K); 199 t2.d[1][1]=t2.d[1][2]=p; 200 //t2.d[2][2]=M1(1,K,K); 201 init_1(t2.d[2][2]); 202 //out((t1*poww(t2,0)).d[1][2]); 203 scanf("%d",&m); 204 while(m--) 205 { 206 scanf("%d%d%d",&u,&v,&d); 207 printf("%d ",solve(u,v,d)); 208 } 209 return 0; 210 }