【BZOJ #4362】Graph

Description

题目链接：

• BZOJ #4362
• 印象里是 FJOI2018 day1T3 的原题，具体题意是否是这样，以及数据范围不太记得了，但记得做法差不多。

给出一张 (n) 个点的有向图 (G(V,E))。对于任意两个点 (u,v)（(u) 可以等于 (v)），(u) 向 (v) 的连边数为

[sum_{i=1}^kout_{u,i} imes in_{v,i} ]

其中 (k) 和数组 (out,in) 均已知，现在给出 (m) 个询问，每次询问给出三个参数 (u, v, d)，你需要回答从节点 (u) 出发，经过不超过 (d) 条边到达节点 (v) 的路径有多少种。答案模 (10^9 + 7)。

(n leq 1000, k leq 20, m leq 50, d < 2^{31})

时空限制：(1 exttt s/256 exttt{MB})

Solution

因为我们要求的是不超过 (d) 条边的方案数，因此新建一个虚拟点 (n+1)，并且令 (out_{n+1,k+1}=in_{n+1,k+1}=1)。

对于每个询问 ((u,v,d))，就临时令 (out_{v,k+1}=1)，其含义是到达了点 (v) 后就多走一步走到 (n+1)，并且 (n+1) 有一条连向自己的自环。这样就将询问转化成，从 (u) 走到 (n+1)，恰好走 (d+1) 步的方案数。

不妨令 (n leftarrow n+1,kleftarrow k+1)，(A,B) 对应如上更改后的 (in) 和 (out) 的矩阵，答案就是

[R=(A imes B^T)^{d+1} ]

这个矩阵的第 (u) 行第 (n) 列，这样直接求是 (mathcal O(mn^3log d)) 的，非常逊色。

注意到 (k leq 20)，以及 (A,B) 分别是 (n imes k) 和 (k imes n) 的矩阵，直接把它们乘起来变成 (n imes n) 的矩阵，然后快速幂显然非常浪费。实际上利用矩阵乘法的结合律，我们可以将 (R) 写成

[R=A(B^T imes A)^dB^T ]

这样中间的矩阵 (B^TA) 就是 (k imes k) 的了，拿去矩阵快速幂就非常合适。

最后因为我们只要求 (R_{u,n})，所以不需要把 (R) 求出来。预处理出 (S=B^TA) 之后，每次询问求出 (S^d)，然后 (mathcal O(k^2)) 枚举一下就可以求出 (R_{u,n})。

每次询问临时令 (out_{v,k+1}=1) 的时候，会更改 (S) 的某一列，这个只要 (mathcal O(k^2)) 直接处理就好了。

总时间复杂度就是 (mathcal Oleft(nk^2+mk^3log d ight))。

#include <bits/stdc++.h>

template <class T>
inline void read(T &x)
{
	static char ch; 
	while (!isdigit(ch = getchar())); 
	x = ch - '0'; 
	while (isdigit(ch = getchar()))
		x = x * 10 + ch - '0'; 
}

template <class T>
inline void putint(T x)
{
	static char buf[15], *tail = buf; 
	if (!x)
		putchar('0'); 
	else
	{
		for (; x; x /= 10) *++tail = x % 10 + '0';
		for (; tail != buf; --tail) putchar(*tail); 
	}
}

typedef long long s64; 

const int mod = 1e9 + 7; 
const s64 MOD = 4611000029LL * mod; 

const int MaxK = 22; 
const int MaxN = 1e3 + 5; 

struct mat
{
	int r, c; 
	int a[22][22]; 

	mat(){}
	mat(int _r, int _c):
		r(_r), c(_c) {memset(a, 0, sizeof(a));}

	inline void init()
	{
		for (int i = 1; i <= r; ++i)
			a[i][i] = 1; 
	}

	inline mat operator * (mat rhs) const
	{
		mat res(r, rhs.c); 
		for (int i = 1; i <= r; ++i)
			for (int j = 1; j <= rhs.c; ++j)
			{
				s64 sum = 0; 
				for (int k = 1; k <= c; ++k)
				{
					sum += 1LL * a[i][k] * rhs.a[k][j]; 
					if (sum >= MOD)
						sum -= MOD; 
				}
				res.a[i][j] = sum % mod; 
			}
		return res; 
	}

	inline mat operator ^ (int p) const
	{
		mat res(r, c), x = *this; 
		res.init(); 
		for (; p; p >>= 1, x = x * x)
			if (p & 1)
				res = res * x; 
		return res; 
	}
}T; 

int n, K, m; 
int prod[MaxK][MaxK]; 
int A[MaxN][MaxK], B[MaxK][MaxN]; 

int main()
{
	read(n), read(K); 
	for (int i = 1; i <= n; ++i)
	{
		for (int j = 1; j <= K; ++j)
		{
			read(A[i][j]); 
			A[i][j] %= mod; 
		}
		for (int j = 1; j <= K; ++j)
		{
			read(B[j][i]); 
			B[j][i] %= mod; 
		}
	}

	++K, ++n; 
	A[n][K] = B[K][n] = 1; 

	for (int i = 1; i <= K; ++i)
		for (int j = 1; j <= K; ++j)
		{
			s64 res = 0; 
			for (int k = 1; k <= n; ++k)
			{
				res += 1LL * B[i][k] * A[k][j]; 
				if (res >= MOD)
					res -= MOD; 
			}
			prod[i][j] = res % mod; 
		}
	T = mat(K, K); 

	read(m); 
	while (m--)
	{
		int u, v, d; 
		read(u), read(v), read(d); 

		A[v][K] = 1; 
		for (int i = 1; i <= K; ++i)
			for (int j = 1; j <= K; ++j)
			{
				T.a[i][j] = prod[i][j]; 
				if (j == K)
					T.a[i][j] = (T.a[i][j] + 1LL * B[i][v] * A[v][j]) % mod; 
			}

		T = T ^ d; 

		int res = 0; 
		for (int i = 1; i <= K; ++i)
		{
			int sum = 0; 
			for (int j = 1; j <= K; ++j)
				sum = (sum + 1LL * T.a[i][j] * B[j][n]) % mod; 
			res = (res + 1LL * A[u][i] * sum) % mod; 
		}
		A[v][K] = 0; 
		
		putint(res), putchar('
'); 
	}

	return 0; 
}