Tag Archives: Mathematics

Combinatorial algorithms

Combinatorial algorithms are algorithms that deal with combinatorial structures, which are sets, ordered n-tuples, and any structures that can be built from them, like graphs.

Combinatorial algorithms include algorithms for:

Recursive integer partitions in C

Here is a recursive algorithm to generate integer partitions in antilexicographic order. The function partitions() takes an integer to partition, and a callback function to call for each partition found.

#include <stdlib.h>

unsigned int min(int a, int b)
{
    return a < b ? a : b;
}

typedef void (*partitionfn)(const unsigned int *, size_t);

static void partitions_recursive(unsigned int n, unsigned int maximum, unsigned int *partition,
        size_t length, partitionfn fun)
{
    unsigned int i;
    if (n == 0) {
        fun(partition, length);
    }
    for (i = min(maximum, n); i >= 1; i--) {
        partition[length] = i;
        partitions_recursive(n - i, i, partition, length + 1, fun);
    }
}

void partitions(unsigned int n, partitionfn fun)
{
    unsigned int *partition = malloc(n * sizeof(unsigned int));
    if (partition) {
        partitions_recursive(n, n, partition, 0, fun);
        free(partition);
    }
}

An example program:

#include <stdio.h>

void print(const unsigned int *partition, size_t length)
{
    unsigned int i;
    for (i = 0; i < length; i++) {
        printf("%u ", partition[i]);
    }
    putchar('\n');
}

int main(void)
{
    partitions(6, print);
    return 0;
}

The output:

6
5 1
4 2
4 1 1
3 3
3 2 1
3 1 1 1
2 2 2
2 2 1 1
2 1 1 1 1
1 1 1 1 1 1

Reference: Partition.java

Permutation cycles in C

A cycle of a permutation is a subset of the elements that replace one another in sequence, until the last element is replaced by the first. For example, consider the permutation below:

\(\sigma=\begin{pmatrix}
0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\
1 & 3 & 7 & 0 & 4 & 8 & 2 & 6 & 5\end{pmatrix}\)

We can find the cycles:
\(0 \rightarrow 1, 1 \rightarrow 3, 3 \rightarrow 0\)
\(2 \rightarrow 7, 7 \rightarrow 6, 6 \rightarrow 2\)
\(4 \rightarrow 4\)
\(5 \rightarrow 8, 8 \rightarrow 5\)

These can be written as:
\((0, 1, 3)(2, 7, 6)(4)(5, 8)\)

It’s customary to omit cycles of length 1, so this would usually be written as
\((0, 1, 3)(2, 7, 6)(5, 8).\)

To find the cycle decomposition of a permutation, we can use an algorithm that is very similar to depth-first search (DFS) on a graph. We begin a new search for each unvisited vertex (number), and visit its neighbour (image in the permutation) until we get back to the first vertex again.

Below is an implementation in C. The function permutation_cycles() takes a permutation in the form of integers starting from 0, its length, and a callback function that is called for each cycle found.

#include <stdlib.h>

typedef void(*cyclefn)(const unsigned int *, size_t);

void permutation_cycles_recursive(const unsigned int *permutation, unsigned int *visited, 
        unsigned int start, unsigned int current, unsigned int *path,
        size_t length, cyclefn fun)
{
    visited[current] = 1;
    path[length] = current;
    if (start == current && length > 0) {
        fun(path, length);
    }
    else {
        permutation_cycles_recursive(permutation, visited, start, permutation[current],
                path, length + 1, fun);
    }
}

void permutation_cycles(const unsigned int *permutation, size_t n, cyclefn fun)
{
    unsigned int i;
    unsigned int *visited = calloc(n, sizeof(unsigned int));
    unsigned int *path = malloc(n * sizeof(unsigned int));
    if (!(visited && path)) {
        free(visited);
        free(path);
        return;
    }
    for (i = 0; i < n; i++) {
        if (!visited[i]) {
            permutation_cycles_recursive(permutation, visited, i, i, 
                    path, 0, fun);
        }
    }
    free(visited);
    free(path);
}

Here is an example program that finds the cycle decomposition of the permutation shown above:

#include <stdio.h>

void print(const unsigned int *cycle, size_t length)
{
    if (length > 1) {
        unsigned int i;
        putchar('(');
        for (i = 0; i < length; i++) {
            printf("%u", cycle[i]);
            if (i < length - 1) {
                printf(", ");
            }
        }
        putchar(')');
    }
}

