### Miracle03's blog

By Miracle03, history, 3 months ago,

Here is the code for undirected case. It's quite similar to directional one, but we only need to keep total degree $\le 2$ for each node. I haven't tested it yet, so I'm not sure if it works well.

To do: 1. do more tests 2. function for hamiltonian cycle 3. When create a cycle, might replace an edge on the cycle with the new edge

#include <bits/stdc++.h>
using namespace std;
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;
}
{
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;
}
}
}
}
vector<vi> used;
unordered_set<int> caneg;
void cut(int a, int b) {
LCT::cut(a, b);
for (int s = 0; s < 2; s++) {
for (int i = 0; i < used[a].size(); i++)
if (used[a][i] == b) {
used[a].erase(used[a].begin() + i);
break;
}
if (used[a].size() == 1) caneg.insert(a);
swap(a, b);
}
}
void link(int a, int b) {
for (int s = 0; s < 2; s++) {
used[a].pb(b);
if (used[a].size() == 2) caneg.erase(a);
swap(a, b);
}
}
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.

LCT::init(n);
if (mx_ch == -1) mx_ch = 1ll * (n + 100) * (n + 50); //default
used.resize(n + 1);
caneg.clear();
for (int i = 1; i <= n; i++) used[i].clear();

vector<vi> edges(n + 1);
for (auto v : eg)
edges[v.fi].pb(v.se),
edges[v.se].pb(v.fi);

for (int i = 1; i <= n; i++)
caneg.insert(i);

int tot = 0;
while (mx_ch >= 0) {
//    cout << tot << ' ' << mx_ch << endl;
vector<pi> eg;
for (auto v : caneg)
for (auto s : edges[v])
eg.pb(mp(v, s));

shuffle(eg.begin(), eg.end(), x);
if (eg.size() == 0) break;
for (auto v : eg) {
mx_ch--;
int a = v.fi, b = v.se;
if (used[a].size() < used[b].size()) swap(a, b);
if (used[b].size() >= 2) continue;
if (x() & 1) continue;
if (LCT::fdr(a) == LCT::fdr(b)) continue;
if (used[a].size() < 2 && used[b].size() < 2)
tot++;
if (used[a].size() == 2) {
int p = used[a][x() % 2];
cut(a, p);
}
}
if (tot == n - 1) {
vi cur;
for (int i = 1; i <= n; i++)
if (used[i].size() <= 1) {
int pl = i, ls = 0;
while (pl) {
cur.pb(pl);
int flag = 0;
for (auto v : used[pl])
if (v != ls) {
ls = pl;
pl = v;
flag = 1;
break;
}
if (!flag) break;
}
break;
}
return cur;
}
}
//failed to find a path
return vi();
}
}

• +108

By Miracle03, history, 3 months ago,

(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;
}
{
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);
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;
}
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.