Tag Archives: Mathematics

Latin squares with backtracking in C

A Latin square is an n × n array filled with n different kinds of object, in which each row and column contains each kind of object only once. For example, below is a 5 × 5 (order 5) Latin square of the integers from 0 to 4:

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

 

Latin squares in which the first row and column are in their natural order as above are called reduced Latin squares. The number of reduced Latin squares of each order goes 1, 1, 1, 4, 56, 9408, 16942080, … (A000315).

Because of the constraints on the structure of Latin squares, there is no simple algorithm for generating them, but backtracking is an effective way of doing so. Below is a backtracking algorithm in C for generating Latin squares, and a program to generate all squares of order 5:

typedef void (*latin_squarefn)(const unsigned int*, unsigned int);

/* Get the cell index at a specified row and column in the Latin square */
unsigned int get_cell(int n, int row, int column)
{
    return (row - 1) * (n - 1) + column - 1;
}

/* Get the row for an array cell */
unsigned int get_row(unsigned int cell, unsigned int n)
{
    return cell / (n - 1) + 1;
}

/* Get the column for an array cell */
unsigned int get_column(unsigned int cell, unsigned int n)
{
    return cell % (n - 1) + 1;
}

/* Check that there are no conflicts */
static unsigned int promising(unsigned int i, const unsigned int *square, unsigned int n)
{
    unsigned int row, column;
    unsigned int ok = 1;
    unsigned int cell;

    /* Check row */
    row = get_row(i, n);
    ok = square[i] != row; /* Row header */
    if (!ok) {
        return 0;
    }
    column = 1;
    cell = get_cell(n, row, column);
    while (cell < i && ok) {
        ok = square[cell] != square[i];
        column++;
        cell = get_cell(n, row, column);
    }
    if (!ok) {
        return 0;
    }
    /* Check column */
    column = get_column(i, n);
    ok = square[i] != column; /* Column header */
    if (!ok) {
        return 0;
    }
    row = 1;
    cell = get_cell(n, row, column);
    while (cell < i && ok) {
        ok = square[cell] != square[i];
        row++;
        cell = get_cell(n, row, column);
    }
    return ok;
}

static void latin_squares_recursive(int i, unsigned int *square, unsigned int n,
        latin_squarefn fun)
{
    const size_t size = (n - 1) * (n - 1);
    if (i < 0 || promising(i, square, n)) {
        if (i == size - 1) {
            fun(square, n);
        }
        else {
            unsigned int val;
            for (val = 0; val < n; val++) {
                square[i + 1] = val;
                latin_squares_recursive(i + 1, square, n, fun);
            }
        }
    }
}

void latin_squares(unsigned int n, latin_squarefn fun)
{
    if (n <= 1) {
        return;
    }
    unsigned int *square = calloc((n - 1) * (n - 1), sizeof(unsigned int));
    if (square == NULL) {
        return;
    }
    latin_squares_recursive(-1, square, n, fun);
    free(square);
}

static void print(const unsigned int *square, unsigned int n)
{
    unsigned int row, col, cell = 0;
    for (col = 0; col < n; col++) {
        printf("%d ", col);
    }
    putchar('\n');
    for (row = 1; row < n; row++) {
        printf("%d ", row);
        for (col = 1; col < n; col++) {
            printf("%d ", square[cell]);
            cell++;
        }
        putchar('\n');
    }
    putchar('\n');
}

int main(void)
{
    latin_squares(5, print);
    return 0;
}

xkcd 287

General solutions get you a 50% tip

I’m not sure what this problem is called – I’m going to call it "multicombination sum" – but I don’t doubt that it is NP-complete, as it’s a variety of knapsack problem in which the values of the items are the same as their weight.

Below are three methods of solving it: a brute force method, using backtracking, and using dynamic programming.

Brute force method

The brute force method is just to construct all of the possible orders that might total $15.05.
The combinatorial algorithm we want is combinations with duplicates, or multicombinations.
Since 3 of the most expensive appetizer – SAMPLER PLATE – exceeds the target, and 7 of the cheapest appetizer – MIXED FRUIT – equals the target (so that’s one of the solutions), we want to iterate over all multicombinations with k ranging from 3 to 7.

unsigned int multiset_sum(const unsigned int *multiset, const unsigned int *values, unsigned int k)
{
    unsigned int i;
    unsigned int sum = 0;
    for (i = 0; i < k; i++) {
        sum += values[multiset[i]];
    }
    return sum;
}

