Make code nice

Signed-off-by: Mattias Andrée <maandree@kth.se>

3 files changed, 506 insertions, 599 deletions

diff --git a/COPYING b/COPYING
index edf2d2f..accda30 100644
--- a/COPYING
+++ b/COPYING
@@ -1,6 +1,6 @@
 𝓞(n³) implementation of the Hungarian algorithm
 
-Copyright (C) 2011, 2014  Mattias Andrée
+Copyright (C) 2011, 2014, 2020  Mattias Andrée
 
 This program is free software. It comes without any warranty, to
 the extent permitted by applicable law. You can redistribute it
diff --git a/Makefile b/Makefile
index 58867b4..2023518 100644
--- a/Makefile
+++ b/Makefile
@@ -8,8 +8,7 @@ WARN = -Wall -Wextra -pedantic -Wdouble-promotion -Wformat=2 -Winit-self -Wmissi
        -Wunsuffixed-float-constants -Wsuggest-attribute=const -Wsuggest-attribute=noreturn       \
        -Wsuggest-attribute=pure -Wsuggest-attribute=format -Wnormalized=nfkc -Wconversion        \
        -fstrict-aliasing -fstrict-overflow -fipa-pure-const -ftree-vrp -fstack-usage             \
-        -funsafe-loop-optimizations
-# excluded: -Wdeclaration-after-statement
+       -funsafe-loop-optimizations  -Wdeclaration-after-statement
 
 
 all: hungarian
@@ -17,14 +16,11 @@ all: hungarian
 hungarian: hungarian.c
 	gcc -std=gnu99 $(OPTIMISE) $(WARN) $(CFLAGS) $(LDFLAGS) $(CPPFLAGS) -o $@ $<
 
-test:
-	./"hungarian"
-
-valgrind:
-	valgrind --tool=memcheck --leak-check=full ./"hungarian"
+test: hungarian
+	./hungarian
 
 clean:
-	-rm hungarian
+	-rm -f -- hungarian
 
 
 .PHONY: all test valgrind clean
diff --git a/hungarian.c b/hungarian.c
index fdc7a7a..38c9f55 100644
--- a/hungarian.c
+++ b/hungarian.c
@@ -1,9 +1,7 @@
-// -*- mode: c, coding: utf-8 -*-
-
 /**
  * 𝓞(n³) implementation of the Hungarian algorithm
  * 
- * Copyright (C) 2011, 2014  Mattias Andrée
+ * Copyright (C) 2011, 2014, 2020  Mattias Andrée
  * 
  * This program is free software. It comes without any warranty, to
  * the extent permitted by applicable law. You can redistribute it
@@ -13,273 +11,189 @@
  */
 
 
+#include <stddef.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
-
-
-#ifndef RANDOM_DEVICE
-#define RANDOM_DEVICE "/dev/urandom"
-#endif
-
-
-#define cell      long
-#define CELL_STR  "%li"
-
-#define llong    int_fast64_t
-#define byte     int_fast8_t
-#define boolean  int_fast8_t
-#define null     0
-#define true     1
-#define false    0
+#include <string.h>
 
 
 
-#ifdef DEBUG
-#  define debug(X) fprintf(stderr, "\033[31m%s\033[m\n", X)
-#else
-#  define debug(X) 
-#endif
-
-
+/**
+ * Cell markings
+ **/
+enum {
+	UNMARKED = 0,
+	MARKED,
+	PRIME
+};
 
 
 /**
- * Cell marking:  none
+ *  Value type for marking
  */
-#define UNMARKED  0L
+typedef int_fast8_t Mark;
 
 /**
- * Cell marking:  marked
+ *  Value type for cells
  */
-#define MARKED    1L
+typedef signed long int Cell;
 
-/**
- * Cell marking:  prime
- */
-#define PRIME     2L
+typedef int_fast8_t Boolean;
 
+typedef int_fast64_t BitSetLimb;
 
 
 /**
  * Bit set, a set of fixed number of bits/booleans
  */
-typedef struct
-{
-    /**
-     * The set of all limbs, a limb consist of 64 bits
-     */
-    llong* limbs;
-    
-    /**
-     * Singleton array with the index of the first non-zero limb
-     */
-    size_t* first;
-    
-    /**
-     * Array the the index of the previous non-zero limb for each limb
-     */
-    size_t* prev;
-    
-    /**
-     * Array the the index of the next non-zero limb for each limb
-     */
-    size_t* next;
-    
+typedef struct {
+	/**
+	 * The set of all limbs, a limb consist of 64 bits
+	 */
+	BitSetLimb *limbs;
+
+	/**
+	 * Singleton array with the index of the first non-zero limb
+	 */
+	size_t first;
+
+	/**
+	 * Array the the index of the previous non-zero limb for each limb
+	 */
+	size_t *prev;
+
+	/**
+	 * Array the the index of the next non-zero limb for each limb
+	 */
+	size_t *next;
+
+	char _buf[];
 } BitSet;
 
 
+typedef struct {
+	size_t row;
+	size_t col;
+} CellPosition;
 
-ssize_t** kuhn_match(cell** table, size_t n, size_t m);
-void kuhn_reduceRows(cell** t, size_t n, size_t m);
-byte** kuhn_mark(cell** t, size_t n, size_t m);
-boolean kuhn_isDone(byte** marks, boolean* colCovered, size_t n, size_t m);
-size_t* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, size_t n, size_t m);
-void kuhn_altMarks(byte** marks, size_t* altRow, size_t* altCol, ssize_t* colMarks, ssize_t* rowPrimes, size_t* prime, size_t n, size_t m);
-void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, size_t n, size_t m);
-ssize_t** kuhn_assign(byte** marks, size_t n, size_t m);
 
