【学习笔记】四处瞎学东西的笔记

大概就是发现杂七杂八学来的东西没有很必要整理的就从略了。

当然有很多是以前学的结果给忘了。

慢慢更,咕咕咕。

Gauss-Jordan 消元法

大概就是考虑普通的高斯消元是把矩阵消成上三角矩阵。这样最后就免不了要比较麻烦地回代。然后高斯-约旦消元的想法是直接消成对角矩阵。具体实现大差不差,只是不再消当前行,且每次将其它行都给消一遍。这样最后出来的就会是对角矩阵。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
db f[N][N] ;

bool gauss_jordan(int L){
for (int i = 1 ; i <= L ; ++ i){
for (int j = i + 1 ; j <= L ; ++ j)
if (f[j][i] > f[i][i]) swap(f[j], f[i]) ;
if (fabs(f[i][i]) < eps) return 0 ;
for (int j = 1 ; j <= n ; ++ j){
if (j == i) continue ;
for (int k = n + 1 ; k >= i + 1 ; -- k)
f[j][k] -= f[j][i] / f[i][i] * f[i][k] ;
}
for (int i = 1 ; i <= n ; ++ i)
ans[i] = f[i][n + 1] / f[i][i] ;
}
return 1 ;
}

链分治维护dp

大概就是所谓的动态 $dp$ 。主要针对于一小部分转移比较方便的 dp 来进行维护。大概操作就是通过轻重链剖分(普通重剖)或者虚实链剖分(动态树)之类的,通过分别维护两类儿子的 $dp$ 值来支持修改。那么不难发现当且仅当这个 dp 信息只需要 up 而不需要 bottom。

然后去用 LCT 写了一发 luogu 的 P4719。之前似乎是做过这题,但当时必然是只会比着抄胡小兔的代码。然后发现轻重链剖是真的难写。LCT的代码要短好多…

考虑大概就是每个点维护两个 dp 值,$g,f$ ,分别表示「只考虑虚儿子」和「考虑了全部的儿子」时的答案。转移比较简单,写成矩阵的形式就是

其中矩阵乘法的定义为 $\max$ 包含 $+$ 。

然后记录点实现上的细节吧:

0、发现 $g$ 不需要用矩阵来维护,直接记录两个值就好了。

1、考虑在 $access$ 的时候动态维护 $g$ 。因为虚实不断变化。方法是加上新的右子树并且减掉割下来的右子树。

2、考虑在 push_up 的时候维护 $f$ 。需要注意的是由于 splay 里面深度单调,而这种转移是需要按照深度来转移的,所以需要 f(rc(x)) * trans(x) * f(lc(x)) 来做。

3、有些实现是需要判左右儿子是否为空。其实可以不用。只要把 $f(0)$ 写成单位矩阵就可以了。

4、这题用 LCT 实现其实是有点大材小用的。因为根本不需要 LCT 的其他什么花里胡哨的操作。

然后去搞了搞加强版,加了 cyj 的 I/O 之后,靠 rp 爆了十几发(可能不到十发,但反正就是很多发) OJ 就过了/cy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
int lans ;
int n, m ;
int base[N] ;

vint E[N] ;

struct mat{
int m[2][2] ;
il mat (int a = -P, int b = -P, int c = -P, int d = -P){
m[0][0] = a ; m[0][1] = b ;
m[1][0] = c ; m[1][1] = d ;
}
il mat friend operator * (const mat & a, const mat & b){
mat c ; //c.reset() ;
c.m[0][1] = max(a.m[0][0] + b.m[0][1], a.m[0][1]+ b.m[1][1]) ;
c.m[1][1] = max(a.m[1][0] + b.m[0][1], a.m[1][1]+ b.m[1][1]) ;
c.m[0][0] = max(a.m[0][0] + b.m[0][0], a.m[0][1]+ b.m[1][0]) ;
c.m[1][0] = max(a.m[1][0] + b.m[0][0], a.m[1][1]+ b.m[1][0]) ;
return c ;
}
} ;

#define f(x) s[x].f
#define fa(x) s[x].fa
#define g_0(x) s[x].g[0]
#define g_1(x) s[x].g[1]
#define lc(x) s[x].son[0]
#define rc(x) s[x].son[1]

struct lct{
mat f ;
int fa ;
int g[2] ;
int son[2] ;
}s[N * 3] ;

il bool w_k(int x){
return (rc(fa(x)) == x) ;
}
il bool notroot(int x){
return (lc(fa(x)) == x || rc(fa(x)) == x) ;
}
il void _up(int x){
f(x) = f(rc(x)) * mat(g_0(x), g_1(x), g_0(x), -P) * f(lc(x)) ; //2
}
il void rotate(int x){
register int da = fa(x) ;
register int dada = fa(da) ;
register bool w = w_k(x), ww = w_k(da) ;
if (notroot(da))
s[dada].son[ww] = x ;
fa(x) = dada ;
fa(s[x].son[w ^ 1]) = da ;
s[da].son[w] = s[x].son[w ^ 1] ;
s[x].son[w ^ 1] = da ; fa(da) = x ;
_up(da) ;
}
void splay(int x){
while (notroot(x)){
if (notroot(fa(x)))
rotate(w_k(fa(x)) == w_k(x) ? fa(x) : x) ;
rotate(x) ;
}
_up(x) ;
}
void access(int x){
register int y = 0 ;
for ( ; x ; x = fa(y = x)){
splay(x) ;
g_1(x) += f(rc(x)).m[0][0] - f(y).m[0][0] ;
g_0(x) -= max(f(y).m[0][0], f(y).m[0][1]) ;
g_0(x) += max(f(rc(x)).m[0][0], f(rc(x)).m[0][1]) ;
rc(x) = y ; _up(x) ;
}
}

void prelude(int x, int dad){
g_1(x) = base[x] ;
for (register auto k : E[x]){
if (k != dad){
fa(k) = x ;
prelude(k, x) ;
g_1(x) += g_0(k) ;
g_0(x) += max(g_0(k), g_1(k)) ;
}
}
f(x) = mat(g_0(x), g_1(x), g_0(x), -P) ;
}

