aboutsummaryrefslogtreecommitdiffstats
path: root/hungarian.c
diff options
context:
space:
mode:
Diffstat (limited to 'hungarian.c')
-rw-r--r--hungarian.c1091
1 files changed, 501 insertions, 590 deletions
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;
+}
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2024-09-19 02:10:57 +0000