typedef void (*multiset1fn)(const unsigned int *, unsigned int);

void multicombination_sum(const unsigned int *items, size_t len, unsigned int target,
        multiset1fn fun)
{
    unsigned int first = target / array_max(items, len);
    unsigned int last = target / array_min(items, len);
    unsigned int *multiset = calloc(last + 1, sizeof(unsigned int));
    unsigned int k;
    if (!multiset) {
        return;
    }
    for (k = first; k <= last; k++) {
        do {
            if (multiset_sum(multiset, items, k) == target) {
                fun(multiset, k);
            }
        } while (next_multicombination(multiset, len, k));
    }
    free(multiset);
}

void order_print(const unsigned int *numbers, unsigned int k)
{
	const char *appetizers[] = {"MIXED FRUIT", "FRENCH FRIES", "SIDE SALAD", 
        "HOT WINGS", "MOZARELLA STICKS", "SAMPLER PLATE"};
    unsigned int i, item, count;
    for (i = 0; i < k; i++) {
        if (i == 0 || numbers[i] != item) {
            if (i > 0) {
                printf("%s (x%d) ", appetizers[item], count);
            }
            count = 1;
            item = numbers[i];
        }
        else {
            count++;
        }
    }
    printf("%s (x%d)\n", appetizers[item], count);
}

int main(void)
{
    unsigned int prices[] = {215, 275, 335, 355, 420, 580};
	const size_t n = sizeof(prices) / sizeof(prices[0]);
    const unsigned int target = 1505;
    multicombination_sum(prices, n, target, (multiset1fn)order_print);

	return 0;
}

Output:

MIXED FRUIT (x1) HOT WINGS (x2) SAMPLER PLATE (x1)
MIXED FRUIT (x7)

This took 1709 iterations to come up with the answer.

Backtracking

We can drastically reduce the search space by building the orders up one item at a time, and backtracking if the target is exceeded.


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

static void multicombination_sum_recursive(int i, unsigned int *combination, 
       const unsigned int *items, size_t len, unsigned int target, multiset2fn fun,
       unsigned int sum)
{
    if (i == (int)len - 1) {
        if (sum == target) {
            fun(combination, items, i);
        }
    }
    else  {
        unsigned int quantity;
        unsigned int max_quantity = (target - sum) / items[i + 1];
        for (quantity = 0; quantity  <= max_quantity; quantity++) {
            combination[i + 1] = quantity;
            multicombination_sum_recursive(i + 1, combination, items, len, target, fun, 
                    sum + quantity * items[i + 1]);
        }
    }
}

void multicombination_sum(const unsigned int *items, size_t len, unsigned int target,
        multiset2fn fun)
{
    unsigned int *combination = malloc(len * sizeof(unsigned int));
    multicombination_sum_recursive(-1, combination, items, len, target, fun, 0);
    free(combination);
}

void order_print(const unsigned int *combination, const unsigned int *items, size_t len)
{
    const char *appetizers[] = {"MIXED FRUIT", "FRENCH FRIES", "SIDE SALAD", 
        "HOT WINGS", "MOZARELLA STICKS", "SAMPLER PLATE"};
    unsigned int i;
    for (i = 0; i <= len; i++) {
        if (combination[i] > 0) {
            printf("%s (x%d) ", appetizers[i], combination[i]);
        }
    }
    putchar('\n');
}

int main(void)
{
    unsigned int prices[] = {215, 275, 335, 355, 420, 580};
    const unsigned int len = sizeof(prices) / sizeof(unsigned int);
    const unsigned int target = 1505;
    multicombination_sum(prices, len, target, (multiset2fn)order_print);
    return 0;
}

Output:

MIXED FRUIT (x1) HOT WINGS (x2) SAMPLER PLATE (x1)
MIXED FRUIT (x7)

This took just 211 iterations.

Dynamic Programming

As I said at the beginning, this problem is a special case of the knapsack problem in which the values of the items are the same as their weight. This means that we can use the classic dynamic programming algorithm for knapsack to solve it. The algorithm works by calculating the most profitable knapsack for each capacity up to the target capacity.

Here is an implementation of the algorithm:

struct knapsack {
    unsigned int profit;
    struct knapsack *prev;
};
typedef struct knapsack knapsack;