int main(){
gi(n), gi(m) ; int x, y ;
for (register int i = 1 ; i <= n ; ++ i) gi(base[i]) ;
for (register int i = 1 ; i < n ; ++ i)
gi(x), gi(y), E[x].p_b(y), E[y].p_b(x) ;
prelude(1, 0) ; f(0) = mat(0, -P, -P, 0) ;
while (m --){
gi(x) ; x ^= lans ;
gi(y) ; access(x) ; splay(x) ;
g_1(x) = g_1(x) + y - base[x] ;
_up(x) ; splay(1) ; base[x] = y ;
lans = max(f(1).m[0][0], f(1).m[0][1]) ;
print(lans), pc('\n');
}
return 0 ;
}

Z-Algorithm

大概就是俩用途:

1、给定一个串 $s$ 。求 $s$ 的所有后缀与 $s$ 的 LCP。

2、给定一个串 $s$ 和一个串 $t$,求 $s$ 的所有后缀与 $t$ 的 LCP。

本质上跟 ManaCher 十分相似。大概就是求的时候维护一段 LCP 最靠右的子段。然后根据相关信息推一下就好了。

1
2
3
4
5
6
7
8
for (int i = 2 ; i <= n ; ++ i){
if (r < i) Z[i] = 1 ;
else Z[i] = min(Z[i - l + 1], r - i + 1) ;
//可以理解为 s[1...i-l+1] = s[l...i],i-l+1向后匹配多少就是i向后可以匹配多少
while (t[Z[i]] == t[i + Z[i] - 1]) ++ Z[i] ; //统计新的贡献
if (i + Z[i] - 1 > r) r = i + Z[i] - 1, l = i ; //更新LCP最靠右的子段
}
Z[1] = n ;

需要注意的是 $Z[1]$ 一般不是良定义的。如果需要最后可以赋值为串长。

当然如果对于第二问,可以采取在 $s$ 后面接一个不属于字符集的分隔符再拼上 $t$ 来做。但是也可以另开一个新的函数来做这个。但注意到这个过程需要提前做一遍 $s$ 的 $Z$ ,因为需要通过 $Z$ 来初始化 $Pz$ 。

这个地方有个小细节 。直接拼一个 # 啥事没有。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
namespace Z_F{
int L ;
int Q ;
int l, r ;
int Z[N] ;
int Pz[N] ;
void gao(int bg, int L){
for (int i = bg ; i < bg + L ; ++ i) Z[i] -- ;
Z[bg] = L ;
}
int get_Z(char *s, int bg, int oo = 0){
L = strlen(s + bg) ; l = bg, r = 0 ;
for (int i = bg + 1 ; i < bg + L ; ++ i){
// cout << i << '\n' ;
if (r < i) Z[i] = 1 ;
else Z[i] = min(Z[i - l + 1], r - i + 1) ;
while (s[Z[i]] == s[i + Z[i] - 1]) ++ Z[i] ;
if (r < i + Z[i] - 1) r = i + Z[i] - 1, l = i ;
}
if (!oo) gao(bg, L) ;
//debug(Z, 1, L) ;
return L ;
}
int exKMP(char *s, int bg, char *t, int gb){
int q = get_Z(s, bg, 1) ; s[q + 1] = '#' ;
L = strlen(t + gb) ; l = gb, r = 0 ;
for (int i = gb ; i < gb + L ; ++ i){
if (r < i) Pz[i] = 1 ;
else Pz[i] = min(Z[i - l + 1], r - i + 1) ;
while (s[Pz[i]] == t[i + Pz[i] - 1]) ++ Pz[i] ;
if (r < i + Pz[i] - 1) r = i + Pz[i] - 1, l = i ;
}
for (int i = gb ; i < gb + L ; ++ i) Pz[i] -- ;
return L ;
}
}

CF126B Password

给一个串,找到一个既是前缀又是后缀并且还是非前后缀子串的最长子串。

$n \leq 10^6$。

被降智了好久…一开始写了一堆奇奇怪怪的东西发现假的很…

大概就是先考虑如何快速找出既是前缀又是后缀的。这个可以 KMP 解决,或者用 Z 来枚举每个后缀的 $i+z_i-1$ 是不是等于 $n$ 。然后他还要求在中间出现过。发现我们只关心最大值且不关心中间的串的具体情况,于是就维护一个到现在为止中间的串能匹配上的最长前缀即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
using namespace Z_F ;

int main(){
scanf("%s", s + 1) ;
n = get_Z(s, 1) ; res = 0 ;
for (int i = 2 ; i <= n ; ++ i){
if (i + Z[i] - 1 == n && Z[i] <= mx)
if (res < Z[i]) res = Z[i], p = i ;
chkmax(mx, Z[i]) ;
}
if (!res)
return puts("Just a legend"), 0 ;
for (int i = p ; i <= p + res - 1 ; ++ i)
putchar(s[i]) ; return puts(""), 0 ;
}

CF432D

给你一个长度为 $n$ 的长字符串,“完美子串”既是它的前缀也是它的后缀,求“完美子串”的个数且统计这些子串的在长字符串中出现的次数。$n\leq 10^6$ 。

先把所有的 border 判出来。统计出现次数可以直接开一个桶记录一下每个 $Z_i$ 值出现了多少次,然后后缀和一下就是答案。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
using namespace Z_F ;

int main(){
scanf("%s", s + 1) ;
n = get_Z(s, 1) ; cnt = 0 ;
for (int i = 1 ; i <= n ; ++ i){
ans[Z[i]] ++ ;
if (i + Z[i] - 1 == n)
res[++ cnt] = Z[i] ;
}
for (int i = n ; i >= 1 ; -- i)
ans[i] += ans[i + 1] ;
cout << cnt << '\n' ;
sort(res + 1, res + cnt + 1) ;
for (int i = 1 ; i <= cnt ; ++ i)
printf("%d %d\n", res[i], ans[res[i]]) ;
}

置换

置换本质上是一个排列到另一个排列的一一映射。置换的乘法就对应于函数复合。

同时一个置换也可以分解成许多轮换,比如 $\begin{pmatrix}1&2&3&4&5\\3&4&5&2&1\end{pmatrix}$ 其中就包含了 $(1\to 3,3\to 5,5\to1)$ 和 $(2\to 4,4\to 2)$ 这两个循环移位。

嗯…以上是抄书。然后考虑 $\sf Burnside$ 引理,即对于一个定义在集合 $\sf G$ 上的置换群 $\rm F$ 而言,设 $f(s)$ 为在群内某个置换 $s$ 本质不变的 $\sf G$ 内元素个数(即「不动点个数」),那么等价类个数为:

