杰杰的女性朋友
时间限制：10s 空间限制：256MB
题目描述
杰杰是魔法界的一名传奇人物。他对魔法具有深刻的洞察力，惊人的领悟力，以及令人叹为观止的创造力。自从他从事魔法竞赛以来，短短几年时间，就已经成为 世界公认的实力最强的魔法选手之一。更让人惊叹的是，他几乎没有借助外界力量，完全凭借自己的努力达到了普通人难以企及的高度。在最近的世界魔法奥林匹克 竞赛上，他使用高超的魔法本领，一路过关斩将，在最后时刻一举击败了前冠军“旅行者”，获得了魔法界最高的荣耀：女神奖杯！女神奖杯可不是一个普通的奖 杯，她能够帮杰杰实现一个愿望。
杰杰本着实事求是的态度，审时度势，向女神奖杯提出了自己的愿望：想要一个女性朋友。杰杰的愿望实现了，可是女性朋友却和他不在一个城市。杰杰想要知道：如果要到达女性朋友的所在城市，有多少种方案供他选择？
杰杰所在的世界有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城市的方案数。
为了温柔的杰杰和他的女性朋友的美好未来，帮助他解答这个问题吧。输入格式
第一行读入两个正整数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城市的方案数。
将每个方案所经过的道路，按顺序写成一个序列（序列可以为空）。两个方案不同，当且仅当他们的道路序列不完全相同。输出格式
对于每个询问，输出一个方案数。由于答案可能太大，输出其除以1000000007后的余数。样例输入
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样例输出
1 51 170107227 271772358 34562176 890241289 8516097 383966304 432287042 326522835提示
数据规模和约定
n<=1000
k<=20
m<=50保证1<=u, v<=n, 其它所有读入为不超过2147483647的非负整数
题目来源
By 佚名提供
FJOI2018一试就是直接使用了这道清华集训原题。
首先列出DP状态转移方程，$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})$$
这样做的复杂度是$O(n^2d)$，而 $d \leqslant 2^{31}-1$，显然无法在时限内出解。
观察这个转移方程，不难看出这是裸的矩阵快速幂，于是可以在$O(n^3 \log d)$时间内出解。
然而这个复杂度仍然不够优，事实上，连FJOI2018现场最低的一档部分分都无法通过。
于是需要进一步观察矩阵的性质：
$f$是一个$1*n$的矩阵，$O$是一个$n*k$的矩阵，$I$是$n*k$的，而$C=OI^T$所以$C$是$n*n$的。由数据范围可知，如果我们能将$n*n$的矩阵乘法优化到$k*k$，那么就可以通过全部数据。
不难发现答案$$f[d]=f[0]*C^d=f[0]*(OI^T)^d=f[0]*O*(I^TO)^{d-1}*I$$而$I^TO$是$k*k$的，所以我们只要求$D=I^TO$就好了。
但是还有一个问题，题目要求的是$$\sum\limits_{i=0}^{d} f[i]$$ 也就是$$f[0]+f[0]*O{(I^{T}O)}^{0}I +f[0]*O(I^TO)^{1}I+ \ldots +f[0]*O(I^TO)^{d-1}I\\=f[0]*O*(\sum\limits_{i=0}^{d-1}D^{i})*I$$这种涉及到幂和的问题就不能直接使用矩阵快速幂解决。
我们可以首先预处理出所有$A_i=D^{2^i} (i=0,1,2,...)$，$B_i=\sum\limits_{j=1}^{2^i} D^{j} (i=0,1,2,...) $，故$B_i=B_{i-1}A_i$
这样我们有$$\sum\limits_{i=0}^{d-1}D^{i}=E+(D+D^{1}+...+D^{d_1})+(D^{d_{1}+1}+D^{d_{1}+2}+...+D^{d_1+d_2})+...$$其中$d_i$为d的二进制第i个1代表的数。$E$为单位矩阵
对于每个括号分别考虑，$$D+D^{1}+...+D^{d_1}=B^{d_1}$$$$D^{d_{1}+1}+D^{d_{1}+2}+...+D^{d_1+d_2}=(B_{d_2}*A_{d_1})$$
以此类推，就可以得到最终的答案。代码片段如下：
rep(i,1,m) rep(j,1,m) B[0][i][j]=A[0][i][j];rep(i,0,L-2){mul(A[i],A[i]);rep(j,1,m) rep(k,1,m) A[i+1][j][k]=C[j][k];rep(j,1,m) up(A[i][j][j],1);mul(B[i],A[i]);rep(j,1,m) rep(k,1,m) B[i+1][j][k]=C[j][k];rep(j,1,m) up(A[i][j][j],P-1);}
1 void cal(int n){ 2 rep(i,1,m) rep(j,1,m) G[i][j]=S[i][j]=0; 3 if(n<0)return; 4 rep(i,1,m) S[i][i]=G[i][i]=1; 5 for(int i=0; i<L; i++) if(n>>i&1){ 6 mul(B[i],G); rep(j,1,m) rep(k,1,m) up(S[j][k],C[j][k]); 7 mul(G,A[i]); rep(j,1,m) rep(k,1,m) G[j][k]=C[j][k]; 8 } 9 }
剩下的只要根据输入数据建矩阵即可
1 #include<cstdio> 2 #include<algorithm> 3 #define rep(i,l,r) for (int i=l; i<=r; i++) 4 using namespace std; 5 6 const int N=1010,K=21,L=31,P=1000000007; 7 int n,m,q,x,y,z,ans,O[N][K],I[N][K],f[N]; 8 int S[K][K],G[K][K],A[L][K][K],B[L][K][K],C[K][K]; 9 10 void up(int &a,int b){ a+=b; if(a>=P)a-=P; } 11 void mul(int a[][K],int b[][K]){ 12 rep(i,1,m) rep(j,1,m) C[i][j]=0; 13 rep(i,1,m) rep(j,1,m) rep(k,1,m) C[i][k]=(C[i][k]+1ll*a[i][j]*b[j][k])%P; 14 } 15 16 void cal(int n){ 17 rep(i,1,m) rep(j,1,m) G[i][j]=S[i][j]=0; 18 if(n<0)return; 19 rep(i,1,m) S[i][i]=G[i][i]=1; 20 for(int i=0; i<L; i++) if(n>>i&1){ 21 mul(B[i],G); rep(j,1,m) rep(k,1,m) up(S[j][k],C[j][k]); 22 mul(G,A[i]); rep(j,1,m) rep(k,1,m) G[j][k]=C[j][k]; 23 } 24 } 25 26 int main(){ 27 freopen("bzoj3583.in","r",stdin); 28 freopen("bzoj3583.out","w",stdout); 29 scanf("%d%d",&n,&m); 30 rep(i,1,n){ 31 rep(j,1,m) scanf("%d",&O[i][j]); 32 rep(j,1,m) scanf("%d",&I[i][j]); 33 } 34 rep(k,1,n) rep(i,1,m) rep(j,1,m) A[0][i][j]=(A[0][i][j]+1ll*I[k][i]*O[k][j])%P; 35 rep(i,1,m) rep(j,1,m) B[0][i][j]=A[0][i][j]; 36 rep(i,0,L-2){ 37 mul(A[i],A[i]); 38 rep(j,1,m) rep(k,1,m) A[i+1][j][k]=C[j][k]; 39 rep(j,1,m) up(A[i][j][j],1); 40 mul(B[i],A[i]); 41 rep(j,1,m) rep(k,1,m) B[i+1][j][k]=C[j][k]; 42 rep(j,1,m) up(A[i][j][j],P-1); 43 } 44 scanf("%d",&q); 45 while (q--){ 46 scanf("%d%d%d",&x,&y,&z); cal(z-1); 47 rep(i,1,m) f[i]=0; ans=0; 48 rep(i,1,m) rep(j,1,m) f[i]=(f[i]+1ll*O[x][j]*S[j][i])%P; 49 rep(i,1,m) ans=(ans+1ll*f[i]*I[y][i])%P; 50 printf("%d\n",(ans+(x==y))%P); 51 } 52 return 0; 53 }
