latin.c

   1 #include <assert.h>
   2 #include <string.h>
   3 #include <stdarg.h>
   4
   5 #include "puzzles.h"
   6 #include "tree234.h"
   7 #include "maxflow.h"
   8
   9 #ifdef STANDALONE_LATIN_TEST
  10 #define STANDALONE_SOLVER
  11 #endif
  12
  13 #include "latin.h"
  14
  15 /* --------------------------------------------------------
  16  * Solver.
  17  */
  18
  19 /*
  20  * Function called when we are certain that a particular square has
  21  * a particular number in it. The y-coordinate passed in here is
  22  * transformed.
  23  */
  24 void latin_solver_place(struct latin_solver *solver, int x, int y, int n)
  25 {
  26     int i, o = solver->o;
  27
  28     assert(n <= o);
  29     assert(cube(x,y,n));
  30
  31     /*
  32      * Rule out all other numbers in this square.
  33      */
  34     for (i = 1; i <= o; i++)
  35         if (i != n)
  36             cube(x,y,i) = FALSE;
  37
  38     /*
  39      * Rule out this number in all other positions in the row.
  40      */
  41     for (i = 0; i < o; i++)
  42         if (i != y)
  43             cube(x,i,n) = FALSE;
  44
  45     /*
  46      * Rule out this number in all other positions in the column.
  47      */
  48     for (i = 0; i < o; i++)
  49         if (i != x)
  50             cube(i,y,n) = FALSE;
  51
  52     /*
  53      * Enter the number in the result grid.
  54      */
  55     solver->grid[YUNTRANS(y)*o+x] = n;
  56
  57     /*
  58      * Cross out this number from the list of numbers left to place
  59      * in its row, its column and its block.
  60      */
  61     solver->row[y*o+n-1] = solver->col[x*o+n-1] = TRUE;
  62 }
  63
  64 int latin_solver_elim(struct latin_solver *solver, int start, int step
  65 #ifdef STANDALONE_SOLVER
  66                       , char *fmt, ...
  67 #endif
  68                       )
  69 {
  70     int o = solver->o;
  71     int fpos, m, i;
  72
  73     /*
  74      * Count the number of set bits within this section of the
  75      * cube.
  76      */
  77     m = 0;
  78     fpos = -1;
  79     for (i = 0; i < o; i++)
  80         if (solver->cube[start+i*step]) {
  81             fpos = start+i*step;
  82             m++;
  83         }
  84
  85     if (m == 1) {
  86         int x, y, n;
  87         assert(fpos >= 0);
  88
  89         n = 1 + fpos % o;
  90         y = fpos / o;
  91         x = y / o;
  92         y %= o;
  93
  94         if (!solver->grid[YUNTRANS(y)*o+x]) {
  95 #ifdef STANDALONE_SOLVER
  96             if (solver_show_working) {
  97                 va_list ap;
  98                 printf("%*s", solver_recurse_depth*4, "");
  99                 va_start(ap, fmt);
 100                 vprintf(fmt, ap);
 101                 va_end(ap);
 102                 printf(":\n%*s  placing %d at (%d,%d)\n",
 103                        solver_recurse_depth*4, "", n, x, YUNTRANS(y));
 104             }
 105 #endif
 106             latin_solver_place(solver, x, y, n);
 107             return +1;
 108         }
 109     } else if (m == 0) {
 110 #ifdef STANDALONE_SOLVER
 111         if (solver_show_working) {
 112             va_list ap;
 113             printf("%*s", solver_recurse_depth*4, "");
 114             va_start(ap, fmt);
 115             vprintf(fmt, ap);
 116             va_end(ap);
 117             printf(":\n%*s  no possibilities available\n",
 118                    solver_recurse_depth*4, "");
 119         }
 120 #endif
 121         return -1;
 122     }
 123
 124     return 0;
 125 }
 126
 127 struct latin_solver_scratch {
 128     unsigned char *grid, *rowidx, *colidx, *set;
 129     int *neighbours, *bfsqueue;
 130 #ifdef STANDALONE_SOLVER
 131     int *bfsprev;
 132 #endif
 133 };
 134
 135 int latin_solver_set(struct latin_solver *solver,
 136                      struct latin_solver_scratch *scratch,
 137                      int start, int step1, int step2
 138 #ifdef STANDALONE_SOLVER
 139                      , char *fmt, ...
 140 #endif
 141                      )
 142 {
 143     int o = solver->o;
 144     int i, j, n, count;
 145     unsigned char *grid = scratch->grid;
 146     unsigned char *rowidx = scratch->rowidx;
 147     unsigned char *colidx = scratch->colidx;
 148     unsigned char *set = scratch->set;
 149
 150     /*
 151      * We are passed a o-by-o matrix of booleans. Our first job
 152      * is to winnow it by finding any definite placements - i.e.
 153      * any row with a solitary 1 - and discarding that row and the
 154      * column containing the 1.
 155      */
 156     memset(rowidx, TRUE, o);
 157     memset(colidx, TRUE, o);
 158     for (i = 0; i < o; i++) {
 159         int count = 0, first = -1;
 160         for (j = 0; j < o; j++)
 161             if (solver->cube[start+i*step1+j*step2])
 162                 first = j, count++;
 163
 164         if (count == 0) return -1;
 165         if (count == 1)
 166             rowidx[i] = colidx[first] = FALSE;
 167     }
 168
 169     /*
 170      * Convert each of rowidx/colidx from a list of 0s and 1s to a
 171      * list of the indices of the 1s.
 172      */
 173     for (i = j = 0; i < o; i++)
 174         if (rowidx[i])
 175             rowidx[j++] = i;
 176     n = j;
 177     for (i = j = 0; i < o; i++)
 178         if (colidx[i])
 179             colidx[j++] = i;
 180     assert(n == j);
 181
 182     /*
 183      * And create the smaller matrix.
 184      */
 185     for (i = 0; i < n; i++)
 186         for (j = 0; j < n; j++)
 187             grid[i*o+j] = solver->cube[start+rowidx[i]*step1+colidx[j]*step2];
 188
 189     /*
 190      * Having done that, we now have a matrix in which every row
 191      * has at least two 1s in. Now we search to see if we can find
 192      * a rectangle of zeroes (in the set-theoretic sense of
 193      * `rectangle', i.e. a subset of rows crossed with a subset of
 194      * columns) whose width and height add up to n.
 195      */
 196
 197     memset(set, 0, n);
 198     count = 0;
 199     while (1) {
 200         /*
 201          * We have a candidate set. If its size is <=1 or >=n-1
 202          * then we move on immediately.
 203          */
 204         if (count > 1 && count < n-1) {
 205             /*
 206              * The number of rows we need is n-count. See if we can
 207              * find that many rows which each have a zero in all
 208              * the positions listed in `set'.
 209              */
 210             int rows = 0;
 211             for (i = 0; i < n; i++) {
 212                 int ok = TRUE;
 213                 for (j = 0; j < n; j++)
 214                     if (set[j] && grid[i*o+j]) {
 215                         ok = FALSE;
 216                         break;
 217                     }
 218                 if (ok)
 219                     rows++;
 220             }
 221
 222             /*
 223              * We expect never to be able to get _more_ than
 224              * n-count suitable rows: this would imply that (for
 225              * example) there are four numbers which between them
 226              * have at most three possible positions, and hence it
 227              * indicates a faulty deduction before this point or
 228              * even a bogus clue.
 229              */
 230             if (rows > n - count) {
 231 #ifdef STANDALONE_SOLVER
 232                 if (solver_show_working) {
 233                     va_list ap;
 234                     printf("%*s", solver_recurse_depth*4,
 235                            "");
 236                     va_start(ap, fmt);
 237                     vprintf(fmt, ap);
 238                     va_end(ap);
 239                     printf(":\n%*s  contradiction reached\n",
 240                            solver_recurse_depth*4, "");
 241                 }
 242 #endif
 243                 return -1;
 244             }
 245
 246             if (rows >= n - count) {
 247                 int progress = FALSE;
 248
 249                 /*
 250                  * We've got one! Now, for each row which _doesn't_
 251                  * satisfy the criterion, eliminate all its set
 252                  * bits in the positions _not_ listed in `set'.
 253                  * Return +1 (meaning progress has been made) if we
 254                  * successfully eliminated anything at all.
 255                  *
 256                  * This involves referring back through
 257                  * rowidx/colidx in order to work out which actual
 258                  * positions in the cube to meddle with.
 259                  */
 260                 for (i = 0; i < n; i++) {
 261                     int ok = TRUE;
 262                     for (j = 0; j < n; j++)
 263                         if (set[j] && grid[i*o+j]) {
 264                             ok = FALSE;
 265                             break;
 266                         }
 267                     if (!ok) {
 268                         for (j = 0; j < n; j++)
 269                             if (!set[j] && grid[i*o+j]) {
 270                                 int fpos = (start+rowidx[i]*step1+
 271                                             colidx[j]*step2);
 272 #ifdef STANDALONE_SOLVER
 273                                 if (solver_show_working) {
 274                                     int px, py, pn;
 275
 276                                     if (!progress) {
 277                                         va_list ap;
 278                                         printf("%*s", solver_recurse_depth*4,
 279                                                "");
 280                                         va_start(ap, fmt);
 281                                         vprintf(fmt, ap);
 282                                         va_end(ap);
 283                                         printf(":\n");
 284                                     }
 285
 286                                     pn = 1 + fpos % o;
 287                                     py = fpos / o;
 288                                     px = py / o;
 289                                     py %= o;
 290
 291                                     printf("%*s  ruling out %d at (%d,%d)\n",
 292                                            solver_recurse_depth*4, "",
 293                                            pn, px, YUNTRANS(py));
 294                                 }
 295 #endif
 296                                 progress = TRUE;
 297                                 solver->cube[fpos] = FALSE;
 298                             }
 299                     }
 300                 }
 301
 302                 if (progress) {
 303                     return +1;
 304                 }
 305             }
 306         }
 307
 308         /*
 309          * Binary increment: change the rightmost 0 to a 1, and
 310          * change all 1s to the right of it to 0s.
 311          */
 312         i = n;
 313         while (i > 0 && set[i-1])
 314             set[--i] = 0, count--;
 315         if (i > 0)
 316             set[--i] = 1, count++;
 317         else
 318             break;                     /* done */
 319     }
 320
 321     return 0;
 322 }
 323
 324 /*
 325  * Look for forcing chains. A forcing chain is a path of
 326  * pairwise-exclusive squares (i.e. each pair of adjacent squares
 327  * in the path are in the same row, column or block) with the
 328  * following properties:
 329  *
 330  *  (a) Each square on the path has precisely two possible numbers.
 331  *
 332  *  (b) Each pair of squares which are adjacent on the path share
 333  *      at least one possible number in common.
 334  *
 335  *  (c) Each square in the middle of the path shares _both_ of its
 336  *      numbers with at least one of its neighbours (not the same
 337  *      one with both neighbours).
 338  *
 339  * These together imply that at least one of the possible number
 340  * choices at one end of the path forces _all_ the rest of the
 341  * numbers along the path. In order to make real use of this, we
 342  * need further properties:
 343  *
 344  *  (c) Ruling out some number N from the square at one end
 345  *      of the path forces the square at the other end to
 346  *      take number N.
 347  *
 348  *  (d) The two end squares are both in line with some third
 349  *      square.
 350  *
 351  *  (e) That third square currently has N as a possibility.
 352  *
 353  * If we can find all of that lot, we can deduce that at least one
 354  * of the two ends of the forcing chain has number N, and that
 355  * therefore the mutually adjacent third square does not.
 356  *
 357  * To find forcing chains, we're going to start a bfs at each
 358  * suitable square, once for each of its two possible numbers.
 359  */
 360 int latin_solver_forcing(struct latin_solver *solver,
 361                          struct latin_solver_scratch *scratch)
 362 {
 363     int o = solver->o;
 364     int *bfsqueue = scratch->bfsqueue;
 365 #ifdef STANDALONE_SOLVER
 366     int *bfsprev = scratch->bfsprev;
 367 #endif
 368     unsigned char *number = scratch->grid;
 369     int *neighbours = scratch->neighbours;
 370     int x, y;
 371
 372     for (y = 0; y < o; y++)
 373         for (x = 0; x < o; x++) {
 374             int count, t, n;
 375
 376             /*
 377              * If this square doesn't have exactly two candidate
 378              * numbers, don't try it.
 379              *
 380              * In this loop we also sum the candidate numbers,
 381              * which is a nasty hack to allow us to quickly find
 382              * `the other one' (since we will shortly know there
 383              * are exactly two).
 384              */
 385             for (count = t = 0, n = 1; n <= o; n++)
 386                 if (cube(x, y, n))
 387                     count++, t += n;
 388             if (count != 2)
 389                 continue;
 390
 391             /*
 392              * Now attempt a bfs for each candidate.
 393              */
 394             for (n = 1; n <= o; n++)
 395                 if (cube(x, y, n)) {
 396                     int orign, currn, head, tail;
 397
 398                     /*
 399                      * Begin a bfs.
 400                      */
 401                     orign = n;
 402
 403                     memset(number, o+1, o*o);
 404                     head = tail = 0;
 405                     bfsqueue[tail++] = y*o+x;
 406 #ifdef STANDALONE_SOLVER
 407                     bfsprev[y*o+x] = -1;
 408 #endif
 409                     number[y*o+x] = t - n;
 410
 411                     while (head < tail) {
 412                         int xx, yy, nneighbours, xt, yt, i;
 413
 414                         xx = bfsqueue[head++];
 415                         yy = xx / o;
 416                         xx %= o;
 417
 418                         currn = number[yy*o+xx];
 419
 420                         /*
 421                          * Find neighbours of yy,xx.
 422                          */
 423                         nneighbours = 0;
 424                         for (yt = 0; yt < o; yt++)
 425                             neighbours[nneighbours++] = yt*o+xx;
 426                         for (xt = 0; xt < o; xt++)
 427                             neighbours[nneighbours++] = yy*o+xt;
 428
 429                         /*
 430                          * Try visiting each of those neighbours.
 431                          */
 432                         for (i = 0; i < nneighbours; i++) {
 433                             int cc, tt, nn;
 434
 435                             xt = neighbours[i] % o;
 436                             yt = neighbours[i] / o;
 437
 438                             /*
 439                              * We need this square to not be
 440                              * already visited, and to include
 441                              * currn as a possible number.
 442                              */
 443                             if (number[yt*o+xt] <= o)
 444                                 continue;
 445                             if (!cube(xt, yt, currn))
 446                                 continue;
 447
 448                             /*
 449                              * Don't visit _this_ square a second
 450                              * time!
 451                              */
 452                             if (xt == xx && yt == yy)
 453                                 continue;
 454
 455                             /*
 456                              * To continue with the bfs, we need
 457                              * this square to have exactly two
 458                              * possible numbers.
 459                              */
 460                             for (cc = tt = 0, nn = 1; nn <= o; nn++)
 461                                 if (cube(xt, yt, nn))
 462                                     cc++, tt += nn;
 463                             if (cc == 2) {
 464                                 bfsqueue[tail++] = yt*o+xt;
 465 #ifdef STANDALONE_SOLVER
 466                                 bfsprev[yt*o+xt] = yy*o+xx;
 467 #endif
 468                                 number[yt*o+xt] = tt - currn;
 469                             }
 470
 471                             /*
 472                              * One other possibility is that this
 473                              * might be the square in which we can
 474                              * make a real deduction: if it's
 475                              * adjacent to x,y, and currn is equal
 476                              * to the original number we ruled out.
 477                              */
 478                             if (currn == orign &&
 479                                 (xt == x || yt == y)) {
 480 #ifdef STANDALONE_SOLVER
 481                                 if (solver_show_working) {
 482                                     char *sep = "";
 483                                     int xl, yl;
 484                                     printf("%*sforcing chain, %d at ends of ",
 485                                            solver_recurse_depth*4, "", orign);
 486                                     xl = xx;
 487                                     yl = yy;
 488                                     while (1) {
 489                                         printf("%s(%d,%d)", sep, xl,
 490                                                YUNTRANS(yl));
 491                                         xl = bfsprev[yl*o+xl];
 492                                         if (xl < 0)
 493                                             break;
 494                                         yl = xl / o;
 495                                         xl %= o;
 496                                         sep = "-";
 497                                     }
 498                                     printf("\n%*s  ruling out %d at (%d,%d)\n",
 499                                            solver_recurse_depth*4, "",
 500                                            orign, xt, YUNTRANS(yt));
 501                                 }
 502 #endif
 503                                 cube(xt, yt, orign) = FALSE;
 504                                 return 1;
 505                             }
 506                         }
 507                     }
 508                 }
 509         }
 510
 511     return 0;
 512 }
 513
 514 struct latin_solver_scratch *latin_solver_new_scratch(struct latin_solver *solver)
 515 {
 516     struct latin_solver_scratch *scratch = snew(struct latin_solver_scratch);
 517     int o = solver->o;
 518     scratch->grid = snewn(o*o, unsigned char);
 519     scratch->rowidx = snewn(o, unsigned char);
 520     scratch->colidx = snewn(o, unsigned char);
 521     scratch->set = snewn(o, unsigned char);
 522     scratch->neighbours = snewn(3*o, int);
 523     scratch->bfsqueue = snewn(o*o, int);
 524 #ifdef STANDALONE_SOLVER
 525     scratch->bfsprev = snewn(o*o, int);
 526 #endif
 527     return scratch;
 528 }
 529
 530 void latin_solver_free_scratch(struct latin_solver_scratch *scratch)
 531 {
 532 #ifdef STANDALONE_SOLVER
 533     sfree(scratch->bfsprev);
 534 #endif
 535     sfree(scratch->bfsqueue);
 536     sfree(scratch->neighbours);
 537     sfree(scratch->set);
 538     sfree(scratch->colidx);
 539     sfree(scratch->rowidx);
 540     sfree(scratch->grid);
 541     sfree(scratch);
 542 }
 543
 544 void latin_solver_alloc(struct latin_solver *solver, digit *grid, int o)
 545 {
 546     int x, y;
 547
 548     solver->o = o;
 549     solver->cube = snewn(o*o*o, unsigned char);
 550     solver->grid = grid;                /* write straight back to the input */
 551     memset(solver->cube, TRUE, o*o*o);
 552
 553     solver->row = snewn(o*o, unsigned char);
 554     solver->col = snewn(o*o, unsigned char);
 555     memset(solver->row, FALSE, o*o);
 556     memset(solver->col, FALSE, o*o);
 557
 558     for (x = 0; x < o; x++)
 559         for (y = 0; y < o; y++)
 560             if (grid[y*o+x])
 561                 latin_solver_place(solver, x, YTRANS(y), grid[y*o+x]);
 562 }
 563
 564 void latin_solver_free(struct latin_solver *solver)
 565 {
 566     sfree(solver->cube);
 567     sfree(solver->row);
 568     sfree(solver->col);
 569 }
 570
 571 int latin_solver_diff_simple(struct latin_solver *solver)
 572 {
 573     int x, y, n, ret, o = solver->o;
 574     /*
 575      * Row-wise positional elimination.
 576      */
 577     for (y = 0; y < o; y++)
 578         for (n = 1; n <= o; n++)
 579             if (!solver->row[y*o+n-1]) {
 580                 ret = latin_solver_elim(solver, cubepos(0,y,n), o*o
 581 #ifdef STANDALONE_SOLVER
 582                                         , "positional elimination,"
 583                                         " %d in row %d", n, YUNTRANS(y)
 584 #endif
 585                                         );
 586                 if (ret != 0) return ret;
 587             }
 588     /*
 589      * Column-wise positional elimination.
 590      */
 591     for (x = 0; x < o; x++)
 592         for (n = 1; n <= o; n++)
 593             if (!solver->col[x*o+n-1]) {
 594                 ret = latin_solver_elim(solver, cubepos(x,0,n), o
 595 #ifdef STANDALONE_SOLVER
 596                                         , "positional elimination,"
 597                                         " %d in column %d", n, x
 598 #endif
 599                                         );
 600                 if (ret != 0) return ret;
 601             }
 602
 603     /*
 604      * Numeric elimination.
 605      */
 606     for (x = 0; x < o; x++)
 607         for (y = 0; y < o; y++)
 608             if (!solver->grid[YUNTRANS(y)*o+x]) {
 609                 ret = latin_solver_elim(solver, cubepos(x,y,1), 1
 610 #ifdef STANDALONE_SOLVER
 611                                         , "numeric elimination at (%d,%d)", x,
 612                                         YUNTRANS(y)
 613 #endif
 614                                         );
 615                 if (ret != 0) return ret;
 616             }
 617     return 0;
 618 }
 619
 620 int latin_solver_diff_set(struct latin_solver *solver,
 621                           struct latin_solver_scratch *scratch,
 622                           int *extreme)
 623 {
 624     int x, y, n, ret, o = solver->o;
 625     /*
 626      * Row-wise set elimination.
 627      */
 628     for (y = 0; y < o; y++) {
 629         ret = latin_solver_set(solver, scratch, cubepos(0,y,1), o*o, 1
 630 #ifdef STANDALONE_SOLVER
 631                                , "set elimination, row %d", YUNTRANS(y)
 632 #endif
 633                                );
 634         if (ret > 0) *extreme = 0;
 635         if (ret != 0) return ret;
 636     }
 637
 638     /*
 639      * Column-wise set elimination.
 640      */
 641     for (x = 0; x < o; x++) {
 642         ret = latin_solver_set(solver, scratch, cubepos(x,0,1), o, 1
 643 #ifdef STANDALONE_SOLVER
 644                                , "set elimination, column %d", x
 645 #endif
 646                                );
 647         if (ret > 0) *extreme = 0;
 648         if (ret != 0) return ret;
 649     }
 650
 651     /*
 652      * Row-vs-column set elimination on a single number.
 653      */
 654     for (n = 1; n <= o; n++) {
 655         ret = latin_solver_set(solver, scratch, cubepos(0,0,n), o*o, o
 656 #ifdef STANDALONE_SOLVER
 657                                , "positional set elimination, number %d", n
 658 #endif
 659                                );
 660         if (ret > 0) *extreme = 1;
 661         if (ret != 0) return ret;
 662     }
 663     return 0;
 664 }
 665
 666 /* This uses our own diff_* internally, but doesn't require callers
 667  * to; this is so it can be used by games that want to rewrite
 668  * the solver so as to use a different set of difficulties.
 669  *
 670  * It returns:
 671  * 0 for 'didn't do anything' implying it was already solved.
 672  * -1 for 'impossible' (no solution)
 673  * 1 for 'single solution'
 674  * >1 for 'multiple solutions' (you don't get to know how many, and
 675  *     the first such solution found will be set.
 676  *
 677  * and this function may well assert if given an impossible board.
 678  */
 679 int latin_solver_recurse(struct latin_solver *solver, int recdiff,
 680                          latin_solver_callback cb, void *ctx)
 681 {
 682     int best, bestcount;
 683     int o = solver->o, x, y, n;
 684
 685     best = -1;
 686     bestcount = o+1;
 687
 688     for (y = 0; y < o; y++)
 689         for (x = 0; x < o; x++)
 690             if (!solver->grid[y*o+x]) {
 691                 int count;
 692
 693                 /*
 694                  * An unfilled square. Count the number of
 695                  * possible digits in it.
 696                  */
 697                 count = 0;
 698                 for (n = 1; n <= o; n++)
 699                     if (cube(x,YTRANS(y),n))
 700                         count++;
 701
 702                 /*
 703                  * We should have found any impossibilities
 704                  * already, so this can safely be an assert.
 705                  */
 706                 assert(count > 1);
 707
 708                 if (count < bestcount) {
 709                     bestcount = count;
 710                     best = y*o+x;
 711                 }
 712             }
 713
 714     if (best == -1)
 715         /* we were complete already. */
 716         return 0;
 717     else {
 718         int i, j;
 719         digit *list, *ingrid, *outgrid;
 720         int diff = diff_impossible;    /* no solution found yet */
 721
 722         /*
 723          * Attempt recursion.
 724          */
 725         y = best / o;
 726         x = best % o;
 727
 728         list = snewn(o, digit);
 729         ingrid = snewn(o*o, digit);
 730         outgrid = snewn(o*o, digit);
 731         memcpy(ingrid, solver->grid, o*o);
 732
 733         /* Make a list of the possible digits. */
 734         for (j = 0, n = 1; n <= o; n++)
 735             if (cube(x,YTRANS(y),n))
 736                 list[j++] = n;
 737
 738 #ifdef STANDALONE_SOLVER
 739         if (solver_show_working) {
 740             char *sep = "";
 741             printf("%*srecursing on (%d,%d) [",
 742                    solver_recurse_depth*4, "", x, y);
 743             for (i = 0; i < j; i++) {
 744                 printf("%s%d", sep, list[i]);
 745                 sep = " or ";
 746             }
 747             printf("]\n");
 748         }
 749 #endif
 750
 751         /*
 752          * And step along the list, recursing back into the
 753          * main solver at every stage.
 754          */
 755         for (i = 0; i < j; i++) {
 756             int ret;
 757
 758             memcpy(outgrid, ingrid, o*o);
 759             outgrid[y*o+x] = list[i];
 760
 761 #ifdef STANDALONE_SOLVER
 762             if (solver_show_working)
 763                 printf("%*sguessing %d at (%d,%d)\n",
 764                        solver_recurse_depth*4, "", list[i], x, y);
 765             solver_recurse_depth++;
 766 #endif
 767
 768             ret = cb(outgrid, o, recdiff, ctx);
 769
 770 #ifdef STANDALONE_SOLVER
 771             solver_recurse_depth--;
 772             if (solver_show_working) {
 773                 printf("%*sretracting %d at (%d,%d)\n",
 774                        solver_recurse_depth*4, "", list[i], x, y);
 775             }
 776 #endif
 777             /* we recurse as deep as we can, so we should never find
 778              * find ourselves giving up on a puzzle without declaring it
 779              * impossible.  */
 780             assert(ret != diff_unfinished);
 781
 782             /*
 783              * If we have our first solution, copy it into the
 784              * grid we will return.
 785              */
 786             if (diff == diff_impossible && ret != diff_impossible)
 787                 memcpy(solver->grid, outgrid, o*o);
 788
 789             if (ret == diff_ambiguous)
 790                 diff = diff_ambiguous;
 791             else if (ret == diff_impossible)
 792                 /* do not change our return value */;
 793             else {
 794                 /* the recursion turned up exactly one solution */
 795                 if (diff == diff_impossible)
 796                     diff = recdiff;
 797                 else
 798                     diff = diff_ambiguous;
 799             }
 800
 801             /*
 802              * As soon as we've found more than one solution,
 803              * give up immediately.
 804              */
 805             if (diff == diff_ambiguous)
 806                 break;
 807         }
 808
 809         sfree(outgrid);
 810         sfree(ingrid);
 811         sfree(list);
 812
 813         if (diff == diff_impossible)
 814             return -1;
 815         else if (diff == diff_ambiguous)
 816             return 2;
 817         else {
 818             assert(diff == recdiff);
 819             return 1;
 820         }
 821     }
 822 }
 823
 824 enum { diff_simple = 1, diff_set, diff_extreme, diff_recursive };
 825
 826 static int latin_solver_sub(struct latin_solver *solver, int maxdiff, void *ctx)
 827 {
 828     struct latin_solver_scratch *scratch = latin_solver_new_scratch(solver);
 829     int ret, diff = diff_simple, extreme;
 830
 831     assert(maxdiff <= diff_recursive);
 832     /*
 833      * Now loop over the grid repeatedly trying all permitted modes
 834      * of reasoning. The loop terminates if we complete an
 835      * iteration without making any progress; we then return
 836      * failure or success depending on whether the grid is full or
 837      * not.
 838      */
 839     while (1) {
 840         /*
 841          * I'd like to write `continue;' inside each of the
 842          * following loops, so that the solver returns here after
 843          * making some progress. However, I can't specify that I
 844          * want to continue an outer loop rather than the innermost
 845          * one, so I'm apologetically resorting to a goto.
 846          */
 847         cont:
 848         latin_solver_debug(solver->cube, solver->o);
 849
 850         ret = latin_solver_diff_simple(solver);
 851         if (ret < 0) {
 852             diff = diff_impossible;
 853             goto got_result;
 854         } else if (ret > 0) {
 855             diff = max(diff, diff_simple);
 856             goto cont;
 857         }
 858
 859         if (maxdiff <= diff_simple)
 860             break;
 861
 862         ret = latin_solver_diff_set(solver, scratch, &extreme);
 863         if (ret < 0) {
 864             diff = diff_impossible;
 865             goto got_result;
 866         } else if (ret > 0) {
 867             diff = max(diff, extreme ? diff_extreme : diff_set);
 868             goto cont;
 869         }
 870
 871         if (maxdiff <= diff_set)
 872             break;
 873
 874         /*
 875          * Forcing chains.
 876          */
 877         if (latin_solver_forcing(solver, scratch)) {
 878             diff = max(diff, diff_extreme);
 879             goto cont;
 880         }
 881
 882         /*
 883          * If we reach here, we have made no deductions in this
 884          * iteration, so the algorithm terminates.
 885          */
 886         break;
 887     }
 888
 889     /*
 890      * Last chance: if we haven't fully solved the puzzle yet, try
 891      * recursing based on guesses for a particular square. We pick
 892      * one of the most constrained empty squares we can find, which
 893      * has the effect of pruning the search tree as much as
 894      * possible.
 895      */
 896     if (maxdiff == diff_recursive) {
 897         int nsol = latin_solver_recurse(solver, diff_recursive, latin_solver, ctx);
 898         if (nsol < 0) diff = diff_impossible;
 899         else if (nsol == 1) diff = diff_recursive;
 900         else if (nsol > 1) diff = diff_ambiguous;
 901         /* if nsol == 0 then we were complete anyway
 902          * (and thus don't need to change diff) */
 903     } else {
 904         /*
 905          * We're forbidden to use recursion, so we just see whether
 906          * our grid is fully solved, and return diff_unfinished
 907          * otherwise.
 908          */
 909         int x, y, o = solver->o;
 910
 911         for (y = 0; y < o; y++)
 912             for (x = 0; x < o; x++)
 913                 if (!solver->grid[y*o+x])
 914                     diff = diff_unfinished;
 915     }
 916
 917     got_result:
 918
 919 #ifdef STANDALONE_SOLVER
 920     if (solver_show_working)
 921         printf("%*s%s found\n",
 922                solver_recurse_depth*4, "",
 923                diff == diff_impossible ? "no solution (impossible)" :
 924                diff == diff_unfinished ? "no solution (unfinished)" :
 925                diff == diff_ambiguous ? "multiple solutions" :
 926                "one solution");
 927 #endif
 928
 929     latin_solver_free_scratch(scratch);
 930
 931     return diff;
 932 }
 933
 934 int latin_solver(digit *grid, int o, int maxdiff, void *ctx)
 935 {
 936     struct latin_solver solver;
 937     int diff;
 938
 939     latin_solver_alloc(&solver, grid, o);
 940     diff = latin_solver_sub(&solver, maxdiff, ctx);
 941     latin_solver_free(&solver);
 942     return diff;
 943 }
 944
 945 void latin_solver_debug(unsigned char *cube, int o)
 946 {
 947 #ifdef STANDALONE_SOLVER
 948     if (solver_show_working) {
 949         struct latin_solver ls, *solver = &ls;
 950         char *dbg;
 951         int x, y, i, c = 0;
 952
 953         ls.cube = cube; ls.o = o; /* for cube() to work */
 954
 955         dbg = snewn(3*o*o*o, unsigned char);
 956         for (y = 0; y < o; y++) {
 957             for (x = 0; x < o; x++) {
 958                 for (i = 1; i <= o; i++) {
 959                     if (cube(x,y,i))
 960                         dbg[c++] = i + '0';
 961                     else
 962                         dbg[c++] = '.';
 963                 }
 964                 dbg[c++] = ' ';
 965             }
 966             dbg[c++] = '\n';
 967         }
 968         dbg[c++] = '\n';
 969         dbg[c++] = '\0';
 970
 971         printf("%s", dbg);
 972         sfree(dbg);
 973     }
 974 #endif
 975 }
 976
 977 void latin_debug(digit *sq, int o)
 978 {
 979 #ifdef STANDALONE_SOLVER
 980     if (solver_show_working) {
 981         int x, y;
 982
 983         for (y = 0; y < o; y++) {
 984             for (x = 0; x < o; x++) {
 985                 printf("%2d ", sq[y*o+x]);
 986             }
 987             printf("\n");
 988         }
 989         printf("\n");
 990     }
 991 #endif
 992 }
 993
 994 /* --------------------------------------------------------
 995  * Generation.
 996  */
 997
 998 digit *latin_generate(int o, random_state *rs)
 999 {
1000     digit *sq;
1001     int *edges, *backedges, *capacity, *flow;
1002     void *scratch;
1003     int ne, scratchsize;
1004     int i, j, k;
1005     digit *row, *col, *numinv, *num;
1006
1007     /*
1008      * To efficiently generate a latin square in such a way that
1009      * all possible squares are possible outputs from the function,
1010      * we make use of a theorem which states that any r x n latin
1011      * rectangle, with r < n, can be extended into an (r+1) x n
1012      * latin rectangle. In other words, we can reliably generate a
1013      * latin square row by row, by at every stage writing down any
1014      * row at all which doesn't conflict with previous rows, and
1015      * the theorem guarantees that we will never have to backtrack.
1016      *
1017      * To find a viable row at each stage, we can make use of the
1018      * support functions in maxflow.c.
1019      */
1020
1021     sq = snewn(o*o, digit);
1022
1023     /*
1024      * In case this method of generation introduces a really subtle
1025      * top-to-bottom directional bias, we'll generate the rows in
1026      * random order.
1027      */
1028     row = snewn(o, digit);
1029     col = snewn(o, digit);
1030     numinv = snewn(o, digit);
1031     num = snewn(o, digit);
1032     for (i = 0; i < o; i++)
1033         row[i] = i;
1034     shuffle(row, i, sizeof(*row), rs);
1035
1036     /*
1037      * Set up the infrastructure for the maxflow algorithm.
1038      */
1039     scratchsize = maxflow_scratch_size(o * 2 + 2);
1040     scratch = smalloc(scratchsize);
1041     backedges = snewn(o*o + 2*o, int);
1042     edges = snewn((o*o + 2*o) * 2, int);
1043     capacity = snewn(o*o + 2*o, int);
1044     flow = snewn(o*o + 2*o, int);
1045     /* Set up the edge array, and the initial capacities. */
1046     ne = 0;
1047     for (i = 0; i < o; i++) {
1048         /* Each LHS vertex is connected to all RHS vertices. */
1049         for (j = 0; j < o; j++) {
1050             edges[ne*2] = i;
1051             edges[ne*2+1] = j+o;
1052             /* capacity for this edge is set later on */
1053             ne++;
1054         }
1055     }
1056     for (i = 0; i < o; i++) {
1057         /* Each RHS vertex is connected to the distinguished sink vertex. */
1058         edges[ne*2] = i+o;
1059         edges[ne*2+1] = o*2+1;
1060         capacity[ne] = 1;
1061         ne++;
1062     }
1063     for (i = 0; i < o; i++) {
1064         /* And the distinguished source vertex connects to each LHS vertex. */
1065         edges[ne*2] = o*2;
1066         edges[ne*2+1] = i;
1067         capacity[ne] = 1;
1068         ne++;
1069     }
1070     assert(ne == o*o + 2*o);
1071     /* Now set up backedges. */
1072     maxflow_setup_backedges(ne, edges, backedges);
1073
1074     /*
1075      * Now generate each row of the latin square.
1076      */
1077     for (i = 0; i < o; i++) {
1078         /*
1079          * To prevent maxflow from behaving deterministically, we
1080          * separately permute the columns and the digits for the
1081          * purposes of the algorithm, differently for every row.
1082          */
1083         for (j = 0; j < o; j++)
1084             col[j] = num[j] = j;
1085         shuffle(col, j, sizeof(*col), rs);
1086         /* We need the num permutation in both forward and inverse forms. */
1087         for (j = 0; j < o; j++)
1088             numinv[num[j]] = j;
1089
1090         /*
1091          * Set up the capacities for the maxflow run, by examining
1092          * the existing latin square.
1093          */
1094         for (j = 0; j < o*o; j++)
1095             capacity[j] = 1;
1096         for (j = 0; j < i; j++)
1097             for (k = 0; k < o; k++) {
1098                 int n = num[sq[row[j]*o + col[k]] - 1];
1099                 capacity[k*o+n] = 0;
1100             }
1101
1102         /*
1103          * Run maxflow.
1104          */
1105         j = maxflow_with_scratch(scratch, o*2+2, 2*o, 2*o+1, ne,
1106                                  edges, backedges, capacity, flow, NULL);
1107         assert(j == o);   /* by the above theorem, this must have succeeded */
1108
1109         /*
1110          * And examine the flow array to pick out the new row of
1111          * the latin square.
1112          */
1113         for (j = 0; j < o; j++) {
1114             for (k = 0; k < o; k++) {
1115                 if (flow[j*o+k])
1116                     break;
1117             }
1118             assert(k < o);
1119             sq[row[i]*o + col[j]] = numinv[k] + 1;
1120         }
1121     }
1122
1123     /*
1124      * Done. Free our internal workspaces...
1125      */
1126     sfree(flow);
1127     sfree(capacity);
1128     sfree(edges);
1129     sfree(backedges);
1130     sfree(scratch);
1131     sfree(numinv);
1132     sfree(num);
1133     sfree(col);
1134     sfree(row);
1135
1136     /*
1137      * ... and return our completed latin square.
1138      */
1139     return sq;
1140 }
1141
1142 /* --------------------------------------------------------
1143  * Checking.
1144  */
1145
1146 typedef struct lcparams {
1147     digit elt;
1148     int count;
1149 } lcparams;
1150
1151 static int latin_check_cmp(void *v1, void *v2)
1152 {
1153     lcparams *lc1 = (lcparams *)v1;
1154     lcparams *lc2 = (lcparams *)v2;
1155
1156     if (lc1->elt < lc2->elt) return -1;
1157     if (lc1->elt > lc2->elt) return 1;
1158     return 0;
1159 }
1160
1161 #define ELT(sq,x,y) (sq[((y)*order)+(x)])
1162
1163 /* returns non-zero if sq is not a latin square. */
1164 int latin_check(digit *sq, int order)
1165 {
1166     tree234 *dict = newtree234(latin_check_cmp);
1167     int c, r;
1168     int ret = 0;
1169     lcparams *lcp, lc;
1170
1171     /* Use a tree234 as a simple hash table, go through the square
1172      * adding elements as we go or incrementing their counts. */
1173     for (c = 0; c < order; c++) {
1174         for (r = 0; r < order; r++) {
1175             lc.elt = ELT(sq, c, r); lc.count = 0;
1176             lcp = find234(dict, &lc, NULL);
1177             if (!lcp) {
1178                 lcp = snew(lcparams);
1179                 lcp->elt = ELT(sq, c, r);
1180                 lcp->count = 1;
1181                 assert(add234(dict, lcp) == lcp);
1182             } else {
1183                 lcp->count++;
1184             }
1185         }
1186     }
1187
1188     /* There should be precisely 'order' letters in the alphabet,
1189      * each occurring 'order' times (making the OxO tree) */
1190     if (count234(dict) != order) ret = 1;
1191     else {
1192         for (c = 0; (lcp = index234(dict, c)) != NULL; c++) {
1193             if (lcp->count != order) ret = 1;
1194         }
1195     }
1196     for (c = 0; (lcp = index234(dict, c)) != NULL; c++)
1197         sfree(lcp);
1198     freetree234(dict);
1199
1200     return ret;
1201 }
1202
1203
1204 /* --------------------------------------------------------
1205  * Testing (and printing).
1206  */
1207
1208 #ifdef STANDALONE_LATIN_TEST
1209
1210 #include <stdio.h>
1211 #include <time.h>
1212
1213 const char *quis;
1214
1215 static void latin_print(digit *sq, int order)
1216 {
1217     int x, y;
1218
1219     for (y = 0; y < order; y++) {
1220         for (x = 0; x < order; x++) {
1221             printf("%2u ", ELT(sq, x, y));
1222         }
1223         printf("\n");
1224     }
1225     printf("\n");
1226 }
1227
1228 static void gen(int order, random_state *rs, int debug)
1229 {
1230     digit *sq;
1231
1232     solver_show_working = debug;
1233
1234     sq = latin_generate(order, rs);
1235     latin_print(sq, order);
1236     if (latin_check(sq, order)) {
1237         fprintf(stderr, "Square is not a latin square!");
1238         exit(1);
1239     }
1240
1241     sfree(sq);
1242 }
1243
1244 void test_soak(int order, random_state *rs)
1245 {
1246     digit *sq;
1247     int n = 0;
1248     time_t tt_start, tt_now, tt_last;
1249
1250     solver_show_working = 0;
1251     tt_now = tt_start = time(NULL);
1252
1253     while(1) {
1254         sq = latin_generate(order, rs);
1255         sfree(sq);
1256         n++;
1257
1258         tt_last = time(NULL);
1259         if (tt_last > tt_now) {
1260             tt_now = tt_last;
1261             printf("%d total, %3.1f/s\n", n,
1262                    (double)n / (double)(tt_now - tt_start));
1263         }
1264     }
1265 }
1266
1267 void usage_exit(const char *msg)
1268 {
1269     if (msg)
1270         fprintf(stderr, "%s: %s\n", quis, msg);
1271     fprintf(stderr, "Usage: %s [--seed SEED] --soak <params> | [game_id [game_id ...]]\n", quis);
1272     exit(1);
1273 }
1274
1275 int main(int argc, char *argv[])
1276 {
1277     int i, soak = 0;
1278     random_state *rs;
1279     time_t seed = time(NULL);
1280
1281     quis = argv[0];
1282     while (--argc > 0) {
1283         const char *p = *++argv;
1284         if (!strcmp(p, "--soak"))
1285             soak = 1;
1286         else if (!strcmp(p, "--seed")) {
1287             if (argc == 0)
1288                 usage_exit("--seed needs an argument");
1289             seed = (time_t)atoi(*++argv);
1290             argc--;
1291         } else if (*p == '-')
1292                 usage_exit("unrecognised option");
1293         else
1294             break; /* finished options */
1295     }
1296
1297     rs = random_new((void*)&seed, sizeof(time_t));
1298
1299     if (soak == 1) {
1300         if (argc != 1) usage_exit("only one argument for --soak");
1301         test_soak(atoi(*argv), rs);
1302     } else {
1303         if (argc > 0) {
1304             for (i = 0; i < argc; i++) {
1305                 gen(atoi(*argv++), rs, 1);
1306             }
1307         } else {
1308             while (1) {
1309                 i = random_upto(rs, 20) + 1;
1310                 gen(i, rs, 0);
1311             }
1312         }
1313     }
1314     random_free(rs);
1315     return 0;
1316 }
1317
1318 #endif
1319
1320 /* vim: set shiftwidth=4 tabstop=8: */