-BitSet new_BitSet(size_t size);
-void BitSet_set(BitSet this, size_t i);
-void BitSet_unset(BitSet this, size_t i);
-ssize_t BitSet_any(BitSet this) __attribute__((pure));
 
-size_t lb(llong x) __attribute__((const));
+/**
+ * Calculates the floored binary logarithm of a positive integer
+ *
+ * @param   value  The integer whose logarithm to calculate
+ * @return         The floored binary logarithm of the integer
+ */
+#if defined(__GNUC__)
+__attribute__((__const__))
+#endif
+static size_t
+lb(BitSetLimb value)
+{
+	size_t rc = 0;
+	BitSetLimb v = value;
 
+	if (v & (int_fast64_t)0xFFFFFFFF00000000LL) { rc |= 32L; v >>= 32; }
+	if (v & (int_fast64_t)0x00000000FFFF0000LL) { rc |= 16L; v >>= 16; }
+	if (v & (int_fast64_t)0x000000000000FF00LL) { rc |=  8L; v >>=  8; }
+	if (v & (int_fast64_t)0x00000000000000F0LL) { rc |=  4L; v >>=  4; }
+	if (v & (int_fast64_t)0x000000000000000CLL) { rc |=  2L; v >>=  2; }
+	if (v & (int_fast64_t)0x0000000000000002LL) { rc |=  1L; }
 
+	return rc;
+}
 
-void print(cell** t, size_t n, size_t m, ssize_t** assignment);
 
-int main(int argc, char** argv)
+/**
+ * Constructor for BitSet
+ *
+ * @param   size  The (fixed) number of bits to bit set should contain
+ * @return        The a unique BitSet instance with the specified size
+ */
+static BitSet *
+bitset_create(size_t size)
 {
-    FILE* urandom = fopen(RANDOM_DEVICE, "r");
-    unsigned int seed;
-    fread(&seed, sizeof(unsigned int), 1, urandom);
-    srand(seed);
-    fclose(urandom);
-    
-    
-    size_t i, j;
-    size_t n = argc < 3 ? 10 : (size_t)atol(*(argv + 1));
-    size_t m = argc < 3 ? 15 : (size_t)atol(*(argv + 2));
-    cell** t = malloc(n * sizeof(cell*));
-    cell** table = malloc(n * sizeof(cell*));
-    if (argc < 3)
-        for (i = 0; i < n; i++)
-	{
-	    *(t + i) = malloc(m * sizeof(cell));
-	    *(table + i) = malloc(m * sizeof(cell));
-	    for (j = 0; j < m; j++)
-	        *(*(table + i) + j) = *(*(t + i) + j) = (cell)(random() & 63);
-	}
-    else
-    {
-        cell x;
-        for (i = 0; i < n; i++)
-	{
-	    *(t + i) = malloc(m * sizeof(cell));
-	    *(table + i) = malloc(m * sizeof(cell));
-	    for (j = 0; j < m; j++)
-	    {
-	        scanf(CELL_STR, &x);
-	        *(*(table + i) + j) = *(*(t + i) + j) = x;
-	    }
-	}
-    }
-    
-    printf("\nInput:\n\n");
-    print(t, n, m, null);
-    
-    ssize_t** assignment = kuhn_match(table, n, m);
-    printf("\nOutput:\n\n");
-    print(t, n, m, assignment);
-    
-    cell sum = 0;
-    for (i = 0; i < n; i++)
-    {
-        sum += *(*(t + *(*(assignment + i) + 0)) + *(*(assignment + i) + 1));
-	free(*(assignment + i));
-	free(*(table + i));
-	free(*(t + i));
-    }
-    free(assignment);
-    free(table);
-    free(t);
-    printf("\n\nSum: %li\n\n", sum);
-    
-    return 0;
+	size_t c     = (size >> 6) + !!(size & 63L);
+	BitSet *this = calloc(1, offsetof(BitSet, _buf) + c * sizeof(BitSetLimb) + 2 * (c + 1) * sizeof(size_t));
+
+	this->limbs =  (BitSetLimb *)&this->_buf[0];
+	this->prev  = (size_t *)&this->_buf[c * sizeof(BitSetLimb)];
+	this->next  = (size_t *)&this->_buf[c * sizeof(BitSetLimb) + c * sizeof(size_t)];
+
+	return this;
 }
 