/* Find the minimum weight item with this profit */
int min_weight_item(unsigned int profit, const unsigned int *weights, const unsigned int *profits,
        size_t len)
{
    int item = -1;
    unsigned int i;
    for (i = 0; i < len; i++) {
        if (profits[i] == profit) {
            if (item == -1 || weights[i] < weights[item]) {
                item = i;
            }
        }
    }
    return item;
}

unsigned int unbounded_knapsack(unsigned int capacity, unsigned int *weights, 
       unsigned int *profits, unsigned int *counts, size_t len)
{
    knapsack *z = malloc((capacity + 1) * sizeof(knapsack));
    unsigned int c, i;
    unsigned int solution, profit;
    z[0].profit = 0;
    z[0].prev = NULL;
    knapsack *current;
    /* Fill in the array */
    for (c = 1; c <= capacity; c++) {
        z.profit = z.profit;
        z.prev = &(z);
        for (i = 0; i < len; i++) {
            if (weights[i] <= c) {
                knapsack *prev = z + (c - weights[i]);
                if (prev->profit + profits[i] > z.profit) {
                    z.profit = prev->profit + profits[i];
                    z.prev = prev;
                }
            }
        }
    }
    /* Read back the best solution */
    for (profit = z[capacity].profit, current = z[capacity].prev;
            current != NULL;
            profit = current->profit, current = current->prev) {
        counts[min_weight_item(profit - current->profit, weights, profits, len)]++;
        
    }
    solution = z[capacity].profit;
    free(z);
    return solution;
}

We just need to use this algorithm and pass the prices of the menu items for both weights and profits:

int main(void)
{
    unsigned int values[] = {215, 275, 335, 355, 420, 580};
    const size_t len = sizeof(values) / sizeof(values[0]);
    unsigned int counts[6] = {0};
    const unsigned int target = 1505;
    unbounded_knapsack(target, values, values, counts, len);
    order_print(counts, len);
    return 0;
}

This produces one of the solutions:

MIXED FRUIT (x7)

We could modify the algorithm to produce all of them.

Functors in C++

A functor in C++ is something that can be called like a function – a function pointer or object of a class with an overloaded operator(). It is the second type that is most commonly being referred to when people talk about functors.

The advantages of having an object of a class with an overloaded operator() over a function pointer are:

  1. The object can maintain state between calls
  2. Calling the overloaded operator() can be more efficient than dereferencing a function pointer

Maintaining state

A class with an overloaded operator() can have member variables, and these allow it to maintain state through function calls. For example, a functor to calculate the Fibonacci numbers could maintain a cache of the ones seen so far to avoid recalculating them.

#include <vector>
#include <iostream>

class Fibonacci
{
public:
    Fibonacci()
    {
        results_.push_back(0);
        results_.push_back(1);
    }
    unsigned int operator()(unsigned int n)
    {
        while (results_.size() < n + 1) {
            results_.push_back(results_[results_.size() - 2] + results_[results_.size() - 1]);
        }
        return results_[n];
    }
private:
    std::vector<unsigned int> results_;
};

int main()
{
    Fibonacci fib;
    for (unsigned int n = 0; n < 10; n++) {
        std::cout << fib(n) << "\n";
    }
}

Efficiency

Functors are commonly used with standard library algorithms like std::for_each() and std::transform(), as in the example below, which fills a vector with the first 100 Fibonaci numbers:

#include <vector>
#include <iostream>
#include <algorithm>
#include <iterator>

int main()
{
    Fibonacci fib;
    std::vector<unsigned int> results(100);
    std::iota(results.begin(), results.end(), 0);
    std::transform(results.begin(), results.end(), results.begin(), fib);
    std::copy(results.begin(), results.end(),
         std::ostream_iterator<unsigned int>(std::cout, " "));
}

When an ordinary function pointer is passed to one of these algorithms, it’s difficult for the compiler to optimise the function call because a function pointer is a variable and so could change to point to a different function. When the pointer is to a virtual member function there is an extra level of indirection to get to the version in the actual class. With an overloaded operator() however, there is no indirection and the function to be called is fixed, so the compiler can optimise the call by inlining the function’s code at the call site. This can make the call much more efficient.

Graph algorithms

Introduction

This is a list of graph algorithms with links to references and implementations.
The graph libraries included are igraph, NetworkX, and Boost Graph Library.

