一种神奇的作差法
暑假时的一些naive想法
黑历史吧
前言：
疫情原因，来不了机房集训，只来得及听一节课的ZJY学长分享网络流……
然后发现自己甚至连费用流板子都没有写过……
学习费用流时我一听说要写 SPFA ，心想：这不是死算法了吗
（然而根据 OIER 业界规则，最大流卡 Dinic，费用流卡 SSP+SPFA 的都是不人道）
在网上翻了一翻，发现真有一种奇怪的方法 Primal Dual（然而是经典算法）能避免 SPFA ，在各种各样形态的图上都较为稳定，于是我就一路递归点起了一颗奇怪的科技树
在写费用流时去复习了一下最短路（没错我最短路都能忘了写），又发现了另一种奇怪的算法 Johnson （然而也是经典算法）
发现了这两种算法运用了同一种思想，于是惊奇地写下这篇博文
首先介绍一下并没有什么关系的势能分析（单纯是因为思想比较像）：
势能分析通常选取一个“势函数”，记为 \(\phi(i)\) 表示一个数据结构在第 i 次操作后的状态，假设对于第 i 次操作，理论时间复杂度为 \(a_i\) ，定义其均摊时间复杂度 \(b_i=a_i+\phi(i)-\phi(i-1)\) ，则总均摊时间复杂度为 \(\sum b_i=\sum a_i+\phi(n)-\phi(0)\) ，如果存在 \(\phi(n)\geq\phi(0)\) ，那么分析出来的时间复杂度就是实际时间复杂度的上界
实际应用中 \(\phi(0)\) 通常是零，而其他势能保持非负。势能分析本质的思想就是形式化地考虑每一次操作对于数据结构的影响
数学不太好，Splay 一堆 \(\log\) 符号分析不动，就来个简单的吧
思路：
如果只是单纯的栈的话直接对于每一个取模后的值 DP 即可，但是现在是双端队列
双端队列 C++ STL::deque 的实现方式其实类似双栈（实际是指针维护的多段数组，思想上也比较类似分块），我们完全可以用双栈的方式对 DP 栈进行维护
直接将队列分成两个栈进行维护，前后直接增删即可
当有一个栈空时，将另一个栈复制过去一半即可
Code:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=50003,P=503;
ll Fd[N][P],Gd[N][P];
int Fu[N],Fv[N],Fs;
int Gu[N],Gv[N],Gs;
int m,p,ids;
void addF(int u,int v){
Fs++;
Fu[Fs]=u; Fv[Fs]=v;
for (int i=0; i<p; i++) Fd[Fs][i]=Fd[Fs-1][i];
for (int i=0; i<p; i++) if (~Fd[Fs-1][i])
Fd[Fs][(i+u)%p]=max(Fd[Fs][(i+u)%p],Fd[Fs-1][i]+v);
}
void addG(int u,int v){
Gs++;
Gu[Gs]=u; Gv[Gs]=v;
for (int i=0; i<p; i++) Gd[Gs][i]=Gd[Gs-1][i];
for (int i=0; i<p; i++) if (~Gd[Gs-1][i])
Gd[Gs][(i+u)%p]=max(Gd[Gs][(i+u)%p],Gd[Gs-1][i]+v);
}
void Move_F_to_G(){
int n=Fs;
Gs=0; Fs=0;
for (int i=(n+1>>1); i; i--) addG(Fu[i],Fv[i]);
for (int i=(n+1>>1)+1; i<=n; i++) addF(Fu[i],Fv[i]);
}
void Move_G_to_F(){
int n=Gs;
Gs=0; Fs=0;
for (int i=(n+1>>1); i; i--) addF(Gu[i],Gv[i]);
for (int i=(n+1>>1)+1; i<=n; i++) addG(Gu[i],Gv[i]);
}
ll st[P][12];
bool flag=1;
int lg[P],bi[12];
void getST(){
for (int i=0; i<p; i++) st[i][0]=Gd[Gs][i];
for (int i=1; i<=lg[p]; i++)
for (int j=0; j+bi[i]<=p; j++)
st[j][i]=max(st[j][i-1],st[j+bi[i-1]][i-1]);
}
ll qu(int l,int r){
int k=lg[r-l+1];
return max(st[l][k],st[r-bi[k]+1][k]);
}
int main(){
scanf("%d",&ids);
scanf("%d%d",&m,&p);
lg[0]=-1; bi[0]=1;
for (int i=1; i<=p; i++) lg[i]=lg[i>>1]+1;
for (int i=1; i<=lg[p]; i++) bi[i]=bi[i-1]*2;
memset(Gd[0],-1,sizeof Gd[0]);
memset(Fd[0],-1,sizeof Fd[0]);
for (int i=0; i<=m; i++) Gd[i][0]=Fd[i][0]=0;
for (int i=1; i<=m; i++){
char s[2];
scanf("%s",s);
if (s[0]=='Q'){
int l,r;
scanf("%d%d",&l,&r);
ll ans=-1;
if (flag) getST();
flag=0;
for (int i=0; i<p; i++) if (~Fd[Fs][i]){
int sl=(l-i+p)%p,sr=(r-i+p)%p; ll t;
if (sl<=sr) t=qu(sl,sr);
else t=max(qu(0,sr),qu(sl,p-1));
if (~t) ans=max(ans,t+Fd[Fs][i]);
}
printf("%lld\n",ans);
}
if (s[0]=='D'){
if (s[1]=='F'){
if (!Fs) Move_G_to_F();
Fs--;
}
if (s[1]=='G'){
if (!Gs) Move_F_to_G();
Gs--;
}
flag=1;
}
if (s[0]=='I'){
int u,v;
scanf("%d%d",&u,&v);
if (s[1]=='F') addF(u,v);
if (s[1]=='G') addG(u,v);
flag=1;
}
}
return 0;
}
这样做复杂度为什么是对的呢？
设势能函数 \(\phi(i)\) 表示两个栈大小差的绝对值
那么当加入一个数时，\(\Delta\phi=\pm1\) ，实际操作复杂度可以看做 \(1\) （不过在本题中并不是，这里仅仅分析双栈），那么均摊复杂度为 0 或 2 。
复制数组时均摊复杂度为 \(\phi(i)-\phi(i)=0\)
综上总时间复杂度在 \(O(n)\) 水平
感觉这道题直接整体分析即可，不需要势能。不过对于 Splay 并查集 等东西的势能分析比较有用
观察势能分析的过程，用到了一个很水的式子
将这个式子用于三角形不等式上，B.Johnson 巨巨在 1977 年提出了一种奇妙的求全源最短路的方法
提到全源最短路，很多人都想到了 Floyd ，但是我当时在想，如果 m 不大的话，直接跑 n 遍单源最短路不是更优吗？
可是我们的 Dijkstra 不支持负权， SPFA 坟头草都长了两年了（话说单源最短路也是这种窘境诶），因此我们有时被迫使用 SPFA （譬如说费用流那种经常有负权的处境）
但是在全源最短路中，SPFA 死透了！（然而还有一口气，为什么待会再说）
如果让 Dijkstra 处理的边权不是负数，怎么做？考虑整体平移，加上某一个定值 \(\delta\)，但是对于两条 \(s \to t\) 的路径，每一条边都加上了 \(\delta\)，如果经过的边数不同，那么最终偏移的量也会不同，造成误判
那么怎样才可以让这个偏移量与经过的边数无关呢
看一看上面那个式子,只要 \(n\) 固定，那么所有的中间偏移量都可以相互抵消，不做考虑
思路来了，令新边权 \(w_e'=w_e+h_u-h_v\) ，路径 \(path (s \to t)\) 总边权 :
在新边权的图上跑 Dijkstra ，再将最后的结果 \(-h_s+h_t\) 即可
问题是如何设置 h 数组使得其能让每一条新边权非负
观察不等式 \(w_e+h_u-h_v\geq 0\) ,变形得 \(h_v \leq h_u+w_e\) ，三角形不等式！！设一个虚拟节点向所有点连边，跑一遍最短路即可
如何处理带负权的最短路？……
来人，快把 SPFA 从坟里刨出来!
由于一定要跑一遍负权图，所以这个方法对于单源最短路是没有优化的
但是全源最短路要跑 n 次，所以可以不计该复杂度
事实上国外是没有 SPFA 这个叫法的，原版 Johnson 是直接跑 Bellman Ford 的（因为这俩的复杂度是一样的）
实在有 PTSD 的可以用 SLF 压压惊
上述过程总复杂度 \(O(nm\log n)\)
Code:
#include <bits/stdc++.h>
using namespace std;
const int _=3003,__=10003;
typedef long long ll;
const ll INF=1000000000;
int n,m;
struct Graph{
int hd[_],nxt[__],ver[__],val[__],tot;
void add(int u,int v,int w){nxt[++tot]=hd[u];hd[u]=tot;ver[tot]=v;val[tot]=w;}
}G;
namespace Johnson_Algorithm{
ll h[_],d[_];
int t[_];
bool v[_];
bool SPFA(){
for (int i=1; i<=n; i++) G.add(0,i,0);
memset(v,0,sizeof v);
memset(t,0,sizeof t);
queue<int> q;
for (int i=1; i<=n; i++) h[i]=INF;
h[0]=0; q.push(0); v[0]=1;
while (!q.empty()){
int x=q.front(); q.pop();
v[x]=0;
for (int i=G.hd[x]; i; i=G.nxt[i]){
int y=G.ver[i],w=G.val[i];
if (h[y]>h[x]+w) {
h[y]=h[x]+w;
if (!v[y]){
v[y]=1;
q.push(y);
t[y]++;
if (t[y]==n+1) return 0;
}
}
}
}
return 1;
}
struct node{
ll dis;
int id;
bool operator <(const node &x)const{return dis>x.dis;}
node(ll di,int i){dis=di;id=i;}
};
ll dijkstra(int s){
priority_queue<node> q;
memset(v,0,sizeof v);
for (int i=1; i<=n; i++) d[i]=INF;
d[s]=0; q.push(node(0,s));
while (!q.empty()){
int x=q.top().id; q.pop();
if (v[x]) continue;
v[x]=1;
for (int i=G.hd[x]; i; i=G.nxt[i]){
int y=G.ver[i],w=G.val[i]+h[x]-h[y];
if (d[y]>d[x]+w){
d[y]=d[x]+w;
if (!v[y]) q.push(node(d[y],y));
}
}
}
ll ans=0;
for (int i=1; i<=n; i++){
if (d[i]==INF) ans+=(ll)i*INF;
else ans+=(ll)i*(d[i]+h[i]-h[s]);
}
return ans;
}
};
using Johnson_Algorithm::SPFA;
using Johnson_Algorithm::dijkstra;
int main(){
G.tot=1;
scanf("%d%d",&n,&m);
for (int i=1; i<=m; i++){
int u,v,w;
scanf("%d%d%d",&u,&v,&w);
G.add(u,v,w);
}
if (!SPFA()) {puts("-1");return 0;}
for (int i=1; i<=n; i++)
printf("%lld\n",dijkstra(i));
return 0;
}
顶标法之所以有优化是因为这里在同一张图上需要多次求最短路可以预处理顶标
还有什么地方需要多次求最短路呢？
费用流！！！
费用流的 SSP 算法每次找出一条费用最小的增广路贪心增广，是伪多项式级别算法
而 Primal Dual 支持多路增广，效率好看很多，相比 zkw （我不会，听博客说的）好在效率稳定，受网络的形态影响不大
费用流原图经常没有负边权，但是人家一增广完就会有一堆负边权出来（因为有反向边）
同样考虑设置顶标 \(h_i\) ，在初始无负权的情况下都置零表示不修改，如果有负权就复活 SPFA (SPFA 反复诈尸
以 \(w_e+h_u-h_v\) 边权跑完最短路后，求出 \(dis\) 数组，将 \(h_i'=h_i+dis_i\) 即可（如果与 i 不连通记为 \(\inf\)）
为什么正确呢？首先对于原来的边有三角形不等式:
在网络流中，由于我们是沿着一颗最短路树（SPT）增广（当然我代码是写的单路增广，因此是链），这些树边满足：
由于反向边 \(w_{\bar e}=-w_e\) ，在增广这些树边时会产生新的反向边，这些边满足：
是非负边，综上在 SSP 的过程中维护顶标，可以用 Dijkstra 代替 SPFA
注意上述过程不是正宗的 Primal Dual，因为还可以优化
应该从汇点出发跑最短路，然后从源点出发，只走 \(dis_v=dis_u+w_e\) 的边，弄出一颗 SPT 多路增广（写过 SPT 或者 Dinic 的应该知道），思想上类似 Dinic 的分层，只不过把 bfs 改成 Dijkstra 了
代码写的是单纯的 SSP+Dijkstra ，额外写了 SPFA ，但这道题没必要，不过想留作模板日后用
卡普通 SSP+SPFA 的题目很少，不知道有没有必须要写 PD 或 zkw 的题
Code:
#include <bits/stdc++.h>
using namespace std;
const int _=5005,__=100005;
int n,m,s,t;
struct CostNetFlow{
int hd[_],nxt[__],ver[__],val[__],flow[__],tot;
void add(int u,int v,int w,int l){
nxt[++tot]=hd[u]; hd[u]=tot; ver[tot]=v; flow[tot]=w; val[tot]=l;
nxt[++tot]=hd[v]; hd[v]=tot; ver[tot]=u; flow[tot]=0; val[tot]=-l;
}
}G;
typedef long long ll;
const int INF=0x3f3f3f3f;
namespace MCMF{
struct node{
int dis,id;
bool operator <(const node &x)const{return dis>x.dis;}
node(int di,int i){dis=di; id=i;}
};
int h[_],d[_];
bool v[_];
int po[_],pe[_];
void SPFA(){
memset(v,0,sizeof v);
memset(h,0x3f,sizeof h);
queue<int> q;
h[s]=0; v[s]=1;
q.push(s);
while (!q.empty()){
int x=q.front(); q.pop();
v[x]=0;
for (int i=G.hd[x]; i; i=G.nxt[i]){
if (!G.flow[i]) continue;
int y=G.ver[i],w=G.val[i];
if (h[y]>h[x]+w){
h[y]=h[x]+w;
if (!v[y]){
v[y]=1;
q.push(y);
}
}
}
}
}
int mxflow,mncost;
bool dijkstra(){
memset(v,0,sizeof v);
memset(d,0x3f,sizeof d);
priority_queue<node> q;
q.push(node(0,s)); d[s]=0;
while (!q.empty()){
int x=q.top().id; q.pop();
if (v[x]) continue;
v[x]=1;
for (int i=G.hd[x]; i; i=G.nxt[i]){
if (!G.flow[i]) continue;
int y=G.ver[i],w=G.val[i]+h[x]-h[y];
if (d[y]>d[x]+w){
d[y]=d[x]+w;
if (!v[y]){
po[y]=x; pe[y]=i;
q.push(node(d[y],y));
}
}
}
}
if (d[t]>=INF) return 0;
for (int i=1; i<=n; i++) h[i]=min(h[i]+d[i],INF);
int nfl=INF;
for (int i=t; i^s; i=po[i])
nfl=min(nfl,G.flow[pe[i]]);
mxflow+=nfl;
mncost+=nfl*h[t];
for (int i=t; i^s; i=po[i])
G.flow[pe[i]]-=nfl,G.flow[pe[i]^1]+=nfl;
return 1;
}
void flow(){
mxflow=0; mncost=0;
memset(h,0,sizeof h);
while (dijkstra());
}
};
int main(){
scanf("%d%d%d%d",&n,&m,&s,&t);
G.tot=1;
for (int i=1; i<=m; i++){
int u,v,w,l;
scanf("%d%d%d%d",&u,&v,&w,&l);
G.add(u,v,w,l);
}
MCMF::flow();
printf("%lld %lld\n",MCMF::mxflow,MCMF::mncost);
return 0;
}
为什么这么多英文？因为我懒得打中文
学术内容参考 OI-wiki，洛谷题解：