但是这个不动点个数有时候会比较难求。而 $\sf Polya$ 定理则在陈述这么一件事:考虑对集合 G 作用的元素进行 $k-$染色,那么某个置换 $s$ 下不动点的个数就是 $k^{\zeta(s)}$,其中 $\zeta(s)$ 是 $s$ 可以分解成的轮换的个数。然后就变成了:

也就是如果「本质不同」的方式使用染色来定义的话,就可以用 $\sf Polya$ 来优化计算过程。

UVA 10294 Arif in Dhaka

求 $n$ 点 $t$ 染色的项链数和手镯数。

项链数比较好求,考虑置换群内的元素分别是顺时针旋转 $0,1,2\cdots n-1$ 个元素,那么假设旋转了 $i$ 个元素,此时每个轮换内会有 $\dfrac{\mathrm{lcm}(n,i)}{i}$ 个元素(总距离除以步长), 那么也就是说总共会有 $\dfrac{n}{\dfrac{\mathrm{lcm}(n,i)}{i}}$ 个轮换,也就是 $\gcd(n,i)$ 个轮换。

手镯数需要加上对称。考虑当 $n$ 为奇数,一共有 $n$ 条对称轴,每条对称轴划分为 $\dfrac{n-1}{2}$ 个长度为 $2$ 的轮换和 $1$ 个长度为 $1$ 的轮换,那么轮换数就是 $\dfrac{n+1}{2}$ ;当 $n$ 为偶数时,会有两种对称方式,分别有 $\dfrac{n+1}{2}$ 和 $\dfrac{n}{2}$ 种轮换。

哦,小细节,由于本质上手镯里面每条对称轴都是一个单独的置换,所以应该乘个系数。

NWERC2006 Leonado’s Notebook

给定一个置换 $A$ ,求是否存在一个置换 $B$ 使得 $B^2=A$ 。

考虑将置换分解成轮换。注意到如果两个轮换彼此之间不交,那么他们之间的乘法是具有交换律的。同时考虑对于一个奇数大小的轮换,将其平方之后依旧是一个轮换,偶数大小的平方之后则可以分解成两个轮换。所以对于给定的置换,长度为奇数的轮换可以不用管,如果长度为偶数的轮换有奇数个就会 gg。

具体证明的话,大概是全部相同大小的置换都可以属于同一个置换群,之后就可以对于每个轮换重标号来证明了。

不过轮换的乘法比较…麻烦。比如 $(1~2~3)(1~2~3)=(1~3~2)$ ,$(1~2~3~4)(1~2~3~4)=(3~4~1~2)=(1,3)(2,4)$ 。感觉自己一碰到这种置换的就 gg。

UVA11077 Find the Permutations

给定一个排列,每次可以交换两个数,求至少需要交换 $k$ 次才能让排列递增的排列数。

发现可以将给定的这个排列看做 $(1,2,3\cdots n)$ 的一个置换。那么考虑将其轮换分解之后,对于一个大小为 $k$ 的轮换,不难发现需要交换 $k-1$ 次。于是可以 dp,$f_{i,j}$ 表示长度为 $i$ 的排列至少交换 $j$ 次的排列方案数。那么每次考虑新加进来的元素是否放到之前的某个轮换里,故有转移

……发现这就是第一类斯特林数。原因是考虑如果有 $p$ 个轮换,那么至少需要交换 $n-p$ 次。而第一类斯特林数的意义就包括 $n$ 元置换可分解为 $k$ 个独立的轮换(轮换)的个数,所以答案就应该是 $s_1(n,n-k)$。而用上面这种推法推出来的东西恰与第一类斯特林数水平对称。

CERC/SWERC2005 Pixel Shuffle

给定一些关于 $n\times n$ 的黑白方格的奇怪操作,求至少多少次操作可以变换会原来的样子。

$n\leq 1024,k\leq10$ 。

发现本质上是给定一个置换 $A$,求一个最小的 $m$ 使得 $A^m=I$ 。然后大概就是考虑轮换的性质。对于 $A$ 中的一个长度为 $p$ 的轮换,如果想要全等,就必须要让 $p|m$ 。所以最后答案就是所有轮换大小的最小公倍数。

找轮换的话大概就是对每个位置编一个号,找就好了。UVA 属实神必 OJ,和过了的程序拍了一万组交上去依旧 WAWAWA。遇到这种情况那必然是不管了。

发现中国 OJ 这么多题,刘汝佳非要找 UVA 写书。真是 mdzz。

不过倒是过了主席的魔改版,这里留个链接:

「My-Code」
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
int T, n, m ;
int tmp[N][N] ;
int base[N][N] ;
int dododo[99] ;
char opt[5000] ;
char tpo[5000] ;
bool vis[N * N] ;
unsigned ll ans ;

