Heuristic algorithm for Hamiltonian path in directed graphs
Difference between en1 and en2, changed 0 character(s)
(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. ↵

History

 
 
 
 
Revisions
 
 
  Rev. Lang. By When Δ Comment
en3 English Rewinding 2021-05-10 04:17:18 49
en2 English Rewinding 2021-05-09 17:16:19 0 (published)
en1 English Rewinding 2021-05-08 17:32:44 7673 Initial revision (saved to drafts)