Contents

Generators

Deterministic

Star

Wikipedia: https://en.wikipedia.org/wiki/Star_%28graph_theory%29

Complete

Wikipedia: https://en.wikipedia.org/wiki/Complete_graph

Note that the igraph function can produce directed complete graphs, as well as those containing 1 or more self-loops.

Tree

Famous

Bull

Wikipedia: https://en.wikipedia.org/wiki/Bull_graph

Chvátal

Wikipedia: https://en.wikipedia.org/wiki/Chv%C3%A1tal_graph

Cubic

Wikipedia: https://en.wikipedia.org/wiki/Cubic_graph

Diamond

Wikipedia: https://en.wikipedia.org/wiki/Diamond_graph

Dodecahedral

Wolfram MathWorld: http://mathworld.wolfram.com/DodecahedralGraph.html

Frucht

Wikipedia: https://en.wikipedia.org/wiki/Frucht_graph

Heawood

Wikipedia: https://en.wikipedia.org/wiki/Heawood_graph

House

Wolfram MathWorld: http://mathworld.wolfram.com/HouseGraph.html

Icosahedral

Wolfram MathWorld: http://mathworld.wolfram.com/IcosahedralGraph.html

Krackhardt Kite

Wolfram MathWorld: http://mathworld.wolfram.com/KrackhardtKite.html

Octahedral

Wolfram MathWorld: http://mathworld.wolfram.com/OctahedralGraph.html

Petersen

Wolfram MathWorld: http://mathworld.wolfram.com/PetersenGraph.html

Tetrahedral

Wolfram MathWorld: http://mathworld.wolfram.com/TetrahedralGraph.html

Tutte

Wikipedia: https://en.wikipedia.org/wiki/Tutte_graph

Random

Erdos-Rényi

Wikipedia: https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model

Watts-Strogatz

Wikipedia: https://en.wikipedia.org/wiki/Watts_and_Strogatz_model

Barabási–Albert

Wikipedia: https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model

Traversal

Lecture notes: http://web.cse.ohio-state.edu/~gurari/course/cis680/cis680Ch14.html

Depth-first search

Wikipedia: https://en.wikipedia.org/wiki/Depth-first_search

Breadth-first search

Wikipedia: https://en.wikipedia.org/wiki/Breadth-first_search

Isomorphism

Wikipedia: https://en.wikipedia.org/wiki/Graph_isomorphism

Connected components

Wikipedia: https://en.wikipedia.org/wiki/Connectivity_%28graph_theory%29

Connectivity

Connectivity testing

Decomposition into connected components

Strong connectivity

Strong connectivity testing

Decomposition into strongly connected components

Weak connectivity

Weak connectivity testing

Decomposition into weakly connected components

Biconnectivity

Biconnected components

Articulation points

Cliques

Wikipedia: https://en.wikipedia.org/wiki/Clique_%28graph_theory%29

Cliques

Maximal cliques

Clique number

Number of maximal cliques

Flows and Cuts

Lecture notes: http://www.cs.princeton.edu/courses/archive/spr04/cos226/lectures/maxflow.4up.pdf

Maximum flow

Maximum flow value

Minimum cut

Minimum cut value

Shortest paths

Wikipedia: https://en.wikipedia.org/wiki/Shortest_path_problem

Shortest paths

Dijkstra

Bellman-Ford

Johnson

A*

Distance measures

Wikipedia: https://en.wikipedia.org/wiki/Distance_%28graph_theory%29

Diameter

Eccentricity

Radius

Centrality measures

Lecture notes: http://cs.brynmawr.edu/Courses/cs380/spring2013/section02/slides/05_Centrality.pdf

Closeness

Betweenness

Google PageRank

Paper: http://ilpubs.stanford.edu:8090/422/1/1999-66.pdf

Directed Acyclic Graphs (DAGs)

Wikipedia: https://en.wikipedia.org/wiki/Directed_acyclic_graph

DAG testing

Topological sort

Transitive closure

Bipartite graphs

Wolfram Mathworld: http://mathworld.wolfram.com/BipartiteGraph.html

Creation

Testing

Creation from adjacency matrix

Note that the igraph documentation calls this an incidence matrix.

Conversion to adjacency matrix

Note that the igraph documentation calls this an incidence matrix.

Projection

Cores

Paper: http://arxiv.org/abs/cs.DS/0310049

