diff options
| -rw-r--r-- | COPYING | 2 | ||||
| -rw-r--r-- | Makefile | 12 | ||||
| -rw-r--r-- | hungarian.c | 1091 |
3 files changed, 506 insertions, 599 deletions
@@ -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 @@ -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; +} |