int main(void)
{
    unsigned int permutation[] = {1, 3, 7, 0, 4, 8, 2, 6, 5};
    const size_t n = sizeof(permutation) / sizeof(permutation[0]);
    permutation_cycles(permutation, n, print);
    putchar('\n');
    return 0;
}

The output:

(0, 1, 3)(2, 7, 6)(5, 8)

Prüfer codes and labeled trees in C

A tree and its Prüfer Code

Arthur Cayley first proved that there are \(n^{n-2}\) labeled trees on \(n\) vertices, and in another proof of the same theorem, Heinz Prüfer showed a bijection between trees on \(n\) vertices and \((n-2)\)-tuples of integers from 1 to \(n – 2\), called Prüfer codes. The proof provides an algorithm that allows us to construct an arbitrary labeled tree from its code.

The algorithm proceeds as follows:

  1. Initialise an array of node degrees for the nodes, which are numbered from 1 to \(n + 2\)
  2. For each occurrence of a node number in the code, increment its degree
  3. Add 1 to all degrees
  4. For each occurrence of a node number in the code:
    1. Find the lowest-numbered node with degree 1
    2. Create an edge between that node and the node in the code
    3. Decrement the degrees of both nodes by 1
  5. There are now 2 nodes left with degree 1, so join them with an edge to complete the tree

Below is an implementation of the algorithm in C. The function prufer_tree() takes a code in the form of an array of unsigned integers and its length, and returns the tree as an array of edges.

#include <stdlib.h>

typedef struct {
    unsigned int first;
    unsigned int second;
} edge;

edge *prufer_tree(const unsigned int *code, size_t len)
{
    unsigned int i, j, e = 0;
    const unsigned int n = len + 2; /* Nodes */
    edge *tree = malloc((n - 1) * sizeof(edge));
    unsigned int *degrees = calloc(n + 1, sizeof(unsigned int)); /* 1-based */
    if (tree == NULL || degrees == NULL) {
        free(tree);
        free(degrees);
        return NULL;
    }
    /* Set the degrees from the occurrences in the code */
    for (i = 0; i < len; i++) {
        degrees[ code[i]]++;
    }
    /* Add 1 to them all */
    for (i = 1; i <= n; i++) {
        degrees[i]++;
    }
    /* Add edges to nodes in the code */
    for (i = 0; i < len; i++) {
        /* Find the lowest-numbered node with degree 1 */
        for (j = 1; degrees[j] != 1; j++);
        /* Add the edge */
        tree[e].first = j;
        tree[e].second = code[i];
        degrees[j]--;
        degrees[ code[i]]--;
        e++;
    }
    /* Find the last 2 degree-1 nodes */
    for (i = 1; degrees[i] != 1; i++);
    for (j = i + 1; degrees[j] != 1; j++);
    /* Add the last edge */
    tree[e].first = i;
    tree[e].second = j;
    free(degrees);

    return tree;
}

Here is an example program that constructs the tree in the image at the top:

#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    unsigned int code[] = {1, 1, 1, 1, 6, 5};
    const size_t len = sizeof(code) / sizeof(unsigned int);
    const unsigned int n = len + 1; /* Edges */
    edge *tree = prufer_tree(code, len);
    unsigned int e;
    for (e = 0; e < n; e++) {
        printf("(%u, %u)\n", tree[e].first, tree[e].second);
    }
    free(tree);
    return 0;
}

The output:

(2, 1)
(3, 1)
(4, 1)
(7, 1)
(1, 6)
(6, 5)
(5, 8)

Partial Latin square completion using backtracking in C

The partial Latin square completion problem is to take a partially-filled Latin square and to determine if it can be completed successfully. For example, below is a 1-based order 5 partial Latin square:

\(
\left( \begin{array}{ccc}
\cdot & 2 & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & \cdot & 3 \\
\cdot & \cdot & 5 & \cdot & \cdot \\
4 & \cdot & \cdot & \cdot & \cdot \\
\cdot & \cdot & \cdot & 2 & \cdot \end{array} \right)
\)

 

The problem is to find out if it can be completed, which in this case it can, to this square:

\(
\left( \begin{array}{ccc}
1 & 2 & 3 & 4 & 5 \\
2 & 1 & 4 & 5 & 3 \\
3 & 4 & 5 & 1 & 2 \\
4 & 5 & 2 & 3 & 1 \\
5 & 3 & 1 & 2 & 4 \end{array} \right)
\)

 

