神秘网络流卡常
UnionRE
·
·
算法·理论
UPD：相比原文章修改了“HLPP 的堆”的错误复杂度。
前言
$n\le 10^6,m\le 3\times 10^6$ 是常见 $O(n^2m)$ 数据范围是吧。
$n\le 10^5$ 是常见 $O(n^2\log n)$ 数据范围是吧。
https://www.luogu.com.cn/training/1087835 这是题单。
## 参考价值
以下内容：
- $\color{green}\blacksquare$ 表示基础。
- $\color{#0afafb} \blacksquare$ 表示推荐。
- $\blue \blacksquare$ 表示在极端环境可以使用。
- $\red \blacksquare$ 表示在超级无敌巨大构式出题人毒瘤条件可以使用。
$$
\gdef \basict{\green{\blacksquare}}
\gdef \good{\color{#0afafb} \blacksquare}
\gdef \extr{\blue \blacksquare}
\gdef \never{\red \blacksquare}
$$
## 网络流！
### $\color{#0afafb} \blacksquare$ ISAP
ISAP 换 dinic。显然能快很多。
因为不仅少 bfs 常数，还有 gap 优化。
### $\color{#0afafb} \blacksquare$ 非递归 ISAP
更快了！！实测快 $30\%$。
参考实现：
```cpp
int isap(){
int flow=0,u=s;
top=0;
memset(buc,0,(t+1)<<2);
bfs();
memcpy(cur,head,(t+1)<<2);
while(lvl[s]<=t){
if(u==t){
int add=INT_MAX;
for(int i=1;i<=top;i++)add=min(add,W[stk[i]]);
for(int i=1;i<=top;i++){
W[stk[i]]-=add;
W[stk[i]^1]+=add;
}
flow+=add;
u=s;
top=0;
continue;
}
bool ok=0;
for(int i=cur[u];i;i=ed[i].nxt){
cur[u]=i;
int v=ed[i].v;
if(lvl[v]==lvl[u]-1&&W[i]&&lvl[v]>=0){
stk[++top]=i;
u=v;
ok=1;
break;
}
}
if(ok)continue;
int minl=t;
for(int i=head[u];i;i=ed[i].nxt){
int v=ed[i].v;
if(W[i]&&lvl[v]>=0)minl=min(minl,lvl[v]);
}
if(!--buc[lvl[u]])break;
lvl[u]=minl+1;
buc[lvl[u]]++;
cur[u]=head[u];
if(u!=s){
int i=stk[top--];
u=ed[i^1].v;
}
}
return flow;
}
```
### $\color{#0afafb} \blacksquare$ 小优化
刚才弹栈时可以记录最小值位置 $I$。然后不直接把 $u$ 设为 $0$ 而是设为 $rev(stk_I)$。$rev$ 表示反边。
可能少一些常数。和 dfs 一样。
### $\color{#0afafb} \blacksquare$ 建模优化
有时出现：
$s\to$ 白点，内部边，黑点 $ \to t$。
考察黑白点数量，如果白点过多，通过交换 $s,t$，尽量让 $s$ 连得少。
### $\color{#0afafb} \blacksquare$ 双向边优化
常见网络流是 $u\xrightarrow{w} v$ 拆成 $u\xrightarrow xv$ 和 $u\xleftarrow 0v$。
但是 $u\xleftrightarrow wv$ 也可以直接两条边完事。
### $\color{#0afafb} \blacksquare$ 建边优化
首先有分段加边。
可是有时候常数还是很大。因此其实可以按边权直接分类。
先建 $\infin$ 边，再建 $w$ 边，最后建 $0$ 边。
### $\color{#0afafb} \blacksquare$ HLPP or ISAP？？
关于速度：稀疏图 ISAP，稠密图 HLPP。
### $\good$ 卡空间？
对于匹配类型，正边边权为 $0$ 的和反边边权为 $1$ 的条数为 $O(n)$。
这样复杂度用 `std::set` 就多个 $\log n$。不过没关系，使用 [spn_set](https://www.luogu.me/paste/pnparijs) 就行了。大概有 $\dfrac1{150}$ 常数。
下面给出一个精简版本：
:::info[code]
```cpp
#include<bits/stdc++.h>
using namespace std;
// by stripe_python
#ifndef SPN_SET_H
#define SPN_SET_H
template <class T, class Cmp, const bool Unique>
class _spn_set_base {
private:
static constexpr size_t _load = 150;
Cmp cmp;
size_t _size;
vector<vector<T>> lists;
vector<T> maxes;
vector<size_t> index;
void build () noexcept ;
pair<size_t, size_t> pos (size_t) noexcept ;
size_t loc (size_t, size_t) noexcept ;
void split (size_t) noexcept ;
void merge (size_t) noexcept ;
void erase (size_t, size_t) ;
public:
_spn_set_base() : cmp(), _size(0) {}
class iterator;
iterator begin ();
iterator end ();
size_t size () const noexcept ;
bool empty () const noexcept ;
void clear () noexcept ;
bool insert (const T&) ;
bool erase (const T&) ;
bool contains (const T&) noexcept ;
};
template <class T, class Cmp, bool Unique>
class _spn_set_base<T, Cmp, Unique>::iterator {
public:
_spn_set_base<T, Cmp, Unique> *set;
size_t pos, idx;
iterator() : set(nullptr), pos(-1), idx(-1) {}
iterator(_spn_set_base<T, Cmp, Unique> *_set, size_t _pos, size_t _idx)
: set(_set), pos(_pos), idx(_idx) {}
operator size_t () {return set->loc(pos, idx);}
T operator* () {return set->lists[pos][idx];}
iterator& operator++ () {
if (++idx >= set->lists[pos].size()) pos++, idx = 0;
return *this;
}
iterator& operator-- () {
if (--idx == (size_t) -1) pos--, idx = set->lists[pos].size() - 1;
return *this;
}
bool operator==(const iterator& other) const {
return set == other.set && pos == other.pos && idx == other.idx;
}
bool operator!=(const iterator& other) const {
return !(*this == other);
}
iterator operator++(int) {
iterator tmp = *this;
++(*this);
return tmp;
}
iterator operator--(int) {
iterator tmp = *this;
--(*this);
return tmp;
}
};
template <class T, class Cmp, bool Unique>
typename _spn_set_base<T, Cmp, Unique>::iterator _spn_set_base<T, Cmp, Unique>::begin() {
return typename _spn_set_base<T, Cmp, Unique>::iterator(this, 0, 0);
}
template <class T, class Cmp, bool Unique>
typename _spn_set_base<T, Cmp, Unique>::iterator _spn_set_base<T, Cmp, Unique>::end() {
return typename _spn_set_base<T, Cmp, Unique>::iterator(this, lists.size(), 0);
}
template <class T, class Cmp, bool Unique>
void _spn_set_base<T, Cmp, Unique>::build() noexcept {
size_t n = lists.size();
index.resize(n, 0);
for (size_t i = 1; i < n; i++) {
index[i] += lists[i - 1].size();
size_t x = i + (i & -i);
if (x < n) index[x] += index[i];
}
}
template <class T, class Cmp, bool Unique>
pair<size_t, size_t> _spn_set_base<T, Cmp, Unique>::pos(size_t idx) noexcept {
if (idx < lists[0].size()) return {0, idx};
if (index.empty()) build();
size_t n = index.size(), p = 0;
for (size_t i = __lg(n); i != (size_t) -1; i--) {
size_t x = p + (1 << i);
if (x < n && idx >= index[x]) idx -= index[x], p = x;
}
return {p, idx};
}
template <class T, class Cmp, bool Unique>
size_t _spn_set_base<T, Cmp, Unique>::loc(size_t pos, size_t idx) noexcept {
if (pos == 0) return idx;
if (index.empty()) build();
while (pos > 0) idx += index[pos], pos &= pos - 1;
return idx;
}
template <class T, class Cmp, bool Unique>
void _spn_set_base<T, Cmp, Unique>::split(size_t pos) noexcept {
if (lists[pos].size() > (_load << 1)) {
vector<T> half(lists[pos].begin() + _load, lists[pos].end());
lists[pos].erase(lists[pos].begin() + _load, lists[pos].end());
maxes[pos] = lists[pos].back();
lists.insert(lists.begin() + pos + 1, half);
maxes.insert(maxes.begin() + pos + 1, half.back());
index.clear();
} else if (!index.empty()) {
size_t n = index.size();
for (size_t i = pos + 1; i < n; i += (i & -i)) index[i]++;
}
}
template <class T, class Cmp, bool Unique>
void _spn_set_base<T, Cmp, Unique>::merge(size_t pos) noexcept {
lists[pos - 1].insert(lists[pos - 1].end(), lists[pos].begin(), lists[pos].end());
maxes[pos - 1] = lists[pos - 1].back();
lists.erase(lists.begin() + pos), maxes.erase(maxes.begin() + pos), index.clear();
}
template <class T, class Cmp, bool Unique>
void _spn_set_base<T, Cmp, Unique>::erase(size_t pos, size_t idx) {
lists[pos].erase(lists[pos].begin() + idx), _size--;
size_t n = lists[pos].size();
if (n > (_load >> 1)) {
maxes[pos] = lists[pos].back();
if (!index.empty()) {
size_t n = index.size();
for (size_t i = pos + 1; i < n; i += (i & -i)) index[i]--;
}
} else if (lists.size() > 1) {
if (pos == 0) pos++;
merge(pos), split(pos - 1);
} else if (n > 0) {
maxes[pos] = lists[pos].back();
} else {
lists.erase(lists.begin() + pos), maxes.erase(maxes.begin() + pos), index.clear();
}
}
template <class T, class Cmp, bool Unique>
size_t _spn_set_base<T, Cmp, Unique>::size() const noexcept {return _size;}
template <class T, class Cmp, bool Unique>
bool _spn_set_base<T, Cmp, Unique>::empty() const noexcept {return maxes.empty();}
template <class T, class Cmp, bool Unique>
void _spn_set_base<T, Cmp, Unique>::clear() noexcept {
_size = 0, lists.clear(), maxes.clear(), index.clear();
}
template <class T, class Cmp, bool Unique>
bool _spn_set_base<T, Cmp, Unique>::insert(const T& val) {
if (maxes.empty()) {
lists.emplace_back(1, val), maxes.emplace_back(val), _size++;
return true;
}
size_t pos = upper_bound(maxes.begin(), maxes.end(), val, cmp) - maxes.begin();
#define _check(val0) do {if (!cmp(val, val0) && !cmp(val0, val)) return false;} while (0)
if (pos == maxes.size()) {
if (Unique) _check(maxes.back());
pos--, lists[pos].emplace_back(val), maxes[pos] = val;
} else {
auto it = upper_bound(lists[pos].begin(), lists[pos].end(), val, cmp);
if (Unique) {
if (it != lists[pos].begin()) _check(*prev(it));
else if (pos > 0) _check(lists[pos - 1].back());
}
lists[pos].insert(it, val);
}
#undef _check
_size++, split(pos);
return true;
}
template <class T, class Cmp, bool Unique>
bool _spn_set_base<T, Cmp, Unique>::erase(const T& val) {
if (maxes.empty()) return false;
size_t pos = lower_bound(maxes.begin(), maxes.end(), val, cmp) - maxes.begin();
if (pos == maxes.size()) return false;
size_t idx = lower_bound(lists[pos].begin(), lists[pos].end(), val, cmp) - lists[pos].begin();
if (cmp(lists[pos][idx], val) || cmp(val, lists[pos][idx])) return false;
erase(pos, idx);
return true;
}
template <class T, class Cmp, bool Unique>
bool _spn_set_base<T, Cmp, Unique>::contains(const T& val) noexcept {
if (maxes.empty()) return false;
size_t pos = lower_bound(maxes.begin(), maxes.end(), val, cmp) - maxes.begin();
if (pos == maxes.size()) return false;
size_t idx = lower_bound(lists[pos].begin(), lists[pos].end(), val, cmp) - lists[pos].begin();
return !cmp(val, lists[pos][idx]) && !cmp(lists[pos][idx], val);
}
template <class T, class Cmp = less<T>>
using spn_set = _spn_set_base<T, Cmp, true>;
template <class T, class Cmp = less<T>>
using spn_multiset = _spn_set_base<T, Cmp, false>;
#endif //SPN_SET_H
```
:::
### $\good$ HLPP 的堆？
狗都不用。
直接用桶维护吧！暴力找最大值，均摊下来 $O(n)$。
:::error[用了就爽变成 $\bm{O(n^2m/w)}$]
```cpp
struct bueue{
int ma=-1,siz=0;
vector<int> t[N];
bitset<N+5> vis;
void push(int x,int y){
if(x>n)return;
siz++;
t[x].push_back(y);
vis[n-x]=1;
ma=max(ma,x);
}
void pop(){
siz--;
t[ma].pop_back();
if(!t[ma].size())vis[n-ma]=0;
ma=n-vis._Find_first();
}
int size(){return siz;}
int top(){return t[ma].back();}
};
```
:::
:::success[用了还是 $\bm{O(n^2\sqrt m)}$]
```cpp
struct bueue{
int ma=-1,siz=0;
vector<int> t[N];
bitset<N> vis;
void push(int x,int y){
if(x>n)return;
siz++;
t[x].push_back(y);
vis[n-x]=1;
ma=max(ma,x);
}
void pop(){
siz--;
t[ma].pop_back();
if(!t[ma].size()){
vis[n-ma]=0;
ma=n-vis._Find_next(n-ma);
}
}
int size(){return siz;}
int top(){return t[ma].back();}
void clear(int n){
siz=0,ma=-1;
vis.reset();
for(int i=0;i<=n;i++)t[i].clear();
}
}pq;
```
:::
### $\good$ `std::vector` 存图与 `reserve(deg[u])`
我们知道 vector 很慢但是其实是因为没有 `reserve(SIZE)`，否则其实和数组一样快。（[参考链接](https://www.luogu.me/article/b34zovu7)）
vector 存图又有吊打链式前向星的内存连续性。
所以我们离线一下图然后 `reserve(deg[u])`。
```cpp
vector<array<int,3>> edges;
edges.reserve(m);
while(m--){
int u,v,w;
u=read(),v=read(),w=read();
deg[u]++;deg[v]++;
edges.push_back({u,v,w});
}
for(int i=1;i<=n;i++)G[i].reserve(deg[i]);
for(auto [u,v,w]:edges){
G[u].push_back({v,w,G[v].size()});
G[v].push_back({u,0,G[u].size()-1});
}
```
### $\good $ 几种不同的 HLPP 实现？
:::success[722ms]
```cpp
int i=cur[u],sz=G[u].size();
while(i<sz){
auto &[v,w,fb]=G[u][i];
if(h[v]==h[u]-1&&w){//推
int use=min(w,ex[u]);
w-=use;
G[v][fb].w+=use;
ex[u]-=use;
ex[v]+=use;
if(ex[v]==use&&v!=s&&v!=t){
cur[v]=0;
pq.push(h[v],v);
}
}
if(!ex[u]){
cur[u]=i;
break;
}
if(w)mi=min(mi,h[v]);
cur[u]=++i;
}
```
:::
:::success[760ms]
```cpp
for(int i=cur[u];i<G[u].size();i++){
auto &[v,w,fb]=G[u][i];
if(h[v]==h[u]-1&&w){//推
int use=min(w,ex[u]);
w-=use;
G[v][fb].w+=use;
ex[u]-=use;
ex[v]+=use;
if(ex[v]==use&&v!=s&&v!=t){
cur[v]=0;
pq.push(h[v],v);
}
}
cur[u]=i;
if(!ex[u])break;
if(w)mi=min(mi,h[v]);
}
```
:::
容易发现快的那个对边的利用率更高。
### $\basict$ 一些减少 HLPP 队列大小常数的？
你应该：
```cpp
if(h[u]<=n)pq.push(h[u],u);
```
而不是直接 push。
### $\good$ 调节重标号次数的参数
可能会有些神秘作用？反正随机调参。
### $\good$ 完全图叠加边权！
如下操作：
```cpp
while(m--){
int u,v,w;
cin>>u>>v>>w;
val[u][v]+=w;
deg[u]++;deg[v]++;
}
for(int i=1;i<=n;i++)G[i].reserve(deg[i]*2);
for(int u=1;u<=n;u++){
for(int v=1;v<=n;v++){
if(val[u][v])adde(u,v,val[u][v]);
}
}
```
另外，似乎**独立地存边权可以更优秀？**
### $\extr$ 网格图优化
不直接建边，而是直接存四个方向的权。
**这样在枚举四周时可以循环展开！**
### $\basict$ 还在用 `memset(A,B,sizeof(A))`？！
尝试 `memset(A,B,n<<3)` 吧。
如果要用 `long long` 的话，就吧 `sizeof(int)` 换成 `sizeof(long long)` 就行了。
### $\basict$ 还在用 `for(int &i=cur[u];i<G[u].size();i++)`？！
不知道啥原理，反正不写 `&i` 而是手动 `cur[u]=i` 会快一些。
## 费用流！
### $\color{#0afafb} \blacksquare$ 原始对偶
虽然是复杂度优化，$O(nmf)\to O(mf\log n)$。
但是确实我一般用来卡常就是了。
### $\extr$ 更快的 dijkstra
用一个很快的堆！比如手写斐波那契堆或者配对堆。可以做到 $O(m+n\log n)$ 的 dijkstra。
总复杂度 $O((m+\log n)f)$。
### $\extr$ 诗人啊？循环展开建图
参考：https://www.luogu.com.cn/record/286358799 我卡了一小会（~~指一个上午~~）。
推荐 wxx 大神的[题解](https://www.luogu.com.cn/article/hzr7vum4)。
### $\extr$ 完全图上正常人就应该用链表优化的 dij！
常数小一点！
```cpp
template <typename T>
struct mylist{
struct data{T val;short nxt,pre;};
data a[N];//next pre
int head=0;
void init(short L,short R){
head=L;
a[L]={0,L+1,-1};
for(int i=L+1;i<=R-1;i++)a[i]={inf,i+1,i-1};
a[R]={inf,-1,R-1};
}
void del(int pos){
if(pos==head)head=a[pos].nxt;
a[a[pos].pre].nxt=a[pos].nxt;
a[a[pos].nxt].pre=a[pos].pre;
}
pair<T,short> qmin(T MAX){
T minn=a[head].val;
short id=head;
for(register int i=a[head].nxt;i!=-1&&minn;i=a[i].nxt){//L 是 s，不是
if(a[i].val<minn){
minn=a[i].val;
id=i;
}
}
return {minn,id};
}
};
mylist<long long> q;
queue<short> p;
inline bool dijkstra(){
for(register int i=s;i<=t;++i)dis[i]=inf;
dis[s]=0;
q.init(s,t);
for(register int i=s;i<=t;i++){
auto [d,u]=q.qmin(inf);
if(u==-1)break;
q.del(u);
if(d==inf)break;
for(int v:G0[u]){
bool flow=G[u][v].w;
long long c=G[u][v].c;
c+=h[u]-h[v];
if(dis[v]>dis[u]+c&&flow){
dis[v]=dis[u]+c;
q.a[v].val=dis[v];
pre[v]=u;
}
}
}
vis.reset();
return dis[t]!=inf;
}
```
### $\never$ 我不明白
https://www.luogu.com.cn/discuss/1295818 为什么 EK 实现的 MCMF 能过啊。
***
:::info[对网络流进行 golf？]{open}
P3376 【模板】网络最大流 AC 代码：
（这还是一个兼顾长度和性能的 ISAP！用时：66ms）
```cpp
#include<bits/stdc++.h>
#define R return
#define E(u,v,w)Y[C]={v,w,H[u]},H[u]=C++
using namespace std;using Z=long;const Z N=410,M=2e5;Z n,m,s,t,L[N],b[N],H[N],c[N],C=2,u,v,w,Q[N],l,r,i,U,F;array<Z,3>Y[M];void B(){b[L[Q[++r]=t]=1]=1;while(l<r){u=Q[++l];for(i=H[u];i;i=Y[i][2])v=Y[i][0],!L[v]&&Y[i^1][1]&&(L[v]=L[u]+1,b[L[v]]++,Q[++r]=v);}}Z D(Z u,Z A){if(u^t){Z F=0;for(Z i=c[u];i&&A;i=Y[i][2]){c[u]=i;auto&[v,w,_]=Y[i];L[v]==L[u]-1&&w&&(U=D(v,min(A,w)),U&&(w-=U,Y[i^1][1]+=U,A-=U,F+=U));}R A?b[L[u]]--,!b[L[u]]&&(L[s]=n+1),L[u]++,b[L[u]]++,F:F;}R A;}main(){for(cin>>n>>m>>s>>t;i<m;i++)cin>>u>>v>>w,E(u,v,w),E(v,u,0);B();while(L[s]<=n)memcpy(c,H,3e3),F+=D(s,1e18);cout<<F;}
```
:::