Code to compute a maximum vertex-weighted matching for a specific graph. Used in an answer here: https://math.stackexchange.com/questions/3138907/.

3138907.c

// This program was used to partially answer https://math.stackexchange.com/questions/3138907. The program can be optimized a little but the main bottleneck is memory accesses when traversing the graph so don't expect any large speedups.

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdint.h>

const uint64_t N = 6;

uint64_t binomial(uint64_t n, uint64_t k) {
    uint64_t r = 1;
    uint64_t d;
    if (k>n) return 0;
    for (d=1; d <= k; d++) {
        r *= n;
        n--;
        r /= d;
    }
    return r;
}

uint8_t next(uint64_t* x, uint64_t m) {
    uint64_t i = m;
    while (i > 0) {
        i--;
        if (x[i] > 0) {
            x[i]--;
            i++;
            break;
        }
    }
    if (i == 0) {
        return 0;
    }
    uint64_t val = x[i-1];
    for(;i<m;i++){
        x[i] = val;
    }
    return 1;
}

uint64_t* create_table(uint64_t n, uint64_t m) {
    uint64_t max;
    // This is a cheap hack
    if (n == m) {
        max = binomial(2*N-1, N);
    } else {
        max = binomial(2*N-2, N-1);
    }
    uint64_t* out = malloc(m * max * sizeof *out);
    uint64_t* temp = malloc(m * sizeof *temp);
    if (out == 0 || temp == 0) {
        printf("could not allocate\n");
    }
    for(uint64_t i = 0; i<m;i++) {
        temp[i] = n-1;
    }
    uint64_t i = 0;
    do {
        memcpy(out+(i*m), temp, m * sizeof *temp);
        i++;
    } while (next(temp, m));
    free(temp);
    return out;
}

int compare_arrays(const void* a, const void* b) {
    uint64_t* p = (uint64_t*) a;
    uint64_t* q = (uint64_t*) b;
    for(uint64_t i = 0; i < N-1; i++) {
        if (p[i] < q[i]) {
            return -1;
        } else if (p[i] > q[i]) {
            return 1;
        }
    }
    return 0;
}

uint64_t row_score(uint64_t* row, uint64_t n) {
    uint64_t out = 1;
    uint64_t streak = 1;
    uint64_t last = row[0];
    for(uint64_t i = 1; i < N; i++) {
        if (row[i] == last) {
            streak++;
        } else {
            last = row[i];
            out *= binomial(n, streak);
            streak = 1;
        }
    }
    out *= binomial(n, streak);
    return out;
}

