X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=blobdiff_plain;f=main.c;h=83728fa7a4aee7cb0168bbdc96d2bd689483669d;hb=3a1fd673483e2e2f17853cc73be3b7f4ffea459c;hp=e246752ef153b4272fb333f13e511fd05aa38252;hpb=34b443b35453ef92904c2ec6d22a998d6078d028;p=matchsticks-search.git diff --git a/main.c b/main.c index e246752..83728fa 100644 --- a/main.c +++ b/main.c @@ -1,27 +1,105 @@ +/* + * 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 #include +#include +#include +#include #include #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 -static int n, m; +static int n, m, maxhamweight; static AdjWord *adjmatrix; static AdjWord adjall; static double best; +static glp_prob *best_prob; +static AdjWord *best_adjmatrix; + +static unsigned printcounter; + +static AdjWord *xalloc_adjmatrix(void) { + return xmalloc(sizeof(*adjmatrix)*n); +} static void prep(void) { adjall = ~((~(AdjWord)0) << m); - adjmatrix = xmalloc(sizeof(*adjmatrix)*n); + adjmatrix = xalloc_adjmatrix(); + glp_term_out(GLP_OFF); } static AdjWord one_adj_bit(int bitnum) { - return (AdjWord)1 << j; + return (AdjWord)1 << bitnum; } static int count_set_adj_bits(AdjWord w) { @@ -31,18 +109,51 @@ static int count_set_adj_bits(AdjWord w) { return total; } -static void optimise(void) { - int i, totalfrags; +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; } \ + }) + retry_with_print: + if (prob) { + glp_delete_prob(prob); + prob = 0; + } + +#define PRINTF(...) if (!doprint) ; else fprintf(stderr, __VA_ARGS__) /* bodgy */ + + PRINTF("%2d ", maxhamweight); + + bool had_max = 0; for (i=0, totalfrags=0; i= 0 */ - glp_set_col_bounds(prob, X_minimum, GLP_LB, 0, 0); + glp_set_col_bnds(prob, X_minimum, GLP_LO, 0, 0); + glp_set_col_name(prob, X_minimum, "minimum"); /* objective is maximising x_minimum */ glp_set_obj_dir(prob, GLP_MAX); glp_set_obj_coef(prob, X_minimum, 1); - for (i=0; i= 0 */ int X_morefrag_i_j = glp_add_cols(prob, 1); glp_set_col_bnds(prob, X_morefrag_i_j, GLP_LO, 0, 0); + if (doprint) { + char buf[255]; + snprintf(buf,sizeof(buf),"mf %d,%d",i,j); + glp_set_col_name(prob, X_morefrag_i_j, buf); + } /* x_total_i += x_morefrag_i_j */ /* x_total_j += x_morefrag_i_j */ @@ -131,37 +247,139 @@ static void optimise(void) { assert(next_matrix_entry == matrix_entries_size); - for (row=1; row<=rows; row++) { - glp_load_matrix(prob, next_matrix_entry, - matrix_entries_XY[1], matrix_entries_XY[0], - matrix_entries); + glp_load_matrix(prob, matrix_entries_size-1, + matrix_entries_XY[1], matrix_entries_XY[0], + matrix_entries); + + int r = glp_simplex(prob, NULL); + PRINTF(" glp=%d", r); + +#define OKERR(e) \ + case e: PRINTF(" " #e ); goto out; +#define BADERR(e) \ + case e: HAVE_PRINTED; printf(" " #e " CRASHING\n"); exit(-1); +#define DEFAULT \ + default: HAVE_PRINTED; printf(" ! CRASHING\n"); exit(-1); + + switch (r) { + OKERR(GLP_ESING); + OKERR(GLP_ECOND); + OKERR(GLP_EBOUND); + OKERR(GLP_EFAIL); + OKERR(GLP_ENOPFS); + OKERR(GLP_ENODFS); + BADERR(GLP_EBADB); + BADERR(GLP_EOBJLL); + BADERR(GLP_EOBJUL); + BADERR(GLP_EITLIM); + BADERR(GLP_ETMLIM); + BADERR(GLP_EINSTAB); + BADERR(GLP_ENOCVG); + case 0: break; + DEFAULT; + } + + r = glp_get_status(prob); + PRINTF(" status=%d", r); + + switch (r) { + OKERR(GLP_NOFEAS); + OKERR(GLP_UNDEF); + BADERR(GLP_FEAS); + BADERR(GLP_INFEAS); + BADERR(GLP_UNBND); + case GLP_OPT: break; + DEFAULT; + } + + double got = glp_get_obj_val(prob); + PRINTF(" %g", got); + if (got <= best) + goto out; + + HAVE_PRINTED; + + best = got; + + if (best_prob) glp_delete_prob(best_prob); + best_prob = prob; + + free(best_adjmatrix); + best_adjmatrix = xalloc_adjmatrix(); + memcpy(best_adjmatrix, adjmatrix, sizeof(*adjmatrix)*n); - printf("nyi\n"); + PRINTF(" BEST \n"); + return; + + } + out: + if (prob) + glp_delete_prob(prob); + if (doprint) { PRINTF(" \r"); fflush(stdout); } } static void iterate_recurse(int i, AdjWord min) { if (i >= n) { - optimise(); + printcounter++; + optimise(!(printcounter & 0xfff)); return; } for (adjmatrix[i] = min; ; adjmatrix[i]++) { + if (count_set_adj_bits(adjmatrix[i]) > maxhamweight) + goto again; + iterate_recurse(i+1, adjmatrix[i]); + + again: if (adjmatrix[i] == adjall) return; } } static void iterate(void) { - iterate_recurse(0, 1); + for (maxhamweight=1; maxhamweight<=m; maxhamweight++) { + double maxminsize = (double)m / maxhamweight; + if (maxminsize <= best) + continue; + + iterate_recurse(0, 1); + } } int main(int argc, char **argv) { + assert(argc==3); n = atoi(argv[1]); m = atoi(argv[2]); prep(); iterate(); + 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; }