Core number

Chordal graphs

Wolfram Mathworld: http://mathworld.wolfram.com/ChordalGraph.html

Chordal testing

Maximum cardinality search

Spanning trees

Wikipedia: https://en.wikipedia.org/wiki/Spanning_tree

Note that in NetworkX they are called arborescences (http://networkx.github.io/documentation/networkx-1.10/reference/algorithms.tree.html#recognition-tests).

Minimum spanning tree

Clustering

Article: http://dollar.biz.uiowa.edu/~street/graphClustering.pdf

Transitivity

Clustering coefficient (local transitivity)

Average clustering (average local transitivity)

Degree sequences

http://mathworld.wolfram.com/DegreeSequence.html

Simple graph testing

Pseudograph testing

Spectral properties

Book chapter: http://www.cs.yale.edu/homes/spielman/PAPERS/SGTChapter.pdf

Laplacian matrix

Reading and writing graphs

Edge list

GraphML

The GraphML File Format: http://graphml.graphdrawing.org/

GML

GML: a portable Graph File Format: http://www.fim.uni-passau.de/fileadmin/files/lehrstuhl/brandenburg/projekte/gml/gml-technical-report.pdf

Pajek

Pajek / Pajek-XXL: http://mrvar.fdv.uni-lj.si/pajek/

GraphViz

GraphViz: http://www.graphviz.org/

K-Permutations

The k-permutations of a set are the permutations of the combinations of size k. They are also known as sequences without repetitions.

There are 24 3-permutations of the 4-set {0, 1, 2, 3}:

[0, 1, 2]
[0, 2, 1]
[1, 0, 2]
[1, 2, 0]
[2, 0, 1]
[2, 1, 0]
[0, 1, 3]
[0, 3, 1]
[1, 0, 3]
[1, 3, 0]
[3, 0, 1]
[3, 1, 0]
[0, 2, 3]
[0, 3, 2]
[2, 0, 3]
[2, 3, 0]
[3, 0, 2]
[3, 2, 0]
[1, 2, 3]
[1, 3, 2]
[2, 1, 3]
[2, 3, 1]
[3, 1, 2]
[3, 2, 1]

The algorithm simply finds the next permutation of the array. If the array is at the last permutation, the next combination from the set is constructed.

Note that this algorithm does not produce the k-permutations in lexicographic order.

unsigned int next_k_permutation(unsigned int *ar, size_t n, unsigned int k)
{
	unsigned int result = next_permutation(ar, k);
	if (result == 0) {
		result = next_combination(ar, n, k);
	}
	return result;
}

Subsets of a multiset

This is the set of all subsets of a multiset, which are themselves multisets, i.e., the power set of a multiset.

For example, the power set of the multiset [a, b, b, c] consists of the sets:

()
(c)
(b)
(b, c)
(b, b)
(b, b, c)
(a)
(a, c)
(a, b)
(a, b, c)
(a, b, b)
(a, b, b, c)

The set of all subsets of a particular size are combinations of a multiset.

The multiset and its subsets are represented as a vector containing, for each element, the count of its occurrences in the multiset. For example, the multiset [a, b, b, c] is represented as [1, 2, 1]. This is similar to the characteristic vector used for power set, but with counts rather than boolean values.

The correspondences for the subsets are then as follows:

[0, 0, 0] = ()
[0, 0, 1] = (c)
[0, 1, 0] = (b)
[0, 1, 1] = (b, c)
[0, 2, 0] = (b, b)
[0, 2, 1] = (b, b, c)
[1, 0, 0] = (a)
[1, 0, 1] = (a, c)
[1, 1, 0] = (a, b)
[1, 1, 1] = (a, b, c)
[1, 2, 0] = (a, b, b)
[1, 2, 1] = (a, b, b, c)

The algorithm can then simply count from [0, 0, 0] to [1, 2, 1], using the values in the source multiset as the upper limit for the value of an element.

Note that this algorithm does not produce the subsets in lexicographic order.

unsigned int next_multiset_subset(const unsigned int *multiset, unsigned int *ar, size_t n)
{
	unsigned int changed = 0;
	int i;

	for (i = n - 1; i >= 0 && !changed; i--) {
		if (ar[i] < multiset[i]) {
			/* Increment */
			ar[i]++;
			changed = 1;
		}
		else {
			/* Roll over */
			ar[i] = 0;
		}
	}
	if (!changed) {
		/* Reset to first combination */
		for (i = 0; i < n; i++) {
			ar[i] = 0;
		}
	}
	return changed;
}

Here is an example program:

#include <stdio.h>

#include <multiset-subset.h>

static void print_list(const unsigned int *ar, size_t len, FILE *fptr)
{
	unsigned int i;
	fputc('(', fptr);
	for (i = 0; i < len; i++) {
		fprintf(fptr, "%d", ar[i]);
		if (i < len - 1) {
			fputs(", ", fptr);
		}
	}
	fputc(')', fptr);
}

int main(void)
{
	unsigned int multiset[] = {1, 2, 1};
	const size_t n = sizeof(multiset) / sizeof(unsigned int);
	unsigned int numbers[] = {0, 0, 0};

	do {
		print_list(numbers, n, stdout);
		putchar('\n');
	} while (next_multiset_subset(multiset, numbers, n));

	return 0;
}

Here is a second example that prints the elements of the multisets:

#include <stdio.h>

#include <multiset-subset.h>

static void print_array(const unsigned int *ar, size_t len, FILE *fptr)
{
	unsigned int i;
	fputc('[', fptr);
	for (i = 0; i < len; i++) {
		fprintf(fptr, "%d", ar[i]);
		if (i < len - 1) {
			fputs(", ", fptr);
		}
	}
	fputc(']', fptr);
}


static void print_multiset_subset(const unsigned int *ar, size_t len, 
		const void **elements, printfn print, FILE *fptr)
{
	unsigned int i, started = 0;
	fputc('(', fptr);
	for (i = 0; i < len; i++) {
		unsigned int j;
		for (j = 0; j < ar[i]; j++) {
			if (started) {
				fputs(", ", fptr);
			}
			print(elements[i], fptr);
			started = 1;
		}
	}
	fputc(')', fptr); 
}


int main(void)
{
	unsigned int multiset[] = {1, 2, 1};
	const size_t n = sizeof(multiset) / sizeof(unsigned int);
	char *elements[] = {"a", "b", "c"};
	unsigned int numbers[] = {0, 0, 0};

	do {
        print_array(numbers, n, stdout);
        printf(" = "); 
		print_multiset_subset(numbers, n, (void*)elements, (printfn)fputs, stdout);
		putchar('\n');
	} while (next_multiset_subset(multiset, numbers, n));

	return 0;
}

Multicombinations

A multicombination is a combination that can contain duplicates (i.e., the combination is a multiset).

Note that the source set does not contain duplicates. For combinations of a set that contains duplicates, see combinations of a multiset.

These are the 3-multicombinations of the 4-set {0, 1, 2, 3}:

[0, 0, 0]
[0, 0, 1]
[0, 0, 2]
[0, 0, 3]
[0, 1, 1]
[0, 1, 2]
[0, 1, 3]
[0, 2, 2]
[0, 2, 3]
[0, 3, 3]
[1, 1, 1]
[1, 1, 2]
[1, 1, 3]
[1, 2, 2]
[1, 2, 3]
[1, 3, 3]
[2, 2, 2]
[2, 2, 3]
[2, 3, 3]
[3, 3, 3]

This algorithm produces the combinations in lexicographic order, with the elements in each combination in increasing order.

To find the next multicombination containing k elements from a set containing n elements, begin with the multicombination containing k zeroes, then at each step:

  1. Find the rightmost element that is less than n – 1
  2. Increment it.
  3. Make the elements after it the same.
unsigned int next_multicombination(unsigned int *ar, size_t n, unsigned int k)
{
	unsigned int changed = 0;
	int i;

	for (i = k - 1; i >= 0 && !changed; i--) {
		if (ar[i] < n - 1) {
			/* Increment this element */
			ar[i]++;
			if (i < k - 1) {
				/* Make the elements after it the same */
				unsigned int j;
				for (j = i + 1; j < k; j++) {
					ar[j] = ar[j - 1];
				}
			}
			changed = 1;
		}
	}
	if (!changed) {
		/* Reset to first combination */
		for (i = 0; i < k; i++) {
			ar[i] = 0;
		}
	}
	return changed;
}

Example program to produce the output above:

void print_array(const unsigned int *ar, size_t len, FILE *fptr)
{
	unsigned int i;
	fputc('[', fptr);
	for (i = 0; i < len; i++) {
		fprintf(fptr, "%d", ar[i]);
		if (i < len - 1) {
			fputs(", ", fptr);
		}
	}
	fputc(']', fptr);
}

int main(void)
{
	unsigned int numbers[3] = {0};
	const size_t k = sizeof(numbers) / sizeof(numbers[0]);
	unsigned int n = 4;

	do {
		print_array(numbers, k, stdout);
		putchar('\n');
	} while (next_multicombination(numbers, n, k));

	return 0;
}

Combinations

Combinations, or k-combinations, are the unordered sets of k elements chosen from a set of size n.

For example, there are 10 3-combinations of the 5-set {0, 1, 2, 3, 4}:

{0, 1, 2}
{0, 1, 3}
{0, 1, 4}
{0, 2, 3}
{0, 2, 4}
{0, 3, 4}
{1, 2, 3}
{1, 2, 4}
{1, 3, 4}
{2, 3, 4}

The set of all combinations 0 ≤ k ≤ n is the power set.

Combinations that can contain duplicates are called multicombinations.

Combinations of a set containing duplicates are combinations of a multiset.

Combinations where order is significant are called k-permutations.

This algorithm produces the combinations in lexicographic order, with the elements in each combination strictly increasing in magnitude.

Begin with the combination containing the the numbers from 0 to k – 1, and at each step:

  1. Find the rightmost element ar[i] that is less than the maximum value it can have (which is (n – 1) – (k – 1) – i)
  2. Increment it
  3. Turn the elements after it into a linear sequence continuing from ar[i]
unsigned int next_combination(unsigned int *ar, size_t n, unsigned int k)
{
	unsigned int finished = 0;
	unsigned int changed = 0;
	unsigned int i;

	if (k > 0) {
		for (i = k - 1; !finished && !changed; i--) {
			if (ar[i] < (n - 1) - (k - 1) + i) {
				/* Increment this element */
				ar[i]++;
				if (i < k - 1) {
					/* Turn the elements after it into a linear sequence */
					unsigned int j;
					for (j = i + 1; j < k; j++) {
						ar[j] = ar[j - 1] + 1;
					}
				}
				changed = 1;
			}
			finished = i == 0;
		}
		if (!changed) {
			/* Reset to first combination */
			for (i = 0; i < k; i++) {
				ar[i] = i;
			}
		}
	}
	return changed;
}

Fibonacci numbers in C

The recurrence formula for the \(n\)th Fibonacci number is:

\(
F_{n} = \begin{cases} F_{n – 1} + F_{n – 2} & \textrm{if}\;n > 1\\ 1 &\textrm{if}\;n = 1\\ 0 &\textrm{if}\;n = 0.\end{cases}
\)

A naive recursive implementation of this would be:

unsigned int fib(unsigned int n)
{
	if (n == 0) return 0;
	if (n == 1) return 1;
	return fib(n - 1) + fib(n - 2);
}

But this runs in exponential time because it keeps recalculating the same results in the recursion.

If we add some instrumentation, we can see what it’s doing:

struct calltree {
    unsigned int value;
    struct calltree *left;
    struct calltree *right;
};
typedef struct calltree calltree;

typedef void (*calltreefn)(unsigned int, unsigned int, void *);

calltree *calltree_create(unsigned int value)
{
    calltree *tree = malloc(sizeof(calltree));
    if (tree) {
        tree->value = value;
        tree->left = NULL;
        tree->right = NULL;
    }
    return tree;
}

void calltree_delete(calltree *tree)
{
    if (tree != NULL) {
        calltree_delete(tree->left);
        calltree_delete(tree->right);
        free(tree);
    }
}

void calltree_for_each(const calltree *tree, unsigned int level, calltreefn fun, void *data)
{
    if (tree != NULL) {
        fun(level, tree->value, data);
        calltree_for_each(tree->left, level + 1, fun, data);
        calltree_for_each(tree->right, level + 1, fun, data);
    }
}

unsigned int fib(unsigned int n, calltree **tree)
{
    *tree = calltree_create(n);
	if (n == 0) return 0;
	if (n == 1) return 1;
	return fib(n - 1, &((*tree)->left)) + fib(n - 2, &((*tree)->right));
}

Here is an example program to calculate \(F_{7}\):

void print_call_text(unsigned int level, unsigned int value, void *data)
{
    unsigned int i;
    for (i = 0; i < level; i++) {
        putchar(' ');
    }
    printf("fib(%d)\n", value);
}

int main(void)
{
    calltree *tree;
    fib(7, &tree);
    calltree_for_each(tree, 0, print_call_text, NULL);
    calltree_delete(tree);
    printf("%u\n", result);
	return 0;
}

The call tree produced:

fib(7)
 fib(6)
  fib(5)
   fib(4)
    fib(3)
     fib(2)
      fib(1)
      fib(0)
     fib(1)
    fib(2)
     fib(1)
     fib(0)
   fib(3)
    fib(2)
     fib(1)
     fib(0)
    fib(1)
  fib(4)
   fib(3)
    fib(2)
     fib(1)
     fib(0)
    fib(1)
   fib(2)
    fib(1)
    fib(0)
 fib(5)
  fib(4)
   fib(3)
    fib(2)
     fib(1)
     fib(0)
    fib(1)
   fib(2)
    fib(1)
    fib(0)
  fib(3)
   fib(2)
    fib(1)
    fib(0)
   fib(1)

We can get rid of all of these recursive calls by caching the results instead:

unsigned int fib(unsigned int n)
{
    unsigned int *numbers;
    unsigned int i;
    unsigned int result;
    if (n == 0) return 0;
    numbers = malloc((n + 1) * sizeof(int));
    if (numbers == NULL) {
        return 0;
    }
    numbers[0] = 0;
    numbers[1] = 1;
    for (i = 2; i <= n; i++) {
        numbers[i] = numbers[i - 2] + numbers[i - 1];
    }
    result = numbers[n];
    free(numbers); 
    return result;
}

If we were going to calculate a lot of Fibonacci numbers we could store the cache between calls rather than creating it anew each time.

Permutations

The permutations of a set are the ways of arranging its elements.

For example, there are 24 permutations of a set of 4 elements:

[0, 1, 2, 3]
[0, 1, 3, 2]
[0, 2, 1, 3]
[0, 2, 3, 1]
[0, 3, 1, 2]
[0, 3, 2, 1]
[1, 0, 2, 3]
[1, 0, 3, 2]
[1, 2, 0, 3]
[1, 2, 3, 0]
[1, 3, 0, 2]
[1, 3, 2, 0]
[2, 0, 1, 3]
[2, 0, 3, 1]
[2, 1, 0, 3]
[2, 1, 3, 0]
[2, 3, 0, 1]
[2, 3, 1, 0]
[3, 0, 1, 2]
[3, 0, 2, 1]
[3, 1, 0, 2]
[3, 1, 2, 0]
[3, 2, 0, 1]
[3, 2, 1, 0]

To find the next permutation of an array ar:

  1. Find the highest index, i1 such that ar[i1] is the first of a pair of elements in ascending order.
    If there isn’t one, the sequence is the highest permutation, so reverse the whole thing to begin again.

  2. Find the highest index i2, such that i2 > i1 and ar[i2] > ar[i1].
  3. Swap ar[i1] and ar[i2].
  4. The elements from ar[i1 + 1] to the end are now in descending order (a later permutation), so reverse them.
static void swap(unsigned int *ar, unsigned int first, unsigned int second)
{
	unsigned int temp = ar[first];
	ar[first] = ar[second];
	ar[second] = temp;
}

static void reverse(unsigned int *ar, size_t len)
{
	unsigned int i, j;

	for (i = 0, j = len - 1; i < j; i++, j--) {
		swap(ar, i, j);
	}
}

unsigned int next_permutation(unsigned int *ar, size_t len)
{
	unsigned int i1, i2;
	unsigned int result = 0;
	
	/* Find the rightmost element that is the first in a pair in ascending order */
	for (i1 = len - 2, i2 = len - 1; ar[i2] <= ar[i1] && i1 != 0; i1--, i2--);
	if (ar[i2] <= ar[i1]) {
		/* If not found, array is highest permutation */
		reverse(ar, len);
	}
	else {
		/* Find the rightmost element to the right of i1 that is greater than ar[i1] */
		for (i2 = len - 1; i2 > i1 && ar[i2] <= ar[i1]; i2--);
		/* Swap it with the first one */
		swap(ar, i1, i2);
		/* Reverse the remainder */
		reverse(ar + i1 + 1, len - i1 - 1);
		result = 1;
	}
	return result;
}