// This algorithm is taken from the thesis 'ALGORITHMS FOR VERTEX-WEIGHTED
// MATCHING IN GRAPHS' by Mahantesh Halappanavar. This is a part of the
// algorithm 'LocalOptimal'. We only need to do one part of it because in this
// application, the "small side"/"right side" is already treated as all zeros.
uint64_t local_optimal(
    uint64_t* right_side,
    uint64_t* weights,
    uint64_t size_left,
    uint64_t size_right
    // If this algorithm would be used on arbitrary graphs, we would need to
    // know the degree of every node to iterate over its neighbours but we know
    // that every node on the smaller side has N neighbours.
    // uint64_t* right_side_lengths
) {
    uint64_t* matching_right = malloc(size_right * sizeof *matching_right);
    uint64_t* matching_left = malloc(size_left * sizeof *matching_left);
    for(uint64_t i = 0; i < size_right; i++) {
        matching_right[i] = size_left;
    }
    for(uint64_t i = 0; i < size_left; i++) {
        matching_left[i] = size_right;
    }


    // Shows which iteration a node was last visited
    uint64_t* visited = calloc(size_left, sizeof *visited);
    uint64_t* visited_right = malloc(size_right * sizeof *visited_right);
    uint64_t visited_right_len = 0;
    uint64_t* visited_left = malloc(size_left * sizeof *visited_right);
    uint64_t visited_left_len = 0;
    uint64_t max_end = 0;
    uint8_t max_end_needs_updating = 1;
    uint64_t* best_path = malloc(size_right * 2 * sizeof *best_path);
    uint64_t* stack = malloc(size_right * 2 * sizeof * stack);
    for(uint64_t rn = 0; rn < size_right; rn++) {
        if (matching_right[rn] != size_left) {
            continue;
        }
        if (max_end_needs_updating) {
            max_end = 0;
            for(uint64_t i = 0; i < size_left; i++) {
                if ((weights[i] > max_end) && (matching_left[i] == size_right)) {
                    max_end = weights[i];
                }
            }
            max_end_needs_updating = 0;
        }
        visited_right_len = 0;
        visited_left_len = 0;
        visited_right[visited_right_len++] = rn;
        uint64_t path_length = 0;
        uint64_t best_score = 0;
        stack[0] = 0;
        uint64_t stack_len = 1;
        uint64_t current_index = N*visited_right[visited_right_len - 1]+stack[stack_len-1];
        while(1) {
            uint64_t nb = right_side[current_index];

            // Check if we've reached this node this iteration
            if(visited[nb] != rn+1) {
                visited[nb] = rn+1;
                uint64_t next = matching_left[nb];
                if (next != size_right) {
                    // 'nb' is part of a matching so we go to its neighbor in
                    // the matching.
                    visited_left[visited_left_len++] = nb;
                    visited_right[visited_right_len++] = next;
                    stack[stack_len++] = 0;
                    current_index = N*next;
                    continue;
                } else if (weights[nb] > best_score) {
                    // We found a new best path!
                    path_length = 2*visited_right_len;
                    best_score = weights[nb];
                    for(uint64_t i = 0; i < path_length/2 - 1;i++) {
                        best_path[2*i] = visited_right[i];
                        best_path[2*i+1] = visited_left[i];
                    }
                    best_path[path_length-2] = visited_right[visited_right_len-1];
                    best_path[path_length-1] = nb;
                    if (max_end == best_score) {
                        max_end_needs_updating = 1;
                        break;
                    }
                }
            }
            stack[stack_len - 1]++;
            current_index++;
            if (stack[stack_len - 1] == N) {
                stack_len--;
                visited_left_len--;
                visited_right_len--;
                if (stack_len == 0) {
                    break;
                }
                current_index = N*visited_right[visited_right_len - 1]+stack[stack_len-1];
            }
        }
        // When we've found the best path, we perform a symmetric difference
        // between the path and our matching. This change our matching so that
        // it includes the node at the end of the path, and all the nodes that
        // the matching previously used on the left side.
        for(uint64_t i = 0; i < path_length; i+=2) {
            matching_right[best_path[i]] = best_path[i+1];
            matching_left[best_path[i+1]] = best_path[i];
        }
    }
    free(stack);
    free(best_path);
    free(visited);
    free(visited_right);
    free(visited_left);
    free(matching_left);

    uint64_t score = 0;
    for(uint64_t i = 0; i < size_right; i++) {
        uint64_t x = matching_right[i];
        if (x < size_left) {
            score += weights[x];
        }
    }
    free(matching_right);
    return score;
}

int main(void) {
    // Generate the graph
    uint64_t* high = create_table(N, N);
    uint64_t* low = create_table(N, N-1);
    uint64_t low_len = binomial(2*N-2, N-1);
    uint64_t high_len = binomial(2*N-1, N);
    printf("low_len: %ld, high_size: %ld\n", low_len, high_len);

    uint64_t* gs = malloc(N * low_len * sizeof *gs);
    uint64_t* gs_lengths = calloc(low_len, sizeof *gs_lengths);
    uint64_t* scratch = malloc((N-1) * sizeof *scratch);
    qsort(low, low_len, sizeof(*low) * (N-1), compare_arrays);
    for(uint64_t i = 0; i < high_len;i++) {
        uint64_t* row = high+(i*N);
        for(uint64_t k = 0; k < N; k++) {
            memcpy(scratch, row, k*sizeof(*row));
            memcpy(scratch+k, row+k+1, sizeof(*row)*(N-k-1));
            uint64_t* pel = bsearch(scratch, low, low_len, (N-1) * sizeof *low, compare_arrays);
            uint64_t el = (pel - &low[0])/(N-1);
            uint8_t found = 0;
            for(uint64_t jj = 0; jj < gs_lengths[el]; jj++) {
                if (gs[el*N+jj] == i) {
                    found = 1;
                }
            }

            if (!found) {
                gs[el*N+gs_lengths[el]] = i;
                gs_lengths[el]++;
            }
        }
    }
    uint64_t* wt = malloc(N * high_len * sizeof(*high));
    for(uint64_t i = 0; i < high_len; i++) {
        wt[i] = row_score(high+(i*N), N);
    }

    // Calculate the total weight of the matching with the highest weight
    uint64_t result = local_optimal(gs, wt, high_len, low_len);
    printf("%ld\n", result);
    return 1;
}