Count number of permutations generated by particular interleaves

permute

#include <vector>
#include <stdint.h>
#include <stdio.h>

static const int NUM_WORDS = 3; // "n" in the article
static const int NUM_ELEMS_PER_WORD = 4; // "k" in the article

static const int NUM_ELEMS = NUM_ELEMS_PER_WORD * NUM_WORDS;

typedef uint32_t PermIndex;
typedef uint8_t Perm[NUM_ELEMS];
typedef uint8_t WordPerm[NUM_WORDS];

static void interleave_step(Perm& pOut, const Perm& pIn)
{
    for (int i = 0; i < NUM_ELEMS_PER_WORD; i++) {
        pOut[i*2 + 0] = pIn[i];
        pOut[i*2 + 1] = pIn[i + NUM_ELEMS_PER_WORD];
    }

    for (int i = NUM_ELEMS_PER_WORD*2; i < NUM_ELEMS; i++)
        pOut[i] = pIn[i];
}

static void permute_words(Perm& pOut, const Perm& pIn, const WordPerm& pWords)
{
    for (int word = 0; word < NUM_WORDS; word++) {
        int iIn = pWords[word] * NUM_ELEMS_PER_WORD;
        int iOut = word * NUM_ELEMS_PER_WORD;
        for (int i = 0; i < NUM_ELEMS_PER_WORD; i++)
            pOut[iOut + i] = pIn[iIn + i];
    }
}

// Permutation ranking / unranking
static void mr_unrank1(PermIndex rank, unsigned int n, uint8_t *vec)
{
    if (n >= 1) {
        std::swap(vec[rank % n], vec[n-1]);
        mr_unrank1(rank / n, n-1, vec);
    }
}

static PermIndex mr_rank1(unsigned int n, uint8_t *vec, uint8_t *inv)
{
    if (n >= 2) {
        PermIndex s = vec[n-1];
        std::swap(vec[n-1], vec[inv[n-1]]);
        std::swap(inv[s], inv[n-1]);
        return s + n * mr_rank1(n-1, vec, inv);
    } else
        return 0;
}

static void get_permutation(PermIndex rank, unsigned int n, uint8_t *vec)
{
    for (unsigned int i = 0; i < n; i++)
        vec[i] = (uint8_t) i;

    mr_unrank1(rank, n, vec);
}

static void print_permutation(unsigned int n, const uint8_t *vec)
{
    for (unsigned int i = 0; i < n; i++)
        fputc("0123456789abcdef"[vec[i]], stdout);
    fputc('\n', stdout);
}

static PermIndex get_rank(const Perm& p)
{
    uint8_t vec[NUM_ELEMS], inv[NUM_ELEMS];
    for (unsigned int i = 0; i < NUM_ELEMS; i++) {
        vec[i] = p[i];
        inv[vec[i]] = i;
    }
    return mr_rank1(NUM_ELEMS, vec, inv);
}

static bool is_set(const std::vector<uint8_t>& v, unsigned int i)
{
    return (v[i/8] & (1 << (i % 8))) != 0;
}

static void set_bit(std::vector<uint8_t>& v, unsigned int i)
{
    v[i/8] |= 1 << (i % 8);
}

static PermIndex factorial(PermIndex n)
{
    return (n < 2) ? 1 : n * factorial(n - 1);
}

int main()
{
    // generate all permutations of NUM_WORDS elements (used later)
    unsigned int nwperms = factorial(NUM_WORDS);
    printf("%d word perms.\n", nwperms);

    // we don't need the last one (identity).
    nwperms--;
    WordPerm *wperms = new WordPerm[nwperms];
    for (unsigned int i = 0; i < nwperms; i++) {
        get_permutation(i, NUM_WORDS, wperms[i]);
        print_permutation(NUM_WORDS, wperms[i]);
    }

    // set up the bit field for found perms
    PermIndex max_count = factorial(NUM_ELEMS);

    size_t num_bytes = (max_count + 7) / 8;
    std::vector<uint8_t> found;
    std::vector<PermIndex> queue;

    found.resize(num_bytes, 0);

    // initialize
    queue.push_back(0);
    set_bit(found, 0); // we found the identity permutation.

    while (!queue.empty()) {
        Perm cur, out;
        get_permutation(queue.back(), NUM_ELEMS, cur);
        queue.pop_back();

        //print_permutation(NUM_ELEMS, cur);

        // at any step, we can either try one of the word permutations or interleave
        for (unsigned int p = 0; p < nwperms; p++) {
            permute_words(out, cur, wperms[p]);
            PermIndex rank = get_rank(out);
            if (!is_set(found, rank)) {
                set_bit(found, rank);
                queue.push_back(rank);
            }
        }

        interleave_step(out, cur);
        PermIndex rank = get_rank(out);
        if (!is_set(found, rank)) {
            set_bit(found, rank);
            queue.push_back(rank);
        }
    }

    // count permutations found
    // quick LUT for counting bits set in a byte (good enough for our purposes)
    uint8_t bitcount[256];
    for (int i = 0; i < 256; i++) {
        int count = 0;
        for (int v = i; v; v &= v - 1)
            count++;

        bitcount[i] = count;
    }

    unsigned long long total = 0;
    for (size_t i = 0; i < num_bytes; i++)
        total += bitcount[found[i]];

    printf("found %lld perms total.\n", total);
    delete[] wperms;

    return 0;
}