Generating Inversion Table

This pro­gram­ming ex­er­cise comes from TAoCP, Vol 3, Section 5.1.1 – 6:

Design an al­go­rithm that com­putes the in­ver­sion table $b_1,b_2\cdots b_n$ corresponding to a given per­mu­ta­tion $a_1{a}_2\cdots a_n$ of $1,2,\cdots ,n$, where the run­ning time is es­sen­tially pro­por­tional to $nlogn$ on typ­i­cal com­put­ers.

Knuth’s so­lu­tion in the book had me stumped for a while. He also men­tions an al­ter­na­tive ap­proach—es­sen­tially a mod­i­fied merge sort. But let’s first un­ravel the bit­wise al­go­rithm.

Implementation

My C++ im­ple­men­ta­tion dif­fers slightly from Knuth’s orig­i­nal: I use 0-indexed ar­rays in­stead of his 1-indexed con­ven­tion. Admittedly, this is­n’t the most el­e­gant ver­sion—my goal was sim­ply to trans­form the pseudocode into some­thing ex­e­cutable.

#include <iostream>
#include <vector>

using namespace std;

int main(int argc, char const* argv[])
{
    // input
    int n; cin >> n;
    vector<int> v(n, 0);
    for (auto& vi: v) cin >> vi;

    // algorithm: TAoCP vol 3. 5.1.1-6
    // init
    vector<int> b(n, 0);
    vector<int> x((n>>1)+1, 0);
    int k = 0, r, s;
    for (int nn=n; nn>1; nn >>= 1) k++; // compute floor(lg(n))

    for (; k >= 0; k--) { // traversal through bits 0 -> \floor(\lg N)
        for (int s = 0; s <= n>>(k+1); s++) // init array x = 0 \forall elements
            x[s] = 0;

        for (int j = 0; j < n; ++j) {
            r = (v[j] >> k) & 1;
            s = v[j] >> (k+1);

            if (r) x[s] += 1;
            else b[v[j]-1] += x[s];
        }
    }

    // output
    for (auto bi: b) cout << bi << " ";
    cout << endl;

    return 0;
}

Analysis

The run­ning time is straight­for­ward to an­a­lyze: the outer loop on $k$ it­er­ates $\lfloor \lg N\rfloor$ times, while the two in­ner loops cost $k+2$ and $N$ op­er­a­tions, re­spec­tively.

But how on earth does this al­go­rithm ac­tu­ally work? The so­lu­tion is too sub­tle to grasp on first en­counter.

Given the in­dex of bit­string $k$, we con­sider 2 strings hav­ing the forms $\overline{s1T}$ and $\overline{s0T}$. Given a fixed $\bar{s}$, $T$ and any $T^{\prime }$, we al­ways know that $\overline{s1T^{\prime }}$ > $\overline{s0T}$. x[] plays a role as a counter which checks how many el­e­ments of $\overline{s1T^{\prime }}$ we have browsed and the in­ver­sion $b[\overline{s0T}]$ up­dates its value by $x[\bar{s}]$.

For ex­am­ple, given the in­put se­quence $5,9,1,8,2,6,4,7,3$, we have the fol­low­ing bi­nary rep­re­sen­ta­tions:

5: 0101
9: 1001
1: 0001
8: 1000
2: 0010
6: 0110
4: 0100
7: 0111
3: 0011

Let’s trace through the al­go­rithm step by step:

1. k=3. $\bar{s}=\bar{0}$, x[0] = 2. b = 1 2 2 2 0 2 2 0 0.
2. k=2. Array x have two items $x[\bar{0}]$, and $x[\bar{1}]$ whose val­ues are 4, 0, re­spec­tively. b = 2 3 6 2 0 2 2 0 0.
3. k=1. There are 3 posi­bil­i­ties of $x[s]$, namely $x[\overline{00}]$, $x[\overline{01}]$, $x[\overline{10}]$ whose val­ues are 2, 2, 0, re­spec­tively. Even­tu­ally, b is 2 3 6 3 0 2 2 0 0.
4. k=0, There are 5 items in x: $x[\overline{000}]=1$, $x[\overline{001}]=1$, $x[\overline{010}]=1$ ,$x[\overline{011}]=1$, $x[\overline{100}]=1$, we get the fi­nal out­put b = 2 3 6 4 0 2 2 1 0.

The merge­sort-based al­go­rithm

Rather than con­struct­ing an in­ver­sion table, this al­go­rithm counts the to­tal num­ber of in­ver­sions in a per­mu­ta­tion by lever­ag­ing the merg­ing pro­ce­dure from the ven­er­a­ble merge­sort.

class inversion {
    private:
        vector<int> x, aux;
        long long int cnt;
    public:
        inversion(const vector<int>& x_): x(x_), aux(x_.size()), cnt(0) {
            merge_sort(x, 0, x.size());
        } //ctor

        void merge_sort(vector<int>& a, int low, int high) {
            if (low >= high-1) return ;

            int mid = low + (high-low)/2;

            merge_sort(a, low, mid);
            merge_sort(a, mid, high);

            merge(a, low, mid, high);

        }

        void merge(vector<int>& a, int low, int mid, int high) {
            int subidx_1  = low;
            int subidx_2  = mid;

            for (int j = low; j < high; ++j) {
                if       (subidx_1 >= subidx_2)       aux[j] = a[subidx_2++];
                else if  (subidx_2 >= high)           aux[j] = a[subidx_1++];
                else if  (a[subidx_1] <= a[subidx_2]) aux[j] = a[subidx_1++];
                else {
                    cnt += (mid-subidx_1); // (1)
                    aux[j] = a[subidx_2++];
                }
            }

            // copy back to a
            copy(aux.begin()+low, aux.begin()+high, a.begin()+low);
        }

        inline long long int count() const { return cnt;}
};

This cor­re­sponds to ex­er­cise 5.2.5 – 21 in TAoCP—unfortunately, Knuth left the so­lu­tion as an ex­er­cise for the reader. The key mod­i­fi­ca­tion is adding line (1) to the merg­ing step, which counts in­ver­sions for a[subidx_2]—the el­e­ment from the sec­ond ar­ray cur­rently be­ing merged.

Here’s the in­sight: el­e­ments from mid to subidx_2 are all less than or equal to a[subidx_2], con­tribut­ing zero in­ver­sions. However, when a[subidx_1] > a[subidx_2], the sorted-ar­ray in­vari­ant tells us that all el­e­ments from subidx_1 to mid must also ex­ceed a[subidx_2]. Hence, we ac­cu­mu­late ex­actly mid - subidx_1 in­ver­sions.

Reference

• gy­wangmtl’s post: This post was in­stru­men­tal in help­ing me un­der­stand the al­go­rithm. And thank you, Google Translate!