chiark / gitweb /
Add -std=gnu99 to Makefile
[matchsticks-search.git] / main.c
1 /*
2  * Searches for "good" ways to divide n matchsticks up and reassemble them
3  * into m matchsticks.  "Good" means the smallest fragment is as big
4  * as possible.
5  *
6  * Invoke as   ./main n m
7  *
8  * The algorithm is faster if the arguments are ordered so that n > m.
9  */
10
11 /*
12  * matchsticks/main.c  Copyright 2014 Ian Jackson
13  *
14  * This program is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * This program is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22  * GNU General Public License for more details.
23  */
24
25 #include <stdio.h>
26 #include <stdint.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <assert.h>
30 #include <unistd.h>
31 #include <stdbool.h>
32 #include <inttypes.h>
33
34 #include <publib.h>
35 #include <glpk.h>
36
37 /*
38  * Algorithm.
39  *
40  * Each input match contributes, or does not contribute, to each
41  * output match; we do not need to consider multiple fragments
42  * relating to the same input/output pair this gives an n*m adjacency
43  * matrix (bitmap).  Given such an adjacency matrix, the problem of
44  * finding the best sizes for the fragments can be expressed as a
45  * linear programming problem.
46  *
47  * We search all possible adjacency matrices, and for each one we run
48  * GLPK's simplex solver.  We represent the adjacency matrix as an
49  * array of bitmaps.
50  *
51  * However, there are a couple of wrinkles:
52  *
53  * To best represent the problem as a standard LP problem, we separate
54  * out the size of each fragment into a common minimum size variable,
55  * plus a fragment-specific extra size variable.  This reduces the LP
56  * problem size at the cost of making the problem construction, and
57  * interpretation of the results, a bit fiddly.
58  *
59  * Many of the adjacency matrices are equivalent.  In particular,
60  * permutations of the columns, or of the rows, do not change the
61  * meaning.  It is only necessasry to consider any one permutation.
62  * We make use of this by considering only adjacency matrices whose
63  * bitmap array contains bitmap words whose numerical values are
64  * nondecreasing in array order.
65  *
66  * Once we have a solution, we also avoid considering any candidate
67  * which involves dividing one of the output sticks into so many
68  * fragment that the smallest fragment would necessarily be no bigger
69  * than our best solution.  That is, we reject candidates where any of
70  * the hamming weights of the adjacency bitmap words are too large.
71  *
72  * And, we want to do the search in order of increasing maximum
73  * hamming weight.  This is because in practice optimal solutions tend
74  * to have low hamming weight, and having found a reasonable solution
75  * early allows us to eliminate a lot of candidates without doing the
76  * full LP.
77  */
78
79 typedef uint32_t AdjWord;
80 #define PRADJ "08"PRIx32
81
82 static int n, m, maxhamweight;
83 static AdjWord *adjmatrix;
84 static AdjWord adjall;
85
86 static double best;
87 static glp_prob *best_prob;
88 static AdjWord *best_adjmatrix;
89
90 static unsigned printcounter;
91
92 static int ncpus = 1;
93
94 static AdjWord *xalloc_adjmatrix(void) {
95   return xmalloc(sizeof(*adjmatrix)*n);
96 }
97
98 static void prep(void) {
99   adjall = ~((~(AdjWord)0) << m);
100   adjmatrix = xalloc_adjmatrix();
101   glp_term_out(GLP_OFF);
102 }
103
104 static AdjWord one_adj_bit(int bitnum) {
105   return (AdjWord)1 << bitnum;
106 }
107
108 static int count_set_adj_bits(AdjWord w) {
109   int j, total;
110   for (j=0, total=0; j<m; j++)
111     total += !!(w & one_adj_bit(j));
112   return total;
113 }
114
115 static void optimise(int doprint) {
116   /* Consider the best answer (if any) for a given adjacency matrix */
117   glp_prob *prob = 0;
118   int i, j, totalfrags;
119
120   /*
121    * Up to a certain point, optimise() can be restarted.  We use this
122    * to go back and print the debugging output if it turns out that we
123    * have an interesting case.  The HAVE_PRINTED macro does this: its
124    * semantics are to go back in time and make sure that we have
125    * printed the description of the search case.
126    */
127 #define HAVE_PRINTED ({                                         \
128       if (!doprint) { doprint = 1; goto retry_with_print; }     \
129     })
130  retry_with_print:
131   if (prob) {
132     glp_delete_prob(prob);
133     prob = 0;
134   }
135
136 #define PRINTF(...) if (!doprint) ; else fprintf(stderr, __VA_ARGS__) /* bodgy */
137
138   PRINTF("%2d ", maxhamweight);
139
140   bool had_max = 0;
141   for (i=0, totalfrags=0; i<n; i++) {
142     int frags = count_set_adj_bits(adjmatrix[i]);
143     had_max += (frags == maxhamweight);
144     totalfrags += frags;
145     PRINTF("%"PRADJ" ", adjmatrix[i]);
146     double maxminsize = (double)m / frags;
147     if (maxminsize <= best) {
148       PRINTF(" too fine");
149       goto out;
150     }
151   }
152   if (!had_max) {
153     /* Skip this candidate as its max hamming weight is lower than
154      * we're currently looking for (which means we must have done it
155      * already).  (The recursive iteration ensures that none of the
156      * words have more than the max hamming weight.) */
157     PRINTF(" nomaxham");
158     goto out;
159   }
160
161   /*
162    * We formulate our problem as an LP problem as follows.
163    * In this file "n" and "m" are the matchstick numbers.
164    *
165    * Each set bit in the adjacency matrix corresponds to taking a
166    * fragment from old match i and making it part of new match j.
167    *
168    * The structural variables (columns) are:
169    *   x_minimum        minimum size of any fragment (bounded below by 0)
170    *   x_morefrag_i_j   the amount by which the size of the fragment
171    *                     i,j exceeds the minimum size (bounded below by 0)
172    *
173    * The auxiliary variables (rows) are:
174    *   x_total_i       total length for each input match (fixed variable)
175    *   x_total_j       total length for each output match (fixed variable)
176    *
177    * The objective function is simply
178    *   maximise x_minimum
179    *
180    * We use X_ and Y_ to refer to GLPK's (1-based) column and row indices.
181    * ME_ refers to entries in the list of constraint matrix elements
182    * which we build up as we go.
183    */
184
185   prob = glp_create_prob();
186
187   int Y_totals_i = glp_add_rows(prob, n);
188   int Y_totals_j = glp_add_rows(prob, m);
189   int X_minimum = glp_add_cols(prob, 1);
190
191   {
192   int next_matrix_entry = 1; /* wtf GLPK! */
193   int matrix_entries_size = next_matrix_entry + n + m + totalfrags*2;
194   double matrix_entries[matrix_entries_size];
195   int matrix_entries_XY[2][matrix_entries_size];
196
197 #define ADD_MATRIX_ENTRY(Y,X) ({                        \
198       assert(next_matrix_entry < matrix_entries_size);  \
199       matrix_entries_XY[0][next_matrix_entry] = (X);    \
200       matrix_entries_XY[1][next_matrix_entry] = (Y);    \
201       matrix_entries[next_matrix_entry] = 0;            \
202       next_matrix_entry++;                              \
203     })
204
205   int ME_totals_i__minimum = next_matrix_entry;
206   for (i=0; i<n; i++) ADD_MATRIX_ENTRY(Y_totals_i+i, X_minimum);
207
208   int ME_totals_j__minimum = next_matrix_entry;
209   for (j=0; j<m; j++) ADD_MATRIX_ENTRY(Y_totals_j+j, X_minimum);
210
211   /* \forall_i x_total_i = m */
212   /* \forall_i x_total_j = n */
213   for (i=0; i<n; i++) glp_set_row_bnds(prob, Y_totals_i+i, GLP_FX, m,m);
214   for (j=0; j<m; j++) glp_set_row_bnds(prob, Y_totals_j+j, GLP_FX, n,n);
215
216   /* x_minimum >= 0 */
217   glp_set_col_bnds(prob, X_minimum, GLP_LO, 0, 0);
218   glp_set_col_name(prob, X_minimum, "minimum");
219
220   /* objective is maximising x_minimum */
221   glp_set_obj_dir(prob, GLP_MAX);
222   glp_set_obj_coef(prob, X_minimum, 1);
223
224   for (i=0; i<n; i++) {
225     for (j=0; j<m; j++) {
226       if (!(adjmatrix[i] & one_adj_bit(j)))
227         continue;
228       /* x_total_i += x_minimum */
229       /* x_total_j += x_minimum */
230       matrix_entries[ ME_totals_i__minimum + i ] ++;
231       matrix_entries[ ME_totals_j__minimum + j ] ++;
232
233       /* x_morefrag_i_j >= 0 */
234       int X_morefrag_i_j = glp_add_cols(prob, 1);
235       glp_set_col_bnds(prob, X_morefrag_i_j, GLP_LO, 0, 0);
236       if (doprint) {
237         char buf[255];
238         snprintf(buf,sizeof(buf),"mf %d,%d",i,j);
239         glp_set_col_name(prob, X_morefrag_i_j, buf);
240       }
241
242       /* x_total_i += x_morefrag_i_j */
243       /* x_total_j += x_morefrag_i_j */
244       int ME_totals_i__mf_i_j = ADD_MATRIX_ENTRY(Y_totals_i+i, X_morefrag_i_j);
245       int ME_totals_j__mf_i_j = ADD_MATRIX_ENTRY(Y_totals_j+j, X_morefrag_i_j);
246       matrix_entries[ME_totals_i__mf_i_j] = 1;
247       matrix_entries[ME_totals_j__mf_i_j] = 1;
248     }
249   }
250
251   assert(next_matrix_entry == matrix_entries_size);
252
253   glp_load_matrix(prob, matrix_entries_size-1,
254                   matrix_entries_XY[1], matrix_entries_XY[0],
255                   matrix_entries);
256
257   int r = glp_simplex(prob, NULL);
258   PRINTF(" glp=%d", r);
259
260 #define OKERR(e) \
261   case e: PRINTF(" " #e ); goto out;
262 #define BADERR(e) \
263   case e: HAVE_PRINTED; printf(" " #e " CRASHING\n"); exit(-1);
264 #define DEFAULT \
265   default: HAVE_PRINTED; printf(" ! CRASHING\n"); exit(-1);
266
267   switch (r) {
268   OKERR(GLP_ESING);
269   OKERR(GLP_ECOND);
270   OKERR(GLP_EBOUND);
271   OKERR(GLP_EFAIL);
272   OKERR(GLP_ENOPFS);
273   OKERR(GLP_ENODFS);
274   BADERR(GLP_EBADB);
275   BADERR(GLP_EOBJLL);
276   BADERR(GLP_EOBJUL);
277   BADERR(GLP_EITLIM);
278   BADERR(GLP_ETMLIM);
279   BADERR(GLP_EINSTAB);
280   BADERR(GLP_ENOCVG);
281   case 0: break;
282   DEFAULT;
283   }
284
285   r = glp_get_status(prob);
286   PRINTF(" status=%d", r);
287
288   switch (r) {
289   OKERR(GLP_NOFEAS);
290   OKERR(GLP_UNDEF);
291   BADERR(GLP_FEAS);
292   BADERR(GLP_INFEAS);
293   BADERR(GLP_UNBND);
294   case GLP_OPT: break;
295   DEFAULT;
296   }
297
298   double got = glp_get_obj_val(prob);
299   PRINTF("  %g", got);
300   if (got <= best)
301     goto out;
302
303   HAVE_PRINTED;
304
305   best = got;
306
307   if (best_prob) glp_delete_prob(best_prob);
308   best_prob = prob;
309
310   free(best_adjmatrix);
311   best_adjmatrix = xalloc_adjmatrix();
312   memcpy(best_adjmatrix, adjmatrix, sizeof(*adjmatrix)*n);
313
314   PRINTF(" BEST        \n");
315   return;
316
317   }
318  out:
319   if (prob)
320     glp_delete_prob(prob);
321   if (doprint) { PRINTF("        \r"); fflush(stdout); }
322 }
323
324 static void iterate_recurse(int i, AdjWord min) {
325   if (i >= n) {
326     printcounter++;
327     optimise(!(printcounter & 0xfff));
328     return;
329   }
330   for (adjmatrix[i] = min;
331        ;
332        adjmatrix[i]++) {
333     if (count_set_adj_bits(adjmatrix[i]) > maxhamweight)
334       goto again;
335     if (i == 0 && (adjmatrix[i] & (1+adjmatrix[i])))
336       goto again;
337
338     iterate_recurse(i+1, adjmatrix[i]);
339
340   again:
341     if (adjmatrix[i] == adjall)
342       return;
343   }
344 }
345
346 static void iterate(void) {
347   for (maxhamweight=1; maxhamweight<=m; maxhamweight++) {
348     double maxminsize = (double)m / maxhamweight;
349     if (maxminsize <= best)
350       continue;
351
352     iterate_recurse(0, 1);
353   }
354 }
355
356 static void report(void) {
357   fprintf(stderr, "\n");
358   if (best_prob) {
359     double min = glp_get_obj_val(best_prob);
360     double a[n][m];
361     int i, j, cols;
362     for (i = 0; i < n; i++)
363       for (j = 0; j < m; j++)
364         a[i][j] = 0;
365     cols = glp_get_num_cols(best_prob);
366     for (i = 1; i <= cols; i++) {
367       int x, y;
368       if (2 != sscanf(glp_get_col_name(best_prob, i), "mf %d,%d", &x, &y))
369         continue;
370       a[x][y] = min + glp_get_col_prim(best_prob, i);
371     }
372     printf("%d into %d: min fragment %g\n", n, m, min);
373     for (i = 0; i < n; i++) {
374       for (j = 0; j < m; j++) {
375         if (a[i][j])
376           printf(" %9.3f", a[i][j]);
377         else
378           printf("          ");
379       }
380       printf("\n");
381     }
382   }
383   if (ferror(stdout) || fclose(stdout)) { perror("stdout"); exit(-1); }
384 }
385  
386 int main(int argc, char **argv) {
387   int opt;
388   while ((opt = getopt(argc,argv,"j:")) >= 0) {
389     switch (opt) {
390     case 'j': ncpus = atoi(optarg); break;
391     case '+': assert(!"bad option");
392     default: abort();
393     }
394   }
395   argc -= optind-1;
396   argv += optind-1;
397   assert(argc==3);
398   n = atoi(argv[1]);
399   m = atoi(argv[2]);
400
401   prep();
402   iterate();
403   report();
404   return 0;
405 }