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

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]\ast in[v][1]+ou[u][2]\ast in[v][2]+...+ou[u][k]\ast 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],\dots ,ou[i][k],in[i][1],in[i][2],\dots ,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]\ast O[k][j]$
那么不妨把$O$换一下变成：
$I[i][j]\ast O[j][k]$
这不就是矩乘的形式吗？

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

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

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

代码

#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;
}