void do_do(int mk){
if (abs(mk) == 1) return ;
if (mk == 2){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
base[i][j] = tmp[j][n - i + 1] ;
}
if (mk == -2){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
base[j][n - i + 1] = tmp[i][j] ;
}
if (abs(mk) == 3){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n / 2 ; ++ j)
base[i][j] = tmp[i][n - j + 1] ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = n / 2 + 1 ; j <= n ; ++ j)
base[i][j] = tmp[i][n - j + 1] ;
}
if (abs(mk) == 4){
for (int i = n / 2 + 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = n / 2 + 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n / 2 ; ++ j)
base[i][j] = tmp[i][n - j + 1] ;
for (int i = n / 2 + 1 ; i <= n ; ++ i)
for (int j = n / 2 + 1 ; j <= n ; ++ j)
base[i][j] = tmp[i][n - j + 1] ;
}
if (abs(mk) == 5){
for (int i = n / 2 + 1 ; i <= n / 4 * 3 ; ++ i)
swap(base[i], base[n - (i - n / 2) + 1]) ;
}
if (mk == 6){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1, j = 1 ; i <= n / 2 ; ++ i, j += 2)
memcpy(base[i], tmp[j], sizeof(tmp[j])) ;
for (int i = n / 2 + 1, j = 2 ; i <= n ; ++ i, j += 2)
memcpy(base[i], tmp[j], sizeof(tmp[j])) ;
}
if (mk == -6){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1, j = 1 ; i <= n / 2 ; ++ i, j += 2)
memcpy(base[j], tmp[i], sizeof(tmp[i])) ;
for (int i = n / 2 + 1, j = 2 ; i <= n ; ++ i, j += 2)
memcpy(base[j], tmp[i], sizeof(tmp[i])) ;
}
if (mk == 7){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1 ; i <= n ; ++ i){
int l = n / 2, r = n / 2 ;
if (i & 1) l = 0, r = 0 ;
for (int j = 1 ; j <= n ; ++ j){
if (i & 1){
if (j & 1) base[i][j] = tmp[i][++ l] ;
else base[i][j] = tmp[i + 1][++ r] ;
}
else{
if (j & 1) base[i][j] = tmp[i - 1][++ l] ;
else base[i][j] = tmp[i][++ r] ;
}
}
}
}
if (mk == -7){
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
tmp[i][j] = base[i][j] ;
for (int i = 1 ; i <= n ; ++ i){
int l = n / 2, r = n / 2 ;
if (i & 1) l = 0, r = 0 ;
for (int j = 1 ; j <= n ; ++ j){
if (i & 1){
if (j & 1) base[i][++ l] = tmp[i][j] ;
else base[i + 1][++ r] = tmp[i][j];
}
else{
if (j & 1) base[i - 1][++ l] = tmp[i][j] ;
else base[i][++ r] = tmp[i][j] ;
}
}
}
}
/*
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
cout << base[i][j] << " \n"[j == n] ;
puts(" - - - - - - - - - - - - - ") ;
*/
}
bool check(){
int owo = 0 ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j)
if (base[i][j] != (++ owo)) return 0 ;
return 1 ;
}
int tot, cnt ;
unsigned ll res[N * N] ;

unsigned ll gcd(unsigned ll x, unsigned ll y){
return y ? gcd(y, x % y) : x ;
}
unsigned ll lcm(unsigned ll x, unsigned ll y){
return x * y / gcd(x, y) ;
}
void record(char *s){
++ cnt ;
if (s[0] == '-')
dododo[cnt] = -1, ++ s ; else dododo[cnt] = 1 ;
if (strcmp(s, "tor") == 0) dododo[cnt] *= 2 ;
else if (strcmp(s, "mys") == 0) dododo[cnt] *= 3 ;
else if (strcmp(s, "myshb") == 0) dododo[cnt] *= 4 ;
else if (strcmp(s, "mysvb") == 0) dododo[cnt] *= 5 ;
else if (strcmp(s, "vid") == 0) dododo[cnt] *= 6 ;
else if (strcmp(s, "xim") == 0) dododo[cnt] *= 7 ;
}
signed main(){
freopen("gen.in", "r", stdin);
freopen("aa.out", "w", stdout);
cin >> T ;
while (T --){
scanf("%d", &n) ; getchar() ;
memset(vis, 0, sizeof(vis)) ;
memset(base, 0, sizeof(base)) ;
ans = 1, tot = cnt = 0 ;
fgets(opt, (sizeof opt / sizeof opt[0]), stdin);
for (int i = 1, k = 0 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j) base[i][j] = ++ k ;
m = strlen(opt) ; int ps = 0 ;
for (int i = m - 1 ; i >= 0 ; -- i)
if (opt[i] != ' ') tpo[ps ++] = opt[i] ;
else tpo[ps] = '\0', record(tpo), ps = 0 ;
tpo[ps] = '\0' ; record(tpo) ;
for (int i = 1 ; i <= cnt ; ++ i) do_do(dododo[i]) ;
for (int i = 1 ; i <= n ; ++ i)
for (int j = 1 ; j <= n ; ++ j){
if (vis[base[i][j]]) continue ;
int x = i, y = j, z, o = 0 ;
while (!vis[base[x][y]]){
vis[base[x][y]] = 1 ; z = base[x][y] ;
x = z / n + (!(bool)(z % n == 0)) ;
y = z % n ; if (!y) y = n ; ++ o ;
// cout << x << " " << y << '\n' ;
}
res[++ tot] = o ;
}
// debug(res, 1, tot) ;
for (int i = 1 ; i <= tot ; ++ i) ans = lcm(ans, res[i]) ;
printf("%llu\n", tot >= 1 ? ans : 0) ;
if (T) puts("") ;
}
}
/*
1
8
3
-2 -6 -7
*/

最大流-ISAP

大概就是设 $d_x$ 表示残量网络上 $x$ 到汇点 $t$ 的距离。那么为了只在残量网络上进行增广,所以每次只能沿着 $d_y+1=d_x$ 的 $(x,y)$ 增广。那么普通的 Dinic 的想法是每次求出一张增广网,然后对这张增广网进行多路增广,为了保证增广路之间不相交所以引入了高度的概念。ISAP 也差不多。只不过是在找增广路的过程中动态地修改标号。即考虑在增广的过程中,可能会有一些边 $(x,y)$ 从残量网络中抹掉,那么如果这条边用来传递标号的边,即 $d_y+1=d_x$,就需要找到另一条边 $(x,z)$ 并给 $x$ 重标号一个高度为 $d_z+1$ 。这个过程可以一开始用一个 bfs 来预处理,不难看出应该倒着预处理,因为在某条边被删掉之后,从 fromto 比从 tofrom 要简单得多,当然这也是为什么 ISAP 对高度的定义和 Dinic 相反的原因(Dinic 其实无所谓怎么定义,因为横竖都要重新生成增广网)。当源点 $s$ 无处增广时算法停止。

那么有一个优化叫做「gap 优化」。大概就是考虑一条增广路上的点的标号必然是单调的。那么如果出现了某个断层,可以知道不再存在增广路。于是可以直接结束。

实现方面,可以记一个当前弧来保证只在残量网络上增广。

最终可以知道 ISAP 本质上是在和 Dinic 做同一件事,因为 ISAP 的重标号一定会是标定最小的那个 $z$ 的 $h_z+1$ ,跟求增广网本质相同。所以复杂度上界也是 $O(n^2m)$ 的。

哦,顺便提一嘴。似乎有证明当容量均为 $1$ 时,复杂度有着更加优秀的上界 $O(m\sqrt[3] {n^2}+m\sqrt m)$ 。

