Bzoj3583 杰杰的女性朋友
Submit: 190 Solved: 98
Description
杰杰是魔法界的一名传奇人物。他对魔法具有深刻的洞察力，惊人的领悟力，以及令人叹为观止的创造力。自从他从事魔法竞赛以来，短短几年时间，就已经成为世界公认的实力最强的魔法选手之一。更让人惊叹的是，他几乎没有借助外界力量，完全凭借自己的努力达到了普通人难以企及的高度。在最近的世界魔法奥林匹克竞赛上，他使用高超的魔法本领，一路过关斩将，在最后时刻一举击败了前冠军“旅行者”，获得了魔法界最高的荣耀：女神奖杯！女神奖杯可不是一个普通的奖杯，她能够帮杰杰实现一个愿望。
杰杰本着实事求是的态度，审时度势，向女神奖杯提出了自己的愿望：想要一个女性朋友。
杰杰的愿望实现了，可是女性朋友却和他不在一个城市。杰杰想要知道：如果要到达女性朋友的所在城市，有多少种方案供他选择？
杰杰所在的世界有n个城市，从1到n进行编号。任意两个城市都通过有向道路连接。每个城市u有k个入点权：in[u][1],in[u][2]...in[u][k]，有k个出点权：ou[u][1],ou[u][2]...ou[u][k]。对于任意两个城市(u,v)（u可以等于v），u到v的道路条数为(ou[u][1]*in[v][1]+ou[u][2]*in[v][2]+...+ou[u][k]*in[v][k])条。杰杰有m次询问，每次询问由三元组(u,v,d)构成，询问从u城市通过不超过d条道路到达v城市的方案数。
为了温柔的杰杰和他的女性朋友的美好未来，帮助他解答这个问题吧。
Input
第一行读入两个正整数n，k，含义如题所示。接下来n行每行2k个整数，第i行代表第i个城市，前k个整数代表i号城市的出点权，后k个整数代表i号城市的入点权：
ou[i][1],ou[i][2],…,ou[i][k],in[i][1],in[i][2],…,in[i][k]
接下来一个整数m，表示m个询问。
接下来m行，每行三个整数:u,v,d，询问从u城市通过不超过d条道路到达v城市的方案数。
将每个方案所经过的道路，按顺序写成一个序列（序列可以为空）。两个方案不同，当且仅当他们的道路序列不完全相同。
Output
对于每个询问，输出一个方案数。由于答案可能太大，输出其除以1000000007后的余数。
Sample Input
5 2
2 5 4 3
7 9 2 4
0 1 5 2
6 3 9 2
2147483647 1000000001 233522 788488
10
1 1 0
2 2 1
2 4 5
4 3 10
3 4 50
1 5 1000
3 5 1000000000
1 2 500000000
4 5 2147483647
3 1 2147483647
Sample Output
51
170107227
271772358
34562176
890241289
8516097
383966304
432287042
326522835
HINT
数据规模和约定
n<=1000
k<=20
m<=50
保证1<=u, v<=n, 其它所有读入为不超过2147483647的非负整数
Source
数学问题 矩阵乘法 脑洞题 卡常数
将In矩阵转置，原来的求和变成了矩阵乘法的形式。
下面的In是原In的转置矩阵
$ (Out*In)^d $ 这是一个n*n的矩阵，看上去肯定TLE得飞起
我们换一个表示方法，$ (O*I) * (O*I) * (O*I)... $这么一串，等于$ O * (I * O)^(d-1) * I $，里面那个是m*m的矩阵，算起来无压力。
对于一个询问u和v，我们只需要先算出$ (IO)^{d-1} $，然后用O里u对应的那一行乘以中间这个矩阵，再乘以I里v对应的那一列，即可出解……吗？并不
注意到题目里问的是不多于d步，我们刚才算的是恰好d步的情况。也就是说中间要乘的其实是 (IO)^1 + (IO)^2 + (IO)^3 + ... +(IO)^d-1
这里需要用到POJ3233 的等比矩阵求和的姿势。
现在真的可以出解了
↓calc里那个快速幂形式的等比矩阵求和是别处找来的板子，有空的时候测试一下它和递归哪个快233
1 #include<iostream> 2 #include<algorithm> 3 #include<cstring> 4 #include<cstdio> 5 #include<cmath> 6 #define LL long long 7 using namespace std; 8 const int mod=1e9+7; 9 int read(){ 10 int x=0,f=1;char ch=getchar(); 11 while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();} 12 while(ch>='0' && ch<='9'){x=x*10-'0'+ch;ch=getchar();} 13 return x*f; 14 } 15 int sz; 16 struct Mat{ 17 LL x[23][23]; 18 void init(int sz){ 19 memset(x,0,sizeof x); 20 for(int i=1;i<=sz;i++)x[i][i]=1; 21 return; 22 } 23 Mat operator * (const Mat &b){ 24 Mat res; 25 for(int i=1;i<=sz;i++) 26 for(int j=1;j<=sz;j++){ 27 res.x[i][j]=0; 28 for(int k=1;k<=sz;k++){ 29 res.x[i][j]+=x[i][k]*b.x[k][j]%mod; 30 if(res.x[i][j]>mod)res.x[i][j]-=mod; 31 } 32 } 33 return res; 34 } 35 }mp,A,B,E; 36 inline void mat_add(Mat &a,const Mat &b){ 37 for(int i=1;i<=sz;i++) 38 for(int j=1;j<=sz;j++){ 39 a.x[i][j]+=b.x[i][j]; 40 if(a.x[i][j]>=mod)a.x[i][j]-=mod; 41 } 42 return; 43 } 44 int n,m,K; 45 LL in[1001][22],out[1001][22]; 46 LL IT[22][1001]; 47 LL ans[22]; 48 LL ANS=0; 49 void calc(int d){ 50 if(d<0){ 51 memset(A.x,0,sizeof A.x); 52 return; 53 } 54 int i,j,k; 55 A.init(K);B.init(K); 56 Mat bas=mp,tmp=mp; 57 // 58 while(d){ 59 if(d&1){ 60 mat_add(A,tmp*B); 61 B=B*bas; 62 } 63 Mat C=bas;mat_add(C,E); 64 tmp=tmp*C; 65 bas=bas*bas; 66 d>>=1; 67 } 68 // 69 return; 70 } 71 int main(){ 72 int i,j; 73 n=read();K=read(); 74 for(i=1;i<=n;i++){ 75 for(j=1;j<=K;j++)out[i][j]=read(); 76 for(j=1;j<=K;j++)in[i][j]=read(); 77 } 78 for(i=1;i<=n;i++) 79 for(j=1;j<=K;j++) 80 IT[j][i]=in[i][j]; 81 for(i=1;i<=K;i++) 82 for(j=1;j<=K;j++) 83 for(int k=1;k<=n;k++) 84 mp.x[i][j]=(mp.x[i][j]+IT[i][k]*out[k][j]%mod)%mod; 85 E.init(K); 86 // 87 sz=K; 88 m=read(); 89 int u,v,d; 90 while(m--){ 91 u=read();v=read();d=read(); 92 calc(d-1); 93 /* for(i=1;i<=sz;i++){ 94 for(j=1;j<=sz;j++){ 95 printf("%d %d :%lld\n",i,j,A.x[i][j]); 96 } 97 puts(""); 98 }*/ 99 ANS=0; 100 for(i=1;i<=sz;i++)ans[i]=0; 101 for(i=1;i<=sz;i++) 102 for(j=1;j<=sz;j++) 103 ans[j]=((LL)ans[j]+out[u][i]*A.x[i][j]%mod)%mod; 104 for(i=1;i<=sz;i++)(ANS+=ans[i]*IT[i][v]%mod)%=mod; 105 if(u==v)ANS=(ANS+1)%mod; 106 if(ANS<0)ANS+=mod; 107 printf("%lld\n",ANS); 108 } 109 return 0; 110 }
