crooks

main.c (18368B)

---

 #include <signal.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
 #include <inttypes.h>
 
 // compilation:
 //   gcc -std=c99 -O2 -Wall main.c
 // usage:
 //   ./a.out N K [SEED]
 //   SEED is used to seed the RNG, otherwise time is used;
 //   used seed is always printed first
 //   N is grid size, K is number of rooks (program will increase K when it
 //   succeeded)
 
 // general idea: always take one of the rooks with most violations and move to a
 // cell with least violations; when stuck for UBOUND steps, do a totally random
 // move
 // in all intermediate configurations there are never two rooks in direct
 // attack, the program tries to minimize the number of cells attacked by wrong
 // number of rooks
 // TODO: optimize
 // TODO: profile
 // TODO: plot asymptotics to make conjectures
 // TODO: optimize UBOUND
 
 // upper bound on grid size (N) for static allocation
 #define MAXN 64
 // number of useless steps to make before doing a random move
 #define UBOUND 10
 
 // C1/C2 is the probability of going on doing an additional random move
 #define C1 3
 #define C2 5
 
 const int debug = 0;
 
 // N: grid size
 // K: number of rooks
 int N, K;
 // stores if a cell is taken by a rook or not
 // coordinates are z, x, y
 // corresponding to height, line, col
 unsigned char G[MAXN][MAXN][MAXN];
 // stores a number indicating by which directions the cell is attacked
 // same coordinates
 // numbers are in [0, 8[:, LSB is z, 2nd LSB is x, 3rd LSB is y
 unsigned char A[MAXN][MAXN][MAXN];
 // stores the coordinates of each rook
 // 0 is z (height), 1 is x (line), 2 is y (col)
 unsigned char R[MAXN*MAXN][3];
 // stores the number of violations (number of attacked cells attacked by two
 // other rooks) for each cell (no matter whether it contains a rook or not)
 unsigned char V[MAXN][MAXN][MAXN];
 
 // proj_xy[i]&(1<<j) is true, iff there is a rook of coordinates
 // (x=i, y=j, z=anything)
 // In other words, it is the projection along the z axis giving a matrice of
 // booleans, with the booleans then grouped in bitfields along the y axis.
 // Similarily for the others.
 uint64_t proj_xy[MAXN];
 uint64_t proj_xz[MAXN];
 uint64_t proj_yx[MAXN];
 uint64_t proj_yz[MAXN];
 uint64_t proj_zx[MAXN];
 uint64_t proj_zy[MAXN];
 // v_zx[k][i] is the number of cases already attacked in both the z and the
 // x direction, in the line with coordinates (z=k, x=i).
 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
 // V[z][x][y] = v_zx[z][x] + v_xy[x][y] + v_zy[z][y]
 uint8_t v_zx[MAXN][MAXN];
 uint8_t v_xy[MAXN][MAXN];
 uint8_t v_zy[MAXN][MAXN];
 
 void print();
 
 int popcount64(uint64_t bitfield) {
 	return __builtin_popcountll((unsigned long long) bitfield);
 }
 
 int getVEmpty(int z, int x, int y) {
 	int result;
 	result = (v_zx[z][x] + v_xy[x][y] + v_zy[z][y]);
 	return result;
 }
 
 // TODO: optimise in position where we know that !G[z][x][y]
 int getV(int z, int x, int y) {
 	int result;
 	result = getVEmpty(z, x, y);
 	if (G[z][x][y])
 		result -= 3; // Do not count a rook as being in violation with itself.
 	return result;
 }
 
 void intHandler() {
   // needed for profiling
   printf("SIGINT caught, exiting\n");
   exit(1);
 }
 
 void print() {
   // print the current grid and also rook positions when debug
   for (int x = 0; x < N; x++) {
     for (int y = 0; y < N; y++) {
       int pos = -1;
       for (int z = 0; z < N; z++) {
         if (G[z][x][y]) {
           pos = z;
           break;
         }
       }
       if (pos < 0)
         printf("* ");
       else {
         if (pos < 9)
           printf("%d ", pos + 1);
         else
           printf("%c ", 'A' + (pos - 9));
       }
     }
     printf("\n");
   }
   if (debug) {
     for (int i = 0; i < K; i++)
     printf("rook %d at %d %d %d\n", i, R[i][1], R[i][2], R[i][0]);
   }
 }
 
 void print_attacks() {
   // print the attack pattern in each cell
   for (int z = 0; z < N; z++) {
     printf("== %d ==\n", z);
     for (int x = 0; x < N; x++) {
       for (int y = 0; y < N; y++) {
         int val = A[z][x][y];
         printf("%d%d%d ", val & (1 << 0), (val & (1 << 1))>>1, (val& (1 << 2))>>2);
       }
       printf("\n");
     }
     printf("\n");
   }
 }
 
 void print_all_violations() {
   // print the number of violation in each cell
   for (int z = 0; z < N; z++) {
     printf("== %d ==\n", z);
     for (int x = 0; x < N; x++) {
       for (int y = 0; y < N; y++) {
         int val = V[z][x][y];
         printf("%d ", val);
       }
       printf("\n");
     }
     printf("\n");
   }
 }
 
 void print_violations() {
   // print the number of violations of the rooks
   int tot = 0;
   for (int i = 0; i < K; i++) {
     int val = V[R[i][0]][R[i][1]][R[i][2]];
     printf("rook %d has %d violations\n", i,
         val);
     tot += val;
   }
   printf("=> %d total violations\n", tot);
 }
 
 int upd_atk_rm(int i) {
   // update A to reflect the fact that rook i is about to be removed
   int z = R[i][0];
   int x = R[i][1];
   int y = R[i][2];
   uint64_t mask_z = 1ull << z;
   uint64_t mask_x = 1ull << x;
   uint64_t mask_y = 1ull << y;
   uint64_t old_vzx1, old_vzx2, old_vzy1, old_vzy2, old_vxy1, old_vxy2;
   proj_xy[x] &= ~mask_y;
   proj_xz[x] &= ~mask_z;
   proj_yx[y] &= ~mask_x;
   proj_yz[y] &= ~mask_z;
   proj_zx[z] &= ~mask_x;
   proj_zy[z] &= ~mask_y;
 
 //printf("rmv rook at %d %d %d\n", z, x, y);
   int new_delta = 0;
   for (int p = 0; p < N; p++) {
     A[p][x][y] &= ~(1 << 0); // p x y attacked by z x y along direction z
     A[z][p][y] &= ~(1 << 1); // z p y attacked by z x y along direction x
     A[z][x][p] &= ~(1 << 2); // z x p attacked by z x y along direction y
 
 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
 // 			   v_zy[z][y] = popcount(proj_yx[y] & proj_zx[z])
 // 			   v_xy[x][y] = popcount(proj_yz[y] & proj_xz[x]);
     // TODO: refactor
     old_vzx1 = v_zx[z][p];
     old_vzx2 = v_zx[p][x];
     old_vzy1 = v_zy[z][p];
     old_vzy2 = v_zy[p][y];
     old_vxy1 = v_xy[x][p];
     old_vxy2 = v_xy[p][y];
     v_zx[z][p] = popcount64(proj_xy[p] & proj_zy[z]);
     v_zx[p][x] = popcount64(proj_xy[x] & proj_zy[p]);
     v_zy[z][p] = popcount64(proj_yx[p] & proj_zx[z]);
     v_zy[p][y] = popcount64(proj_yx[y] & proj_zx[p]);
     v_xy[x][p] = popcount64(proj_yz[p] & proj_xz[x]);
     v_xy[p][y] = popcount64(proj_yz[y] & proj_xz[p]);
     if (proj_xz[p] & mask_z)
 	    new_delta += v_zx[z][p] - old_vzx1;
     if (proj_zx[p] & mask_x)
 	    new_delta += v_zx[p][x] - old_vzx2;
     if (proj_yz[p] & mask_z)
 	    new_delta += v_zy[z][p] - old_vzy1;
     if (proj_zy[p] & mask_y)
 	    new_delta += v_zy[p][y] - old_vzy2;
     if (proj_yx[p] & mask_x)
 	    new_delta += v_xy[x][p] - old_vxy1;
     if (proj_xy[p] & mask_y)
 	    new_delta += v_xy[p][y] - old_vxy2;
   }
   new_delta -= getV(z, x, y);
   return new_delta;
 }
 
 int upd_atk_add(int i) {
   // update A to reflect the fact that rook i was just added
   int z = R[i][0];
   int x = R[i][1];
   int y = R[i][2];
   uint64_t mask_z = 1ull << z;
   uint64_t mask_x = 1ull << x;
   uint64_t mask_y = 1ull << y;
   uint64_t old_vzx1, old_vzx2, old_vzy1, old_vzy2, old_vxy1, old_vxy2;
   proj_xy[x] |= mask_y;
   proj_xz[x] |= mask_z;
   proj_yx[y] |= mask_x;
   proj_yz[y] |= mask_z;
   proj_zx[z] |= mask_x;
   proj_zy[z] |= mask_y;
   //printf("added rook at %d %d %d\n", z, x, y);
   int new_delta = 0;
   for (int p = 0; p < N; p++) {
     A[p][x][y] |= 1 << 0;
     A[z][p][y] |= 1 << 1;
     A[z][x][p] |= 1 << 2;
 
 // Interesting invariants: v_zx[z][x] = popcount(proj_xy[x] & proj_zy[z])
 // 			   v_zy[z][y] = popcount(proj_yx[y] & proj_zx[z])
 // 			   v_xy[x][y] = popcount(proj_yz[y] & proj_xz[x]);
     // TODO: refactor
     old_vzx1 = v_zx[z][p];
     old_vzx2 = v_zx[p][x];
     old_vzy1 = v_zy[z][p];
     old_vzy2 = v_zy[p][y];
     old_vxy1 = v_xy[x][p];
     old_vxy2 = v_xy[p][y];
     v_zx[z][p] = popcount64(proj_xy[p] & proj_zy[z]);
     v_zx[p][x] = popcount64(proj_xy[x] & proj_zy[p]);
     v_zy[z][p] = popcount64(proj_yx[p] & proj_zx[z]);
     v_zy[p][y] = popcount64(proj_yx[y] & proj_zx[p]);
     v_xy[x][p] = popcount64(proj_yz[p] & proj_xz[x]);
     v_xy[p][y] = popcount64(proj_yz[y] & proj_xz[p]);
     if (p != z) {
 	    if (proj_zx[p] & mask_x)
 		    new_delta += v_zx[p][x] - old_vzx2;
 	    if (proj_zy[p] & mask_y)
 		    new_delta += v_zy[p][y] - old_vzy2;
     }
     if (p != x) {
 	    if (proj_xz[p] & mask_z)
 		    new_delta += v_zx[z][p] - old_vzx1;
 	    if (proj_xy[p] & mask_y)
 		    new_delta += v_xy[p][y] - old_vxy2;
     }
     if (p != y) {
 	    if (proj_yz[p] & mask_z)
 		    new_delta += v_zy[z][p] - old_vzy1;
 	    if (proj_yx[p] & mask_x)
 		    new_delta += v_xy[x][p] - old_vxy1;
     }
   }
   new_delta += getV(z, x, y);
   return new_delta;
 }
 
 void num_attacks() {
   // generate A from scratch
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++)
         A[i][j][k] = 0;
   for (int i = 0; i < K; i++) {
     if (R[i][0] < MAXN)
       upd_atk_add(i);
   }
 }
 
 void check_attacks() {
   // debug: regenerate A from scratch
   // and check that the previous version was already correct
   int backup[MAXN][MAXN][MAXN];
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++)
         backup[i][j][k] = A[i][j][k];
   num_attacks();
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++)
         if (backup[i][j][k] != A[i][j][k]) {
           print_attacks();
           printf("ERROR at %d %d %d\n", i, j, k);
           exit(1);
         }
 }
 
 int num_violations() {
   // regenerate V from scratch
   // returns total number of violations
   // TODO: find a way to update this rather than recompute
   int tot = 0;
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++)
         V[i][j][k] = 0; // clear
   for (int z = 0; z < N; z++)
     for (int x = 0; x < N; x++)
       for (int y = 0; y < N; y++) {
         for (int p = 0; p < N; p++) {
           // must count cells that have already two attack directions and are
           // missing the attack direction from the current cell
           if (p != z && (A[p][x][y] | (1 << 0)) == 7)
           {V[z][x][y]++;} // putting rook in z x y would attack p x y thrice
           if (p != x && (A[z][p][y] | (1 << 1)) == 7)
           {V[z][x][y]++;} // putting rook in z x y would attack z p y thrice
           if (p != y && (A[z][x][p] | (1 << 2)) == 7)
           {V[z][x][y]++;} // putting rook in z x y would attack z x p thrice
         }
         if (G[z][x][y]) {
           tot += V[z][x][y];
           //printf("add %d to num_violations for %d %d %d\n", V[z][x][y], x, y, z);
         }
       }
   return tot;
 }
 
 void check_increment(int ntot) {
   //printf("CHECK: before\n");
   //print_all_violations();
   // debug: regenerate A, V from scratch
   // and check that the previous version was already correct
   int backupA[MAXN][MAXN][MAXN];
   int backupV[MAXN][MAXN][MAXN];
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++) {
         backupA[i][j][k] = A[i][j][k];
         backupV[i][j][k] = V[i][j][k];
       }
   num_attacks();
   int rntot = num_violations();
   //printf("CHECK: after\n");
   //print_all_violations();
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++) {
         if (backupA[i][j][k] != A[i][j][k]) {
           print_attacks();
           printf("ERROR attacks at %d %d %d\n", i, j, k);
           exit(1);
         }
         if (backupV[i][j][k] != V[i][j][k]) {
           print();
           print_attacks();
           print_all_violations();
           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]);
           exit(1);
         }
       }
   if (ntot != rntot) {
     print_all_violations();
     print_violations();
     printf("total violation miscount, had %d computed %d\n", ntot, rntot);
     exit(1);
   }
 }
 
 
 void init_pos() {
   // create initial position
   for (int i = 0; i < N; i++)
     for (int j = 0; j < N; j++)
       for (int k = 0; k < N; k++)
         G[i][j][k] = 0;
   int remaining = K;
   int height = 0, pos = 0;
   while (remaining > 0) {
     int z = height;
     int x = pos;
     int y = (pos + height) % N;
     G[z][x][y] = 1;
     remaining--;
     R[remaining][0] = z;
     R[remaining][1] = x;
     R[remaining][2] = y;
     pos++;
     if (pos >= N) {
       pos = 0;
       height++;
     }
   }
 }
 
 int choose_rook(int prev) {
   // choose the next rook to move: a random one among the ones with maximal
   // violations, forbiding prev (the previously moved rook)
   // TODO: better notion including "age" of rook to favor moving old rooks (?)
   // and small probability to move rooks with not exactly maximal num of viols
   // TODO: make it faster (priority queue?)
   int best = -1;
   int nvs[MAXN*MAXN];
   int nnvs = 0;
   int nv;
   for (int i = 0; i < K; i++) {
     if (i == prev)
       continue;
     if ((nv = getV(R[i][0], R[i][1], R[i][2])) > best) {
       best = nv;
       nvs[0] = i;
       nnvs = 1;
     } else if (nv == best) {
       nvs[nnvs++] = i;
     }
   }
   // choose random rook among the optimal ones
   return nvs[rand() % nnvs];
 }
 
 int random_move() {
   // make a totally random move (when we are stuck)
   int rook = (rand() % K);
   //printf("move %d at random\n", rook);
   int z = R[rook][0];
   int x = R[rook][1];
   int y = R[rook][2];
   int delta = 0;
   G[z][x][y] = 0;
   delta += upd_atk_rm(rook);
   int nx, ny, nz;
   uint64_t mask_x, mask_y;
   do {
     nx = (rand() % N);
     ny = (rand() % N);
     nz = (rand() % N);
     mask_x = 1ull << x;
     mask_y = 1ull << y;
     if (!A[nz][nx][ny] && !G[nz][nx][ny]) {
 #if 0
 	    printf("proj_zy[nz]:%lx, mask_y: %lx:%lx -- ", proj_zy[nz], mask_y, proj_zy[nz] & mask_y);
 	    printf("proj_xy[nx]:%lx, mask_y: %lx:%lx -- ", proj_xy[nx], mask_y, proj_xy[nx] & mask_y);
 	    printf("proj_zx[nz]:%lx, mask_x: %lx:%lx\n", proj_zx[nz], mask_x, proj_zx[nz] & mask_x);
 #endif
       break; // found a spot that's both empty and not in direct attack
     }
     // TODO: nothing tests if this loop stands a chance of terminating, it would
     // be neater to exit the program then
   } while (1); // ((proj_zy[nz] & mask_y) || (proj_xy[nx] & mask_y) || (proj_zx[nz] & mask_x));
   G[nz][nx][ny] = 1;
   R[rook][0] = nz;
   R[rook][1] = nx;
   R[rook][2] = ny;
   delta += upd_atk_add(rook);
   return delta;
 }
 
 int position(int i) {
 	// position rook i that was just removed
 	// TODO: more efficient (pqueue?)
 	// TODO more clever (chance to select cell with non-optimal number of
 	// violations)
 	int best = N*N*N+1;
 	int bs[MAXN*MAXN*MAXN][3];
 	int nbs = 0;
 	int v_zxy;
 	uint64_t mask_x, mask_y;
 	for (int z = 0; z < N; z++) {
 		for (int x = 0; x < N; x++) {
 			mask_x = 1ull << x;
 			if((proj_zx[z] & mask_x) || (v_zx[z][x] > best))
 				{continue;}
 			for (int y = 0; y < N; y++) {
 				mask_y = 1ull << y;
 				if ((proj_zy[z] | proj_xy[x]) & mask_y)
 					{continue;}
 				// getVEmpty, because we know that the cell is
 				// not attacked, thus empty.
 				v_zxy = getVEmpty(z, x, y);
 				if (v_zxy < best) {
 				    	best = v_zxy;
 					bs[0][0] = z;
 					bs[0][1] = x;
 					bs[0][2] = y;
 				    	nbs = 1;
 					if (debug)
 						printf("%d %d %d is a fitting position\n", z, x, y);
 				}
 				else if (v_zxy == best) {
 				    	bs[nbs][0] = z;
 				    	bs[nbs][1] = x;
 				    	bs[nbs][2] = y;
 				    	nbs++;
 				    	if (debug)
 						printf("%d %d %d is a fitting position\n", z, x, y);
 				}
 			}
 		}
 	}
 	// choose random cell among optimal ones
 	int npos = rand() % nbs;
 	int bz = bs[npos][0];
 	int bx = bs[npos][1];
 	int by = bs[npos][2];
 	if(debug) printf("position %d to %d %d %d with viol %d\n", i, bx, by, bz, getV(bz, by, bx));
 	G[bz][bx][by] = 1;
 	R[i][0] = bz;
 	R[i][1] = bx;
 	R[i][2] = by;
 	return upd_atk_add(i);
 }
 
 int main(int argc, char **argv) {
   int oK;
   int stopK;
   int first = 1;
   int seed;
   N = atoi(argv[1]);
   oK = atoi(argv[2]);
   stopK = atoi(argv[3]);
   if (argc == 5)
     seed = atoi(argv[4]);
   else {
     seed = time(NULL);
   }
   printf("using random seed %d\n", seed);
   srand(seed);
   signal(SIGINT, intHandler);
 
   int tot, otot;
   for (K = oK; K < stopK; K++) {
     printf("try %d rooks\n", K);
     // starting position if we just started, otherwise add the rook to the
     // existing position
     if (first) {
       init_pos();
       first = 0;
       num_attacks();
       tot = num_violations();
     } else {
       tot += position(K-1);
     }
 
     int useless = 0;
     int prev_rook = -1;
 
     // improve while violations remain
     while (tot > 0) {
       if (debug) print_violations();
 
       // choose a rook
       int rook = choose_rook(prev_rook);
       prev_rook = rook;
 
       // when doing a useless step, offset the choice
       //rook = (rook + useless1) % K;
       if(debug) printf("chosen rook %d\n", rook);
       if(debug) print();
       if(debug) print_attacks();
       int z = R[rook][0];
       int x = R[rook][1];
       int y = R[rook][2];
 
       // remove the rook
       //printf("tot is %d will rm %d\n", tot, rook);
       G[z][x][y] = 0;
       otot = tot;
       tot += upd_atk_rm(rook);
       //printf("tot is %d after rm\n", tot);
       R[rook][0] = MAXN;
       R[rook][1] = MAXN;
       R[rook][2] = MAXN;
 
       if(debug) print_attacks();
       //check_increment();
 
       // recompute V without the rook
       //num_violations();
       //printf("check_incr\n");
       //check_increment(tot);
 
       // choose next position
       tot += position(rook);
 
       int pred = V[R[rook][0]][R[rook][1]][R[rook][2]];
       if(debug) print();
       if(debug) print_all_violations();
 
       // recompute number of violations
       //tot = num_violations();
 
       if(debug) print_all_violations();
       int npred = V[R[rook][0]][R[rook][1]][R[rook][2]];
       if (npred != pred) {
         printf("violation mispredict: thought %d and had %d!\n", pred, npred);
         exit(1);
       }
       if(debug) print_violations();
       if(debug) print_attacks();
       if(debug) print();
       //printf("total viol %d\n", tot);
 
       if (tot >= otot) {
         // no improvement at this round
         useless++;
       } else {
         useless = 0;
       }
       if (useless > UBOUND) {
         // too many useless steps, do a random move
         useless = 0;
         do {
           tot += random_move();
         } while ((rand() % C2) < C1);
       }
     }
     // we have found a configuration for K rooks; display it
     printf("!!! DONE %d rooks\n", K);
     print();
     fflush(stdout);
   }
   return 0;
 }