X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?p=matchsticks-search.git;a=blobdiff_plain;f=main.c;h=6f0df725949c9b26fccb613768f74c283bc9e591;hp=538073724fa159d3734140bb3ed336be43a1990e;hb=0169095a9c345a65ee9e091df63a0e52d3910d1d;hpb=2abc59d8ece71404651d3a8d1701ddf6fa024b5b diff --git a/main.c b/main.c index 5380737..6f0df72 100644 --- a/main.c +++ b/main.c @@ -1,3 +1,26 @@ +/* + * Searches for "good" ways to divide n matchsticks up and reassemble them + * into m matchsticks. "Good" means the smallest fragment is as big + * as possible. + * + * Invoke as ./main n m + * + * The algorithm is faster if the arguments are ordered so that n > m. + */ + +/* + * matchsticks/main.c Copyright 2014 Ian Jackson + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + */ #include #include @@ -10,6 +33,48 @@ #include #include +/* + * Algorithm. + * + * Each input match contributes, or does not contribute, to each + * output match; we do not need to consider multiple fragments + * relating to the same input/output pair this gives an n*m adjacency + * matrix (bitmap). Given such an adjacency matrix, the problem of + * finding the best sizes for the fragments can be expressed as a + * linear programming problem. + * + * We search all possible adjacency matrices, and for each one we run + * GLPK's simplex solver. We represent the adjacency matrix as an + * array of bitmaps. + * + * However, there are a couple of wrinkles: + * + * To best represent the problem as a standard LP problem, we separate + * out the size of each fragment into a common minimum size variable, + * plus a fragment-specific extra size variable. This reduces the LP + * problem size at the cost of making the problem construction, and + * interpretation of the results, a bit fiddly. + * + * Many of the adjacency matrices are equivalent. In particular, + * permutations of the columns, or of the rows, do not change the + * meaning. It is only necessasry to consider any one permutation. + * We make use of this by considering only adjacency matrices whose + * bitmap array contains bitmap words whose numerical values are + * nondecreasing in array order. + * + * Once we have a solution, we also avoid considering any candidate + * which involves dividing one of the output sticks into so many + * fragment that the smallest fragment would necessarily be no bigger + * than our best solution. That is, we reject candidates where any of + * the hamming weights of the adjacency bitmap words are too large. + * + * And, we want to do the search in order of increasing maximum + * hamming weight. This is because in practice optimal solutions tend + * to have low hamming weight, and having found a reasonable solution + * early allows us to eliminate a lot of candidates without doing the + * full LP. + */ + typedef uint32_t AdjWord; #define PRADJ "08"PRIx32 @@ -45,9 +110,17 @@ static int count_set_adj_bits(AdjWord w) { } static void optimise(int doprint) { + /* Consider the best answer (if any) for a given adjacency matrix */ glp_prob *prob = 0; int i, j, totalfrags; + /* + * Up to a certain point, optimise() can be restarted. We use this + * to go back and print the debugging output if it turns out that we + * have an interesting case. The HAVE_PRINTED macro does this: its + * semantics are to go back in time and make sure that we have + * printed the description of the search case. + */ #define HAVE_PRINTED ({ \ if (!doprint) { doprint = 1; goto retry_with_print; } \ }) @@ -57,7 +130,7 @@ static void optimise(int doprint) { prob = 0; } -#define PRINTF if (!doprint) ; else printf /* bodgy */ +#define PRINTF(...) if (!doprint) ; else fprintf(stderr, __VA_ARGS__) /* bodgy */ PRINTF("%2d ", maxhamweight); @@ -74,6 +147,10 @@ static void optimise(int doprint) { } } if (!had_max) { + /* Skip this candidate as its max hamming weight is lower than + * we're currently looking for (which means we must have done it + * already). (The recursive iteration ensures that none of the + * words have more than the max hamming weight.) */ PRINTF(" nomaxham"); goto out; } @@ -128,8 +205,8 @@ static void optimise(int doprint) { int ME_totals_j__minimum = next_matrix_entry; for (j=0; j maxhamweight) goto again; + if (i == 0 && (adjmatrix[i] & (1+adjmatrix[i]))) + goto again; iterate_recurse(i+1, adjmatrix[i]); @@ -277,9 +356,32 @@ int main(int argc, char **argv) { m = atoi(argv[2]); prep(); iterate(); - printf("\n"); - if (best_prob) - glp_print_sol(best_prob,"/dev/stdout"); + fprintf(stderr, "\n"); + if (best_prob) { + double min = glp_get_obj_val(best_prob); + double a[n][m]; + int i, j, cols; + for (i = 0; i < n; i++) + for (j = 0; j < m; j++) + a[i][j] = 0; + cols = glp_get_num_cols(best_prob); + for (i = 1; i <= cols; i++) { + int x, y; + if (2 != sscanf(glp_get_col_name(best_prob, i), "mf %d,%d", &x, &y)) + continue; + a[x][y] = min + glp_get_col_prim(best_prob, i); + } + printf("%d into %d: min fragment %g\n", n, m, min); + for (i = 0; i < n; i++) { + for (j = 0; j < m; j++) { + if (a[i][j]) + printf(" %9.3f", a[i][j]); + else + printf(" "); + } + printf("\n"); + } + } if (ferror(stdout) || fclose(stdout)) { perror("stdout"); exit(-1); } return 0; }