-void print(cell** t, size_t n, size_t m, ssize_t** assignment)
+
+/**
+ * Gets the index of any set bit in a bit set
+ * 
+ * @param   this  The bit set
+ * @return        The index of any set bit
+ */
+#if defined(__GNUC__)
+__attribute__((__pure__))
+#endif
+static ssize_t
+bitset_any(BitSet *this)
 {
-    size_t i, j;
-    
-    ssize_t** assigned = malloc(n * sizeof(ssize_t*));
-    for (i = 0; i < n; i++)
-    {
-        *(assigned + i) = malloc(m * sizeof(ssize_t));
-	for (j = 0; j < m; j++)
-	    *(*(assigned + i) + j) = 0;
-    }
-    if (assignment != null)
-        for (i = 0; i < n; i++)
-	    (*(*(assigned + **(assignment + i)) + *(*(assignment + i) + 1)))++;
-    
-    for (i = 0; i < n; i++)
-    {
-	printf("    ");
-	for (j = 0; j < m; j++)
-	{
-	    if (*(*(assigned + i) + j))
-	      printf("\033[%im", (int)(30 + *(*(assigned + i) + j)));
-	    printf("%5li%s\033[m   ", (cell)(*(*(t + i) + j)), (*(*(assigned + i) + j) ? "^" : " "));
-        }
-	printf("\n\n");
-	
-	free(*(assigned + i));
-    }
-    
-    free(assigned);
+	size_t i;
+
+	if (!this->first)
+		return -1;
+
+	i = this->first - 1;
+	return (ssize_t)(lb(this->limbs[i] & -this->limbs[i]) + (i << 6));
 }
 
 
+/**
+ * Turns off a bit in a bit set
+ * 
+ * @param  this  The bit set
+ * @param  i     The index of the bit to turn off
+ */
+static void
+bitset_unset(BitSet *this, size_t i)
+{
+	size_t p, n, j = i >> 6;
+	BitSetLimb old = this->limbs[j];
+
+	this->limbs[j] &= ~(1LL << (i & 63L));
+
+	if (!this->limbs[j] ^ !old) {
+		j++;
+		p = this->prev[j];
+		n = this->next[j];
+		this->prev[n] = p;
+		this->next[p] = n;
+		if (this->first == j)
+			this->first = n;
+	}
+}
+
 
 /**
- * Calculates an optimal bipartite minimum weight matching using an
- * O(n³)-time implementation of The Hungarian Algorithm, also known
- * as Kuhn's Algorithm.
+ * Turns on a bit in a bit set
  * 
- * @param   table  The table in which to perform the matching
- * @param   n      The height of the table
- * @param   m      The width of the table
- * @return         The optimal assignment, an array of row–coloumn pairs
+ * @param  this  The bit set
+ * @param  i     The index of the bit to turn on
  */