The partial Latin square completion problem is NP-complete, but can be solved using backtracking.

Here is an implementation in C. It takes a partial Latin square in the form of a 1-dimensional array, with 0s for the unfilled cells, and a callback function which is called with the solution, if found.

#include <stdlib.h>

static void get_row_column(unsigned int cell, unsigned int *row, unsigned int *column, size_t size)
{
    *row = cell / size;
    *column =  cell % size;
}

static unsigned int get_value(const unsigned int *square, size_t size, unsigned int row,
        unsigned int column)
{
    return square[row * size + column];
}

static unsigned int unique_in_column(const unsigned int *square, size_t size, unsigned int cell)
{
    unsigned int row;
    unsigned int column;
    get_row_column(cell, &row, &column, size);
    unsigned int r;
    unsigned int unique = 1;
    for (r = 0; r < size && unique; r++) {
        unique = r == row || get_value(square, size, r, column) != square[cell];
    } 
    return unique;
}

static unsigned int unique_in_row(const unsigned int *square, size_t size, unsigned int cell)
{
    unsigned int row;
    unsigned int column;
    get_row_column(cell, &row, &column, size);
    unsigned int c;
    unsigned int unique = 1;
    for (c = 0; c < size && unique; c++) {
        unique = c == column || get_value(square, size, row, c) != square[cell];
    } 
    return unique;
}

static unsigned int promising(const unsigned int *square, size_t size, unsigned int cell)
{
    return square[cell] == 0
        || (unique_in_row(square, size, cell)
                && unique_in_column(square, size, cell));
}

static unsigned int get_next_cell(const unsigned int *square, size_t size, unsigned int *cell)
{
    unsigned int found = 0;
    unsigned int i;
    for (i = *cell; i < size * size && !found; i++) {
        if (square[i] == 0) {
            found = 1;
            *cell = i;
        }
    }
    return found;
}

static void latin_square_completion_recursive(unsigned int *square, size_t size, unsigned int cell,
        unsigned int *succeeded, latin_square_completionfn fun)
{
    if (!get_next_cell(square, size, &cell)) {
        *succeeded = 1;
        fun(square, size);
    }
    else {
        unsigned int i;
        for (i = 1; i <= size && !*succeeded; i++) {
            unsigned int temp = square[cell];
            square[cell] = i;
            if (promising(square, size, cell)) {
                latin_square_completion_recursive(square, size, cell + 1, succeeded, fun);
            }
            square[cell] = temp;
        }
    }
}

void latin_square_completion(unsigned int *square, size_t size,
        latin_square_completionfn fun)
{
    unsigned int cell = 0;
    unsigned int succeeded = 0;
    latin_square_completion_recursive(square, size, cell, &succeeded, fun);
}

Here is a program to complete the partial Latin square shown at the beginning:

#include <stdio.h>

static void print(const unsigned int *square, size_t size)
{
    unsigned int i;
    for (i = 0; i < size * size; i++) {
        printf("%d ", square[i]);
        if (i % size == size - 1) {
            printf("\n");
        }
    }
    printf("\n");
}

int main(void)
{
    unsigned int square[] = {0, 2, 0, 0, 0,
                             0, 0, 0, 0, 3,
                             0, 0, 5, 0, 0,
                             4, 0, 0, 0, 0,
                             0, 0, 0, 2, 0};
    latin_square_completion(square, 5, print);
    return 0;
}

The output:

1 2 3 4 5
2 1 4 5 3
3 4 5 1 2
4 5 2 3 1
5 3 1 2 4

Integer partitions in Python

These functions produce integer partitions in reverse lexicographic and lexicographic order respectively.

def revlex_partitions(n):
    """
    Generate all integer partitions of n
    in reverse lexicographic order
    """
    if n == 0:
        yield []
        return
    for p in revlex_partitions(n - 1):
        if len(p) == 1 or (len(p) > 1 and p[-1] < p[-2]):
            p[-1] += 1
            yield p
            p[-1] -= 1
        p.append(1)
        yield p
        p.pop()

def lex_partitions(n):
    """
    Generate all integer partitions of n
    in lexicographic order
    """
    if n == 0:
        yield []
        return
    for p in lex_partitions(n - 1):
        p.append(1)
        yield p
        p.pop()
        if len(p) == 1 or (len(p) > 1 and p[-1] < p[-2]):
            p[-1] += 1
            yield p
            p[-1] -= 1

Example program:

