the code

1 files changed, 706 insertions, 0 deletions

diff --git a/hungarian.c b/hungarian.c
new file mode 100644
index 0000000..7171cf3
--- /dev/null
+++ b/hungarian.c
@@ -0,0 +1,706 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+
+#define cell      long
+#define CELL_STR  "%li"
+
+#define llong    long long
+#define byte     char
+#define boolean  long
+#define null     0
+#define true     1
+#define false    0
+
+
+
+//new cell[X]
+#define new_cells(X)  new_longs(X)
+
+//new boolean[X]
+#define new_booleans(X)  new_longs(X)
+
+//new byte[X]
+#define new_bytes(X)  malloc(X)
+
+//new llong[X]
+#define new_llongs(X)  malloc((X) << 3)
+
+//new long[X]
+#if !(defined __LP64__ || defined __LLP64__)
+   #define new_longs(X)  malloc((X) << 2) /*32-bit*/
+#else
+  #define new_longs(X)  malloc((X) << 3) /*64-bit*/
+#endif
+
+//new float[X]
+#define new_floats(X)  malloc((X) << 2)
+
+//new double[X]
+#define new_doubles(X)  malloc((X) << 3)
+
+//new ?[][X]
+#define new_arrays(X)  new_longs(X)
+
+
+#ifdef DEBUG
+  #define debug(X) fprintf(stderr, "\e[31m%s\e[m\n", X)
+#else
+  #define debug(X) 
+#endif
+
+
+
+
+/**
+ * Cell marking:  none
+ */
+#define UNMARKED  0L
+
+/**
+ * Cell marking:  marked
+ */
+#define MARKED    1L
+
+/**
+ * Cell marking:  prime
+ */
+#define PRIME     2L
+
+
+
+/**
+ * Bit set, a set of fixed number of bits/booleans
+ */
+typedef struct
+{
+    /**
+     * The set of all limbs, a limb consist of 64 bits
+     */
+    cell* limbs;
+
+    /**
+     * Singleton array with the index of the first non-zero limb
+     */
+    long* first;
+
+    /**
+     * Array the the index of the previous non-zero limb for echo limb
+     */
+    long* prev;
+
+    /**
+     * Array the the index of the next non-zero limb for echo limb
+     */
+    long* next;
+
+} BitSet;
+
+
+
+long** kuhn_match(cell** table, long n, long m);
+void kuhn_reduceRows(cell** t, long n, long m);
+byte** kuhn_mark(cell** t, long n, long m);
+boolean kuhn_isDone(byte** marks, boolean* colCovered, long n, long m);
+long* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, long n, long m);
+void kuhn_altMarks(byte** marks, long* altRow, long* altCol, long* colMarks, long* rowPrimes, long* prime, long n, long m);
+void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, long n, long m);
+long** kuhn_assign(byte** marks, long n, long m);
+
+BitSet new_BitSet(long size);
+void BitSet_set(BitSet this, long i);
+void BitSet_unset(BitSet this, long i);
+long BitSet_any(BitSet this);
+
+long lb(llong x);
+
+
+
+void print(cell** t, long n, long m, long** assignment);
+
+int main(int argc, char** argv)
+{
+    unsigned a, d;
+    asm("cpuid");
+    asm volatile("rdtsc" : "=a" (a), "=d" (d));
+    srand(((llong)a) | (((llong)d) << 32LL));
+
+  
+    long n = 10, m = 15;
+    
+    long i;
+    cell** t = new_arrays(n);
+    cell** table = new_arrays(n);
+    long j;
+    if (argc < 2)
+        for (i = 0; i < n; i++)
+	{
+	    *(t + i) = new_llongs(m);
+	    *(table + i) = new_llongs(m);
+	    for (j = 0; j < m; j++)
+	        *(*(table + i) + j) = *(*(t + i) + j) = (cell)(random() & 63);
+	}
+    else
+    {
+        long x;
+	scanf("%li", &n);
+	scanf("%li", &m);
+        t = new_arrays(n);
+        table = new_arrays(n);
+        for (i = 0; i < n; i++)
+	{
+	    *(t + i) = new_llongs(m);
+	    *(table + i) = new_llongs(m);
+	    for (j = 0; j < m; j++)
+	    {
+	        scanf(CELL_STR, &x);
+	        *(*(table + i) + j) = *(*(t + i) + j) = x;
+	    }
+	}
+    }
+    
+    //kuhn_match(table, n, m);
+    
+    printf("\nInput:\n\n");
+    print(t, n, m, null);
+    
+    long** 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));
+    printf("\n\nSum: %li\n\n", sum);
+
+    return 0;
+}
+
+void print(cell** t, long n, long m, long** assignment)
+{
+    long i, j;
+
+    long** assigned = new_arrays(n);
+    for (i = 0; i < n; i++)
+    {
+        *(assigned + i) = new_longs(m);
+	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("\e[%lim", 30 + *(*(assigned + i) + j));
+	    printf("%5li%s\e[m   ", (cell)(*(*(t + i) + j)), (*(*(assigned + i) + j) ? "^" : " "));
+        }
+	printf("\n\n");
+    }
+}
+
+
+
+/**
+ * Calculates an optimal bipartite minimum weight matching using an
+ * O(n³)-time implementation of The Hungarian Algorithm, also known
+ * as Kuhn's Algorithm. This implemention is restricted to square
+ * tables.
+ * 
+ * @param   table  The table in which to perform the matching
+ * @param   n      The dimension of the table
+ * @return         The optimal assignment, an array of row–coloumn pairs
+ */
+long** kuhn_match(cell** table, long n, long m)
+{
+    long 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 = new_booleans(n);
+    boolean* colCovered = new_booleans(m);
+    for (i = 0; i < n; i++)
+    {
+        *(rowCovered + i) = false;
+        *(colCovered + i) = false;
+    }
+    for (i = n; i < m; i++)
+        *(colCovered + i) = false;
+
+    long* altRow = new_longs(n * m);
+    long* altCol = new_longs(n * m);
+    
+    long* rowPrimes = new_longs(n);
+    long* colMarks  = new_longs(m);
+
+    long* 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;
+		break;
+	    }
+	    kuhn_addAndSubtract(table, rowCovered, colCovered, n, m);
+	}
+    }
+    
+    return kuhn_assign(marks, n, m);
+}
+
+
+/**
+ * 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  t  The table in which to perform the reduction
+ * @param  n  The table's height
+ * @param  m  The table's width
+ */
+void kuhn_reduceRows(cell** t, long n, long m)
+{
+    long 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);
+
+	for (j = 0; j < m; j++)
+	    *(ti + j) -= min;
+    }
+}
+
+
+/**
+ * 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   t  The table in which to perform the reduction
+ * @param   n  The table's height
+ * @param   m  The table's width
+ * @return     A matrix of markings as described in the summary
+ */
+byte** kuhn_mark(cell** t, long n, long m)
+{
+    long i, j;
+    byte** marks = new_arrays(n);
+    byte* marksi;
+    for (i = 0; i < n; i++)
+    {
+	marksi = *(marks + i) = new_bytes(m);
+        for (j = 0; j < m; j++)
+ 	    *(marksi + j) = UNMARKED;
+    }
+
+    boolean* rowCovered = new_booleans(n);
+    boolean* colCovered = new_booleans(m);
+    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;
+	    }
+    
+    return marks;
+}
+
+
+/**
+ * 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, long n, long m)
+{
+    long i, j;
+    for (j = 0; j < m; j++)
+        for (i = 0; i < n; i++)
+	    if (*(*(marks + i) + j) == MARKED)
+	    {
+	        *(colCovered + j) = true;
+		break;
+	    }
+    
+    long count = 0;
+    for (j = 0; j < m; j++)
+        if (*(colCovered + j))
+	    count++;
+        
+    return count == n;
+}
+
+
+/**
+ * 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
+ */
+long* kuhn_findPrime(cell** t, byte** marks, boolean* rowCovered, boolean* colCovered, long n, long m)
+{
+    long i, j;
+    BitSet zeroes = new_BitSet(n * n);
+    
+    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);
+    
+    long p, row, col;
+    boolean markInRow;
+    
+    for (;;)
+    {
+        p = BitSet_any(zeroes);
+	if (p < 0)
+	    return null;
+
+	row = p / m;
+	col = p % m;
+	
+	*(*(marks + row) + col) = PRIME;
+	
+	markInRow = false;
+	for (j = 0; j < m; j++)
+	    if (*(*(marks + row) + j) == MARKED)
+	    {
+		markInRow = true;
+		col = j;
+	    }
+            
+	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);
+		}
+                
+	    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 ((!*(rowCovered + row)) && (!*(colCovered + col)))
+	        BitSet_set(zeroes, row * m + col);
+	    else
+	        BitSet_unset(zeroes, row * m + col);
+	}
+	else
+	{
+	    long* rc = new_longs(2);
+	    *rc = row;
+	    *(rc + 1) = col;
+	    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
+ */
+void kuhn_altMarks(byte** marks, long* altRow, long* altCol, long* colMarks, long* rowPrimes, long* prime, long n, long m)
+{
+    long 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) = i;
+	    else if (*(*(marks + i) + j) == PRIME)
+	        *(rowPrimes + i) = j;
+    
+    long row, col;
+    for (;;)
+    {
+        row = *(colMarks + *(altCol + index));
+	if (row < 0)
+	    break;
+	
+	index++;
+	*(altRow + index) = row;
+	*(altCol + index) = *(altCol + index - 1);
+	
+	col = *(rowPrimes + *(altRow + index));
+	
+	index++;
+	*(altRow + index) = *(altRow + index - 1);
+	*(altCol + index) = 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;
+    }
+}
+
+
+/**
+ * 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  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
+ */
+void kuhn_addAndSubtract(cell** t, boolean* rowCovered, boolean* colCovered, long n, long m)
+{
+    long 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;
+	}
+}
+
+
+/**
+ * Creates a list of the assignment cells
+ * 
+ * @param   marks  Matrix markings
+ * @param   n      The table's height
+ * @param   m      The table's width
+ * @return         The assignment, an array of row–coloumn pairs
+ */
+long** kuhn_assign(byte** marks, long n, long m)
+{
+    long** assignment = new_arrays(n);
+
+    long i, j;
+    for (i = 0; i < n; i++)
+    {
+        *(assignment + i) = new_longs(2);
+        for (j = 0; j < m; j++)
+	    if (*(*(marks + i) + j) == MARKED)
+	    {
+	        **(assignment + i) = i;
+		*(*(assignment + i) + 1) = j;
+	    }
+    }
+
+    return assignment;
+}
+
+
+/**
+ * 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(long size)
+{
+    BitSet this;
+    
+    long c = size >> 6L;
+    if (size & 63L)
+        c++;
+    
+    this.limbs = new_llongs(c);
+    this.prev = new_longs(c + 1L);
+    this.next = new_longs(c + 1L);
+    *(this.first = new_longs(1)) = 0;
+    
+    long 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;
+}
+
+/**
+ * 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, long i)
+{
+    long 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
+ * 
+ * @param  this  The bit set
+ * @param  i     The index of the bit to turn off
+ */
+void BitSet_unset(BitSet this, long i)
+{
+    long j = i >> 6L;
+    llong old = *(this.limbs + j);
+    
+    *(this.limbs + j) &= ~(1LL << (llong)(i & 63L));
+  
+    if ((!*(this.limbs + j)) ^ (!old))
+    {
+	j++;
+	long p = *(this.prev + j);
+	long n = *(this.next + j);
+	*(this.prev + n) = p;
+	*(this.next + p) = n;
+	if (*(this.first) == j)
+	    *(this.first) = *(this.next + j);
+    }
+}
+
+/**
+ * Gets the index of any set bit in a bit set
+ * 
+ * @param   this  The bit set
+ * @return        The index of any set bit
+ */
+long BitSet_any(BitSet this)
+{
+    if (*(this.first) == 0L)
+        return -1;
+            
+    long i = *(this.first) - 1L;
+    return lb(*(this.limbs + i) & -*(this.limbs + i)) + (i << 6L);
+}
+
+
+/**
+ * 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
+ */
+long lb(llong value)
+{
+    long rc = 0L;
+    llong v = value;
+            
+    if (v & 0xFFFFFFFF00000000LL)  {  rc |= 32L;  v >>= 32LL;  }
+    if (v & 0x00000000FFFF0000LL)  {  rc |= 16L;  v >>= 16LL;  }
+    if (v & 0x000000000000FF00LL)  {  rc |=  8L;  v >>=  8LL;  }
+    if (v & 0x00000000000000F0LL)  {  rc |=  4L;  v >>=  4LL;  }
+    if (v & 0x000000000000000CLL)  {  rc |=  2L;  v >>=  2LL;  }
+    if (v & 0x0000000000000002LL)     rc |=  1L;
+    
+    return rc;
+}
+