-ssize_t** kuhn_match(cell** table, size_t n, size_t m)
+static void
+bitset_set(BitSet *this, size_t i)
 {
-    size_t i;
-    
-    /* not copying table since it will only be used once */
-    
-    kuhn_reduceRows(table, n, m);
-    byte** marks = kuhn_mark(table, n, m);
-    
-    boolean* rowCovered = malloc(n * sizeof(boolean));
-    boolean* colCovered = malloc(m * sizeof(boolean));
-    for (i = 0; i < n; i++)
-    {
-        *(rowCovered + i) = false;
-        *(colCovered + i) = false;
-    }
-    for (i = n; i < m; i++)
-        *(colCovered + i) = false;
-    
-    size_t* altRow = malloc(n * m * sizeof(ssize_t));
-    size_t* altCol = malloc(n * m * sizeof(ssize_t));
-    
-    ssize_t* rowPrimes = malloc(n * sizeof(ssize_t));
-    ssize_t* colMarks  = malloc(m * sizeof(ssize_t));
-    
-    size_t* prime;
-    
-    for (;;)
-    {
-	if (kuhn_isDone(marks, colCovered, n, m))
-	    break;
-	
-        for (;;)
-	{
-	    prime = kuhn_findPrime(table, marks, rowCovered, colCovered, n, m);
-	    if (prime != null)
-	    {
-		kuhn_altMarks(marks, altRow, altCol, colMarks, rowPrimes, prime, n, m);
-		for (i = 0; i < n; i++)
-		{
-		    *(rowCovered + i) = false;
-		    *(colCovered + i) = false;
-		}
-		for (i = n; i < m; i++)
-		    *(colCovered + i) = false;
-		free(prime);
-		break;
-	    }
-	    kuhn_addAndSubtract(table, rowCovered, colCovered, n, m);
+	size_t j = i >> 6;
+	BitSetLimb old = this->limbs[j];
+
+	this->limbs[j] |= 1LL << (i & 63L);
+
+	if (!this->limbs[j] ^ !old) {
+		j++;
+		this->prev[this->first] = j;
+		this->prev[j] = 0;
+		this->next[j] = this->first;
+		this->first = j;
 	}
-    }
-    
-    free(rowCovered);
-    free(colCovered);
-    free(altRow);
-    free(altCol);
-    free(rowPrimes);
-    free(colMarks);
-    
-    ssize_t** rc = kuhn_assign(marks, n, m);
-    
-    for (i = 0; i < n; i++)
-        free(*(marks + i));
-    free(marks);
-    
-    return rc;
 }
 
 
@@ -288,26 +202,58 @@ ssize_t** kuhn_match(cell** table, size_t n, size_t m)
  * lowest cells value is zero, and all cells' values is decrease
  * with the same value [the minium value in the row].
  * 
- * @param  t  The table in which to perform the reduction
  * @param  n  The table's height
  * @param  m  The table's width
+ * @param  t  The table in which to perform the reduction
  */
-void kuhn_reduceRows(cell** t, size_t n, size_t m)
+static void
+kuhn_reduce_rows(size_t n, size_t m, Cell **t)
 {
-    size_t i, j;
-    cell min;
-    cell* ti;
-    for (i = 0; i < n; i++)
-    {
-        ti = *(t + i);
-        min = *ti;
-	for (j = 1; j < m; j++)
-	    if (min > *(ti + j))
-	        min = *(ti + j);
-	
+	size_t i, j;
+	Cell min, *ti;
+
+	for (i = 0; i < n; i++) {
+		ti = t[i];
+		min = *ti;
+		for (j = 1; j < m; j++)
+			if (min > ti[j])
+				min = ti[j];
+		for (j = 0; j < m; j++)
+			ti[j] -= min;
+	}
+}
+
+
+/**
+ * Determines whether the marking is complete, that is
+ * if each row has a marking which is on a unique column.
+ *
+ * @param   n            The table's height
+ * @param   m            The table's width
+ * @param   marks        The marking matrix
+ * @param   col_covered  Column cover array
+ * @return               Whether the marking is complete
+ */
+static Boolean
+kuhn_is_done(size_t n, size_t m, Mark **marks, Boolean col_covered[m])
+{
+	size_t i, j, count = 0;
+
+	memset(col_covered, 0, m * sizeof(*col_covered));
+
+	for (j = 0; j < m; j++) {
+		for (i = 0; i < n; i++) {
+			if (marks[i][j] == MARKED) {
+				col_covered[j] = 1;
+				break;
+			}
+		}
+	}
+
 	for (j = 0; j < m; j++)
-	    *(ti + j) -= min;
-    }
+		count += (size_t)col_covered[j];
+
+	return count == n;
 }
 
 
@@ -316,241 +262,182 @@ void kuhn_reduceRows(cell** t, size_t n, size_t m)
  * value is zero [minimal for the row]. Each marking will
  * be on an unique row and an unique column.
  * 
- * @param   t  The table in which to perform the reduction
  * @param   n  The table's height
  * @param   m  The table's width
+ * @param   t  The table in which to perform the reduction
  * @return     A matrix of markings as described in the summary
  */
-byte** kuhn_mark(cell** t, size_t n, size_t m)
+static Mark **
+kuhn_mark(size_t n, size_t m, Cell **t)
 {
-    size_t i, j;
-    byte** marks = malloc(n * sizeof(byte*));
-    byte* marksi;
-    for (i = 0; i < n; i++)
-    {
-      marksi = *(marks + i) = malloc(m * sizeof(byte));
-        for (j = 0; j < m; j++)
- 	    *(marksi + j) = UNMARKED;
-    }
-    
-    boolean* rowCovered = malloc(n * sizeof(boolean));
-    boolean* colCovered = malloc(m * sizeof(boolean));
-    for (i = 0; i < n; i++)
-    {
-        *(rowCovered + i) = false;
-        *(colCovered + i) = false;
-    }
-    for (i = 0; i < m; i++)
-        *(colCovered + i) = false;
-    
-    for (i = 0; i < n; i++)
-        for (j = 0; j < m; j++)
-	    if ((!*(rowCovered + i)) && (!*(colCovered + j)) && (*(*(t + i) + j) == 0))
-	    {
-	        *(*(marks + i) + j) = MARKED;
-		*(rowCovered + i) = true;
-		*(colCovered + j) = true;
-	    }
-    
-    free(rowCovered);
-    free(colCovered);
-    return marks;
-}
-
+	size_t i, j;
+	Mark **marks;
+	Boolean *row_covered, *col_covered;
+
+	marks = malloc(n * sizeof(Mark *));
+	for (i = 0; i < n; i++)
+		marks[i] = calloc(m, sizeof(Mark)); /* UNMARKED == 0 */
+
+	row_covered = calloc(n, sizeof(Boolean));
+	col_covered = calloc(m, sizeof(Boolean));
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < m; j++) {
+			if (!row_covered[i] && !col_covered[j] && !t[i][j]) {
+				marks[i][j] = MARKED;
+				row_covered[i] = 1;
+				col_covered[j] = 1;
+			}
+		}
+	}
 
-/**
- * Determines whether the marking is complete, that is
- * if each row has a marking which is on a unique column.
- *
- * @param   marks       The marking matrix
- * @param   colCovered  An array which tells whether a column is covered
- * @param   n           The table's height
- * @param   m           The table's width
- * @return              Whether the marking is complete
- */
-boolean kuhn_isDone(byte** marks, boolean* colCovered, size_t n, size_t m)
-{
-    size_t i, j;
-    for (j = 0; j < m; j++)
-        for (i = 0; i < n; i++)
-	    if (*(*(marks + i) + j) == MARKED)
-	    {
-	        *(colCovered + j) = true;
-		break;
-	    }
-    
-    size_t count = 0;
-    for (j = 0; j < m; j++)
-        if (*(colCovered + j))
-	    count++;
-    
-    return count == n;
+	free(row_covered);
+	free(col_covered);
+	return marks;
 }
 
 
 /**
  * Finds a prime
  * 
- * @param   t           The table
- * @param   marks       The marking matrix
- * @param   rowCovered  Row cover array
- * @param   colCovered  Column cover array
- * @param   n           The table's height
- * @param   m           The table's width
- * @return              The row and column of the found print, <code>null</code> will be returned if none can be found
+ * @param   n            The table's height
+ * @param   m            The table's width
+ * @param   t            The table
+ * @param   marks        The marking matrix
+ * @param   row_covered  Row cover array
+ * @param   col_covered  Column cover array
+ * @param   primep       Output parameter for the row and column of the found prime
+ * @return               1 if a prime was found, 0 otherwise
  */