「My-Code」
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
namespace ISAP{
struct edge{
int from ;
int flow ;
int cap ;
int to ;
edge (int a = 0, int b = 0, int c = 0){
from = a, to = b, cap = c, flow = 0 ;
}
} e[N] ;
int n ;
int cnt ;
int ans ;
int S, T ;
int h, t ;
int que[N] ;
int pre[N] ;
int hgt[N] ;
int cur[N] ;
int gap[N] ;
vector <int> E[N] ;
void add(int x, int y, int c){
e[++ cnt] = edge(x, y, c) ; E[x].push_back(cnt) ;
e[++ cnt] = edge(y, x, 0) ; E[y].push_back(cnt) ;
}
void reset(int x, int y, int z){
memset(cur, 0, sizeof(cur)) ;
memset(hgt, 0, sizeof(hgt)) ;
memset(gap, 0, sizeof(gap)) ;
S = x, T = y, n = z ; cnt = -1 ;
for (int i = 1 ; i <= n ; ++ i) E[i].clear() ;
}
void preBFS(){
h = 1, t = 0 ;
que[++ t] = T ;
while (h <= t){
int x = que[h ++] ;
for (auto k : E[x])
if (!hgt[e[k].to] && e[k].to != T)
hgt[que[++ t] = e[k].to] = hgt[x] + 1 ;
}
for (int i = 1 ; i <= n ; ++ i) gap[hgt[i]] ++ ;
// debug(hgt, 1, n) ;
}
void Augment(){
int x = T ;
int f = Inf ;
while (x != S){
// debug(x) ;
f = min(f, e[pre[x]].cap - e[pre[x]].flow) ;
x = e[pre[x]].from ;
}
ans += f ; x = T ;
while (x != S){
e[pre[x]].flow += f ;
e[pre[x] ^ 1].flow -= f ;
x = e[pre[x]].from ;//1
}
}
bool Advance(int x, int &y){
bool ret = 0 ;
for (int k = cur[x] ; k < E[x].size() ; ++ k){
int p = E[x][k] ;
// debug(e[p].to) ;//0
if (e[p].flow >= e[p].cap) continue ;
if (hgt[e[p].to] + 1 != hgt[x]) continue ;
pre[e[p].to] = p ; cur[x] = k ;
y = e[p].to ; ret = 1 ; break ;
}
return ret ;
}
void isap(){
preBFS() ;
int x = S ;
while (hgt[S] < n){
// debug(x) ;
if (x == T)
Augment(), x = S ;
if (!Advance(x, x)){
int h = n - 1 ;
for (auto k : E[x]){
if (e[k].flow < e[k].cap)
h = min(h, hgt[e[k].to]) ;
}
if (!-- gap[hgt[x]]) break ;
gap[hgt[x] = h + 1] ++ ; cur[x] = 0 ;
if (x != S) x = e[pre[x]].from ;
}
// break ;
}
}
}

int n, m, s, t ;

KM 算法

呃,这算法真的是学了很久。以下默认求的是最大权匹配。

是用来求「完美匹配」基础上的「最佳权值匹配」。「完美匹配」大概就是说对于一张二分图 $G=\{V_1,V_2,E\}$,设 $|V_1|\lt |V_2|$,对于某个匹配 $M’$ 如果可以使得 $V_1$ 中的点均可以被匹配到,那么就称 $M’$ 为 $G$ 的一个完美匹配。

然后 KM 算法大概是对于每个点设置一个顶标 $p_i$ 或者 $q_i$ (区分左右部,左部记为 $p$,右部记为 $q$),保证 $p_i+q_j\geq w(i,j)$ 。定义 $G$ 的一个相等子图 $G_0=\{V_1,V_2,E_0\}$,需要有 $\forall e(u,v)\in E_0,p_u+q_v=w(u,v)$。那么不难知道如果相等子图 $G_0$ 存在完美匹配 $M_0$,那么 $M_0$ 就是最大权完美匹配。这一点不难证明,只需要注意 $G_0$ 同样包含所有 $G$ 中的顶点即可。

那么KM 算法本质上就是在不断扩大相等子图的范围实现增广,直到存在完美匹配。假设当前 $G_0$ 中的点集分别为 $\rm \{S,T\}$,那么考虑为了扩大相等子图,需要修改顶标。

具体的,让 $\rm S$ 中的所有点顶标增加一个权值 $a$,让 $T$ 中所有点的顶标减小一个权值 $a$ ,那么可知 $G_0$ 中原来的点依旧属于 $G_0$ ,同时对于某个点对 $(u,v),u\in \mathrm{S},v\in V_2\setminus T$,可以知道会有 $p_u+q_v$ 变小,那么就会出现可能把 $(u,v)$ 也加入相等子图的情况。那么不难知道为了保证顶标合法且有新的边加入相等子图,$a$ 需要满足 $a=\min\{p_u+q_v-w(u,v)\},u\in \mathrm{S},v\in V_2\setminus T$ 。考虑加入了这条边之后,要么 $v$ 是未盖点,要么 $v$ 是匹配点。就可以直接进行类似匈牙利算法的过程。考虑对于每个点,至多会增广 $n$ 次,每次增广是 $O(n)$ 的,修改顶标可以 $O(n)$ 查找,那么复杂度就应该是 $O(n^3)$(虽然我一开始写的是暴力找 $a$ 的 $O(n^4)$,但实际上可以记一个 slack)解决这个问题。

看上去很棒?实际上大多数人实现的 dfs 都是 $O(n^4)$ 的。原因在于每次找增广路并不是 $O(n)$ 而是 $O(m)$ 的。为了解决这个问题就要避免类似 dfs 的过程里反复遍历交错树的情况。于是考虑 bfs ,这样聚合分析一下每次找增广路就会是严格 $O(m)$ 的。

