Nima Aghdaii


Calculating Mahonian Numbers

27 Dec 2020 » codechef, math, algorithm, problems

Problem Inversions

Find the number of permutations of \({1,2,...,N}\) with exactly \(K\) inversions (modulus \(2\)).


Summary of my approach

Step 1

Figure it is called “mahonian triangle numbers”, either by using OEIS or by reading at the bottom of this wikipedia page inversions.

Step 2

Find the relationship with “Truncated Euler Product Triangle” from this paper.

You can basically write mahonian numbers as: \(M(r, c) = \sum_{k=0}^{c} {r+c-k \choose r} P(r, k)\)

Step 3

Realize that if \(c < r\) then \(P(r, k)\) becomes only a function of \(k\) (Every row is a prefix of the next row!). The paper also points out that: \(P(r, c) = p(c)\) for \(c < r + 2\).

Step 4

Since we want the answer modulus 2, \({N \choose K} \bmod 2\) can be done in \(O(1)\). Read my other post to see how.

Step 5

Finally, if you generate p(c) for c up to 100, you realize that it is very sparse!

1 -1 -1 0 0 1 0 1 0 0 0 0 -1 0 0 -1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ...

Doing some research, it turns out the non-zero indices are a series called Generalized Pentagonal Numbers and there is an efficient way to generate them one by one. Using Pentagonal Number Theorem, you can generate the non-zero indices. The k-th index \(G_k\) is calculated as: \(G_{k} = \frac{k(3k − 1)}{2}\) for k = 1, −1, 2, −2, 3, -3, …


Code and Time Complexity

There are \(\sqrt{N}\) non-zero terms and for each term we do \(O(1)\) calculation. So the total time complexity is \(O(\sqrt{N})\).

Here is the code I submitted during the contest and below is a cleaned up version of it.

#include <iostream>

using namespace std;

typedef long long ll;

ll comb2(ll n, ll p) {
    return (p & (n - p)) ? 0 : 1;
}

ll g(ll n) {
    return (3 * n * n - n) / 2;
}

ll calc(ll r, ll c, ll n, ll num) {
    return (2 + comb2(r + c - n, r) * (num % 2 == 0 ? 1 : -1) ) % 2;
}

ll mahonian(ll r, ll c) {
    ll res = calc(r, c, 0, 0);
    for (ll num = 1; ; num++) {
        for (ll mul: {1, -1}) {
            ll nxt = g(num * mul);
            if (nxt >= c + 1) return res;
            res = (res + calc(r, c, nxt, num)) % 2;
        }
    }
}

int main() {
    ll tc, n, k;
    cin >> tc;
    while (tc--) {
        cin >> n >> k;
        cout << mahonian(n-1, k) << endl;
    }
    return 0;
}