-size_t* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, size_t n, size_t m)
+static Boolean
+kuhn_find_prime(size_t n, size_t m, Cell **t, Mark **marks, Boolean row_covered[n], Boolean col_covered[m], CellPosition *primep)
 {
-    size_t i, j;
-    BitSet zeroes = new_BitSet(n * m);
-    
-    for (i = 0; i < n; i++)
-        if (!*(rowCovered + i))
-	    for (j = 0; j < m; j++)
-	        if ((!*(colCovered + j)) && (*(*(t + i) + j) == 0))
-		  BitSet_set(zeroes, i * m + j);
-    
-    ssize_t p;
-    size_t row, col;
-    boolean markInRow;
-    
-    for (;;)
-    {
-        p = BitSet_any(zeroes);
-	if (p < 0)
-        {
-	    free(zeroes.limbs);
-	    free(zeroes.first);
-	    free(zeroes.next);
-	    free(zeroes.prev);
-	    return null;
-	}
-	
-	row = (size_t)p / m;
-	col = (size_t)p % m;
-	
-	*(*(marks + row) + col) = PRIME;
+	size_t i, j, row, col;
+	ssize_t p;
+	Boolean mark_in_row;
+	BitSet *zeroes = bitset_create(n * m);
+
+	for (i = 0; i < n; i++)
+		if (!row_covered[i])
+			for (j = 0; j < m; j++)
+				if (!col_covered[j] && !t[i][j])
+					bitset_set(zeroes, i * m + j);
+
+	for (;;) {
+		p = bitset_any(zeroes);
+		if (p < 0) {
+			free(zeroes);
+			return 0;
+		}
+
+		row = (size_t)p / m;
+		col = (size_t)p % m;
 	
-	markInRow = false;
-	for (j = 0; j < m; j++)
-	    if (*(*(marks + row) + j) == MARKED)
-	    {
-		markInRow = true;
-		col = j;
-	    }
+		marks[row][col] = PRIME;
 	
-	if (markInRow)
-	{
-	    *(rowCovered + row) = true;
-	    *(colCovered + col) = false;
-	    
-	    for (i = 0; i < n; i++)
-	        if ((*(*(t + i) + col) == 0) && (row != i))
-		{
-		    if ((!*(rowCovered + i)) && (!*(colCovered + col)))
-		        BitSet_set(zeroes, i * m + col);
-		    else
-		        BitSet_unset(zeroes, i * m + col);
+		mark_in_row = 0;
+		for (j = 0; j < m; j++) {
+			if (marks[row][j] == MARKED) {
+				mark_in_row = 1;
+				col = j;
+			}
 		}
-	    
-	    for (j = 0; j < m; j++)
-	        if ((*(*(t + row) + j) == 0) && (col != j))
-		{
-		    if ((!*(rowCovered + row)) && (!*(colCovered + j)))
-		        BitSet_set(zeroes, row * m + j);
-		    else
-		        BitSet_unset(zeroes, row * m + j);
+
+		if (mark_in_row) {
+			row_covered[row] = 1;
+			col_covered[col] = 0;
+
+			for (i = 0; i < n; i++) {
+				if (!t[i][col] && row != i) {
+					if (!row_covered[i] && !col_covered[col])
+						bitset_set(zeroes, i * m + col);
+					else
+						bitset_unset(zeroes, i * m + col);
+				}
+			}
+
+			for (j = 0; j < m; j++) {
+				if (!t[row][j] && col != j) {
+					if (!row_covered[row] && !col_covered[j])
+						bitset_set(zeroes, row * m + j);
+					else
+						bitset_unset(zeroes, row * m + j);
+				}
+			}
+
+			if (!row_covered[row] && !col_covered[col])
+				bitset_set(zeroes, row * m + col);
+			else
+				bitset_unset(zeroes, row * m + col);
+		} else {
+			free(zeroes);
+			primep->row = row;
+			primep->col = col;
+			return 1;
 		}
-	    
-	    if ((!*(rowCovered + row)) && (!*(colCovered + col)))
-	        BitSet_set(zeroes, row * m + col);
-	    else
-	        BitSet_unset(zeroes, row * m + col);
-	}
-	else
-	{
-	    size_t* rc = malloc(2 * sizeof(size_t));
-	    *rc = row;
-	    *(rc + 1) = col;
-	    free(zeroes.limbs);
-	    free(zeroes.first);
-	    free(zeroes.next);
-	    free(zeroes.prev);
-	    return rc;
 	}
-    }
 }
 
 
 /**
  * Removes all prime marks and modifies the marking
  *
- * @param  marks      The marking matrix
- * @param  altRow     Marking modification path rows
- * @param  altCol     Marking modification path columns
- * @param  colMarks   Markings in the columns
- * @param  rowPrimes  Primes in the rows
- * @param  prime      The last found prime
- * @param  n          The table's height
- * @param  m          The table's width
+ * @param  n           The table's height
+ * @param  m           The table's width
+ * @param  marks       The marking matrix
+ * @param  alt         Marking modification paths
+ * @param  col_marks   Markings in the columns
+ * @param  row_primes  Primes in the rows
+ * @param  prime       The last found prime
  */
-void kuhn_altMarks(byte** marks, size_t* altRow, size_t* altCol, ssize_t* colMarks, ssize_t* rowPrimes, size_t* prime, size_t n, size_t m)
+static void
+kuhn_alt_marks(size_t n, size_t m, Mark **marks, CellPosition alt[n * m],
+               ssize_t col_marks[m], ssize_t row_primes[n], const CellPosition *prime)
 {
-    size_t index = 0, i, j;
-    *altRow = *prime;
-    *altCol = *(prime + 1);
-    
-    for (i = 0; i < n; i++)
-    {
-        *(rowPrimes + i) = -1;
-        *(colMarks + i) = -1;
-    }
-    for (i = n; i < m; i++)
-        *(colMarks + i) = -1;
-    
-    for (i = 0; i < n; i++)
-        for (j = 0; j < m; j++)
-	    if (*(*(marks + i) + j) == MARKED)
-	        *(colMarks + j) = (ssize_t)i;
-	    else if (*(*(marks + i) + j) == PRIME)
-	        *(rowPrimes + i) = (ssize_t)j;
-    
-    ssize_t row, col;
-    for (;;)
-    {
-        row = *(colMarks + *(altCol + index));
-	if (row < 0)
-	    break;
-	
-	index++;
-	*(altRow + index) = (size_t)row;
-	*(altCol + index) = *(altCol + index - 1);
-	
-	col = *(rowPrimes + *(altRow + index));
-	
-	index++;
-	*(altRow + index) = *(altRow + index - 1);
-	*(altCol + index) = (size_t)col;
-    }
-    
-    byte* markx;
-    for (i = 0; i <= index; i++)
-    {
-        markx = *(marks + *(altRow + i)) + *(altCol + i);
-        if (*markx == MARKED)
-	    *markx = UNMARKED;
-	else
-	    *markx = MARKED;
-    }
-    
-    byte* marksi;
-    for (i = 0; i < n; i++)
-    {
-        marksi = *(marks + i);
-        for (j = 0; j < m; j++)
-	    if (*(marksi + j) == PRIME)
-	        *(marksi + j) = UNMARKED;
-    }
+	size_t i, j, index = 0;
+	ssize_t row, col;
+	Mark *markx, *marksi;
+
+	alt[0].row = prime->row;
+	alt[0].col = prime->col;
+
+	for (i = 0; i < n; i++)
+		row_primes[i] = -1;
+
+	for (i = 0; i < m; i++)
+		col_marks[i] = -1;
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < m; j++) {
+			if (marks[i][j] == MARKED)
+				col_marks[j] = (ssize_t)i;
+			else if (marks[i][j] == PRIME)
+				row_primes[i] = (ssize_t)j;
+		}
+	}
+
+	while ((row = col_marks[alt[index].col]) >= 0) {
+		index++;
+		alt[index].row = (size_t)row;
+		alt[index].col = alt[index - 1].col;
+
+		col = row_primes[alt[index].row];
+		index++;
+		alt[index].row = alt[index - 1].row;
+		alt[index].col = (size_t)col;
+	}
+
+	for (i = 0; i <= index; i++) {
+		markx = &marks[alt[i].row][alt[i].col];
+		*markx = *markx == MARKED ? UNMARKED : MARKED;
+	}
+
+	for (i = 0; i < n; i++) {
+		marksi = marks[i];
+		for (j = 0; j < m; j++)
+			if (marksi[j] == PRIME)
+				marksi[j] = UNMARKED;
+	}
 }
 
 
@@ -559,29 +446,31 @@ void kuhn_altMarks(byte** marks, size_t* altRow, size_t* altCol, ssize_t* colMar
  * the the minimum value in the table is added, subtracted or
  * neither from the cells.
  *
- * @param  t           The table to manipulate
- * @param  rowCovered  Array that tell whether the rows are covered
- * @param  colCovered  Array that tell whether the columns are covered
- * @param  n           The table's height
- * @param  m           The table's width
+ * @param  n            The table's height
+ * @param  m            The table's width
+ * @param  t            The table to manipulate
+ * @param  row_covered  Array that tell whether the rows are covered
+ * @param  col_covered  Array that tell whether the columns are covered
  */
-void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, size_t n, size_t m)
+static void
+kuhn_add_and_subtract(size_t n, size_t m, Cell **t, Boolean row_covered[n], Boolean col_covered[m])
 {
-    size_t i, j;
-    cell min = 0x7FFFffffL;
-    for (i = 0; i < n; i++)
-        if (!*(rowCovered + i))
-	    for (j = 0; j < m; j++)
-	        if ((!*(colCovered + j)) && (min > *(*(t + i) + j)))
-		    min = *(*(t + i) + j);
-    
-    for (i = 0; i < n; i++)
-        for (j = 0; j < m; j++)
-	{
-	    if (*(rowCovered + i))
-	        *(*(t + i) + j) += min;
-	    if (*(colCovered + j) == false)
-	        *(*(t + i) + j) -= min;
+	size_t i, j;
+	Cell min = 0x7FFFFFFFL;
+
+	for (i = 0; i < n; i++)
+		if (!row_covered[i])
+			for (j = 0; j < m; j++)
+				if (!col_covered[j] && min > t[i][j])
+					min = t[i][j];
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < m; j++) {
+			if (row_covered[i])
+				t[i][j] += min;
+			if (!col_covered[j])
+				t[i][j] -= min;
+		}
 	}
 }
 
