[BZOJ3583]杰杰的女性朋友(矩阵快速幂)

杰杰的女性朋友

时间限制：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]=sumlimits_{j=1}^{n} (f[t-1][j]*sumlimits_{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$就好了。

但是还有一个问题，题目要求的是$$sumlimits_{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*(sumlimits_{i=0}^{d-1}D^{i})*I$$这种涉及到幂和的问题就不能直接使用矩阵快速幂解决。

我们可以首先预处理出所有$A_i=D^{2^i} (i=0,1,2,...)$，$B_i=sumlimits_{j=1}^{2^i} D^{j} (i=0,1,2,...) $，故$B_i=B_{i-1}A_i$

这样我们有$$sumlimits_{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);
    }  

void cal(int n){
    rep(i,1,m) rep(j,1,m) G[i][j]=S[i][j]=0;
    if(n<0)return;
    rep(i,1,m) S[i][i]=G[i][i]=1;
    for(int i=0; i<L; i++) if(n>>i&1){
        mul(B[i],G); rep(j,1,m) rep(k,1,m) up(S[j][k],C[j][k]);
        mul(G,A[i]); rep(j,1,m) rep(k,1,m) G[j][k]=C[j][k];
    }
}

剩下的只要根据输入数据建矩阵即可

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=l; i<=r; i++)
using namespace std;

const int N=1010,K=21,L=31,P=1000000007;
int n,m,q,x,y,z,ans,O[N][K],I[N][K],f[N];
int S[K][K],G[K][K],A[L][K][K],B[L][K][K],C[K][K];

void up(int &a,int b){ a+=b; if(a>=P)a-=P; }
void mul(int a[][K],int b[][K]){
    rep(i,1,m) rep(j,1,m) C[i][j]=0;
    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;
}

void cal(int n){
    rep(i,1,m) rep(j,1,m) G[i][j]=S[i][j]=0;
    if(n<0)return;
    rep(i,1,m) S[i][i]=G[i][i]=1;
    for(int i=0; i<L; i++) if(n>>i&1){
        mul(B[i],G); rep(j,1,m) rep(k,1,m) up(S[j][k],C[j][k]);
        mul(G,A[i]); rep(j,1,m) rep(k,1,m) G[j][k]=C[j][k];
    }
}

int main(){
    freopen("bzoj3583.in","r",stdin);
    freopen("bzoj3583.out","w",stdout);
    scanf("%d%d",&n,&m);
    rep(i,1,n){
        rep(j,1,m) scanf("%d",&O[i][j]);
        rep(j,1,m) scanf("%d",&I[i][j]);
    }
    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;
    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);
    }
    scanf("%d",&q);
    while (q--){
        scanf("%d%d%d",&x,&y,&z); cal(z-1);
        rep(i,1,m) f[i]=0; ans=0;
        rep(i,1,m) rep(j,1,m) f[i]=(f[i]+1ll*O[x][j]*S[j][i])%P;
        rep(i,1,m) ans=(ans+1ll*f[i]*I[y][i])%P;
        printf("%d
",(ans+(x==y))%P);
    }
    return 0;
}