jzoj3545. 【清华集训2014】杰杰的女性朋友

Description

杰杰是魔法界的一名传奇人物。他对魔法具有深刻的洞察力，惊人的领悟力，以及令人叹为观止的创造力。自从他从事魔法竞赛以来，短短几年时间，就已经成为世界公认的实力最强的魔法选手之一。更让人惊叹的是，他几乎没有借助外界力量，完全凭借自己的努力达到了普通人难以企及的高度。在最近的世界魔法奥林匹克竞赛上，他使用高超的魔法本领，一路过关斩将，在最后时刻一举击败了前冠军“旅行者”，获得了魔法界最高的荣耀：女神奖杯！女神奖杯可不是一个普通的奖杯，她能够帮杰杰实现一个愿望。

杰杰本着实事求是的态度，审时度势，向女神奖杯提出了自己的愿望：想要一个女性朋友。

杰杰的愿望实现了，可是女性朋友却和他不在一个城市。杰杰想要知道：如果要到达女性朋友的所在城市，有多少种方案供他选择？

杰杰所在的世界有n个城市，从1到n进行编号。任意两个城市都通过有向道路连接。每个城市u有k个入点权：in[u][1],in[u][2]...in[u][k]in[u][1],in[u][2]...in[u][k]in[u][1],in[u][2]...in[u][k]，有k个出点权：ou[u][1],ou[u][2]...ou[u][k]ou[u][1],ou[u][2]...ou[u][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])(ou[u][1]*in[v][1]+ou[u][2]*in[v][2]+...+ou[u][k]*in[v][k])(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]ou[i][1],ou[i][2],…,ou[i][k],in[i][1],in[i][2],…,in[i][k]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

1
51
170107227
271772358
34562176
890241289
8516097
383966304
432287042
326522835

Data Constraint

赛时

终于期末考试爆炸完毕，来信息学划水了。

这题比赛是没有多想，然后就水了个暴力上去，想拿30。
然后爆0。
QWQ

题解

题解还是蛮妙的。
首先，我们考虑答案计算。
我们发现，答案计算柿子长这样：
I[i][j]∗O[k][j]I[i][j]*O[k][j]I[i][j]∗O[k][j]
那么不妨把OOO换一下变成：
I[i][j]∗O[j][k]I[i][j]*O[j][k]I[i][j]∗O[j][k]
这不就是矩乘的形式吗？

那么走刚好ddd步答案即为：(I∗O)d(I*O)^d(I∗O)d
矩阵快速幂？
没戳，但是问题是我们的I∗OI*OI∗O是一个1000∗10001000*10001000∗1000大小的矩阵，所以不可接受。
然鹅交换一下就是O∗IO*IO∗I是一个20∗2020*2020∗20的矩阵，可以接受！
但是矩阵不满足交换律，但是满足结合律。如果我们展开：
I∗O∗I∗O∗……∗I∗O∗I∗OI*O*I*O*……*I*O*I*OI∗O∗I∗O∗……∗I∗O∗I∗O
结合一下
I∗(O∗I∗O∗……∗I∗O∗I)∗O=I∗(OI)d−1∗OI*(O*I*O*……*I*O*I)*O=I*(OI)^{d-1}*OI∗(O∗I∗O∗……∗I∗O∗I)∗O=I∗(OI)d−1∗O
那么中间就可以快速rush了，然后外面就再1000∗20∗201000*20*201000∗20∗20就可以做完。

然鹅这样算只能算出刚好ddd步的答案。
然后接下来利用一个倍增的奇妙思想即可做出1−d1-d1−d步的答案。

设f[i]f[i]f[i]表示(OI)2i(OI)^{2^i}(OI)2i设g[i]g[i]g[i]表示∑j=12i(OI)j\sum_{j=1}^{2^i}(OI)^j∑j=12i(OI)j
然后利用倍增的思想就可以求出g[i]g[i]g[i]
然后对于任意数我们可以把它拆成若干个2x2^x2x的形式，用ggg来计算即可。

代码