@@ -589,143 +478,165 @@ void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, siz
 /**
  * Creates a list of the assignment cells
  * 
- * @param   marks  Matrix markings
  * @param   n      The table's height
  * @param   m      The table's width
+ * @param   marks  Matrix markings
  * @return         The assignment, an array of row–coloumn pairs
  */
-ssize_t** kuhn_assign(byte** marks, size_t n, size_t m)
+static CellPosition *
+kuhn_assign(size_t n, size_t m, Mark **marks)
 {
-    ssize_t** assignment = malloc(n * sizeof(ssize_t*));
-    
-    size_t i, j;
-    for (i = 0; i < n; i++)
-    {
-        *(assignment + i) = malloc(2 * sizeof(ssize_t));
-        for (j = 0; j < m; j++)
-	    if (*(*(marks + i) + j) == MARKED)
-	    {
-		**(assignment + i) = (ssize_t)i;
-		*(*(assignment + i) + 1) = (ssize_t)j;
-	    }
-    }
-    
-    return assignment;
-}
-
+	CellPosition *assignment = malloc(n * sizeof(CellPosition));
+	size_t i, j;
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < m; j++) {
+			if (marks[i][j] == MARKED) {
+				assignment[i].row = i;
+				assignment[i].col = j;
+			}
+		}
+	}
 