「My-Code」
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
//Km3=KM of O(n^3), Km4=KM of O(n^4)
namespace KM{
const int N = 510 ;
int n ;
ll ans ;
ll lack ;
ll lv[N] ;
ll rv[N] ;
int h, t ;
int q[N] ;
int pre[N] ;
int matl[N] ;
int matr[N] ;
ll slack[N] ;
bool visl[N] ;
bool visr[N] ;
ll val[N][N] ;
void reset(int x){
n = x ;
for (int i = 1 ; i <= n ; ++ i){
matl[i] = matr[i] = 0 ;
pre[i] = lv[i] = rv[i] = 0 ;
for (int j = 1 ; j <= n ; ++ j)
lv[i] = max(lv[i], val[i][j]) ;
}
}
void upd(){
for (int i = 1 ; i <= n ; ++ i){
if (visl[i]) lv[i] -= lack ;
if (visr[i]) rv[i] += lack ;
else slack[i] -= lack ;
}
}
bool do_match4(int x){
visl[x] = 1 ;
for (int y = 1 ; y <= n ; ++ y)
if (!visr[y]){
ll t = lv[x] + rv[y] - val[x][y] ;
if (!t){
visr[y] = 1 ;
if (!matl[y] || do_match4(matl[y])){
matl[y] = x ; matr[x] = y ; return 1 ;
}
}
else slack[y] = min(slack[y], t) ;
}
return 0 ;
}
void Km4(){
for (int i = 1 ; i <= n ; ++ i){
for (int j = 1 ; j <= n ; ++ j)
slack[j] = 1ll << 60 ;
while (1){
lack = 1ll << 60 ;
fill(visl + 1, visl + n + 1, 0) ;
fill(visr + 1, visr + n + 1, 0) ;
if (do_match4(i)) break ;
for (int i = 1 ; i <= n ; ++ i)
if (!visr[i]) lack = min(lack, slack[i]) ;
upd() ;
}
}
for (int i = 1 ; i <= n ; ++ i)
ans += lv[i], ans += rv[i] ;
}
void match_back(int x){
int z ;
while (x){
z = x ;
matl[x] = pre[x] ;
x = matr[pre[x]] ;//2
matr[pre[z]] = z ;
}
}
void do_match3(int x){
q[++ t] = x ;
while (1){
while (h <= t){
int z = q[h ++] ;
for (int y = 1 ; y <= n ; ++ y){
if (visr[y]) continue ;
ll v = lv[z] + rv[y] - val[z][y] ;//3
if (v > slack[y]) continue ; pre[y] = z ;
if (!v){
visr[y] = 1 ;
if (!matl[y]) return match_back(y) ;
else visl[q[++ t] = matl[y]] = 1 ;//5
}
else slack[y] = min(slack[y], v) ;
}
}
lack = 1ll << 60 ; int z = 0 ;
for (int i = 1 ; i <= n ; ++ i)
if (!visr[i])
if (lack > slack[i])
lack = slack[i], z = i ;
upd() ; //4
//debug(z) ;
if (!matl[z]) return match_back(z) ;
visr[z] = visl[matl[z]] = 1 ; q[++ t] = matl[z] ;
}
}
void Km3(){
for (int i = 1 ; i <= n ; ++ i){
memset(visl, 0, sizeof(visl)) ;
memset(visr, 0, sizeof(visr)) ;
memset(slack, 127, sizeof(slack)) ;
h = 1 ; t = 0 ; visl[i] = 1 ; do_match3(i) ;
}
for (int i = 1 ; i <= n ; ++ i)
ans += lv[i], ans += rv[i] ;
}
}

int nl, nr, n, m ;

int main(){
int x, y, z ;
cin >> n >> m ;
memset(KM :: val, -127, sizeof(KM :: val)) ;//1
for (int i = 1 ; i <= m ; ++ i){
x = qr(), y = qr(), z = qr() ;
KM :: val[x][y] = z ;
}
KM :: reset(n) ;
KM :: Km3() ; cout << KM :: ans << '\n' ;
for (int i = 1 ; i <= n ; ++ i)
cout << KM :: matl[i] << " " ; puts("") ;
return 0 ;
}

upd: 重写了一遍,易错的地方都在板子里标出来了。

ZKW 费用流

比较神奇的费用流。考虑一般费用流的做法是改良版的EK,即每次对着残量网络跑一遍最短路,然后对着一条边都满足 $d_v=d_u+val(u,v)$ 的路径进行增广。单路增广必然会很慢,于是来考虑一下为什么必须要单路增广:考虑最短路算法保证在算法结束后,对于残量网络上的每个节点 $x$ 的所有入点 $y$ 有 $d_y+val(y,x)\geq d_x$,且至少存在一个 $y$ 可以使等号成立。那么在残量网络上找出一条增广路之后,与普通的最大流(边权为 $1$)不同的是,剩余的参量网络对于每个 $x$ 不一定会继续存在一个 $y$ 使得等号成立,所以这就需要重新找增广路。

而 $zkw$ 费用流的思想则是如何不用每次跑一遍最短路来增广,需要引入顶标概念,此处的顶标就是从 $\rm S$ 开始的最短路长度 $d_i$,即有

那么发现跟 KM 的本质极其相似:只增广等号成立的点,那么为了可以使得等号成立就需要不断扩大相等子图,扩大的方式也如出一辙,只需要让所有增广过的点的顶标减去一个 $ \rm lack$ ,而这个 $\rm lack$ 就理应是

这样就可以继续增广了。

然后写法方面需要注意一些细节。我所知道的应该是有两种写法:

1、只依赖自动调整。这个写法有个很关键的依赖,就是 $t_{\rm S}$ 单调不增。换句话说,每次调整 $t_{\rm S}$ 必须要满足存在某条边可以增广,但是此时 $t_{\rm S}$ 并不是单纯地跟某条边的边权挂钩,它只是一个在不断减小的量。所以不难理解,如果存在某个时刻增广到了 $\rm T$ ,那么就一定有 $\sum c=-t_{\rm S}$ ,因为无论什么时候,$\rm S$ 都会在扩展完毕的集合里。

于是就可以知道,$t_{\rm S}$ 的初值应该赋为 $\rm S$ 所有出边里边权最大的那一条。或者可以直接赋值为 $-\rm Inf$ ,让其自己调整。

2、大概是个优化?可以考虑在最开始的时候跑一遍最短路。设增广路上的点按顺序分别是 $\rm S,u_1,u_2\cdots u_m,T$ 那么此时每条增广路的 $cost$ 之和就是

化简一下可以得到就是 $t_{\rm T}-t_{\rm S}$ 。

其实第一个 $case$ 最后的结果本质上也是 $t_T-t_S$ ,然而因为 $T$ 的 $t$ 始终为 $0$,所以可以忽略。看怎么理解了,也可以直接推出 $-t_S$ 。

事实证明第二种确实比第一种要快。

「My-Code」
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
void chkmin(int &x, int y){
x = x > y ? y : x ;
}

