目录
壹 —— 矩阵
“矩阵就是映射!!!”
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
\qquad
———— 《程序员的数学》
表面上看,矩阵无非就是把一堆数整整齐齐地排列在一起,像这样
[ 1 5 1 8 1 1 0 6 4 4 0 7 ] \begin{bmatrix} 1 &5 &1 &8 \\ 1 &1 &0 &6 \\ 4 &4 &0 &7 \end{bmatrix} 114514100867
有数字的地方自然就有运算
矩阵的运算:
加、减
要求两个矩阵形态相同,结果为对应位置上的数字的代数和 / 差
乘
与加、减不同
[ 1 2 1 1 3 4 ] × [ 2 2 2 3 1 1 ] = [ 1 × 2 + 2 × 2 + 1 × 1 1 × 2 + 2 × 3 + 1 × 1 1 × 2 + 3 × 2 + 4 × 1 1 × 2 + 3 × 3 + 4 × 1 ] \begin{bmatrix} 1 &2 &1 \\ 1 &3 &4 \end{bmatrix} \times \begin{bmatrix} 2 &2 \\ 2 &3 \\ 1 &1 \end{bmatrix}=\begin{bmatrix} 1\times2+2\times2+1\times1 &1\times2+2\times3+1\times1 \\ 1\times2+3\times2+4\times1 &1\times2+3\times3+4\times1 \end{bmatrix} [112314]× 221231 =[1×2+2×2+1×11×2+3×2+4×11×2+2×3+1×11×2+3×3+4×1]
从上述过程,我们可以发现一些性质:
- n × k n\times k n×k 的矩阵乘上 k × m k\times m k×m 的矩阵能够得到一个 n × m n\times m n×m 的矩阵
- 矩阵乘法满足结合率
- 代码实现两个矩阵相乘的复杂度为 O ( n 3 ) O(n^3) O(n3)
其中性质 2、3 是矩阵优化的关键
从以上三种运算使我们发现一些特殊的矩阵:
特殊矩阵
- 单位矩阵,记为 I I I,形如: [ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ] \begin{bmatrix} 1 &0 &0 &0 \\ 0 & 1 & 0 &0 \\ 0 & 0 &1 &0 \\ 0 & 0 & 0 &1 \end{bmatrix} 1000010000100001
的 方阵
为什么叫单位矩阵?我们会发现 任何矩阵与其相乘都会得到本身(前提是能够相乘)
单位——可以理解为整数中的 “1”
-
零矩阵,记为 o o o ,顾名思义: [ 0 0 0 0 0 0 0 0 0 ] \begin{bmatrix} 0& 0 &0 \\ 0 & 0 & 0\\ 0& 0 &0 \end{bmatrix} 000000000
任何矩阵与其相乘结果都是零矩阵 -
转置矩阵, A A A 的转置矩阵记为 A T A^T AT [ 1 3 5 2 4 6 ] ⟶ [ 1 2 3 4 5 6 ] \begin{bmatrix} 1 &3 &5 \\ 2 &4 &6 \end{bmatrix}\longrightarrow\begin{bmatrix} 1 &2 \\ 3 & 4\\ 5 &6 \end{bmatrix} [123456]⟶ 135246
将 ( i , j ) (i,j) (i,j) 位置的数字变为 ( j , i ) (j,i) (j,i) ,就得到了转置矩阵
转置矩阵 与 原矩阵 的行列式相等
- 分块矩阵 —— 矩阵套矩阵
作矩阵乘法时,有时我们会感到整体作很困难,那么此时我们就可以将其 分块 [ 1 2 3 1 0 0 ] ⟶ [ 1 A 1 o ] ,其中 A = [ 2 3 ] \begin{bmatrix} 1 & 2 &3 \\ 1 & 0 &0 \end{bmatrix}\longrightarrow \begin{bmatrix} 1 & A\\ 1 & o \end{bmatrix},其中\ A=\begin{bmatrix} 2 &3 \end{bmatrix} [112030]⟶[11Ao],其中 A=[23]
分块矩阵整体 与 另一个矩阵 作乘法,与将其拆开后的 结果相同
除
矩阵没有直接的除法,但是有求逆
A A A 的逆矩阵,即 使得 A × X = I A\times X=I A×X=I 的矩阵 X X X,记为 A − 1 A^{-1} A−1
其存在的充要条件时 行列式不为零
后面可能还会用到?
乘方
矩阵乘方 与 数的乘方 定义相同,都是若干次乘法运算
但矩阵乘方可以形象的理解为:映射的叠加
矩阵的应用:
本质上讲,矩阵优化的是 具有线性关系的转移过程
简单的例子1:
求 斐波那契数列 的第 K K K 项对 m m m 取余的结果,其中 K ≤ 1 e 18 K\leq 1e18 K≤1e18
分析
常规的递推显然无法解决,我们从矩阵的角度考虑:
首先建出初始矩阵,令 F = [ f 1 f 0 ] F=\begin{bmatrix} f_1 &f_0 \end{bmatrix} F=[f1f0]
然后是伴随矩阵(用方阵比较简便),要使得 F F F 乘上伴随矩阵后得到递推的下一层,令 T = [ 1 1 1 0 ] T=\begin{bmatrix} 1 &1\\ 1 &0 \end{bmatrix} T=[1110]
那么 F × T = [ f 2 f 1 ] F\times T =[f_2\ \ \ \ f_1] F×T=[f2 f1],同理 F × T K = [ f K + 1 f K ] F\times T^K = [f_{K+1} \ \ \ f_K] F×TK=[fK+1 fK]
因为矩阵有结合律,我们对后面的 T K T^K TK 整体处理,用类似整数快速幂的方式对其进行优化
void Mul_fib()
{
int tmp[2] = {} ;
for(int i = 0 ; i < 2 ; i ++ ) {
for(int k = 0 ; k < 2 ; k ++ ) {
tmp[i] += ( f[k]*t[k][i] ) % mod ;
tmp[i] %= mod ;
}
}
f[0] = tmp[0] , f[1] = tmp[1] ;
}
void Mul_t()
{
int tmp[2][2] = {} ;
for(int i = 0 ; i < 2 ; i ++ ) {
for(int j = 0 ; j < 2 ; j ++ ) {
for(int k = 0 ; k < 2 ; k ++ ) {
tmp[i][j] += ( t[i][k] * t[k][j] ) % mod ;
tmp[i][j] %= mod ;
}
}
}
t[0][0] = tmp[0][0] , t[0][1] = tmp[0][1] ;
t[1][0] = tmp[1][0] , t[1][1] = tmp[1][1] ;
}
int main()
{
while( scanf("%d" , &n ) && n != -1 ) {
f[0] = 1 , f[1] = 0 ;
t[0][0] = 1 , t[0][1] = 1 , t[1][0] = 1 , t[1][1] = 0 ;
while( n ) {
if( n&1 ) Mul_fib() ;
n = n>>1 ;
Mul_t() ;
}
printf("%d\n" , f[1] ) ;
}
return 0 ;
}
这样就得到了 O ( n 3 l o g K ) O(n^3logK) O(n3logK) 的算法,其中 n n n 为矩阵边长, K K K 为操作次数
矩阵可以维护什么样的递推关系呢?
我们再来看一个例子:
例 2 维护特殊的递推式子
求 T ( n ) = ∑ i = 1 n i × f i T(n)=\sum_{i=1}^{n}i\times f_i T(n)=∑i=1ni×fi 的值,其中 f i f_i fi 表示斐波那契第 i i i 项
首先考虑能否在矩阵中维护 i f i if_i ifi,令 g i = i f i g_i=if_i gi=ifi
有:
g
i
=
i
×
f
i
=
i
(
f
i
−
1
+
f
i
−
2
)
=
(
i
−
1
)
×
f
i
−
1
+
(
i
−
2
)
×
f
i
−
2
+
f
i
−
1
+
2
f
i
−
2
g_i=i\times f_i=i\ (f_{i-1}+f_{i-2})=(i-1)\times f_{i-1}+(i-2)\times f_{i-2}+f_{i-1}+2f_{i-2}
gi=i×fi=i (fi−1+fi−2)=(i−1)×fi−1+(i−2)×fi−2+fi−1+2fi−2
=
g
i
−
1
+
g
i
−
2
+
f
i
−
1
+
2
f
i
−
2
\ \quad=g_{i-1}+g_{i-2}+f_{i-1}+2f_{i-2}
=gi−1+gi−2+fi−1+2fi−2
这样经过变形,我们得到了 线性的递推关系
那么只需要在 F F F 矩阵中维护 [ g i g i − 1 f i f i − 1 ] [g_i\ \ g_{i-1} \ \ f_i \ \ f_{i-1}] [gi gi−1 fi fi−1] 四个量
引入前缀和,只需额外维护 S i − 1 S_{i-1} Si−1,每次把 g i g_i gi 加上去即可
总结来说: [ g i g i − 1 f i f i − 1 S i − 1 ] × [ 1 1 0 0 1 1 0 0 0 0 1 0 1 1 0 2 0 1 0 0 0 0 0 0 1 ] = [ g i + 1 g i f i + 1 f i S i ] \begin{bmatrix} g_i &g_{i-1} &f_i &f_{i-1} &S_{i-1} \end{bmatrix}\ \ \times \ \ \begin{bmatrix} 1 &1 &0 &0 &1 \\ 1 &0 &0 &0 &0 \\ 1 &0 &1 &1 &0 \\ 2 &0 &1 &0 &0 \\ 0 &0 &0 &0 &1 \end{bmatrix}=\begin{bmatrix} g_{i+1} &g_i &f_{i+1} &f_i &S_i \end{bmatrix} [gigi−1fifi−1Si−1] × 1112010000001100010010001 =[gi+1gifi+1fiSi]
由上面的例子我们可以发现,矩阵可以维护 线性递推关系
如果遇到 多个线性递推相乘,可以根据每个递推式的性质 向前化简( i i i 化成 i − 1 i-1 i−1 ),只要得到的还是线性式子,就能用矩阵优化
例 3 分块矩阵应用
本题是分块矩阵的一个应用,将
A
A
A 看做一个数字正常转移即可
下面来看 实际DP 过程中矩阵的优化
例 4 技巧:二维转一维 / 周期处理
本题的关键在于:二维的方格很难整体维护
于是将二维转化为一维,给每个节点编上号,转化成 一维向量,正常转移
然后操作序列长度比较小,容易想到用最小公倍数作为一个周期
求出一个周期整体的矩阵,作快速幂,最后 for 一遍补上周期余数的贡献
#include<bits/stdc++.h>
using namespace std ;
#define int long long
typedef long long LL ;
const int N = 11 ;
// 矩阵乘法 适于优化 一维递推关系,二维时考虑 转化成一维
int n , m , mp[N][N] , len[N] , Lcm = 1 , tot ;
int dx[305] , dy[305] ;
char opt[N][6*N] ;
LL T , f[N*N] , t[N*N][N*N] , tp[6*N][N*N][N*N] , ans ;
void Mul_tp_t( int op ) // t = t * tp
{
LL tmp[N*N][N*N] = {} ;
for(int i = 0 ; i <= tot ; i ++ ) {
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[i][j] = ( tmp[i][j] + t[i][k]*tp[op][k][j] ) ;
}
}
}
memcpy( t , tmp , sizeof t ) ;
}
void pre_square()
{
for(int i = 0 ; i <= tot ; i ++ ) t[i][i] = 1 ;
for(int op = 1 ; op <= Lcm ; op ++ ) { // 先计算一个 周期 内的情况
tp[op][0][0] = 1 ;
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= m ; j ++ ) {
int id = (i-1)*m + j , typ = mp[i][j] ;
char p = opt[typ][(op-1)%len[typ]] ; // tp[op] : op-1 -> op
if( p >= '0' && p <= '9' ) {
tp[op][0][id] = (p^48) ;
tp[op][id][id] = 1 ;
}
else if( p == 'D' ) {
tp[op][id][id] = 0 ;
}
else {
int to = (i+dx[p]-1)*m + (j+dy[p]) ;
if( i+dx[p] >= 1 && i+dx[p] <= n && j+dy[p] >= 1 && j+dy[p] <= m ) {
tp[op][id][to] = 1 ;
}
}
}
}
Mul_tp_t( op ) ;
}
}
void Mul_f()
{
LL tmp[N*N] = {} ;
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[j] = ( tmp[j] + f[k]*t[k][j] ) ;
}
}
memcpy( f , tmp , sizeof f ) ;
}
void Mul_t()
{
LL tmp[N*N][N*N] = {} ;
for(int i = 0 ; i <= tot ; i ++ ) {
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[i][j] = ( tmp[i][j] + t[i][k]*t[k][j] ) ;
}
}
}
memcpy( t , tmp , sizeof t ) ;
}
void Mul_tp_f( int op )
{
LL tmp[N*N] = {} ;
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[j] = ( tmp[j] + f[k]*tp[op][k][j] ) ;
}
}
memcpy( f , tmp , sizeof f ) ;
}
signed main()
{
dx['N'] = -1 ; dx['S'] = 1 ; dy['E'] = 1 ; dy['W'] = -1 ;
int act ;
scanf("%lld%lld%lld%lld" , &n , &m , &T , &act ) ;
tot = n*m ;
char c ;
for(int i = 1 ; i <= n ; i ++ ) {
scanf("\n") ;
for(int j = 1 ; j <= m ; j ++ ) {
c = getchar() ;
mp[i][j] = c^48 ;
}
}
for(int i = 0 ; i < act ; i ++ ) {
scanf("\n%s" , opt[i] ) ;
len[i] = strlen( opt[i] ) ;
Lcm = Lcm*len[i]/__gcd( Lcm , len[i] ) ;
}
pre_square() ;
f[0] = 1 ;
LL R = T / Lcm ;
while( R ) {
if( R & 1 ) {
Mul_f() ;
}
R = R >> 1 ;
Mul_t() ;
}
T = T % Lcm ;
for(int op = 1 ; op <= T ; op ++ ) {
Mul_tp_f( op ) ;
}
for(int i = 1 ; i <= tot ; i ++ ) ans = max( ans , f[i] ) ;
printf("%lld" , ans ) ;
return 0 ;
}
常用套路是,先以操作次数(通常都很大)为阶段,写出转移式,然后将一个阶段整体视为一个向量,根据转移方程构造伴随矩阵
例 5 结合较复杂的 DP
本题前置:修复DNA
“不包含” 的限制 常规方法很难维护,字符串间不知道怎么转移
于是想到 在 AC 自动机上作 DP ,简单说就是以 AC 自动机上每个节点为一个 DP 状态,将不合法状态赋值 inf 。每增加一个长度时,考虑 节点 j j j 转移到 节点 t r [ j ] [ c ] tr[j][c] tr[j][c]
对于本题,我们发现这个转移式子符合 阶段 i i i 的所有状态都由 i − 1 i-1 i−1 转来,且转移方程符合线性,那么可以矩阵乘法优化
#include<bits/stdc++.h>
using namespace std ;
typedef long long LL ;
const LL mod = 998244353 ;
LL n , m ;
char ch[8] ;
int tr[1<<8][2] , tot , fail[1<<8] ;
bool ed[1<<8] ;
void Insert()
{
int len = strlen( ch+1 ) , p = 0 ;
for(int i = 1 ; i <= len ; i ++ ) {
int c = ch[i] - 'a' ;
if( !tr[p][c] ) tr[p][c] = ++tot ;
p = tr[p][c] ;
}
ed[p] = 1 ;
}
void build_ac()
{
queue<int> q ;
if( tr[0][0] ) q.push( tr[0][0] ) ;
if( tr[0][1] ) q.push( tr[0][1] ) ;
while( !q.empty() ) {
int x = q.front() ; q.pop() ;
if( !x ) continue ;
if( ed[fail[x]] ) ed[x] = 1 ; // 传标记
for(int i = 0 ; i < 2 ; i ++ ) {
if( tr[x][i] ) {
if( ed[x] ) ed[tr[x][i]] = 1 ; // 下传标记
fail[tr[x][i]] = tr[fail[x]][i] , q.push( tr[x][i] ) ;
}
else {
tr[x][i] = tr[fail[x]][i] ;
}
}
}
}
LL f[(1<<8)] , t[(1<<8)][(1<<8)] ;
void Mul_f()
{
LL tmp[1<<8] = {} ;
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[j] = ( tmp[j] + f[k]*t[k][j]%mod ) % mod ;
}
}
memcpy( f , tmp , sizeof f ) ;
}
void Mul_t()
{
LL tmp[1<<8][1<<8] = {} ;
for(int i = 0 ; i <= tot ; i ++ ) {
for(int j = 0 ; j <= tot ; j ++ ) {
for(int k = 0 ; k <= tot ; k ++ ) {
tmp[i][j] = ( tmp[i][j] + t[i][k]*t[k][j]%mod ) % mod ;
}
}
}
memcpy( t , tmp , sizeof t ) ;
}
int main()
{
scanf("%lld%lld" , &n , &m ) ;
for(int i = 1 ; i <= m ; i ++ ) {
scanf("\n%s" , ch+1 ) ;
Insert() ;
}
build_ac() ;
for(int p = 0 ; p <= tot ; p ++ ) {
for(int c = 0 ; c < 2 ; c ++ ) {
if( !ed[tr[p][c]] && !ed[p] ) t[p][tr[p][c]] = 1 ;
}
}
f[0] = 1 ;
while( n ) {
if( n & 1 ) Mul_f() ;
n = n >> 1 ;
Mul_t() ;
}
LL ans = 0 ;
for(int p = 0 ; p <= tot ; p ++ ) if( !ed[p] ) ans = ( ans + f[p] ) % mod ;
printf("%lld" , ans ) ;
return 0 ;
}
状态机类DP
Tip : 如何根据原来的转移方程构造伴随矩阵呢?实际上,只需要对 阶段 i i i 的内层循环做一次,标记上 i − 1 i-1 i−1 的每个状态对 i i i 的每个状态贡献是多少,伴随矩阵就构造成了
总结一下:伴随矩阵 T T T 的 第 i i i 行 j j j 列表示:上一个阶段的 f [ i ] f[i] f[i] 对 这个阶段 f [ j ] f[j] f[j] 的贡献
以后会反复用到这一点
与本题像,但切入点不同的 [POI2010] CHO-Hamsters
矩乘的思想是一样的,不同在于这道题对字符串的限制没有那么强,用不着 AC 自动机(因此数据范围也就开到了很大)
矩阵加速图上问题
万变不离其宗,图上问题还是要通过分析,写出节点间的 线性转移关系,但可能会结合很多很多技巧
邻接矩阵的乘方
我们来思考邻接矩阵乘方的意义:
A
=
[
0
1
0
0
0
1
0
1
0
0
0
1
0
0
1
1
1
1
0
0
0
0
0
1
0
]
A=\begin{bmatrix} 0 & 1 & 0 & 0 &0 \\ 1 & 0 & 1& 0 &0 \\ 0 & 1 &0 &0 & 1\\ 1 & 1 & 1 & 0&0 \\ 0 & 0& 0& 1 &0 \end{bmatrix}
A=
0101010110010100000100100
( i , j ) = 1 (i,j)=1 (i,j)=1 表示 i → j i \to j i→j 有边,也可以理解为 f [ i ] → f [ j ] f[i] \to f[j] f[i]→f[j] 的贡献 系数 为 1 1 1
那么当 A 1 × A 2 A_1\times A_2 A1×A2 , A 1 ( i , k ) × A 2 ( k , j ) A_1(i,k)\times A_2(k,j) A1(i,k)×A2(k,j) 的结果就贡献到 A 2 A^2 A2 的 ( i , j ) (i,j) (i,j) 上
那么我们是否可以将其理解为 对于任意 k k k ,走一步 i → k i\to k i→k 的方案数 乘上 k → j k\to j k→j 的方案数
得到 走两步, i → j i\to j i→j 的方案总数
那么推广:邻接矩阵 A A A 的 k k k 次幂就代表走 k k k 步, i → j i\to j i→j 的方案数
类似的,我们可以依据这个原理实现很多 矩乘的黑科技!!!
技巧 1:拆点
带边权,转移的式子就变成了 i − v i-v i−v 阶段转到 i i i ,不符合矩阵优化的条件
又由于边权很小,考虑拆点:
将 i i i 号点拆成 i 0 , i 1 , i 2 . . . i 9 i_0,i_1,i_2...i_9 i0,i1,i2...i9 ,分别表示 “还差多少步就到 i i i 号点了”
连边 ( x , y , v ) (x,y,v) (x,y,v) 即为 x 0 → y v − 1 x_0 \to y_{v-1} x0→yv−1,同时每个点内部要连 i k → i k − 1 i_k\to i_{k-1} ik→ik−1
那么这样就做到了 i − 1 i-1 i−1 状态 转到 i i i,直接矩乘
技巧 2:点边互换
严格的说,这不能算是矩乘中的技巧,而是图论题的一种处理方式
本题的难点在于限制 “一条边不能连着走两遍”
表面上看不好标记走过的边,深入分析会发现限制等价于 走完一条边后不能立刻走它的反边
那么 点边互换 ,以边为状态,进行转移
正边对其反边的贡献为 0
#include<bits/stdc++.h>
using namespace std ;
typedef long long LL ;
const int N = 110 ;
const LL mod = 45989 ;
// 一条边不能连着走两次,转化一下等价于不能待在一条边不走,那么把边 按照方向 看成两个节点,转化为上一题
int n , m , T , st , ed ;
int fr[2*N] , to[2*N] ;
LL as[2*N][2*N] , t[2*N][2*N] , ans ;
void Mul_s()
{
LL tmp[2*N][2*N] = {} ;
for(int i = 0 ; i < 2*m ; i ++ ) {
for(int j = 0 ; j < 2*m ; j ++ ) {
for(int k = 0 ; k < 2*m ; k ++ ) {
tmp[i][j] = ( tmp[i][j] + as[i][k]*t[k][j]%mod ) % mod ;
}
}
}
memcpy( as , tmp , sizeof as ) ;
}
void Mul_t()
{
LL tmp[2*N][2*N] = {} ;
for(int i = 0 ; i < 2*m ; i ++ ) {
for(int j = 0 ; j < 2*m ; j ++ ) {
for(int k = 0 ; k < 2*m ; k ++ ) {
tmp[i][j] = ( tmp[i][j] + t[i][k]*t[k][j]%mod ) % mod ;
}
}
}
memcpy( t , tmp , sizeof t ) ;
}
int main()
{
scanf("%d%d%d%d%d" , &n , &m , &T , &st , &ed ) ;
int x , y ;
for(int i = 1 ; i <= m ; i ++ ) {
scanf("%d%d" , &x , &y ) ;
fr[2*i-2] = to[2*i-1] = x ;
to[2*i-2] = fr[2*i-1] = y ;
}
for(int i = 0 ; i < 2*m ; i ++ ) {
for(int j = 0 ; j < 2*m ; j ++ ) {
if( i == (j^1) ) continue ; // 不能走 自己的反边
if( to[i] == fr[j] ) t[i][j] = 1 ;
}
}
T -- ;
for(int i = 0 ; i < 2*m ; i ++ ) as[i][i] = 1 ;
while( T ) {
if( T & 1 ) Mul_s() ;
T = T >> 1 ;
Mul_t() ;
}
for(int i = 0 ; i < 2*m ; i ++ ) {
for(int j = 0 ; j < 2*m ; j ++ ) {
if( fr[i] == st && to[j] == ed ) ans = ( ans + as[i][j] ) % mod ;
}
}
printf("%lld" , ans ) ;
return 0 ;
}
综合题 / bitset 优化矩乘
首先容易想到,按
d
i
d_i
di 排序,分段作 矩乘,这样保证 每一段的伴随矩阵相同,可以用快速幂加速
考虑怎么更新答案,可以想到最终答案一定是 前面的整段次数+当前段内到 n n n 的次数
对于前者就是 d i d_i di。后者我们来思考怎么维护
不妨钦定一个起点 j j j (技巧),表示第 d i d_i di 次到达的节点,那么只需要求出 只使用当前段内的边, j j j 到 n n n 的最少次数
那我们就可以想到直接跑最短路了(因为边集确定)
如何维护 j j j —— 第 d i d_i di 次能到达的点集呢? 这就要靠矩乘了。 ( i , j ) (i,j) (i,j) 表示 K K K 次, i i i 能否到 j j j(01)
可是这样来算一下复杂度: O ( m l o g d i n 3 ) O(mlogd_in^3) O(mlogdin3),爆掉了
但冷静一下发现:我们维护的矩阵元素可全是 01
由 01 想到 —— b i t s e t bitset bitset !!!那么复杂度 O ( m l o g d i n 3 w ) O(\large\frac{m logd_in^3}{w}) O(wmlogdin3),可过
#include<bits/stdc++.h>// 好题
using namespace std ;
typedef long long LL ;
const int N = 151 ;
// 基于一个性质:若当前边的集合确定,就可以 快速地求最短路
int n , m ;
struct nod { int x , y , d ; } a[N] ;
bool cmp( nod x , nod y ) { return x.d < y.d ; }
bitset<151> egh[N] , egl[N] , f , th[N] , tl[N] ;
int ans = 2e9 , dis[N][N][N] ; // 加入第 k 条边后 , i->j 的距离
void get_dis()
{
memset( dis , 0x3f , sizeof dis ) ;
for(int k = 1 ; k <= m ; k ++ ) {
for(int i = 1 ; i <= n ; i ++ ) dis[k][i][i] = 0 ;
dis[k][a[k].x][a[k].y] = 1 ;
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= n ; j ++ ) {
dis[k][i][j] = min( min( dis[k][i][j] , dis[k-1][i][j] ) , dis[k-1][i][a[k].x]+1+dis[k-1][a[k].y][j] ) ;
}
}
}
}
void mul_f()
{
bitset<151> tmp ;
for(int j = 1 ; j <= n ; j ++ ) {
if( (f&tl[j]).count() ) tmp[j] = 1 ;
}
f = tmp ;
}
void mul_t()
{
bitset<151> tmph[N] , tmpl[N] ;
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= n ; j ++ ) {
if( (th[i]&tl[j]).count() ) tmph[i][j] = 1 , tmpl[j][i] = 1 ;
}
}
for(int i = 1 ; i <= n ; i ++ ) tl[i] = tmpl[i] , th[i] = tmph[i] ;
}
bool fg ;
int main()
{
scanf("%d%d" , &n , &m ) ;
for(int i = 1 ; i <= m ; i ++ ) {
scanf("%d%d%d" , &a[i].x , &a[i].y , &a[i].d ) ;
if( !a[i].d ) fg = 1 ;
}
if( !fg ) {
printf("Impossible") ;
return 0 ;
}
sort( a+1 , a+m+1 , cmp ) ;
get_dis() ;
f[1] = 1 ;
for(int i = 1 ; i <= m ; i ++ ) {
int T = a[i].d - a[i-1].d ;
for(int i = 1 ; i <= n ; i ++ ) tl[i] = egl[i] , th[i] = egh[i] ;
while( T ) {
if( T&1 ) mul_f() ;
T = T >> 1 ;
mul_t() ;
}
for(int j = 1 ; j <= n ; j ++ ) {
if( f[j] ) {
ans = min( ans , a[i].d + dis[i][j][n] ) ;
}
}
egh[a[i].x][a[i].y] = 1 ;
egl[a[i].y][a[i].x] = 1 ;
}
if( ans < 2e9 ) printf("%d" , ans ) ;
else printf("-1") ;
return 0 ;
}
关于 bitset 优化矩乘的 mul 函数,本人的想法是省掉 k k k 维,这样的话就分别需要 T T T 矩阵行和列的信息了,因此开两个分别维护 行 / 列
然鹅还有一种常数更小的做法,留在这里,以后要是被卡常了再来学~~
void self_mul(){//自己乘
bitset< N > Tmp[N];
for(int i = 1; i <= n; i++)
for(int k = 1; k <= n; k++)
if(Mt[i][k]) Tmp[i] |= Mt[k];
for(int i = 1; i <= n; i++) Mt[i] = Tmp[i];
}
上面我们讨论的都是用矩乘维护比较简单的信息,方案数(矩阵定义) / 01
下面来看一些 对矩阵乘法进行重定义,使之维护更复杂的信息(如 Min , 路径长… )
技巧 3、4 :矩乘重定义 / 二进制预处理
首先忽略多次询问,我们思考如何维护第
K
K
K 天的 魔法值集合
从矩阵本质出发: ( i , j ) (i,j) (i,j) 表示 i i i 号点 对 下一阶段 j j j号点 的贡献
对于 T 1 T^1 T1 ,若 i , j i,j i,j 直接相连,赋值为 1 ,否则为 0
矩阵定义完了,我们来定义矩乘的规则:
第一天 ( i , k ) (i,k) (i,k) 与 第二天 ( k , j ) (k,j) (k,j) 叠加 ,得到 前两天 ( i , j ) (i,j) (i,j) (这里用 “天” 是为了便于理解,事实上等价于任意两个相邻的阶段)
显然只有 (i,k) 与 (k,j) 都为 1 时,i 的值才能顺利贡献到 j 上
贡献的含义?对本题来说,不就是让 f [ j ] f[j] f[j] 异或 上 f [ i ] f[i] f[i] 吗?
那么由于是异或,贡献肯定要么为 1 1 1 要么为 0 0 0 ,只需统计满足条件的 k k k 个数的奇偶性即可
01 矩阵,常规套路 bitset 优化后,复杂度 O ( q l o g a i n 3 w ) \large O(\frac{qloga_in^3}{w} ) O(wqlogain3)
其实已经可过了
#include<bits/stdc++.h>
using namespace std ;
#define LL int
const int N = 101 ;
int n , m , Q ;
LL a[N] , f[N] ;
// 复杂度问题,用 bitset 优化
// 复杂度 qn^3logK/W
bitset<N> tl[N] , th[N] , tml[N] , tmh[N] ;
void Mul_t()
{
bitset<N> tmph[N] , tmpl[N] ;
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= n ; j ++ ) {
if( ((th[i]&tl[j]).count())&1 ) tmph[i][j] = tmpl[j][i] = 1 ; //求合法k的个数奇偶性
}
}
for(int i = 1 ; i <= n ; i ++ ) {
tl[i] = tmpl[i] , th[i] = tmph[i] ;
}
}
void Mul_f()
{
LL tmp[N] = {} ;
for(int j = 1 ; j <= n ; j ++ ) {
for(int k = 1 ; k <= n ; k ++ ) {
if( th[k][j] ) tmp[j] ^= f[k] ;
}
}
memcpy( f , tmp , sizeof f ) ;
}
int main()
{
scanf("%d%d%d" , &n , &m , &Q ) ;
for(int i = 1 ; i <= n ; i ++ ) {
scanf("%d" , &a[i] ) ;
}
int x , y ;
for(int i = 1 ; i <= m ; i ++ ) {
scanf("%d%d" , &x , &y ) ;
th[x][y] = th[y][x] = 1 , tl[y][x] = tl[x][y] = 1 ;
}
for(int i = 1 ; i <= n ; i ++ ) tml[i] = tl[i] , tmh[i] = th[i] ;
LL K ;
while( Q -- ) {
scanf("%d" , &K ) ;
for(int i = 1 ; i <= n ; i ++ ) f[i] = a[i] ;
while( K ) {
if( K & 1 ) Mul_f() ;
K = K >> 1 ;
Mul_t() ;
}
printf("%d\n" , f[1] ) ;
for(int i = 1 ; i <= n ; i ++ ) tl[i] = tml[i] , th[i] = tmh[i] ;
}
return 0 ;
}
然鹅本题还有更快的做法
我们会发现上面的代码中很多 Mul_t 操作是重复的,许多 n 3 n^3 n3 在被浪费
所以我们考虑 预处理所有 T 2 k T^{2^k} T2k 的伴随矩阵,计算 F F F 时只需要做 l o g log log 次 n 2 n^2 n2 的矩乘运算就好啦
预处理一次 n 3 l o g K m n^3logK_m n3logKm
查询时等价于 1 × n 1\times n 1×n 与 n × n n\times n n×n 矩阵相乘 , 复杂度 q n 2 l o g K qn^2logK qn2logK
总体复杂度就是 O ( n 3 l o g K m + q n 2 l o g K ) \large O(n^3logK_m+qn^2logK) O(n3logKm+qn2logK),十分优秀
有了上面的这些技巧,我们就可以轻松切掉 NOI 2020 美食家
这道题中体现了 分段、二进制预处理、拆点、点边转换、重定义等技巧,绝对是一道好题
技巧 5 :DP转移 与 当前阶段本身有关 情况的处理
传送门
本题第一个难点在于:用什么划分阶段?
容易想到以 走的每一步,但这样与本题的 K K K 的联系不大,没法统计次数
再想,我们当然要以 施展魔法的次数 划分,那么有转移
d p [ i ] [ j ] = min k → j ( d p [ i − 1 ] [ k ] − E [ k ] [ j ] , d p [ i ] [ k ] + E [ k ] [ j ] ) dp[i][j]=\min_{k\to j}(\ \ dp[i-1][k]-E[k][j]\ \ ,\ dp[i][k]+E[k][j]\ \ ) dp[i][j]=k→jmin( dp[i−1][k]−E[k][j] , dp[i][k]+E[k][j] )
观察这个式子,我们发现考虑一个点对 ( k , j ) (k,j) (k,j), i i i 阶段的值 不全由 i − 1 i-1 i−1 阶段转来,不满足矩乘应用的条件
那么我们尝试将 同一阶段 的转移放在一起,也就是说 i i i 阶段先由 i − 1 i-1 i−1 更新,再由 i i i 自身更新,我们只需构造两个矩阵即可
又因为矩阵的映射是 可叠加 的,我们直接将这两个矩阵相乘得到的就是最终的 伴随矩阵
LL dis[N][N] , t[N][N] , f[N] ;
void floyd()
{
for(int i = 1 ; i <= n ; i ++ ) dis[i][i] = 0 , t[i][i] = 0 ; // dis 是因为 同一次内随便跑最短路,t 是因为可以空用一次
for(int k = 1 ; k <= n ; k ++ ) {
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= n ; j ++ ) {
if( i == j ) continue ;
dis[i][j] = min( dis[i][j] , dis[i][k]+dis[k][j] ) ;
}
}
}
for(int i = 1 ; i <= n ; i ++ ) f[i] = dis[1][i] ;
//mul: dis -> t 叠加映射
LL tmp[N][N] ; memset( tmp , 0x3f , sizeof tmp ) ;
for(int i = 1 ; i <= n ; i ++ ) {
for(int j = 1 ; j <= n ; j ++ ) {
for(int k = 1 ; k <= n ; k ++ ) {
tmp[i][j] = min( tmp[i][j] , t[i][k]+dis[k][j] ) ;
}
}
}
memcpy( t , tmp , sizeof t ) ;
}
int main()
{
scanf("%d%d%d" , &n , &m , &K ) ;
memset( dis , 0x3f , sizeof dis ) ;
memset( t , 0x3f , sizeof t ) ;
int x , y , v ;
for(int i = 1 ; i <= m ; i ++ ) {
scanf("%d%d%d" , &x , &y , &v ) ;
dis[x][y] = v ;
t[x][y] = -v ;
}
floyd() ;
}
通过本题我们可以得到到一个技巧:当 i i i 阶段的DP值 不仅与 i − 1 i-1 i−1 有关,还与自身阶段有关时,考虑构造多个伴随矩阵,然后 叠加映射
隐隐感觉这个技巧可以处理更复杂的情况