-/**
- * Constructor for BitSet
- *
- * @param   size  The (fixed) number of bits to bit set should contain
- * @return        The a unique BitSet instance with the specified size
- */
-BitSet new_BitSet(size_t size)
-{
-    BitSet this;
-    
-    size_t c = size >> 6L;
-    if (size & 63L)
-        c++;
-    
-    this.limbs = malloc(c * sizeof(llong));
-    this.prev = malloc((c + 1) * sizeof(size_t));
-    this.next = malloc((c + 1) * sizeof(size_t));
-    *(this.first = malloc(sizeof(size_t))) = 0;
-    
-    size_t i;
-    for (i = 0; i < c; i++)
-    {
-        *(this.limbs + i) = 0LL;
-        *(this.prev + i) = *(this.next + i) = 0L;
-    }
-    *(this.prev + c) = *(this.next + c) = 0L;
-    
-    return this;
+	return assignment;
 }
 
-/**
- * Turns on a bit in a bit set
- * 
- * @param  this  The bit set
- * @param  i     The index of the bit to turn on
- */
-void BitSet_set(BitSet this, size_t i)
-{
-    size_t j = i >> 6L;
-    llong old = *(this.limbs + j);
-    
-    *(this.limbs + j) |= 1LL << (llong)(i & 63L);
-    
-    if ((!*(this.limbs + j)) ^ (!old))
-    {
-        j++;
-	*(this.prev + *(this.first)) = j;
-	*(this.prev + j) = 0;
-	*(this.next + j) = *(this.first);
-	*(this.first) = j;
-    }
-}
 
 /**
- * Turns off a bit in a bit set
+ * Calculates an optimal bipartite minimum weight matching using an
+ * O(n³)-time implementation of The Hungarian Algorithm, also known
+ * as Kuhn's Algorithm.
  * 
- * @param  this  The bit set
- * @param  i     The index of the bit to turn off
+ * @param   n      The height of the table
+ * @param   m      The width of the table
+ * @param   table  The table in which to perform the matching
+ * @return         The optimal assignment, an array of row–coloumn pairs
  */
