crooks

combinatorial rooks
git clone https://a3nm.net/git/crooks/
Log | Files | Refs

main.c (18368B)

```      1 #include <signal.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <time.h>
5 #include <inttypes.h>
6
7 // compilation:
8 //   gcc -std=c99 -O2 -Wall main.c
9 // usage:
10 //   ./a.out N K [SEED]
11 //   SEED is used to seed the RNG, otherwise time is used;
12 //   used seed is always printed first
13 //   N is grid size, K is number of rooks (program will increase K when it
14 //   succeeded)
15
16 // general idea: always take one of the rooks with most violations and move to a
17 // cell with least violations; when stuck for UBOUND steps, do a totally random
18 // move
19 // in all intermediate configurations there are never two rooks in direct
20 // attack, the program tries to minimize the number of cells attacked by wrong
21 // number of rooks
22 // TODO: optimize
23 // TODO: profile
24 // TODO: plot asymptotics to make conjectures
25 // TODO: optimize UBOUND
26
27 // upper bound on grid size (N) for static allocation
28 #define MAXN 64
29 // number of useless steps to make before doing a random move
30 #define UBOUND 10
31
32 // C1/C2 is the probability of going on doing an additional random move
33 #define C1 3
34 #define C2 5
35
36 const int debug = 0;
37
38 // N: grid size
39 // K: number of rooks
40 int N, K;
41 // stores if a cell is taken by a rook or not
42 // coordinates are z, x, y
43 // corresponding to height, line, col
44 unsigned char G[MAXN][MAXN][MAXN];
45 // stores a number indicating by which directions the cell is attacked
46 // same coordinates
47 // numbers are in [0, 8[:, LSB is z, 2nd LSB is x, 3rd LSB is y
48 unsigned char A[MAXN][MAXN][MAXN];
49 // stores the coordinates of each rook
50 // 0 is z (height), 1 is x (line), 2 is y (col)
51 unsigned char R[MAXN*MAXN][3];
52 // stores the number of violations (number of attacked cells attacked by two
53 // other rooks) for each cell (no matter whether it contains a rook or not)
54 unsigned char V[MAXN][MAXN][MAXN];
55
56 // proj_xy[i]&(1<<j) is true, iff there is a rook of coordinates
57 // (x=i, y=j, z=anything)
58 // In other words, it is the projection along the z axis giving a matrice of
59 // booleans, with the booleans then grouped in bitfields along the y axis.
60 // Similarily for the others.
61 uint64_t proj_xy[MAXN];
62 uint64_t proj_xz[MAXN];
63 uint64_t proj_yx[MAXN];
64 uint64_t proj_yz[MAXN];
65 uint64_t proj_zx[MAXN];
66 uint64_t proj_zy[MAXN];
67 // v_zx[k][i] is the number of cases already attacked in both the z and the
68 // x direction, in the line with coordinates (z=k, x=i).
69 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
70 // V[z][x][y] = v_zx[z][x] + v_xy[x][y] + v_zy[z][y]
71 uint8_t v_zx[MAXN][MAXN];
72 uint8_t v_xy[MAXN][MAXN];
73 uint8_t v_zy[MAXN][MAXN];
74
75 void print();
76
77 int popcount64(uint64_t bitfield) {
78 	return __builtin_popcountll((unsigned long long) bitfield);
79 }
80
81 int getVEmpty(int z, int x, int y) {
82 	int result;
83 	result = (v_zx[z][x] + v_xy[x][y] + v_zy[z][y]);
84 	return result;
85 }
86
87 // TODO: optimise in position where we know that !G[z][x][y]
88 int getV(int z, int x, int y) {
89 	int result;
90 	result = getVEmpty(z, x, y);
91 	if (G[z][x][y])
92 		result -= 3; // Do not count a rook as being in violation with itself.
93 	return result;
94 }
95
96 void intHandler() {
97   // needed for profiling
98   printf("SIGINT caught, exiting\n");
99   exit(1);
100 }
101
102 void print() {
103   // print the current grid and also rook positions when debug
104   for (int x = 0; x < N; x++) {
105     for (int y = 0; y < N; y++) {
106       int pos = -1;
107       for (int z = 0; z < N; z++) {
108         if (G[z][x][y]) {
109           pos = z;
110           break;
111         }
112       }
113       if (pos < 0)
114         printf("* ");
115       else {
116         if (pos < 9)
117           printf("%d ", pos + 1);
118         else
119           printf("%c ", 'A' + (pos - 9));
120       }
121     }
122     printf("\n");
123   }
124   if (debug) {
125     for (int i = 0; i < K; i++)
126     printf("rook %d at %d %d %d\n", i, R[i][1], R[i][2], R[i][0]);
127   }
128 }
129
130 void print_attacks() {
131   // print the attack pattern in each cell
132   for (int z = 0; z < N; z++) {
133     printf("== %d ==\n", z);
134     for (int x = 0; x < N; x++) {
135       for (int y = 0; y < N; y++) {
136         int val = A[z][x][y];
137         printf("%d%d%d ", val & (1 << 0), (val & (1 << 1))>>1, (val& (1 << 2))>>2);
138       }
139       printf("\n");
140     }
141     printf("\n");
142   }
143 }
144
145 void print_all_violations() {
146   // print the number of violation in each cell
147   for (int z = 0; z < N; z++) {
148     printf("== %d ==\n", z);
149     for (int x = 0; x < N; x++) {
150       for (int y = 0; y < N; y++) {
151         int val = V[z][x][y];
152         printf("%d ", val);
153       }
154       printf("\n");
155     }
156     printf("\n");
157   }
158 }
159
160 void print_violations() {
161   // print the number of violations of the rooks
162   int tot = 0;
163   for (int i = 0; i < K; i++) {
164     int val = V[R[i][0]][R[i][1]][R[i][2]];
165     printf("rook %d has %d violations\n", i,
166         val);
167     tot += val;
168   }
169   printf("=> %d total violations\n", tot);
170 }
171
172 int upd_atk_rm(int i) {
173   // update A to reflect the fact that rook i is about to be removed
174   int z = R[i][0];
175   int x = R[i][1];
176   int y = R[i][2];
177   uint64_t mask_z = 1ull << z;
178   uint64_t mask_x = 1ull << x;
179   uint64_t mask_y = 1ull << y;
180   uint64_t old_vzx1, old_vzx2, old_vzy1, old_vzy2, old_vxy1, old_vxy2;
187
188 //printf("rmv rook at %d %d %d\n", z, x, y);
189   int new_delta = 0;
190   for (int p = 0; p < N; p++) {
191     A[p][x][y] &= ~(1 << 0); // p x y attacked by z x y along direction z
192     A[z][p][y] &= ~(1 << 1); // z p y attacked by z x y along direction x
193     A[z][x][p] &= ~(1 << 2); // z x p attacked by z x y along direction y
194
195 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
196 // 			   v_zy[z][y] = popcount(proj_yx[y] & proj_zx[z])
197 // 			   v_xy[x][y] = popcount(proj_yz[y] & proj_xz[x]);
198     // TODO: refactor
199     old_vzx1 = v_zx[z][p];
200     old_vzx2 = v_zx[p][x];
201     old_vzy1 = v_zy[z][p];
202     old_vzy2 = v_zy[p][y];
203     old_vxy1 = v_xy[x][p];
204     old_vxy2 = v_xy[p][y];
205     v_zx[z][p] = popcount64(proj_xy[p] & proj_zy[z]);
206     v_zx[p][x] = popcount64(proj_xy[x] & proj_zy[p]);
207     v_zy[z][p] = popcount64(proj_yx[p] & proj_zx[z]);
208     v_zy[p][y] = popcount64(proj_yx[y] & proj_zx[p]);
209     v_xy[x][p] = popcount64(proj_yz[p] & proj_xz[x]);
210     v_xy[p][y] = popcount64(proj_yz[y] & proj_xz[p]);
212 	    new_delta += v_zx[z][p] - old_vzx1;
214 	    new_delta += v_zx[p][x] - old_vzx2;
216 	    new_delta += v_zy[z][p] - old_vzy1;
218 	    new_delta += v_zy[p][y] - old_vzy2;
220 	    new_delta += v_xy[x][p] - old_vxy1;
222 	    new_delta += v_xy[p][y] - old_vxy2;
223   }
224   new_delta -= getV(z, x, y);
225   return new_delta;
226 }
227
229   // update A to reflect the fact that rook i was just added
230   int z = R[i][0];
231   int x = R[i][1];
232   int y = R[i][2];
233   uint64_t mask_z = 1ull << z;
234   uint64_t mask_x = 1ull << x;
235   uint64_t mask_y = 1ull << y;
236   uint64_t old_vzx1, old_vzx2, old_vzy1, old_vzy2, old_vxy1, old_vxy2;
243   //printf("added rook at %d %d %d\n", z, x, y);
244   int new_delta = 0;
245   for (int p = 0; p < N; p++) {
246     A[p][x][y] |= 1 << 0;
247     A[z][p][y] |= 1 << 1;
248     A[z][x][p] |= 1 << 2;
249
250 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
251 // 			   v_zy[z][y] = popcount(proj_yx[y] & proj_zx[z])
252 // 			   v_xy[x][y] = popcount(proj_yz[y] & proj_xz[x]);
253     // TODO: refactor
254     old_vzx1 = v_zx[z][p];
255     old_vzx2 = v_zx[p][x];
256     old_vzy1 = v_zy[z][p];
257     old_vzy2 = v_zy[p][y];
258     old_vxy1 = v_xy[x][p];
259     old_vxy2 = v_xy[p][y];
260     v_zx[z][p] = popcount64(proj_xy[p] & proj_zy[z]);
261     v_zx[p][x] = popcount64(proj_xy[x] & proj_zy[p]);
262     v_zy[z][p] = popcount64(proj_yx[p] & proj_zx[z]);
263     v_zy[p][y] = popcount64(proj_yx[y] & proj_zx[p]);
264     v_xy[x][p] = popcount64(proj_yz[p] & proj_xz[x]);
265     v_xy[p][y] = popcount64(proj_yz[y] & proj_xz[p]);
266     if (p != z) {
268 		    new_delta += v_zx[p][x] - old_vzx2;
270 		    new_delta += v_zy[p][y] - old_vzy2;
271     }
272     if (p != x) {
274 		    new_delta += v_zx[z][p] - old_vzx1;
276 		    new_delta += v_xy[p][y] - old_vxy2;
277     }
278     if (p != y) {
280 		    new_delta += v_zy[z][p] - old_vzy1;
282 		    new_delta += v_xy[x][p] - old_vxy1;
283     }
284   }
285   new_delta += getV(z, x, y);
286   return new_delta;
287 }
288
289 void num_attacks() {
290   // generate A from scratch
291   for (int i = 0; i < N; i++)
292     for (int j = 0; j < N; j++)
293       for (int k = 0; k < N; k++)
294         A[i][j][k] = 0;
295   for (int i = 0; i < K; i++) {
296     if (R[i][0] < MAXN)
298   }
299 }
300
301 void check_attacks() {
302   // debug: regenerate A from scratch
303   // and check that the previous version was already correct
304   int backup[MAXN][MAXN][MAXN];
305   for (int i = 0; i < N; i++)
306     for (int j = 0; j < N; j++)
307       for (int k = 0; k < N; k++)
308         backup[i][j][k] = A[i][j][k];
309   num_attacks();
310   for (int i = 0; i < N; i++)
311     for (int j = 0; j < N; j++)
312       for (int k = 0; k < N; k++)
313         if (backup[i][j][k] != A[i][j][k]) {
314           print_attacks();
315           printf("ERROR at %d %d %d\n", i, j, k);
316           exit(1);
317         }
318 }
319
320 int num_violations() {
321   // regenerate V from scratch
322   // returns total number of violations
323   // TODO: find a way to update this rather than recompute
324   int tot = 0;
325   for (int i = 0; i < N; i++)
326     for (int j = 0; j < N; j++)
327       for (int k = 0; k < N; k++)
328         V[i][j][k] = 0; // clear
329   for (int z = 0; z < N; z++)
330     for (int x = 0; x < N; x++)
331       for (int y = 0; y < N; y++) {
332         for (int p = 0; p < N; p++) {
333           // must count cells that have already two attack directions and are
334           // missing the attack direction from the current cell
335           if (p != z && (A[p][x][y] | (1 << 0)) == 7)
336           {V[z][x][y]++;} // putting rook in z x y would attack p x y thrice
337           if (p != x && (A[z][p][y] | (1 << 1)) == 7)
338           {V[z][x][y]++;} // putting rook in z x y would attack z p y thrice
339           if (p != y && (A[z][x][p] | (1 << 2)) == 7)
340           {V[z][x][y]++;} // putting rook in z x y would attack z x p thrice
341         }
342         if (G[z][x][y]) {
343           tot += V[z][x][y];
344           //printf("add %d to num_violations for %d %d %d\n", V[z][x][y], x, y, z);
345         }
346       }
348 }
349
350 void check_increment(int ntot) {
351   //printf("CHECK: before\n");
352   //print_all_violations();
353   // debug: regenerate A, V from scratch
354   // and check that the previous version was already correct
355   int backupA[MAXN][MAXN][MAXN];
356   int backupV[MAXN][MAXN][MAXN];
357   for (int i = 0; i < N; i++)
358     for (int j = 0; j < N; j++)
359       for (int k = 0; k < N; k++) {
360         backupA[i][j][k] = A[i][j][k];
361         backupV[i][j][k] = V[i][j][k];
362       }
363   num_attacks();
364   int rntot = num_violations();
365   //printf("CHECK: after\n");
366   //print_all_violations();
367   for (int i = 0; i < N; i++)
368     for (int j = 0; j < N; j++)
369       for (int k = 0; k < N; k++) {
370         if (backupA[i][j][k] != A[i][j][k]) {
371           print_attacks();
372           printf("ERROR attacks at %d %d %d\n", i, j, k);
373           exit(1);
374         }
375         if (backupV[i][j][k] != V[i][j][k]) {
376           print();
377           print_attacks();
378           print_all_violations();
379           printf("ERROR violations at z=%d x=%d y=%d, had %d computed %d\n", i, j, k, backupV[i][j][k], V[i][j][k]);
380           exit(1);
381         }
382       }
383   if (ntot != rntot) {
384     print_all_violations();
385     print_violations();
386     printf("total violation miscount, had %d computed %d\n", ntot, rntot);
387     exit(1);
388   }
389 }
390
391
392 void init_pos() {
393   // create initial position
394   for (int i = 0; i < N; i++)
395     for (int j = 0; j < N; j++)
396       for (int k = 0; k < N; k++)
397         G[i][j][k] = 0;
398   int remaining = K;
399   int height = 0, pos = 0;
400   while (remaining > 0) {
401     int z = height;
402     int x = pos;
403     int y = (pos + height) % N;
404     G[z][x][y] = 1;
405     remaining--;
406     R[remaining][0] = z;
407     R[remaining][1] = x;
408     R[remaining][2] = y;
409     pos++;
410     if (pos >= N) {
411       pos = 0;
412       height++;
413     }
414   }
415 }
416
417 int choose_rook(int prev) {
418   // choose the next rook to move: a random one among the ones with maximal
419   // violations, forbiding prev (the previously moved rook)
420   // TODO: better notion including "age" of rook to favor moving old rooks (?)
421   // and small probability to move rooks with not exactly maximal num of viols
422   // TODO: make it faster (priority queue?)
423   int best = -1;
424   int nvs[MAXN*MAXN];
425   int nnvs = 0;
426   int nv;
427   for (int i = 0; i < K; i++) {
428     if (i == prev)
429       continue;
430     if ((nv = getV(R[i][0], R[i][1], R[i][2])) > best) {
431       best = nv;
432       nvs[0] = i;
433       nnvs = 1;
434     } else if (nv == best) {
435       nvs[nnvs++] = i;
436     }
437   }
438   // choose random rook among the optimal ones
439   return nvs[rand() % nnvs];
440 }
441
442 int random_move() {
443   // make a totally random move (when we are stuck)
444   int rook = (rand() % K);
445   //printf("move %d at random\n", rook);
446   int z = R[rook][0];
447   int x = R[rook][1];
448   int y = R[rook][2];
449   int delta = 0;
450   G[z][x][y] = 0;
451   delta += upd_atk_rm(rook);
452   int nx, ny, nz;
454   do {
455     nx = (rand() % N);
456     ny = (rand() % N);
457     nz = (rand() % N);
458     mask_x = 1ull << x;
459     mask_y = 1ull << y;
460     if (!A[nz][nx][ny] && !G[nz][nx][ny]) {
461 #if 0
465 #endif
466       break; // found a spot that's both empty and not in direct attack
467     }
468     // TODO: nothing tests if this loop stands a chance of terminating, it would
469     // be neater to exit the program then
471   G[nz][nx][ny] = 1;
472   R[rook][0] = nz;
473   R[rook][1] = nx;
474   R[rook][2] = ny;
476   return delta;
477 }
478
479 int position(int i) {
480 	// position rook i that was just removed
481 	// TODO: more efficient (pqueue?)
482 	// TODO more clever (chance to select cell with non-optimal number of
483 	// violations)
484 	int best = N*N*N+1;
485 	int bs[MAXN*MAXN*MAXN][3];
486 	int nbs = 0;
487 	int v_zxy;
489 	for (int z = 0; z < N; z++) {
490 		for (int x = 0; x < N; x++) {
491 			mask_x = 1ull << x;
492 			if((proj_zx[z] & mask_x) || (v_zx[z][x] > best))
493 				{continue;}
494 			for (int y = 0; y < N; y++) {
495 				mask_y = 1ull << y;
496 				if ((proj_zy[z] | proj_xy[x]) & mask_y)
497 					{continue;}
498 				// getVEmpty, because we know that the cell is
499 				// not attacked, thus empty.
500 				v_zxy = getVEmpty(z, x, y);
501 				if (v_zxy < best) {
502 				    	best = v_zxy;
503 					bs[0][0] = z;
504 					bs[0][1] = x;
505 					bs[0][2] = y;
506 				    	nbs = 1;
507 					if (debug)
508 						printf("%d %d %d is a fitting position\n", z, x, y);
509 				}
510 				else if (v_zxy == best) {
511 				    	bs[nbs][0] = z;
512 				    	bs[nbs][1] = x;
513 				    	bs[nbs][2] = y;
514 				    	nbs++;
515 				    	if (debug)
516 						printf("%d %d %d is a fitting position\n", z, x, y);
517 				}
518 			}
519 		}
520 	}
521 	// choose random cell among optimal ones
522 	int npos = rand() % nbs;
523 	int bz = bs[npos][0];
524 	int bx = bs[npos][1];
525 	int by = bs[npos][2];
526 	if(debug) printf("position %d to %d %d %d with viol %d\n", i, bx, by, bz, getV(bz, by, bx));
527 	G[bz][bx][by] = 1;
528 	R[i][0] = bz;
529 	R[i][1] = bx;
530 	R[i][2] = by;
532 }
533
534 int main(int argc, char **argv) {
535   int oK;
536   int stopK;
537   int first = 1;
538   int seed;
539   N = atoi(argv[1]);
540   oK = atoi(argv[2]);
541   stopK = atoi(argv[3]);
542   if (argc == 5)
543     seed = atoi(argv[4]);
544   else {
545     seed = time(NULL);
546   }
547   printf("using random seed %d\n", seed);
548   srand(seed);
549   signal(SIGINT, intHandler);
550
551   int tot, otot;
552   for (K = oK; K < stopK; K++) {
553     printf("try %d rooks\n", K);
554     // starting position if we just started, otherwise add the rook to the
555     // existing position
556     if (first) {
557       init_pos();
558       first = 0;
559       num_attacks();
560       tot = num_violations();
561     } else {
562       tot += position(K-1);
563     }
564
565     int useless = 0;
566     int prev_rook = -1;
567
568     // improve while violations remain
569     while (tot > 0) {
570       if (debug) print_violations();
571
572       // choose a rook
573       int rook = choose_rook(prev_rook);
574       prev_rook = rook;
575
576       // when doing a useless step, offset the choice
577       //rook = (rook + useless1) % K;
578       if(debug) printf("chosen rook %d\n", rook);
579       if(debug) print();
580       if(debug) print_attacks();
581       int z = R[rook][0];
582       int x = R[rook][1];
583       int y = R[rook][2];
584
585       // remove the rook
586       //printf("tot is %d will rm %d\n", tot, rook);
587       G[z][x][y] = 0;
588       otot = tot;
589       tot += upd_atk_rm(rook);
590       //printf("tot is %d after rm\n", tot);
591       R[rook][0] = MAXN;
592       R[rook][1] = MAXN;
593       R[rook][2] = MAXN;
594
595       if(debug) print_attacks();
596       //check_increment();
597
598       // recompute V without the rook
599       //num_violations();
600       //printf("check_incr\n");
601       //check_increment(tot);
602
603       // choose next position
604       tot += position(rook);
605
606       int pred = V[R[rook][0]][R[rook][1]][R[rook][2]];
607       if(debug) print();
608       if(debug) print_all_violations();
609
610       // recompute number of violations
611       //tot = num_violations();
612
613       if(debug) print_all_violations();
614       int npred = V[R[rook][0]][R[rook][1]][R[rook][2]];
615       if (npred != pred) {
616         printf("violation mispredict: thought %d and had %d!\n", pred, npred);
617         exit(1);
618       }
619       if(debug) print_violations();
620       if(debug) print_attacks();
621       if(debug) print();
622       //printf("total viol %d\n", tot);
623
624       if (tot >= otot) {
625         // no improvement at this round
626         useless++;
627       } else {
628         useless = 0;
629       }
630       if (useless > UBOUND) {
631         // too many useless steps, do a random move
632         useless = 0;
633         do {
634           tot += random_move();
635         } while ((rand() % C2) < C1);
636       }
637     }
638     // we have found a configuration for K rooks; display it
639     printf("!!! DONE %d rooks\n", K);
640     print();
641     fflush(stdout);
642   }
643   return 0;
644 }
645
646
```