def main():
    for p in revlex_partitions(5):
        print(p)
    print()
    for p in lex_partitions(5):
        print(p)

if __name__ == '__main__':
    main()

Output:

[5]
[4, 1]
[3, 2]
[3, 1, 1]
[2, 2, 1]
[2, 1, 1, 1]
[1, 1, 1, 1, 1]

[1, 1, 1, 1, 1]
[2, 1, 1, 1]
[2, 2, 1]
[3, 1, 1]
[3, 2]
[4, 1]
[5]

Reference: David Eppstein, IntegerPartitions.py

Partition function using dynamic programming in C

A partition of an integer is a way of writing it as a sum of positive integers. For example, partitions of 5 are:

5
4 + 1
3 + 2
3 + 1 + 1
2 + 2 + 1
2 + 1 + 1 + 1
1 + 1 + 1 + 1 + 1

The order of numbers in a partition doesn’t matter. The number of partitions of a number \(n\) is given by the partition function \(p(n)\).

The first few values of \(p(n)\) are 1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627, 792, 1002, 1255, 1575, 1958, 2436, 3010, 3718, 4565, 5604, … (A000041).

The partition function can be computed using a dynamic programming table-filling algorithm. Here it is in C:

#include <stdlib.h>

unsigned int partition(unsigned int n)
{
    unsigned int i, j;
    unsigned int result;
    unsigned int **table;

    if (n == 0) {
        return 1;
    }
   
    table = malloc((n + 1) * sizeof(unsigned int *));
    for (i = 0; i <= n; i++) {
        table[i] = malloc((n + 1) * sizeof(unsigned int));
    }

    for (i = 0;i <= n; i++) {
        table[i][0]=0;
    }
    for (i = 1; i <= n; i++) {
        table[0][i] = 1;
    }

    for (i = 1; i <= n; i++) {
        for (j = 1; j <= n; j++) {
            if (j > i) {
                table[i][j] = table[i][j - 1];
            }
            else {
                table[i][j] = table[i][j - 1] + table[i - j][j];
            }
        }
    }

    result = table[n][n];
    for (i = 0; i <= n; i++) {
        free(table[i]);
    }
    free(table);

    return result;
}

Example program:

#include <stdio.h>

int main(void) 
{
    unsigned int i;
    for (i = 0; i < 31; i++) {
        printf("%d\n",partition(i));
    }

    return 0;
}

Output:

1
1
2
3
5
7
11
15
22
30
42
56
77
101
135
176
231
297
385
490
627
792
1002
1255
1575
1958
2436
3010
3718
4565
5604

Greedy Set Cover in Python

In the Set Cover problem, we are given a universal set \(U\) and a family of subsets \(S_1, \ldots, S_k \subseteq U\). A set cover is a collection of subsets \(C\) from \(S_1, \ldots, S_k\) whose union is the universal set \(U\). We would like to minimise \(\left\vert{C}\right\vert\).

The problem of finding the optimum \(C\) is NP-Complete, but a greedy algorithm can give an \(O(log_e n)\) approximation to optimal solution.

The greedy algorithm selects the set \(S_i\) containing the largest number of uncovered points at each step, until all of the points have been covered.

Below is an implementation in Python:

def set_cover(universe, subsets):
    """Find a family of subsets that covers the universal set"""
    elements = set(e for s in subsets for e in s)
    # Check the subsets cover the universe
    if elements != universe:
        return None
    covered = set()
    cover = []
    # Greedily add the subsets with the most uncovered points
    while covered != elements:
        subset = max(subsets, key=lambda s: len(s - covered))
        cover.append(subset)
        covered |= subset

    return cover

def main():
    universe = set(range(1, 11))
    subsets = [set([1, 2, 3, 8, 9, 10]),
        set([1, 2, 3, 4, 5]),
        set([4, 5, 7]),
        set([5, 6, 7]),
        set([6, 7, 8, 9, 10])]
    cover = set_cover(universe, subsets)
    print(cover)

if __name__ == '__main__':
    main()

Output:

[set([1, 2, 3, 8, 9, 10]), set([4, 5, 7]), set([5, 6, 7])]

Longest Increasing Subsequence using dynamic programming in C

A longest increasing subsequence (LIS) of a sequence is a maximal subsequence which is not necessarily contiguous, and is monotonically increasing across its length, and is as long as the longest such subsequence for that sequence. A given sequence may have more than one longest increasing subsequence.

