【BZOJ #4362】Graph

Description

题目链接:

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

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

\[\sum_{i=1}^kout_{u,i}\times 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\texttt s/256 \texttt{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,k\leftarrow k+1\)\(A,B\) 对应如上更改后的 \(in\)\(out\) 的矩阵,答案就是

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

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

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

\[R=A(B^T\times A)^dB^T \]

这样中间的矩阵 \(B^TA\) 就是 \(k\times 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 O\left(nk^2+mk^3\log d\right)\)

#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('\n'); 
	}

	return 0; 
}
posted @ 2020-02-17 18:07  changle_cyx  阅读(203)  评论(0编辑  收藏  举报