/** * š¯“˛(nĀ³) implementation of the Hungarian algorithm * * 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 * and/or modify it under the terms of the Do What The Fuck You Want * To Public License, Version 2, as published by Sam Hocevar. See * http://sam.zoy.org/wtfpl/COPYING for more details. */ #include #include #include #include #include #include /** * Cell markings **/ enum { UNMARKED = 0, MARKED, PRIME }; /** * Value type for marking */ typedef int_fast8_t Mark; /** * Value type for cells */ typedef signed long int Cell; 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 */ 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; /** * 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; } /** * 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) { 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; } /** * 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; 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; } } /** * Turns on a bit in a bit set * * @param this The bit set * @param i The index of the bit to turn on */ static void bitset_set(BitSet *this, size_t i) { 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; } } /** * Reduces the values on each rows so that, for each row, the * lowest cells value is zero, and all cells' values is decrease * with the same value [the minium value in the row]. * * @param n The table's height * @param m The table's width * @param t The table in which to perform the reduction */ static void kuhn_reduce_rows(size_t n, size_t m, Cell **t) { 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++) count += (size_t)col_covered[j]; return count == n; } /** * Create a matrix with marking of cells in the table whose * value is zero [minimal for the row]. Each marking will * be on an unique row and an unique column. * * @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 */ static Mark ** kuhn_mark(size_t n, size_t m, Cell **t) { 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; } } } free(row_covered); free(col_covered); return marks; } /** * Finds a prime * * @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 */ 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, 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; marks[row][col] = PRIME; mark_in_row = 0; for (j = 0; j < m; j++) { if (marks[row][j] == MARKED) { mark_in_row = 1; col = 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; } } } /** * Removes all prime marks and modifies the marking * * @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 */ 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 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; } } /** * Depending on whether the cells' rows and columns are covered, * the the minimum value in the table is added, subtracted or * neither from the cells. * * @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 */ 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 (!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; } } } /** * Creates a list of the assignment cells * * @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 */ static CellPosition * kuhn_assign(size_t n, size_t m, Mark **marks) { 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; } } } return assignment; } /** * Calculates an optimal bipartite minimum weight matching using an * O(nĀ³)-time implementation of The Hungarian Algorithm, also known * as Kuhn's Algorithm. * * @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 */ static CellPosition * kuhn_match(size_t n, size_t m, Cell **table) { 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; } static void print(size_t n, size_t m, Cell **t, CellPosition assignment[n]) { 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); } int main(int argc, char *argv[]) { 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; }