namespace zkw{
#define to(k) e[k].to
#define cap(k) e[k].cap
#define cost(k) e[k].cost
#define from(k) e[k].from
#define flow(k) e[k].flow
const int N = 200010 ;

int cnt ;
struct edge{
int from ;
int flow ;
int cost ;
int cap ;
int to ;
edge (int a = 0, int b = 0, int c = 0, int d = 0){
from = a, to = b, cap = c, cost = d, flow = 0 ;
}
} e[N] ;
int mflow ;
int mcost ;
int val[N] ;
bool vis[N] ;
int n, S, T ;
vector <int> E[N] ;
void reset(int x, int y, int z){
memset(val, 0, sizeof(val)) ;
S = x, T = y, n = z ; cnt = -1 ;
for (int i = 1 ; i <= n ; ++ i) E[i].clear() ;
}
void add(int x, int y, int c, int d){
// debug(x, ' '), debug(y, ' '), debug(c, ' '), debug(d) ;
e[++ cnt] = edge(x, y, c, d) ; E[x].push_back(cnt) ;
e[++ cnt] = edge(y, x, 0, -d) ; E[y].push_back(cnt) ;
}
int h, t ;
int que[N] ;
void spfa(){
memset(val, 127, sizeof(val)) ;
val[que[h = t = 1] = S] = 0, vis[S] = 1 ;
while (h <= t){
int x = que[h ++] ;
vis[x] = 0 ; //debug(x) ;
for (auto k : E[x]){
if (!cap(k)) continue ;
if (val[to(k)] > val[x] + cost(k)){
val[to(k)] = val[x] + cost(k) ;
if (!vis[to(k)])
vis[que[++ t] = to(k)] = 1 ;
}
}
}
}
int advance(int x, int f){
vis[x] = 1 ;
if (x == T){
// debug(f, ' '), debug(- val[S]) ;
mcost += (val[T] - val[S]) * f ;
mflow += f ; return f ;
}
int fl = 0 ;
for (auto k : E[x])
if (flow(k) < cap(k))
if (!vis[to(k)] && val[to(k)] == val[x] + cost(k))
if ((fl = advance(to(k), min(f, cap(k) - flow(k)))))
return flow(k) += fl, flow(k ^ 1) -= fl, fl ;
return 0 ;
}
bool modify(){
int lack = 1 << 30 ;
if (vis[T] == 1) return 1 ;
for (int i = 0 ; i <= cnt ; ++ i)
if (flow(i) < cap(i) && vis[from(i)] && !vis[to(i)])
chkmin(lack, val[from(i)] - val[to(i)] + cost(i)) ;
if (lack == 1 << 30) return 0 ;
for (int i = 1 ; i <= n ; ++ i)
if (vis[i]) val[i] -= lack ;
return 1 ;
}
void solve(){
spfa() ;
// debug(val, 1, n) ;
mflow = 0 ; mcost = 0 ;
do{
memset(vis, 0, sizeof(vis)) ;
advance(S, 1 << 30) ;
}while (modify()) ;
cout << mflow << " " << mcost << '\n' ;
}
}

int n, m, s, t ;

int main(){
cin >> n >> m >> s >> t ;
zkw :: reset(s, t, n) ;
int x, y, f, c ;
for (int i = 1 ; i <= m ; ++ i){
scanf("%d%d%d%d", &x, &y, &f, &c) ;
zkw :: add(x, y, f, c) ;
}
zkw :: solve() ;
return 0 ;
}

上下界网络流

这一部分是跟 __stdcall 学的。感觉他讲的真的是最详细的了。这个地方安利一下:

以下默认用 $(u,v,\mathrm{lower},\mathrm{upper})$ 来描述一条带有上下界的弧,用 $(u,v,{\rm cap})$ 表示正常流图里的一条弧。

无源汇上下界可行流

建模

大概是考虑先建立一个虚的源和汇 $\rm S, \rm T$。然后考虑对于一条弧 $(u,v,c_1,c_2)$ 将其拆成 $({\rm S}, v,c_1),(u, {\rm T},c_1),(u,v,c_2-c_1)$ 三条正常流图的弧(称前两条弧为附加弧),跑最大流。之后对于每条非附加弧,流量加上这条弧原本该有的下界 $c_2$ 就是每条弧的实际流量。

原理

大概就是所有附加弧保证了一定会流满下界。非附加弧承载的则是「自由流量」,是没有下界可以随便流的弧。不难发现如果所有附加弧都满载,最后流的一定会是一个可行流。

有源汇上下界可行流

建模

考虑在上一个模型的基础上稍加修改。设原图的源和汇分别是 $s,t$ (下同),那么添加一条 $(t,s,+\infty)$ 的弧即可。这样最后跑出来的就是原图的一个可行流,同时弧 $(t,s)$ 的流量就是原图的总流量。

原理

考虑有源汇和无源汇本质上区别就在于,有源汇的图里,$s,t$ 是不满足流量平衡的。所以为了维护流量平衡,就需要建一条 $(t,s,+\infty)$ 的弧来疏通流量。

有源汇上下界最大流

建模

考虑先跑一下可行流,如果有可行流,那么就在跑完的残量网络上再跑一次从原图的源到原图的汇的最大流,即 $s\to t$ 的最大流。那么第二次最大流的结果减就是所求的最大流。

原理

考虑在一次可行流之后,原图的一定是满足流量下界的。同时因为 $t\to s$ 弧的存在,在新的流图上 $t\to s$ 积攒在反向边 $(s,t)$ 中的流量会直接流给 $t$ 。同时只需要让剩下的边流到不能流就好了。

有源汇上下界最小流

建模 1

首先按照无源汇可行流的方式建图。跑一下。之后再把 $(t,s,+\infty)$ 建回去,跑一遍从 $\rm S\to T$ 的最大流,然后去 check 是否合法,如果合法即为答案。

原理 1

考虑本质上是要让可行流最小。根据上文可以知道一组可能流的流量就是 $t\to s$ 的流量,那么为了让可以不流过 $t\to s$ 的流量尽量不要流过这条弧,可以先不考虑这条弧尽量流一波,之后再跑包含这条弧的最大流。这样一方面是可行流,另一方面也会让可行流最小。

建模 2

首先按照有源汇可行流的方式建图。跑一 z下,记此时的可行流为 $f_1$ 。之后从 $\rm T\to S$ 跑一遍最大流,记此时的流为 $f_2$。那么答案就是 $f_1-f_2$ 。

原理 2

