Loading [MathJax]/extensions/TeX/newcommand.js
\newcommand{\ord}[1]{\mathcal{O}\left(#1\right)} \newcommand{\abs}[1]{\lvert #1 \rvert} \newcommand{\floor}[1]{\lfloor #1 \rfloor} \newcommand{\ceil}[1]{\lceil #1 \rceil} \newcommand{\opord}{\operatorname{\mathcal{O}}} \newcommand{\argmax}{\operatorname{arg\,max}} \newcommand{\str}[1]{\texttt{"#1"}}

2015年10月14日 星期三

[ Augmenting Path Algorithm ] 二分圖匹配增廣路徑演算法

就是一直找增廣路徑,找到不能再找就是最大匹配惹
複雜度\ord{n1*(m+n2+n1)},用dinic會比較快但是code比較長
模板:
#define MAXN1 505
#define MAXN2 505
int n1,n2;/*n1個點連向n2個點*/
int match[MAXN2];/*每個屬於n2的點匹配了哪個點*/
vector<int > g[MAXN1];/*圖*/
bool vis[MAXN2];/*是否走訪過*/
bool dfs(int u){
for(size_t i=0;i<g[u].size();++i){
int v=g[u][i];
if(vis[v])continue;
vis[v]=1;
if(match[v]==-1||dfs(match[v])){
match[v]=u;
return 1;
}
}
return 0;
}
inline int max_match(){
int ans=0;
memset(match,-1,sizeof(int)*n2);
for(int i=0;i<n1;++i){
memset(vis,0,sizeof(bool)*n2);
if(dfs(i))++ans;
}
return ans;
}

[ Improved Shortest Augmenting Paths,ISAP ,Flow ] 網路流ISAP算法

使用gap優化(間隙優化)
使用當前弧優化

這裡提供遞迴與非遞迴兩種實現
基本上兩個版本效能差不多
遞迴版本:
#include<string.h>
#include<limits.h>
#include<vector>
#include<algorithm>
template<typename T>
struct ISAP{
static const int MAXN=105;
static const T INF=INT_MAX;
int n;//點數
int d[MAXN],gap[MAXN],cur[MAXN];
//層次、gap[i]=層次為i的點之個數、當前弧優化
struct edge{
int v,pre;
T cap,flow,r;
edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
};
int g[MAXN];
std::vector<edge> e;
void init(int _n){
memset(g,-1,sizeof(int)*((n=_n)+1));
e.clear();
}
void add_edge(int u,int v,T cap,bool directed=false){
e.push_back(edge(v,g[u],cap));
g[u]=e.size()-1;
e.push_back(edge(u,g[v],directed?0:cap));
g[v]=e.size()-1;
}
T dfs(int u,int s,int t,T cur_flow=INF){
if(u==t)return cur_flow;
T tf=cur_flow,df;
for(int &i=cur[u];~i;i=e[i].pre){
if(e[i].r&&d[u]==d[e[i].v]+1){
df=dfs(e[i].v,s,t,std::min(tf,e[i].r));
e[i].flow+=df;
e[i^1].flow-=df;
e[i].r-=df;
e[i^1].r+=df;
if(!(tf-=df)||d[s]==n)return cur_flow-tf;
}
}
int minh=n;
for(int i=cur[u]=g[u];~i;i=e[i].pre){
if(e[i].r&&d[e[i].v]<minh)minh=d[e[i].v];
}
if(!--gap[d[u]])d[s]=n;
else ++gap[d[u]=++minh];
return cur_flow-tf;
}
T isap(int s,int t,bool clean=true){
memset(d,0,sizeof(int)*(n+1));
memset(gap,0,sizeof(int)*(n+1));
memcpy(cur,g,sizeof(int)*(n+1));
if(clean){
for(size_t i=0;i<e.size();++i){
e[i].flow=0;
e[i].r=e[i].cap;
}
}
T max_flow=0;
for(gap[0]=n;d[s]<n;)max_flow+=dfs(s,s,t);
return max_flow;
}
std::vector<int> cut_e;//最小割邊集
bool vis[MAXN];
void dfs_cut(int u){
vis[u]=1;//表示u屬於source的最小割集
for(int i=g[u];~i;i=e[i].pre){
if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
}
}
T min_cut(int s,int t){
T ans=isap(s,t);
memset(vis,0,sizeof(bool)*(n+1));
dfs_cut(s),cut_e.clear();
for(int u=0;u<=n;++u){
if(vis[u])for(int i=g[u];~i;i=e[i].pre){
if(!vis[e[i].v])cut_e.push_back(i);
}
}
return ans;
}
};
view raw isap_dfs.cpp hosted with ❤ by GitHub

非遞迴版本:
#include<string.h>
#include<limits.h>
#include<vector>
#include<algorithm>
template<typename T>
struct ISAP{
static const int MAXN=105;
static const T INF=INT_MAX;
int n;//點數
int d[MAXN],gap[MAXN],pre[MAXN],cur[MAXN];
//層次、gap[i]=層次為i的點之個數、前一個節點、當前弧優化
struct edge{
int v,pre;
T cap,flow,r;
edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
};
int g[MAXN];
std::vector<edge> e;
void init(int _n){
memset(g,-1,sizeof(int)*((n=_n)+1));
e.clear();
}
void add_edge(int u,int v,T cap,bool directed=false){
e.push_back(edge(v,g[u],cap));
g[u]=e.size()-1;
e.push_back(edge(u,g[v],directed?0:cap));
g[v]=e.size()-1;
}
T isap(int s,int t,bool clean=true){
T cur_flow,max_flow=0;
int u=s;
memset(d,0,sizeof(int)*(n+1));
memset(gap,0,sizeof(int)*(n+1));
memcpy(cur,g,sizeof(int)*(n+1));
if(clean){
for(size_t i=0;i<e.size();++i){
e[i].flow=0;
e[i].r=e[i].cap;
}
}
for(gap[0]=n;d[s]<n;){
if(u==t){
cur_flow=INF;
for(int i=s;i!=t;i=e[cur[i]].v){
if(cur_flow>e[cur[i]].r)cur_flow=e[cur[i]].r;
}
for(int i=s;i!=t;i=e[cur[i]].v){
e[cur[i]].flow+=cur_flow;
e[cur[i]^1].flow-=cur_flow;
e[cur[i]].r-=cur_flow;
e[cur[i]^1].r+=cur_flow;
}
max_flow+=cur_flow;
u=s;
}
int i;
for(i=cur[u];~i;i=e[i].pre){
if(e[i].r&&d[u]==d[e[i].v]+1)break;
}
if(~i){
cur[u]=i;
pre[e[i].v]=u;
u=e[i].v;
}else{
if(!--gap[d[u]])break;
for(d[u]=n,i=cur[u]=g[u];~i;i=e[i].pre){
if(e[i].r&&d[e[i].v]<d[u])d[u]=d[e[i].v];
}
++gap[++d[u]];
if(u!=s)u=pre[u];
}
}
return max_flow;
}
std::vector<int> cut_e;//最小割邊集
bool vis[MAXN];
void dfs_cut(int u){
vis[u]=1;//表示u屬於source的最小割集
for(int i=g[u];~i;i=e[i].pre){
if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
}
}
T min_cut(int s,int t){
T ans=isap(s,t);
memset(vis,0,sizeof(bool)*(n+1));
dfs_cut(s),cut_e.clear();
for(int u=0;u<=n;++u){
if(vis[u])for(int i=g[u];~i;i=e[i].pre){
if(!vis[e[i].v])cut_e.push_back(i);
}
}
return ans;
}
};
view raw isap.cpp hosted with ❤ by GitHub

[ Dinic's algorithm ,Flow ] 網路流Dinic模板

使用當前弧優化

多路增廣模板:
#include<string.h>
#include<limits.h>
#include<vector>
#include<queue>
#include<algorithm>
template<typename T>
struct DINIC{
static const int MAXN=105;
static const T INF=INT_MAX;
int n;/*點數*/
int level[MAXN],cur[MAXN];/*層次、當前弧優化*/
struct edge{
int v,pre;
T cap,flow,r;
edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
};
int g[MAXN];
std::vector<edge> e;
void init(int _n){
memset(g,-1,sizeof(int)*((n=_n)+1));
e.clear();
}
void add_edge(int u,int v,T cap,bool directed=false){
e.push_back(edge(v,g[u],cap));
g[u]=e.size()-1;
e.push_back(edge(u,g[v],directed?0:cap));
g[v]=e.size()-1;
}
int bfs(int s,int t){
memset(level,0,sizeof(int)*(n+1));
memcpy(cur,g,sizeof(int)*(n+1));
std::queue<int >q;
q.push(s);
level[s]=1;
while(q.size()){
int u=q.front();q.pop();
for(int i=g[u];~i;i=e[i].pre){
if(!level[e[i].v]&&e[i].r){
level[e[i].v]=level[u]+1;
q.push(e[i].v);
if(e[i].v==t)return 1;
}
}
}
return 0;
}
T dfs(int u,int t,T cur_flow=INF){
if(u==t||!cur_flow)return cur_flow;
T df,tf=0;
for(int &i=cur[u];~i;i=e[i].pre){
if(level[e[i].v]==level[u]+1&&e[i].r){
if(df=dfs(e[i].v,t,std::min(cur_flow,e[i].r))){
e[i].flow+=df;
e[i^1].flow-=df;
e[i].r-=df;
e[i^1].r+=df;
tf+=df;
if(!(cur_flow-=df))break;
}
}
}
if(!tf)level[u]=0;
return tf;
}
T dinic(int s,int t,bool clean=true){
if(clean){
for(size_t i=0;i<e.size();++i){
e[i].flow=0;
e[i].r=e[i].cap;
}
}
T ans=0;
while(bfs(s,t))ans+=dfs(s,t);
return ans;
}
std::vector<int> cut_e;//最小割邊集
bool vis[MAXN];
void dfs_cut(int u){
vis[u]=1;//表示u屬於source的最小割集
for(int i=g[u];~i;i=e[i].pre){
if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
}
}
T min_cut(int s,int t){
T ans=dinic(s,t);
memset(vis,0,sizeof(bool)*(n+1));
dfs_cut(s),cut_e.clear();
for(int u=0;u<=n;++u){
if(vis[u])for(int i=g[u];~i;i=e[i].pre){
if(!vis[e[i].v])cut_e.push_back(i);
}
}
return ans;
}
};
view raw dinic.cpp hosted with ❤ by GitHub

單路增廣模板:
#include<string.h>
#include<limits.h>
#include<vector>
#include<queue>
#include<algorithm>
template<typename T>
struct DINIC{
static const int MAXN=105;
static const T INF=INT_MAX;
int n;/*點數*/
int level[MAXN],cur[MAXN];/*層次、當前弧優化*/
struct edge{
int v,pre;
T cap,flow,r;
edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),flow(0),r(cap){}
};
int g[MAXN];
std::vector<edge> e;
void init(int _n){
memset(g,-1,sizeof(int)*((n=_n)+1));
e.clear();
}
void add_edge(int u,int v,T cap,bool directed=false){
e.push_back(edge(v,g[u],cap));
g[u]=e.size()-1;
e.push_back(edge(u,g[v],directed?0:cap));
g[v]=e.size()-1;
}
int bfs(int s,int t){
memset(level,0,sizeof(int)*(n+1));
memcpy(cur,g,sizeof(int)*(n+1));
std::queue<int >q;
q.push(s);
level[s]=1;
while(q.size()){
int u=q.front();q.pop();
for(int i=g[u];~i;i=e[i].pre){
if(!level[e[i].v]&&e[i].r){
level[e[i].v]=level[u]+1;
q.push(e[i].v);
if(e[i].v==t)return 1;
}
}
}
return 0;
}
T dfs(int u,int t,T cur_flow=INF){
if(u==t)return cur_flow;
T df;
for(int &i=cur[u];~i;i=e[i].pre){
if(level[e[i].v]==level[u]+1&&e[i].r){
if(df=dfs(e[i].v,t,std::min(cur_flow,e[i].r))){
e[i].flow+=df;
e[i^1].flow-=df;
e[i].r-=df;
e[i^1].r+=df;
return df;
}
}
}
return level[u]=0;
}
T dinic(int s,int t,bool clean=true){
if(clean){
for(size_t i=0;i<e.size();++i){
e[i].flow=0;
e[i].r=e[i].cap;
}
}
T ans=0,mf=0;
while(bfs(s,t))while(mf=dfs(s,t))ans+=mf;
return ans;
}
std::vector<int> cut_e;//最小割邊集
bool vis[MAXN];
void dfs_cut(int u){
vis[u]=1;//表示u屬於source的最小割集
for(int i=g[u];~i;i=e[i].pre){
if(e[i].flow<e[i].cap&&!vis[e[i].v])dfs_cut(e[i].v);
}
}
T min_cut(int s,int t){
T ans=dinic(s,t);
memset(vis,0,sizeof(bool)*(n+1));
dfs_cut(s),cut_e.clear();
for(int u=0;u<=n;++u){
if(vis[u])for(int i=g[u];~i;i=e[i].pre){
if(!vis[e[i].v])cut_e.push_back(i);
}
}
return ans;
}
};
view raw dinic2.cpp hosted with ❤ by GitHub

[ Edmonds–Karp algorithm ,Flow] 網路流Edmonds–Karp算法

這是最基本的多項式時間網路流演算法,因此在這裡我多了一些變化,首先是一般用鄰接矩陣,複雜度\ord{\abs{V}^5}的做法:
#include<string.h>
#include<limits.h>
#include<queue>
#include<algorithm>
#define MAXN 100
int n,m;/*點數、邊數*/
int g[MAXN+5][MAXN+5],f[MAXN+5][MAXN+5],r[MAXN+5][MAXN+5];
/*容量上限、流量、剩餘容量*/
int v[MAXN+5];/*經過的點*/
int pa[MAXN+5];/*紀錄路徑*/
int d[MAXN+5];/*紀錄流量瓶頸*/
inline int bfs(int s,int t){
memset(v,0,sizeof(v));
std::queue<int >q;
q.push(s);
v[s]=1;
pa[s]=-1;
d[s]=INT_MAX;
while(q.size()){
int x=q.front();q.pop();
for(int i=1;i<=n;++i){
if(v[i]||!r[x][i])continue;
v[i]=1;
pa[i]=x;
d[i]=std::min(d[x],r[x][i]);
q.push(i);
if(i==t)return d[t];
}
}
return 0;
}
inline int edmonds_karp(int s,int t){
memset(f,0,sizeof(f));
memcpy(r,g,sizeof(g));
int ans=0,tmd;
for(;(tmd=bfs(s,t));ans+=tmd){
for(int i=t;~pa[i];i=pa[i]){
f[pa[i]][i]+=tmd;
f[i][pa[i]]-=tmd;
r[i][pa[i]]=g[i][pa[i]]-f[i][pa[i]];
r[pa[i]][i]=g[pa[i]][i]-f[pa[i]][i];
}
}
return ans;
}
接著是使用鄰接串列,複雜度是\ord{\abs{E}^2\abs{V}}的作法:
#include<cstring>
#include<climits>
#include<queue>
#include<algorithm>
template<typename T>
struct Edmonds_Karp{
static const int MAXN=105;
static const T INF=INT_MAX;
struct edge{
int v,pre;
T cap,r;
edge(int v,int pre,T cap):v(v),pre(pre),cap(cap),r(cap){}
};
int g[MAXN],n;
bool vis[MAXN];
int pa[MAXN];
T d[MAXN];
std::vector<edge> e;
void init(int _n){
memset(g,-1,sizeof(int)*((n=_n)+1));
e.clear();
}
void add_edge(int u,int v,T cap,bool directed=false){
e.push_back(edge(v,g[u],cap));
g[u]=e.size()-1;
e.push_back(edge(u,g[v],directed?0:cap));
g[v]=e.size()-1;
}
T bfs(int s,int t){
memset(vis,0,sizeof(bool)*(n+1));
std::queue<int> q;
q.push(s);
vis[s]=1;
pa[s]=-1;
d[s]=INF;
while(q.size()){
int u=q.front();q.pop();
for(int i=g[u];~i;i=e[i].pre){
if(vis[e[i].v]||e[i].r==0)continue;
vis[e[i].v]=1;
pa[e[i].v]=i;
d[e[i].v]=std::min(d[u],e[i].r);
q.push(e[i].v);
if(e[i].v==t)return d[t];
}
}
return 0;
}
T edmonds_karp(int s,int t){
T ans=0,tmd;
for(;(tmd=bfs(s,t));ans+=tmd){
for(int u=t;~pa[u];){
int i=pa[u];
e[i].r-=tmd;
e[i^1].r+=tmd;
u=e[i^1].v;
}
}
return ans;
}
};
最近幫學長寫作業的時候,又發現了一個不需要多紀錄反邊的做法,但是每個點要同時記錄入邊和出邊,這樣可以用較少的記憶體完成演算法:
#include<cstring>
#include<climits>
#include<queue>
#include<algorithm>
template<typename T>
struct Edmonds_Karp{
static const int MAXN=105;
static const T INF=INT_MAX;
struct Edge{
int u,v;
T cap,r;
Edge(int u,int v,const T &cap):u(u),v(v),cap(cap),r(cap){}
};
std::vector<Edge> edge;
std::vector<int> in[MAXN];
std::vector<int> out[MAXN];
bool vis[MAXN];
int pa[MAXN];
T d[MAXN];
int n;
void init(int _n){
edge.clear();
for(n=_n++;_n--;){
in[_n].clear();
out[_n].clear();
}
}
void add_edge(int u,int v,const T &cap){
in[u].emplace_back(edge.size());
out[v].emplace_back(edge.size());
edge.emplace_back(u,v,cap);
}
T bfs(int s,int t){
memset(vis,0,sizeof(bool)*(n+1));
std::queue<int> q;
q.push(s);
vis[s]=1;
pa[s]=-1;
d[s]=INF;
while(q.size()){
int u=q.front();q.pop();
for(auto i:in[u]){
const auto &e=edge[i];
if(vis[e.v]||e.r==0)continue;
vis[e.v]=1;
pa[e.v]=i;
d[e.v]=std::min(d[u],e.r);
q.push(e.v);
if(e.v==t)return d[t];
}
for(auto i:out[u]){
const auto &e=edge[i];
if(vis[e.u]||e.r==e.cap)continue;
vis[e.u]=1;
pa[e.u]=i;
d[e.u]=std::min(d[u],e.cap-e.r);
q.push(e.u);
if(e.u==t)return d[t];
}
}
return 0;
}
T edmonds_karp(int s,int t){
T ans=0,tmd;
for(;(tmd=bfs(s,t));ans+=tmd){
for(int u=t;~pa[u];){
auto &e=edge[pa[u]];
if(e.v==u){
e.r-=tmd;
u=e.u;
}else{
e.r+=tmd;
u=e.v;
}
}
}
return ans;
}
};

[ Travelling Salesman Problem, TSP ] 旅行推銷員問題

假設圖是完全圖的情況下,從原點出發經過所有點一次並回到原點的最短路徑問題
以下為模板:
#include<limits.h>
#include<algorithm>
#define MAXN 16
int n;/*點數*/
int dp[MAXN][1<<MAXN];/*DP陣列*/
int g[MAXN][MAXN];/*圖*/
int dfs(int s,int d){
if(dp[d][s])return dp[d][s];
if((1<<d)==s)return g[0][d];
dp[d][s]=INT_MAX;
for(int i=0;i<n;++i){
if(((1<<i)&s)&&i!=d){
dp[d][s]=std::min(dp[d][s],dfs(s^(1<<d),i)+g[i][d]);
}
}
return dp[d][s];
}
view raw TSP.cpp hosted with ❤ by GitHub

2015年10月7日 星期三

[ print itself ] 印出自己的code

一個會輸出自己的code:
#include<stdio.h>
#define s "#include<stdio.h>%c#define s %c%s%c%cint main(){%c%cprintf(s,10,34,s,34,10,10,9,10,9,10);%c%creturn 0;%c}"
int main(){
printf(s,10,34,s,34,10,10,9,10,9,10);
return 0;
}
view raw print self.cpp hosted with ❤ by GitHub

可以看一下輸出:http://codepad.org/ZFGeRjkn

2015年10月6日 星期二

[ Rabin-Karp rolling hash ] Rabin-Karp 字串hash演算法

Michael O. Rabin和Richard M. Karp在1987年提出一個想法,即可以對模式串進行哈希運算並將其哈希值與文本中子串的哈希值進行比對。
設字串陣列為S[n],字元編號為0n-1
首先建立陣列Hash[n+1]Hash[i]表示
(S[0]*P^{i-1}+S[1]*P^{i-2}+...+S[i-1])\%PM,Hash[0]=0,其中PPM是質數
假設要取編號L到編號R左閉右開區間的hash值,則其為(Hash[R]-Hash[L]*P^{R-L})\%PM
如果要用double hash則用不同質數計算兩種不同的hash值,用pair拼起來即可
以下提供模板:
template <unsigned P, unsigned PM> class RollingHash {
static vector<unsigned> Base;
vector<unsigned> Hash;
using LL = long long;
public:
RollingHash(const string &S) : Hash(S.size() + 1) {
if (Base.empty())
Base.resize(1, 1);
for (size_t i = 1; i <= S.size(); ++i)
Hash[i] = (LL(Hash[i - 1]) * P + S[i - 1]) % PM;
while (Base.size() < Hash.size())
Base.emplace_back((LL(Base.back()) * P) % PM);
}
size_t size() const { return Hash.size() - 1; }
const vector<unsigned> &getHash() const { return Hash; }
unsigned hash() const { return Hash.back(); }
unsigned hash(unsigned L, unsigned R) const {
return (Hash[R] + PM - (LL(Hash[L]) * Base[R - L]) % PM) % PM;
}
};
template <unsigned P, unsigned PM> vector<unsigned> RollingHash<P, PM>::Base;