(It's a translated version of my blog at https://peehs-moorhsum.blog.uoj.ac/blog/6384. )↵
↵
Finding a Hamiltonian path in a directed graph is a well-known NP problem. ↵
↵
However, about a year ago, I came up with the following heuristic algorithm which has GREAT performance on random graphs(by first generating a hamiltonian path, adding random edges, then randomly permuting indices) and many CP problems. ↵
↵
Here is the idea of the algorithm:↵
↵
We maintain a subset $S$ of edges, ensuring there's no cycle and each vertex has in-degree and out-degree at most 1. Whenever $|S|$ reaches $n-1$, we find a hamiltonian path. ↵
↵
We start with $S$ empty. Each time we pick a random edge $e \notin S$. If we can add $e$ to $S$ without violating any rule, we simply add it to $S$. Otherwise, if $e$ won't create a cycle and violate with exactly one edge e', we replace e' with $e$ with 50% probability. Repeat until $|S|$ reaches $n-1$. ↵
↵
We can use LCT to check cycles. Here's my code, which takes less than 1 second on most random graphs with $|V|=100000, |E|=500000$:↵
↵
~~~~~↵
namespace hamil {↵
template <typename T> bool chkmax(T &x,T y){return x<y?x=y,true:false;}↵
template <typename T> bool chkmin(T &x,T y){return x>y?x=y,true:false;}↵
#define vi vector<int>↵
#define pb push_back↵
#define mp make_pair↵
#define pi pair<int, int>↵
#define fi first↵
#define se second↵
#define ll long long↵
namespace LCT {↵
vector<vi> ch;↵
vi fa, rev;↵
void init(int n) {↵
ch.resize(n + 1);↵
fa.resize(n + 1);↵
rev.resize(n + 1);↵
for (int i = 0; i <= n; i++)↵
ch[i].resize(2), ↵
ch[i][0] = ch[i][1] = fa[i] = rev[i] = 0;↵
}↵
bool isr(int a)↵
{↵
return !(ch[fa[a]][0] == a || ch[fa[a]][1] == a);↵
} ↵
void pushdown(int a)↵
{↵
if(rev[a])↵
{↵
rev[ch[a][0]] ^= 1, rev[ch[a][1]] ^= 1;↵
swap(ch[a][0], ch[a][1]);↵
rev[a] = 0;↵
}↵
}↵
void push(int a)↵
{↵
if(!isr(a)) push(fa[a]);↵
pushdown(a); ↵
}↵
void rotate(int a)↵
{↵
int f = fa[a], gf = fa[f];↵
int tp = ch[f][1] == a;↵
int son = ch[a][tp ^ 1];↵
if(!isr(f)) ↵
ch[gf][ch[gf][1] == f] = a; ↵
fa[a] = gf;↵
↵
ch[f][tp] = son;↵
if(son) fa[son] = f;↵
↵
ch[a][tp ^ 1] = f, fa[f] = a;↵
}↵
void splay(int a)↵
{↵
push(a);↵
while(!isr(a))↵
{↵
int f = fa[a], gf = fa[f];↵
if(isr(f)) rotate(a);↵
else↵
{↵
int t1 = ch[gf][1] == f, t2 = ch[f][1] == a;↵
if(t1 == t2) rotate(f), rotate(a);↵
else rotate(a), rotate(a); ↵
} ↵
} ↵
}↵
void access(int a)↵
{↵
int pr = a;↵
splay(a);↵
ch[a][1] = 0;↵
while(1)↵
{↵
if(!fa[a]) break; ↵
int u = fa[a];↵
splay(u);↵
ch[u][1] = a;↵
a = u;↵
}↵
splay(pr);↵
}↵
void makeroot(int a)↵
{↵
access(a);↵
rev[a] ^= 1;↵
}↵
void link(int a, int b)↵
{↵
makeroot(a);↵
fa[a] = b;↵
}↵
void cut(int a, int b)↵
{↵
makeroot(a);↵
access(b);↵
fa[a] = 0, ch[b][0] = 0;↵
}↵
int fdr(int a)↵
{↵
access(a);↵
while(1)↵
{↵
pushdown(a);↵
if (ch[a][0]) a = ch[a][0];↵
else {↵
splay(a);↵
return a;↵
}↵
}↵
}↵
}↵
vi out, in;↵
vi work(int n, vector<pi> eg, ll mx_ch = -1) { ↵
// mx_ch : max number of adding/replacing default is (n + 100) * (n + 50) ↵
// n : number of vertices. 1-indexed. ↵
// eg: vector<pair<int, int> > storing all the edges. ↵
// return a vector<int> consists of all indices of vertices on the path. return empty list if failed to find one. ↵
out.resize(n + 1), in.resize(n + 1);↵
LCT::init(n);↵
for (int i = 0; i <= n; i++) in[i] = out[i] = 0;↵
if (mx_ch == -1) mx_ch = 1ll * (n + 100) * (n + 50); //default↵
vector<vi> from(n + 1), to(n + 1);↵
for (auto v : eg)↵
from[v.fi].pb(v.se), ↵
to[v.se].pb(v.fi);↵
unordered_set<int> canin, canout;↵
for (int i = 1; i <= n; i++)↵
canin.insert(i), ↵
canout.insert(i); ↵
mt19937 x(time(0));↵
int tot = 0;↵
while (mx_ch >= 0) {↵
// cout << tot << ' ' << mx_ch << endl;↵
vector<pi> eg;↵
for (auto v : canout)↵
for (auto s : from[v])↵
if (in[s] == 0) {↵
assert(canin.count(s));↵
continue;↵
}↵
else eg.pb(mp(v, s));↵
for (auto v : canin)↵
for (auto s : to[v])↵
eg.pb(mp(s, v));↵
shuffle(eg.begin(), eg.end(), x);↵
if (eg.size() == 0) break;↵
for (auto v : eg) {↵
mx_ch--;↵
if (in[v.se] && out[v.fi]) continue;↵
if (LCT::fdr(v.fi) == LCT::fdr(v.se)) continue;↵
if (in[v.se] || out[v.fi]) ↵
if (x() & 1) continue;↵
if (!in[v.se] && !out[v.fi]) ↵
tot++;↵
if (in[v.se]) {↵
LCT::cut(in[v.se], v.se);↵
canin.insert(v.se);↵
canout.insert(in[v.se]);↵
out[in[v.se]] = 0;↵
in[v.se] = 0;↵
}↵
if (out[v.fi]) {↵
LCT::cut(v.fi, out[v.fi]);↵
canin.insert(out[v.fi]);↵
canout.insert(v.fi);↵
in[out[v.fi]] = 0;↵
out[v.fi] = 0;↵
}↵
LCT::link(v.fi, v.se);↵
canin.erase(v.se);↵
canout.erase(v.fi);↵
in[v.se] = v.fi;↵
out[v.fi] = v.se;↵
}↵
if (tot == n - 1) {↵
vi cur;↵
for (int i = 1; i <= n; i++) ↵
if (!in[i]) {↵
int pl = i;↵
while (pl) {↵
cur.pb(pl), ↵
pl = out[pl];↵
}↵
break;↵
} ↵
return cur;↵
}↵
}↵
//failed to find a path↵
return vi();↵
}↵
}↵
~~~~~↵
↵
Note that if the start vertex $S$ and/or end vertex $T$ of the path is given, we simply create two new vertices $A, B$ and add edges from $A$ to $S$, from $T$ to $B$. Any Hamiltonian path in the new graph automatically has $A\to S$ as the first edge and $T\to B$ as the last edge. If we want to find a Hamiltonian cycle, we simply enumerate an edge, then convert the problem to finding a path given start/end vertices. For bidirectional case, we can simply add each edge twice. ↵
↵
Finding a Hamiltonian path in a directed graph is a well-known NP problem. ↵
↵
However, about a year ago, I came up with the following heuristic algorithm which has GREAT performance on random graphs(by first generating a hamiltonian path, adding random edges, then randomly permuting indices) and many CP problems. ↵
↵
Here is the idea of the algorithm:↵
↵
We maintain a subset $S$ of edges, ensuring there's no cycle and each vertex has in-degree and out-degree at most 1. Whenever $|S|$ reaches $n-1$, we find a hamiltonian path. ↵
↵
We start with $S$ empty. Each time we pick a random edge $e \notin S$. If we can add $e$ to $S$ without violating any rule, we simply add it to $S$. Otherwise, if $e$ won't create a cycle and violate with exactly one edge e', we replace e' with $e$ with 50% probability. Repeat until $|S|$ reaches $n-1$. ↵
↵
We can use LCT to check cycles. Here's my code, which takes less than 1 second on most random graphs with $|V|=100000, |E|=500000$:↵
↵
~~~~~↵
namespace hamil {↵
template <typename T> bool chkmax(T &x,T y){return x<y?x=y,true:false;}↵
template <typename T> bool chkmin(T &x,T y){return x>y?x=y,true:false;}↵
#define vi vector<int>↵
#define pb push_back↵
#define mp make_pair↵
#define pi pair<int, int>↵
#define fi first↵
#define se second↵
#define ll long long↵
namespace LCT {↵
vector<vi> ch;↵
vi fa, rev;↵
void init(int n) {↵
ch.resize(n + 1);↵
fa.resize(n + 1);↵
rev.resize(n + 1);↵
for (int i = 0; i <= n; i++)↵
ch[i].resize(2), ↵
ch[i][0] = ch[i][1] = fa[i] = rev[i] = 0;↵
}↵
bool isr(int a)↵
{↵
return !(ch[fa[a]][0] == a || ch[fa[a]][1] == a);↵
} ↵
void pushdown(int a)↵
{↵
if(rev[a])↵
{↵
rev[ch[a][0]] ^= 1, rev[ch[a][1]] ^= 1;↵
swap(ch[a][0], ch[a][1]);↵
rev[a] = 0;↵
}↵
}↵
void push(int a)↵
{↵
if(!isr(a)) push(fa[a]);↵
pushdown(a); ↵
}↵
void rotate(int a)↵
{↵
int f = fa[a], gf = fa[f];↵
int tp = ch[f][1] == a;↵
int son = ch[a][tp ^ 1];↵
if(!isr(f)) ↵
ch[gf][ch[gf][1] == f] = a; ↵
fa[a] = gf;↵
↵
ch[f][tp] = son;↵
if(son) fa[son] = f;↵
↵
ch[a][tp ^ 1] = f, fa[f] = a;↵
}↵
void splay(int a)↵
{↵
push(a);↵
while(!isr(a))↵
{↵
int f = fa[a], gf = fa[f];↵
if(isr(f)) rotate(a);↵
else↵
{↵
int t1 = ch[gf][1] == f, t2 = ch[f][1] == a;↵
if(t1 == t2) rotate(f), rotate(a);↵
else rotate(a), rotate(a); ↵
} ↵
} ↵
}↵
void access(int a)↵
{↵
int pr = a;↵
splay(a);↵
ch[a][1] = 0;↵
while(1)↵
{↵
if(!fa[a]) break; ↵
int u = fa[a];↵
splay(u);↵
ch[u][1] = a;↵
a = u;↵
}↵
splay(pr);↵
}↵
void makeroot(int a)↵
{↵
access(a);↵
rev[a] ^= 1;↵
}↵
void link(int a, int b)↵
{↵
makeroot(a);↵
fa[a] = b;↵
}↵
void cut(int a, int b)↵
{↵
makeroot(a);↵
access(b);↵
fa[a] = 0, ch[b][0] = 0;↵
}↵
int fdr(int a)↵
{↵
access(a);↵
while(1)↵
{↵
pushdown(a);↵
if (ch[a][0]) a = ch[a][0];↵
else {↵
splay(a);↵
return a;↵
}↵
}↵
}↵
}↵
vi out, in;↵
vi work(int n, vector<pi> eg, ll mx_ch = -1) { ↵
// mx_ch : max number of adding/replacing default is (n + 100) * (n + 50) ↵
// n : number of vertices. 1-indexed. ↵
// eg: vector<pair<int, int> > storing all the edges. ↵
// return a vector<int> consists of all indices of vertices on the path. return empty list if failed to find one. ↵
out.resize(n + 1), in.resize(n + 1);↵
LCT::init(n);↵
for (int i = 0; i <= n; i++) in[i] = out[i] = 0;↵
if (mx_ch == -1) mx_ch = 1ll * (n + 100) * (n + 50); //default↵
vector<vi> from(n + 1), to(n + 1);↵
for (auto v : eg)↵
from[v.fi].pb(v.se), ↵
to[v.se].pb(v.fi);↵
unordered_set<int> canin, canout;↵
for (int i = 1; i <= n; i++)↵
canin.insert(i), ↵
canout.insert(i); ↵
mt19937 x(time(0));↵
int tot = 0;↵
while (mx_ch >= 0) {↵
// cout << tot << ' ' << mx_ch << endl;↵
vector<pi> eg;↵
for (auto v : canout)↵
for (auto s : from[v])↵
if (in[s] == 0) {↵
assert(canin.count(s));↵
continue;↵
}↵
else eg.pb(mp(v, s));↵
for (auto v : canin)↵
for (auto s : to[v])↵
eg.pb(mp(s, v));↵
shuffle(eg.begin(), eg.end(), x);↵
if (eg.size() == 0) break;↵
for (auto v : eg) {↵
mx_ch--;↵
if (in[v.se] && out[v.fi]) continue;↵
if (LCT::fdr(v.fi) == LCT::fdr(v.se)) continue;↵
if (in[v.se] || out[v.fi]) ↵
if (x() & 1) continue;↵
if (!in[v.se] && !out[v.fi]) ↵
tot++;↵
if (in[v.se]) {↵
LCT::cut(in[v.se], v.se);↵
canin.insert(v.se);↵
canout.insert(in[v.se]);↵
out[in[v.se]] = 0;↵
in[v.se] = 0;↵
}↵
if (out[v.fi]) {↵
LCT::cut(v.fi, out[v.fi]);↵
canin.insert(out[v.fi]);↵
canout.insert(v.fi);↵
in[out[v.fi]] = 0;↵
out[v.fi] = 0;↵
}↵
LCT::link(v.fi, v.se);↵
canin.erase(v.se);↵
canout.erase(v.fi);↵
in[v.se] = v.fi;↵
out[v.fi] = v.se;↵
}↵
if (tot == n - 1) {↵
vi cur;↵
for (int i = 1; i <= n; i++) ↵
if (!in[i]) {↵
int pl = i;↵
while (pl) {↵
cur.pb(pl), ↵
pl = out[pl];↵
}↵
break;↵
} ↵
return cur;↵
}↵
}↵
//failed to find a path↵
return vi();↵
}↵
}↵
~~~~~↵
↵
Note that if the start vertex $S$ and/or end vertex $T$ of the path is given, we simply create two new vertices $A, B$ and add edges from $A$ to $S$, from $T$ to $B$. Any Hamiltonian path in the new graph automatically has $A\to S$ as the first edge and $T\to B$ as the last edge. If we want to find a Hamiltonian cycle, we simply enumerate an edge, then convert the problem to finding a path given start/end vertices. For bidirectional case, we can simply add each edge twice. ↵