考虑正向边流量和反向边流量是对称的,那么考虑原图中的一个可行流中,某些弧的流量可以变小,那么如果使反向边流量尽量大,就可以使正向边流量尽量小。同时可以知道这样是一定合法的,因为附加弧不会被流,会被流的只有非附加弧,同时非附加弧不可能出现流量为负。所以不难知道这样是对的。

「My-Code」
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
namespace ISAP{
#define to(k) e[k].to
#define cap(k) e[k].cap
#define cost(k) e[k].cost
#define from(k) e[k].from
#define flow(k) e[k].flow
struct edge{
int from ;
ll flow ;
ll cap ;
int to ;
edge (int a = 0, int b = 0, ll c = 0){
from = a, to = b, cap = c, flow = 0 ;
}
} e[N] ;
int n ;
ll ans ;
int cnt ;
int S, T ;
int h, t ;
int que[N] ;
int pre[N] ;
int hgt[N] ;
int cur[N] ;
int gap[N] ;
vector <int> E[N] ;
void add(int x, int y, ll c){
e[++ cnt] = edge(x, y, c) ; E[x].push_back(cnt) ;
e[++ cnt] = edge(y, x, 0) ; E[y].push_back(cnt) ;
}
void reset(int x, int y, int z){
memset(cur, 0, sizeof(cur)) ;
memset(hgt, 0, sizeof(hgt)) ;
memset(gap, 0, sizeof(gap)) ;
S = x, T = y, n = z ; cnt = -1 ;
for (int i = 1 ; i <= n ; ++ i) E[i].clear() ;
}
void prebfs(){
que[h = t = 1] = T ;
while (h <= t){
int x = que[h ++] ;
for (auto k : E[x])
if (!hgt[to(k)] && to(k) != T)
hgt[que[++ t] = to(k)] = hgt[x] + 1 ;
}
for (int i = 1 ; i <= n ; ++ i) gap[hgt[i]] ++ ;
}
void Aug(){
int x = T ;
ll z = 1ll << 60 ;
while (x != S){
chkmin(z, cap(pre[x]) - flow(pre[x])) ;
x = from(pre[x]) ;
}
ans += z ; x = T ;
while (x != S){
flow(pre[x]) += z ;
flow(pre[x] ^ 1) -= z ;
x = from(pre[x]) ;
}
}
bool Advance(int x, int &y){
bool ret = 0 ;
for (int d = cur[x] ; d < E[x].size() ; ++ d){
int k = E[x][d] ;
if (flow(k) >= cap(k)) continue ;
if (hgt[x] != hgt[to(k)] + 1) continue ;
cur[x] = d ; ret = 1 ;
pre[y = to(k)] = k ; break ;
}
return ret ;
}
void isap(){
prebfs() ; int x = S ;
while (hgt[S] < n){
// cout << x << '\n' ;
if (x == T)
Aug(), x = S ;
if (!Advance(x, x)){
int minh = n - 1 ;
for (auto k : E[x])
if (flow(k) < cap(k))
chkmin(minh, hgt[to(k)]) ;
if (!-- gap[hgt[x]]) return ;
gap[hgt[x] = minh + 1] ++ ; cur[x] = 0 ;
if (x != S) x = from(pre[x]) ;
}
}
}
}

namespace NO_Lower_Upper_Feasible{
using namespace :: ISAP ;
void pre(int n, int m){
int x, y, c1, c2 ;
reset(n + 1, n + 2, n + 2) ;
for (int i = 1 ; i <= m ; ++ i){
cin >> x >> y >> c1 >> c2 ;
add(S, y, c1) ;
add(x, T, c1) ;
add(x, y, c2 - c1) ;
}
}
void solve(){
isap() ;
for (int i = 1 ; i <= cnt ; ++ i){
if (from(i) == S && flow(i) != cap(i)) goto fuck ;
else if (to(i) == T && flow(i) != cap(i)) goto fuck ;
}
cout << "YES" << "\n" ;
for (int i = 4 ; i <= cnt ; i += 6){
// cout << from(i) << " " << to(i) << '\n' ;
cout << flow(i) + cap(i - 2) << "\n" ;
}
return ;
fuck : cout << "NO" << '\n' ;
}
}

namespace YES_Lower_Upper_Maximum{
using namespace :: ISAP ;
void pre(int n, int m, int s, int t){
int x, y, c1, c2 ;
reset(n + 1, n + 2, n + 2) ;
for (int i = 1 ; i <= m ; ++ i){
cin >> x >> y >> c1 >> c2 ;
add(S, y, c1) ;
add(x, T, c1) ;
add(x, y, c2 - c1) ;
}
add(t, s, 1ll << 60) ;
}
void solve(int s, int t){
isap() ; //debug(ans) ;
for (int i = 1 ; i <= cnt ; ++ i){
if (from(i) == S && flow(i) != cap(i)) goto fuck ;
else if (to(i) == T && flow(i) != cap(i)) goto fuck ;
}
S = s, T = t ;
ans = 0 ; isap() ;
cout << ans << "\n" ;
return ;
fuck : cout << "please go home to sleep" << '\n' ;
}
}
namespace YES_Lower_Upper_Minimum{
using namespace :: ISAP ;
void pre(int n, int m, int s, int t){
int x, y ; ll c1, c2 ;
reset(n + 1, n + 2, n + 2) ;
for (int i = 1 ; i <= m ; ++ i){
x = qr() ;
y = qr() ;
c1 = qr() ;
c2 = qr() ;
add(S, y, c1) ;
add(x, T, c1) ;
add(x, y, c2 - c1) ;
}
}
void solve(int s, int t){
isap() ; //debug(ans) ;
add(t, T, 0) ;
add(S, s, 0) ;
add(t, s, 1ll << 60) ;
ans = 0 ; isap() ;
for (int i = 1 ; i <= cnt ; ++ i){
if (from(i) == S && flow(i) != cap(i)) goto fuck ;
else if (to(i) == T && flow(i) != cap(i)) goto fuck ;
}
cout << ans << "\n" ;
return ;
fuck : cout << "please go home to sleep" << '\n' ;
}
}

上下界费用流

和最大流没什么本质区别。唯一的不同在于,比如求最小费用最大流的时候第二遍只清空 mflow 而不清空 mcost。因为有一部分流是通过 $t\to s$ 这个不定费用弧流过来的,这时是不计算流量的。感觉其他的也大差不差?稍微改一波就可以了。