#include <cstdio>
#include <iostream>
#include <cmath>
#include <cstring>
using namespace std;
const long long mo=1000000007;int n,m,x,y,z;
long long ans,ru[1010][25],cu[1010][25],mi[32];struct matrix1{long long a[1010][1010];int n,m;
}I,O;struct matrix{long long a[25][25];int n,m;
}f[32],g[32];inline int read() {int x = 0, f = 0; char c = getchar();while (c < '0' || c > '9') f = (c == '-') ? 1 : f, c = getchar();while (c >= '0' && c <= '9') x = (x << 1) + (x << 3) + (c ^ 48), c = getchar();return f ? -x : x;
}matrix1 cheng1(matrix1 a,matrix1 b)
{matrix1 c;return c;
}matrix cheng(matrix a,matrix b)
{matrix c;for (int i=1;i<=a.n;i++){for (int j=1;j<=b.m;j++){c.a[i][j]=0;for (int k=1;k<=a.m;k++){c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mo)%mo;}}}c.n=a.n;c.m=b.m;return c;
}matrix qsm(matrix a,long long b)
{matrix t;matrix y;for (int i=1;i<=a.n;i++){for (int j=1;j<=a.m;j++){t.a[i][j]=0;y.a[i][j]=a.a[i][j];}t.a[i][i]=1;}t.n=y.n=a.n;t.m=y.m=a.m;while (b>0){if (b%2==1) t=cheng(t,y);y=cheng(y,y);b=b/2;}return t;
}matrix jia(matrix a,matrix b)
{matrix c;for (int i=1;i<=a.n;i++){for (int j=1;j<=a.m;j++){c.a[i][j]=(a.a[i][j]+b.a[i][j])%mo;}}return c;
}int main()
{scanf("%d%d",&n,&m);for (int i=1;i<=n;i++){for (int j=1;j<=m;j++){scanf("%lld",&cu[i][j]);cu[i][j]=cu[i][j]%mo;O.a[i][j]=cu[i][j];O.n=n;O.m=m;}for (int j=1;j<=m;j++){scanf("%lld",&ru[i][j]);ru[i][j]=ru[i][j]%mo;I.a[j][i]=ru[i][j];I.n=m;I.m=n;}}matrix A;for (int i=1;i<=m;i++){for (int j=1;j<=m;j++){A.a[i][j]=0;for (int k=1;k<=n;k++){A.a[i][j]=(A.a[i][j]+I.a[i][k]*O.a[k][j]%mo)%mo;}}}A.n=m;A.m=m;mi[0]=1;for (int i=0;i<=30;i++){if (i<31) mi[i+1]=mi[i]*2;f[i].n=f[i].m=m;f[i]=qsm(A,mi[i]);}g[0]=A;for (int i=1;i<=30;i++){g[i]=jia(g[i-1],cheng(g[i-1],f[i-1]));}int q;scanf("%d",&q);while (q>0){q--;scanf("%d%d%d",&x,&y,&z);if (z==0){printf("%d\n",(x==y?1:0));continue;}z--;ans=(x==y?1:0);matrix qzh;matrix qzj;qzh.n=qzh.m=qzj.n=qzj.m=m;memset(qzh.a,0,sizeof(qzh.a));memset(qzj.a,0,sizeof(qzj.a));for (int i=1;i<=m;i++){qzh.a[i][i]=qzj.a[i][i]=1;}for (int i=30;i>=0;i--){if (mi[i]<=z){qzh=jia(qzh,cheng(qzj,g[i]));qzj=cheng(qzj,f[i]);z-=mi[i];}}matrix1 jl;for (int i=1;i<=O.n;i++){for (int j=1;j<=qzh.m;j++){jl.a[i][j]=0;for (int k=1;k<=O.m;k++){jl.a[i][j]=(jl.a[i][j]+O.a[i][k]*qzh.a[k][j]%mo)%mo;}}}for (int i=1;i<=m;i++) ans=(ans+jl.a[x][i]*I.a[i][y]%mo)%mo;printf("%lld\n",ans);}return 0;
}