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     if (!extreme) {
 627         /*
 628          * Row-wise set elimination.
 629          */
 630         for (y = 0; y < o; y++) {
 631             ret = latin_solver_set(solver, scratch, cubepos(0,y,1), o*o, 1
 632 #ifdef STANDALONE_SOLVER
 633                                    , "set elimination, row %d", YUNTRANS(y)
 634 #endif
 635                                   );
 636             if (ret != 0) return ret;
 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) return ret;
 648         }
 649     } else {
 650         /*
 651          * Row-vs-column set elimination on a single number
 652          * (much tricker for a human to do!)
 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) return ret;
 661         }
 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;
 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, 0);
 863         if (ret < 0) {
 864             diff = diff_impossible;
 865             goto got_result;
 866         } else if (ret > 0) {
 867             diff = max(diff, diff_set);
 868             goto cont;
 869         }
 870
 871         if (maxdiff <= diff_set)
 872             break;
 873
 874         ret = latin_solver_diff_set(solver, scratch, 1);
 875         if (ret < 0) {
 876             diff = diff_impossible;
 877             goto got_result;
 878         } else if (ret > 0) {
 879             diff = max(diff, diff_extreme);
 880             goto cont;
 881         }
 882
 883         /*
 884          * Forcing chains.
 885          */
 886         if (latin_solver_forcing(solver, scratch)) {
 887             diff = max(diff, diff_extreme);
 888             goto cont;
 889         }
 890
 891         /*
 892          * If we reach here, we have made no deductions in this
 893          * iteration, so the algorithm terminates.
 894          */
 895         break;
 896     }
 897
 898     /*
 899      * Last chance: if we haven't fully solved the puzzle yet, try
 900      * recursing based on guesses for a particular square. We pick
 901      * one of the most constrained empty squares we can find, which
 902      * has the effect of pruning the search tree as much as
 903      * possible.
 904      */
 905     if (maxdiff == diff_recursive) {
 906         int nsol = latin_solver_recurse(solver, diff_recursive, latin_solver, ctx);
 907         if (nsol < 0) diff = diff_impossible;
 908         else if (nsol == 1) diff = diff_recursive;
 909         else if (nsol > 1) diff = diff_ambiguous;
 910         /* if nsol == 0 then we were complete anyway
 911          * (and thus don't need to change diff) */
 912     } else {
 913         /*
 914          * We're forbidden to use recursion, so we just see whether
 915          * our grid is fully solved, and return diff_unfinished
 916          * otherwise.
 917          */
 918         int x, y, o = solver->o;
 919
 920         for (y = 0; y < o; y++)
 921             for (x = 0; x < o; x++)
 922                 if (!solver->grid[y*o+x])
 923                     diff = diff_unfinished;
 924     }
 925
 926     got_result:
 927
 928 #ifdef STANDALONE_SOLVER
 929     if (solver_show_working)
 930         printf("%*s%s found\n",
 931                solver_recurse_depth*4, "",
 932                diff == diff_impossible ? "no solution (impossible)" :
 933                diff == diff_unfinished ? "no solution (unfinished)" :
 934                diff == diff_ambiguous ? "multiple solutions" :
 935                "one solution");
 936 #endif
 937
 938     latin_solver_free_scratch(scratch);
 939
 940     return diff;
 941 }
 942
 943 int latin_solver(digit *grid, int o, int maxdiff, void *ctx)
 944 {
 945     struct latin_solver solver;
 946     int diff;
 947
 948     latin_solver_alloc(&solver, grid, o);
 949     diff = latin_solver_sub(&solver, maxdiff, ctx);
 950     latin_solver_free(&solver);
 951     return diff;
 952 }
 953
 954 void latin_solver_debug(unsigned char *cube, int o)
 955 {
 956 #ifdef STANDALONE_SOLVER
 957     if (solver_show_working) {
 958         struct latin_solver ls, *solver = &ls;
 959         char *dbg;
 960         int x, y, i, c = 0;
 961
 962         ls.cube = cube; ls.o = o; /* for cube() to work */
 963
 964         dbg = snewn(3*o*o*o, unsigned char);
 965         for (y = 0; y < o; y++) {
 966             for (x = 0; x < o; x++) {
 967                 for (i = 1; i <= o; i++) {
 968                     if (cube(x,y,i))
 969                         dbg[c++] = i + '0';
 970                     else
 971                         dbg[c++] = '.';
 972                 }
 973                 dbg[c++] = ' ';
 974             }
 975             dbg[c++] = '\n';
 976         }
 977         dbg[c++] = '\n';
 978         dbg[c++] = '\0';
 979
 980         printf("%s", dbg);
 981         sfree(dbg);
 982     }
 983 #endif
 984 }
 985
 986 void latin_debug(digit *sq, int o)
 987 {
 988 #ifdef STANDALONE_SOLVER
 989     if (solver_show_working) {
 990         int x, y;
 991
 992         for (y = 0; y < o; y++) {
 993             for (x = 0; x < o; x++) {
 994                 printf("%2d ", sq[y*o+x]);
 995             }
 996             printf("\n");
 997         }
 998         printf("\n");
 999     }
1000 #endif
1001 }
1002
1003 /* --------------------------------------------------------
1004  * Generation.
1005  */
1006
1007 digit *latin_generate(int o, random_state *rs)
1008 {
1009     digit *sq;
1010     int *edges, *backedges, *capacity, *flow;
1011     void *scratch;
1012     int ne, scratchsize;
1013     int i, j, k;
1014     digit *row, *col, *numinv, *num;
1015
1016     /*
1017      * To efficiently generate a latin square in such a way that
1018      * all possible squares are possible outputs from the function,
1019      * we make use of a theorem which states that any r x n latin
1020      * rectangle, with r < n, can be extended into an (r+1) x n
1021      * latin rectangle. In other words, we can reliably generate a
1022      * latin square row by row, by at every stage writing down any
1023      * row at all which doesn't conflict with previous rows, and
1024      * the theorem guarantees that we will never have to backtrack.
1025      *
1026      * To find a viable row at each stage, we can make use of the
1027      * support functions in maxflow.c.
1028      */
1029
1030     sq = snewn(o*o, digit);
1031
1032     /*
1033      * In case this method of generation introduces a really subtle
1034      * top-to-bottom directional bias, we'll generate the rows in
1035      * random order.
1036      */
1037     row = snewn(o, digit);
1038     col = snewn(o, digit);
1039     numinv = snewn(o, digit);
1040     num = snewn(o, digit);
1041     for (i = 0; i < o; i++)
1042         row[i] = i;
1043     shuffle(row, i, sizeof(*row), rs);
1044
1045     /*
1046      * Set up the infrastructure for the maxflow algorithm.
1047      */
1048     scratchsize = maxflow_scratch_size(o * 2 + 2);
1049     scratch = smalloc(scratchsize);
1050     backedges = snewn(o*o + 2*o, int);
1051     edges = snewn((o*o + 2*o) * 2, int);
1052     capacity = snewn(o*o + 2*o, int);
1053     flow = snewn(o*o + 2*o, int);
1054     /* Set up the edge array, and the initial capacities. */
1055     ne = 0;
1056     for (i = 0; i < o; i++) {
1057         /* Each LHS vertex is connected to all RHS vertices. */
1058         for (j = 0; j < o; j++) {
1059             edges[ne*2] = i;
1060             edges[ne*2+1] = j+o;
1061             /* capacity for this edge is set later on */
1062             ne++;
1063         }
1064     }
1065     for (i = 0; i < o; i++) {
1066         /* Each RHS vertex is connected to the distinguished sink vertex. */
1067         edges[ne*2] = i+o;
1068         edges[ne*2+1] = o*2+1;
1069         capacity[ne] = 1;
1070         ne++;
1071     }
1072     for (i = 0; i < o; i++) {
1073         /* And the distinguished source vertex connects to each LHS vertex. */
1074         edges[ne*2] = o*2;
1075         edges[ne*2+1] = i;
1076         capacity[ne] = 1;
1077         ne++;
1078     }
1079     assert(ne == o*o + 2*o);
1080     /* Now set up backedges. */
1081     maxflow_setup_backedges(ne, edges, backedges);
1082
1083     /*
1084      * Now generate each row of the latin square.
1085      */
1086     for (i = 0; i < o; i++) {
1087         /*
1088          * To prevent maxflow from behaving deterministically, we
1089          * separately permute the columns and the digits for the
1090          * purposes of the algorithm, differently for every row.
1091          */
1092         for (j = 0; j < o; j++)
1093             col[j] = num[j] = j;
1094         shuffle(col, j, sizeof(*col), rs);
1095         shuffle(num, j, sizeof(*num), rs);
1096         /* We need the num permutation in both forward and inverse forms. */
1097         for (j = 0; j < o; j++)
1098             numinv[num[j]] = j;
1099
1100         /*
1101          * Set up the capacities for the maxflow run, by examining
1102          * the existing latin square.
1103          */
1104         for (j = 0; j < o*o; j++)
1105             capacity[j] = 1;
1106         for (j = 0; j < i; j++)
1107             for (k = 0; k < o; k++) {
1108                 int n = num[sq[row[j]*o + col[k]] - 1];
1109                 capacity[k*o+n] = 0;
1110             }
1111
1112         /*
1113          * Run maxflow.
1114          */
1115         j = maxflow_with_scratch(scratch, o*2+2, 2*o, 2*o+1, ne,
1116                                  edges, backedges, capacity, flow, NULL);
1117         assert(j == o);   /* by the above theorem, this must have succeeded */
1118
1119         /*
1120          * And examine the flow array to pick out the new row of
1121          * the latin square.
1122          */
1123         for (j = 0; j < o; j++) {
1124             for (k = 0; k < o; k++) {
1125                 if (flow[j*o+k])
1126                     break;
1127             }
1128             assert(k < o);
1129             sq[row[i]*o + col[j]] = numinv[k] + 1;
1130         }
1131     }
1132
1133     /*
1134      * Done. Free our internal workspaces...
1135      */
1136     sfree(flow);
1137     sfree(capacity);
1138     sfree(edges);
1139     sfree(backedges);
1140     sfree(scratch);
1141     sfree(numinv);
1142     sfree(num);
1143     sfree(col);
1144     sfree(row);
1145
1146     /*
1147      * ... and return our completed latin square.
1148      */
1149     return sq;
1150 }
1151
1152 /* --------------------------------------------------------
1153  * Checking.
1154  */
1155
1156 typedef struct lcparams {
1157     digit elt;
1158     int count;
1159 } lcparams;
1160
1161 static int latin_check_cmp(void *v1, void *v2)
1162 {
1163     lcparams *lc1 = (lcparams *)v1;
1164     lcparams *lc2 = (lcparams *)v2;
1165
1166     if (lc1->elt < lc2->elt) return -1;
1167     if (lc1->elt > lc2->elt) return 1;
1168     return 0;
1169 }
1170
1171 #define ELT(sq,x,y) (sq[((y)*order)+(x)])
1172
1173 /* returns non-zero if sq is not a latin square. */
1174 int latin_check(digit *sq, int order)
1175 {
1176     tree234 *dict = newtree234(latin_check_cmp);
1177     int c, r;
1178     int ret = 0;
1179     lcparams *lcp, lc, *aret;
1180
1181     /* Use a tree234 as a simple hash table, go through the square
1182      * adding elements as we go or incrementing their counts. */
1183     for (c = 0; c < order; c++) {
1184         for (r = 0; r < order; r++) {
1185             lc.elt = ELT(sq, c, r); lc.count = 0;
1186             lcp = find234(dict, &lc, NULL);
1187             if (!lcp) {
1188                 lcp = snew(lcparams);
1189                 lcp->elt = ELT(sq, c, r);
1190                 lcp->count = 1;
1191                 aret = add234(dict, lcp);
1192                 assert(aret == lcp);
1193             } else {
1194                 lcp->count++;
1195             }
1196         }
1197     }
1198
1199     /* There should be precisely 'order' letters in the alphabet,
1200      * each occurring 'order' times (making the OxO tree) */
1201     if (count234(dict) != order) ret = 1;
1202     else {
1203         for (c = 0; (lcp = index234(dict, c)) != NULL; c++) {
1204             if (lcp->count != order) ret = 1;
1205         }
1206     }
1207     for (c = 0; (lcp = index234(dict, c)) != NULL; c++)
1208         sfree(lcp);
1209     freetree234(dict);
1210
1211     return ret;
1212 }
1213
1214
1215 /* --------------------------------------------------------
1216  * Testing (and printing).
1217  */
1218
1219 #ifdef STANDALONE_LATIN_TEST
1220
1221 #include <stdio.h>
1222 #include <time.h>
1223
1224 const char *quis;
1225
1226 static void latin_print(digit *sq, int order)
1227 {
1228     int x, y;
1229
1230     for (y = 0; y < order; y++) {
1231         for (x = 0; x < order; x++) {
1232             printf("%2u ", ELT(sq, x, y));
1233         }
1234         printf("\n");
1235     }
1236     printf("\n");
1237 }
1238
1239 static void gen(int order, random_state *rs, int debug)
1240 {
1241     digit *sq;
1242
1243     solver_show_working = debug;
1244
1245     sq = latin_generate(order, rs);
1246     latin_print(sq, order);
1247     if (latin_check(sq, order)) {
1248         fprintf(stderr, "Square is not a latin square!");
1249         exit(1);
1250     }
1251
1252     sfree(sq);
1253 }
1254
1255 void test_soak(int order, random_state *rs)
1256 {
1257     digit *sq;
1258     int n = 0;
1259     time_t tt_start, tt_now, tt_last;
1260
1261     solver_show_working = 0;
1262     tt_now = tt_start = time(NULL);
1263
1264     while(1) {
1265         sq = latin_generate(order, rs);
1266         sfree(sq);
1267         n++;
1268
1269         tt_last = time(NULL);
1270         if (tt_last > tt_now) {
1271             tt_now = tt_last;
1272             printf("%d total, %3.1f/s\n", n,
1273                    (double)n / (double)(tt_now - tt_start));
1274         }
1275     }
1276 }
1277
1278 void usage_exit(const char *msg)
1279 {
1280     if (msg)
1281         fprintf(stderr, "%s: %s\n", quis, msg);
1282     fprintf(stderr, "Usage: %s [--seed SEED] --soak <params> | [game_id [game_id ...]]\n", quis);
1283     exit(1);
1284 }
1285
1286 int main(int argc, char *argv[])
1287 {
1288     int i, soak = 0;
1289     random_state *rs;
1290     time_t seed = time(NULL);
1291
1292     quis = argv[0];
1293     while (--argc > 0) {
1294         const char *p = *++argv;
1295         if (!strcmp(p, "--soak"))
1296             soak = 1;
1297         else if (!strcmp(p, "--seed")) {
1298             if (argc == 0)
1299                 usage_exit("--seed needs an argument");
1300             seed = (time_t)atoi(*++argv);
1301             argc--;
1302         } else if (*p == '-')
1303                 usage_exit("unrecognised option");
1304         else
1305             break; /* finished options */
1306     }
1307
1308     rs = random_new((void*)&seed, sizeof(time_t));
1309
1310     if (soak == 1) {
1311         if (argc != 1) usage_exit("only one argument for --soak");
1312         test_soak(atoi(*argv), rs);
1313     } else {
1314         if (argc > 0) {
1315             for (i = 0; i < argc; i++) {
1316                 gen(atoi(*argv++), rs, 1);
1317             }
1318         } else {
1319             while (1) {
1320                 i = random_upto(rs, 20) + 1;
1321                 gen(i, rs, 0);
1322             }
1323         }
1324     }
1325     random_free(rs);
1326     return 0;
1327 }
1328
1329 #endif
1330
1331 /* vim: set shiftwidth=4 tabstop=8: */