-void BitSet_unset(BitSet this, size_t i)
+static CellPosition *
+kuhn_match(size_t n, size_t m, Cell **table)
 {
-    size_t j = i >> 6L;
-    llong old = *(this.limbs + j);
-    
-    *(this.limbs + j) &= ~(1LL << (llong)(i & 63L));
-    
-    if ((!*(this.limbs + j)) ^ (!old))
-    {
-        j++;
-	size_t p = *(this.prev + j);
-	size_t n = *(this.next + j);
-	*(this.prev + n) = p;
-	*(this.next + p) = n;
-	if (*(this.first) == j)
-	    *(this.first) = n;
-    }
+	size_t i;
+	ssize_t *row_primes, *col_marks;
+	Mark **marks;
+	Boolean *row_covered, *col_covered;
+	CellPosition *ret, prime, *alt;
+
+	/* Not copying table since it will only be used once. */
+
+	row_covered = calloc(n, sizeof(Boolean));
+	col_covered = calloc(m, sizeof(Boolean));
+
+	row_primes = malloc(n * sizeof(ssize_t));
+	col_marks  = malloc(m * sizeof(ssize_t));
+
+	alt = malloc(n * m * sizeof(CellPosition));
+
+	kuhn_reduce_rows(n, m, table);
+	marks = kuhn_mark(n, m, table);
+
+	while (!kuhn_is_done(n, m, marks, col_covered)) {
+		while (!kuhn_find_prime(n, m, table, marks, row_covered, col_covered, &prime))
+			kuhn_add_and_subtract(n, m, table, row_covered, col_covered);
+		kuhn_alt_marks(n, m, marks, alt, col_marks, row_primes, &prime);
+		memset(row_covered, 0, n * sizeof(*row_covered));
+		memset(col_covered, 0, m * sizeof(*col_covered));
+	}
+
+	free(row_covered);
+	free(col_covered);
+	free(alt);
+	free(row_primes);
+	free(col_marks);
+
+	ret = kuhn_assign(n, m, marks);
+
+	for (i = 0; i < n; i++)
+		free(marks[i]);
+	free(marks);
+
+	return ret;
 }
 
-/**
- * Gets the index of any set bit in a bit set
- * 
- * @param   this  The bit set
- * @return        The index of any set bit
- */
-ssize_t BitSet_any(BitSet this)
+
+
+static void
+print(size_t n, size_t m, Cell **t, CellPosition assignment[n])
 {
-    if (*(this.first) == 0L)
-        return -1;
-    
-    size_t i = *(this.first) - 1;
-    return (ssize_t)(lb(*(this.limbs + i) & -*(this.limbs + i)) + (i << 6L));
+	size_t i, j, (*assigned)[n][m];
+
+	assigned = calloc(1, sizeof(ssize_t [n][m]));
+
+	if (assignment)
+		for (i = 0; i < n; i++)
+			(*assigned)[assignment[i].row][assignment[i].col] += 1;
+
+	for (i = 0; i < n; i++) {
+		printf("    ");
+		for (j = 0; j < m; j++) {
+			if ((*assigned)[i][j])
+				printf("\033[%im", (int)(30 + (*assigned)[i][j]));
+			printf("%5li%s\033[m   ", (Cell)t[i][j], (*assigned)[i][j] ? "^" : " ");
+		}
+		printf("\n\n");
+	}
+
+	free(assigned);
 }
 
 
-/**
- * Calculates the floored binary logarithm of a positive integer
- *
- * @param   value  The integer whose logarithm to calculate
- * @return         The floored binary logarithm of the integer
- */
-size_t lb(llong value)
+int
+main(int argc, char *argv[])
 {
-    size_t rc = 0;
-    llong v = value;
-    
-    if (v & (int_fast64_t)0xFFFFFFFF00000000LL)  {  rc |= 32L;  v >>= 32LL;  }
-    if (v & (int_fast64_t)0x00000000FFFF0000LL)  {  rc |= 16L;  v >>= 16LL;  }
-    if (v & (int_fast64_t)0x000000000000FF00LL)  {  rc |=  8L;  v >>=  8LL;  }
-    if (v & (int_fast64_t)0x00000000000000F0LL)  {  rc |=  4L;  v >>=  4LL;  }
-    if (v & (int_fast64_t)0x000000000000000CLL)  {  rc |=  2L;  v >>=  2LL;  }
-    if (v & (int_fast64_t)0x0000000000000002LL)     rc |=  1L;
-    
-    return rc;
-}
+	FILE *urandom;
+	unsigned int seed;
+	size_t i, j, n, m;
+	Cell **t, **table, x, sum = 0;
+	CellPosition *assignment;
+
+	urandom = fopen("/dev/urandom", "r");
+	fread(&seed, sizeof(unsigned int), 1, urandom);
+	srand(seed);
+	fclose(urandom);
+
+	n     = argc < 3 ? 10 : (size_t)atol(argv[1]);
+	m     = argc < 3 ? 15 : (size_t)atol(argv[2]);
+	t     = malloc(n * sizeof(Cell *));
+	table = malloc(n * sizeof(Cell *));
+
+	if (argc < 3) {
+		for (i = 0; i < n; i++) {
+			t[i]     = malloc(m * sizeof(Cell));
+			table[i] = malloc(m * sizeof(Cell));
+			for (j = 0; j < m; j++)
+				table[i][j] = t[i][j] = (Cell)(random() & 63);
+		}
+	} else {
+		for (i = 0; i < n; i++) {
+			t[i]     = malloc(m * sizeof(Cell));
+			table[i] = malloc(m * sizeof(Cell));
+			for (j = 0; j < m; j++) {
+				scanf("%li", &x);
+				table[i][j] = t[i][j] = x;
+			}
+		}
+	}
+
+	printf("\nInput:\n\n");
+	print(n, m, t, NULL);
+
+	assignment = kuhn_match(n, m, table);
+	printf("\nOutput:\n\n");
+	print(n, m, t, assignment);
 
+	for (i = 0; i < n; i++) {
+		sum += t[assignment[i].row][assignment[i].col];
+		free(table[i]);
+		free(t[i]);
+	}
+	free(assignment);
+	free(table);
+	free(t);
+	printf("\n\nSum: %li\n\n", sum);
+
+	return 0;
+}