This is a dynamic programming algorithm for finding an LIS of a sequence of integers. It works by building chains of increasing subsequences as it scans the main sequence, and keeping track of a maximum one. At the end of the scan the maximum subsequence chain is read back and copied to an output array.

struct lis_entry {
    int value;
    int score;
    struct lis_entry *prev;
};
typedef struct lis_entry lis_entry;

unsigned int longest_increasing_subsequence(const unsigned int *arr, size_t len,
        unsigned int **result)
{
    lis_entry *lis;
    unsigned int i, j, max = 0;
    lis_entry *head;

    lis = malloc(len * sizeof(lis_entry));
    if (lis == NULL) {
        return 0;
    }

    /* Initialise entries */
    for (i = 0; i < len; i++ ) {
        lis[i].value = arr[i];
        lis[i].score = 1;
        lis[i].prev = NULL;
    }
 
    /* Calculate LIS scores and build chains */
    for (i = 1; i < len; i++ ) {
        for (j = 0; j < i; j++ ) {
            if (arr[i] > arr[j] && lis[i].score < lis[j].score + 1) {
                lis[i].score = lis[j].score + 1;
                lis[i].prev = &lis[j];
                if (lis[i].score > max) {
                    max = lis[i].score;
                    head = &lis[i];
                }
            }
        }
    }

    /* Allocate the output array and copy the max chain */
    *result = malloc(max * sizeof(int));
    if (*result != NULL) {
        for (i = max - 1; head != NULL; i--) {
            (*(result))[i] = head->value;
            head = head->prev;
        }
    }
    else {
        max = 0;
    }
 
    free(lis);
 
    return max;
}

An example program using the first few terms of the Van der Corput Sequence:

int main(void)
{
    unsigned int arr[] = {0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15};
    const size_t len = sizeof(arr)/sizeof(arr[0]);
    unsigned int *result;
    unsigned int i;
    unsigned int lis_length = longest_increasing_subsequence(arr, len, &result);
    printf("Length of LIS is %d\n", lis_length);
    for (i = 0; i < lis_length; i++) {
        printf("%u\n", result[i]);
    }
    free(result);
    return 0;
}

Output:

Length of LIS is 6
0
4
6
9
13
15

Subset-Sum with backtracking in C

The subset-sum problem is to find a subset of a set of integers that sums to a given value. The decision problem of finding out if such a subset exists is NP-complete. One way of solving the problem is to use backtracking. This is a backtracking solution in C that finds all of the subsets that sum to the target value.

typedef void(*subset_sumfn)(const unsigned int *, size_t);

static unsigned int promising(int i, size_t len, unsigned int weight, unsigned int total, 
        unsigned int target, const unsigned int *weights)
{
    return (weight + total >= target) && (weight == target || weight + weights[i + 1] <= target);
}

static unsigned int sum(const unsigned int *weights, size_t len)
{
    unsigned int total = 0;
    unsigned int i;
    for (i = 0; i < len; i++) {
        total += weights[i];
    }
    return total;
}

static void subset_sum_recursive(const unsigned int *weights, size_t len, unsigned int target, 
        int i, unsigned int weight, unsigned int total, unsigned int *include, subset_sumfn fun)
{
    if (promising(i, len, weight, total, target, weights)) {
        if (weight == target) {
            fun(include, i + 1);
        }
        else if (i < (int)len - 1){
            include[i + 1] = 1;
            subset_sum_recursive(weights, len, target, i + 1, weight + weights[i + 1], 
                   total - weights[i + 1], include, fun);
            include[i + 1] = 0;
            subset_sum_recursive(weights, len, target, i + 1, weight,
                    total - weights[i + 1], include, fun);
        }
    }
}

void subset_sum(const unsigned int *weights, size_t len, unsigned int target, subset_sumfn fun)
{
    const unsigned int total = sum(weights, len);
    unsigned int *include = calloc(len, sizeof(unsigned int));
    if (include == NULL) {
        return;
    }
    subset_sum_recursive(weights, len, target, -1, 0, total, include, fun);
    free(include);
}

int main(void)
{
    unsigned int weights[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
    const unsigned int len = sizeof(weights) / sizeof(unsigned int);
    const unsigned int target = 7;
    subset_sum(weights, len, target, print_vector);
    return 0;
}

The output is in the form of binary strings showing which elements of the original set are in the subset, so for example the first binary string corresponds to 1 + 2 + 4 = 7.

1 1 0 1
1 0 0 0 0 1
0 1 0 0 1
0 0 1 1
0 0 0 0 0 0 1