## Unbiased grid generation: B2347 method

Programs which generate, solve, and analyze Sudoku puzzles

### Unbiased grid generation: B2347 method

About time I released the code, I guess.

Hidden Text: Show
Code: Select all
`/* urg2 : fast uniform random grids * * This uses the MTwist package available at http://www.cs.hmc.edu/~geoff/mtwist.html * * Compilation: *   \mingw\bin\mingw32-gcc -o urg2 urg2.c mtwist-0.6\mtwist.c mtwist-0.6\randistrs.c -Wall -O9 */ #include <stdio.h>#include <stdlib.h>#include <mem.h>#include <time.h>#include "mtwist-0.6/mtwist.h"#include "mtwist-0.6/randistrs.h"typedef unsigned long long u64;typedef unsigned int u32;typedef struct {  u32   b347;      /* max value is <280^3 */  u64   cdf;       /* max value is 510645067241472 */  short ext[16];   /* max value is 11664 */} rep_t;static rep_t rep[2865];static int   mrow[280][3];static u32 tmpl[10][16];static int ntmpl[10];static mt_state state;   /* rng state *//* ---------------------------------------------------------------------------------- * randomly pick a B5689 solution to a B12347 partial grid * ---------------------------------------------------------------------------------- */static void solve_b12347(int (*soln)[9], int segnr, int gridnr) {  int xok[6], yok[6], next;  int goodx[2][2], goody[2][2];  int ii, jj;  int i, j;    for (i=3; i<9; i++) yok[i-3] = 0x3fe ^ (1<<soln[i][0]) ^ (1<<soln[i][1]) ^ (1<<soln[i][2]);  for (j=3; j<9; j++) xok[j-3] = 0x3fe ^ (1<<soln[0][j]) ^ (1<<soln[1][j]) ^ (1<<soln[2][j]);        for (i=3,j=0; i<9; i++) if (yok[i-3]&2) { goody[j/2][j%2] = i; j++; }  for (j=3,i=0; j<9; j++) if (xok[j-3]&2) { goodx[i/2][i%2] = j; i++; }    for (next=segnr; next<segnr+1; next++) {      /* part 0: set 1s     */    soln[ goody[0][  (next>>0)&1] ][ goodx[0][  (next>>2)&1] ] = 1;    soln[ goody[0][1^(next>>0)&1] ][ goodx[1][  (next>>3)&1] ] = 1;    soln[ goody[1][  (next>>1)&1] ][ goodx[0][1^(next>>2)&1] ] = 1;    soln[ goody[1][1^(next>>1)&1] ][ goodx[1][1^(next>>3)&1] ] = 1;    /* part 1: get u32 templates for digits 2...9     * (needs to be fast, so we didn't use recursion)     */    {      int pos[6][6];      int digit;      for (ii=i=0; ii<6; ii++)        for (jj=0; jj<6; jj++)          if (soln[3+ii][3+jj] == 0)            pos[ii][jj] = i++;      for (digit=2; digit<=8; digit++) {        int iposs[4], jposs[4];        int i1, i2, j1, j2;        for (ii=i=0; ii<6; ii++) if (yok[ii]&(1<<digit)) iposs[i++] = ii;        for (jj=i=0; jj<6; jj++) if (xok[jj]&(1<<digit)) jposs[i++] = jj;        ntmpl[digit] = 0;        for (i1=0; i1<2; i1++) {          for (j1=0; j1<2; j1++) {            if (soln[3+iposs[i1]][3+jposs[j1]] == 0) {              for (i2=2; i2<4; i2++) {                if (soln[3+iposs[i2]][3+jposs[j1^1]] == 0) {                  for (j2=2; j2<4; j2++) {                    if (  soln[3+iposs[i1^1]][3+jposs[j2  ]]==0 \                       && soln[3+iposs[i2^1]][3+jposs[j2^1]]==0 ) {                      tmpl[digit][ntmpl[digit]++] =                   \                            (((u32)1)<<pos[iposs[i1  ]][jposs[j1  ]]) \                          | (((u32)1)<<pos[iposs[i1^1]][jposs[j2  ]]) \                          | (((u32)1)<<pos[iposs[i2  ]][jposs[j1^1]]) \                          | (((u32)1)<<pos[iposs[i2^1]][jposs[j2^1]]) ;                    }                  } /* j2 */                }              } /* i2 */            }          } /* j1 */        }  /* i1 */      } /* digit */    }  /* end of template-packing block */    /* part 2: enumerate solutions     * (quicker not to use compression for this tiny problem)     */    {      u32 mask;      int i2, i3, i4, i5, i6, i7, i8;      int ct = 0;      for (i2=0; i2<ntmpl[2]; i2++) {        mask = tmpl[2][i2];        for (i3=0; i3<ntmpl[3]; i3++) {          if (0 == (mask&tmpl[3][i3])) {            mask |= tmpl[3][i3];            for (i4=0; i4<ntmpl[4]; i4++) {              if (0 == (mask&tmpl[4][i4])) {                mask |= tmpl[4][i4];                for (i5=0; i5<ntmpl[5]; i5++) {                  if (0 == (mask&tmpl[5][i5])) {                    mask |= tmpl[5][i5];                    for (i6=0; i6<ntmpl[6]; i6++) {                      if (0 == (mask&tmpl[6][i6])) {                        mask |= tmpl[6][i6];                        for (i7=0; i7<ntmpl[7]; i7++) {                          if (0 == (mask&tmpl[7][i7])) {                            mask |= tmpl[7][i7];                            for (i8=0; i8<ntmpl[8]; i8++) {                              if (0==(mask&tmpl[8][i8]) && ct++==gridnr) {#define SETT(D)                                       \  {                                                   \    int i, j, k;                                      \    for (i=3,k=0; i<9; i++) for (j=3; j<9; j++) {     \      if (soln[i][j]!=1 && (tmpl[D][i##D]&(1<<k++)) ) \        soln[i][j] = D;                               \    }                                                 \  }SETT(2) SETT(3) SETT(4) SETT(5) SETT(6) SETT(7) SETT(8)for (i=0; i<81; i++) if (((int*)soln)[i] == 0) ((int*)soln)[i] = 9;return;#undef SETT                              }                            } /* i8 */                            mask ^= tmpl[7][i7];                          }                        } /* i7 */                        mask ^= tmpl[6][i6];                      }                    } /* i6 */                    mask ^= tmpl[5][i5];                  }                }  /* i5 */                mask ^= tmpl[4][i4];              }            }  /* i4 */            mask ^= tmpl[3][i3];          }        } /* i3 */      } /* i2 */    }  /* end of part 2 */    /* part 3: reset 1s     */    fprintf(stderr,"DOOM --- fell out of part 2 of the solver!\n"); exit(1);    soln[ goody[0][  (next>>0)&1] ][ goodx[0][  (next>>2)&1] ] = 0;    soln[ goody[0][1^(next>>0)&1] ][ goodx[1][  (next>>3)&1] ] = 0;    soln[ goody[1][  (next>>1)&1] ][ goodx[0][1^(next>>2)&1] ] = 0;    soln[ goody[1][1^(next>>1)&1] ][ goodx[1][1^(next>>3)&1] ] = 0;  } /* next */  return;}/* ---------------------------------------------------------------------------------- * randomly map a B2347 representative to a solved grid * ---------------------------------------------------------------------------------- */#define PLOOP(P) for (P=&perm[0]; P<&perm[18]; P+=3)static void solve_rep(rep_t *rep, int segnr, int gridnr) {  int perm[18] = { 0,1,2, 0,2,1, 1,0,2, 1,2,0, 2,0,1, 2,1,0 };  int grid[9][9], soln[9][9];  int b2[36][3][3], nb2, *b2mcp[2], b3mc[3];  int b4[36][3][3], nb4, *b4mrp[2], b7mr[3];  int i;  for (i=0; i<3; i++) {    grid[0][3+i] =  mrow[279][i] / 100;    grid[1][3+i] = (mrow[279][i]/10) % 10;    grid[2][3+i] =  mrow[279][i] % 10;    grid[0][6+i] =  mrow[ rep->b347/78400][i] / 100;    grid[1][6+i] = (mrow[ rep->b347/78400][i]/10) % 10;    grid[2][6+i] =  mrow[ rep->b347/78400][i] % 10;    grid[3+i][0] =  mrow[(rep->b347/280)%280][i] / 100;    grid[3+i][1] = (mrow[(rep->b347/280)%280][i]/10) % 10;    grid[3+i][2] =  mrow[(rep->b347/280)%280][i] % 10;    grid[6+i][0] =  mrow[ rep->b347%280][i] / 100;    grid[6+i][1] = (mrow[ rep->b347%280][i]/10) % 10;    grid[6+i][2] =  mrow[ rep->b347%280][i] % 10;  }  restart:   /* reset the solution grid */    memset(soln,0,81*sizeof(int));    /* get a list of good resolutions of B2 w.r.t. B3 */    b3mc[0] = (1<<grid[0][6]) + (1<<grid[1][6]) + (1<<grid[2][6]);  b3mc[1] = (1<<grid[0][7]) + (1<<grid[1][7]) + (1<<grid[2][7]);  b3mc[2] = (1<<grid[0][8]) + (1<<grid[1][8]) + (1<<grid[2][8]);  nb2=0; PLOOP(b2mcp[0]) PLOOP(b2mcp[1]) {    int xx[3][3], m;    xx[0][0] = grid[b2mcp[0][0]][3]; xx[0][1] = grid[b2mcp[1][0]][4]; xx[0][2] = grid[0][5];    xx[1][0] = grid[b2mcp[0][1]][3]; xx[1][1] = grid[b2mcp[1][1]][4]; xx[1][2] = grid[1][5];    xx[2][0] = grid[b2mcp[0][2]][3]; xx[2][1] = grid[b2mcp[1][2]][4]; xx[2][2] = grid[2][5];    if ( (m=(1<<xx[0][0])+(1<<xx[0][1])+(1<<xx[0][2]))==b3mc[0] || m==b3mc[1] || m==b3mc[2] ) continue;    if ( (m=(1<<xx[1][0])+(1<<xx[1][1])+(1<<xx[1][2]))==b3mc[0] || m==b3mc[1] || m==b3mc[2] ) continue;    if ( (m=(1<<xx[2][0])+(1<<xx[2][1])+(1<<xx[2][2]))==b3mc[0] || m==b3mc[1] || m==b3mc[2] ) continue;    memcpy(b2+nb2++,xx,9*sizeof(int));  }  /* pick a random B2 */    nb2 = rds_iuniform(&state,0,nb2);  for (i=0; i<9; i++) soln[i/3][3+i%3] = b2[nb2][i/3][i%3];    /* get a list of good resolutions of B4 w.r.t. B7 */    b7mr[0] = (1<<grid[6][0]) + (1<<grid[6][1]) + (1<<grid[6][2]);  b7mr[1] = (1<<grid[7][0]) + (1<<grid[7][1]) + (1<<grid[7][2]);  b7mr[2] = (1<<grid[8][0]) + (1<<grid[8][1]) + (1<<grid[8][2]);  nb4=0; PLOOP(b4mrp[0]) PLOOP(b4mrp[1]) {    int xx[3][3], m;    xx[0][0] = grid[3][b4mrp[0][0]]; xx[0][1] = grid[3][b4mrp[0][1]]; xx[0][2] = grid[3][b4mrp[0][2]];    xx[1][0] = grid[4][b4mrp[1][0]]; xx[1][1] = grid[4][b4mrp[1][1]]; xx[1][2] = grid[4][b4mrp[1][2]];    xx[2][0] = grid[5][         0 ]; xx[2][1] = grid[5][         1 ]; xx[2][2] = grid[5][         2 ];    if ( (m=(1<<xx[0][0])+(1<<xx[1][0])+(1<<xx[2][0]))==b7mr[0] || m==b7mr[1] || m==b7mr[2] ) continue;    if ( (m=(1<<xx[0][1])+(1<<xx[1][1])+(1<<xx[2][1]))==b7mr[0] || m==b7mr[1] || m==b7mr[2] ) continue;    if ( (m=(1<<xx[0][2])+(1<<xx[1][2])+(1<<xx[2][2]))==b7mr[0] || m==b7mr[1] || m==b7mr[2] ) continue;    memcpy(b4+nb4++,xx,9*sizeof(int));  }  /* pick a random B4 */    nb4 = rds_iuniform(&state,0,nb4);  for (i=0; i<9; i++) soln[3+i/3][i%3] = b4[nb4][i/3][i%3];    /* create a random one of the eight B3 resolutions */  {    int b2mr[3], *b3mcp[3], ct;    b2mr[0] = (1<<soln[0][3]) + (1<<soln[0][4]) + (1<<soln[0][5]);    b2mr[1] = (1<<soln[1][3]) + (1<<soln[1][4]) + (1<<soln[1][5]);    b2mr[2] = (1<<soln[2][3]) + (1<<soln[2][4]) + (1<<soln[2][5]);    ct=rds_iuniform(&state,1,9); PLOOP(b3mcp[0]) PLOOP(b3mcp[1]) PLOOP(b3mcp[2]) {#define B3(X,Y) soln[X][6+Y]      B3(0,0) = grid[b3mcp[0][0]][6]; B3(0,1) = grid[b3mcp[1][0]][7]; B3(0,2) = grid[b3mcp[2][0]][8];      B3(1,0) = grid[b3mcp[0][1]][6]; B3(1,1) = grid[b3mcp[1][1]][7]; B3(1,2) = grid[b3mcp[2][1]][8];      B3(2,0) = grid[b3mcp[0][2]][6]; B3(2,1) = grid[b3mcp[1][2]][7]; B3(2,2) = grid[b3mcp[2][2]][8];      if ( b2mr[0] & ((1<<B3(0,0))+(1<<B3(0,1))+(1<<B3(0,2))) ) continue;      if ( b2mr[1] & ((1<<B3(1,0))+(1<<B3(1,1))+(1<<B3(1,2))) ) continue;      if ( b2mr[2] & ((1<<B3(2,0))+(1<<B3(2,1))+(1<<B3(2,2))) ) continue;#undef B3      if (!--ct) goto got_b3;    }  }got_b3:    /* create a random one of the eight B7 resolutions */  {    int b4mc[3], *b7mrp[3], ct;    b4mc[0] = (1<<soln[3][0]) + (1<<soln[4][0]) + (1<<soln[5][0]);    b4mc[1] = (1<<soln[3][1]) + (1<<soln[4][1]) + (1<<soln[5][1]);    b4mc[2] = (1<<soln[3][2]) + (1<<soln[4][2]) + (1<<soln[5][2]);    ct=rds_iuniform(&state,1,9); PLOOP(b7mrp[0]) PLOOP(b7mrp[1]) PLOOP(b7mrp[2]) {#define B7(X,Y) soln[6+X][Y]      B7(0,0) = grid[6][b7mrp[0][0]]; B7(0,1) = grid[6][b7mrp[0][1]]; B7(0,2) = grid[6][b7mrp[0][2]];      B7(1,0) = grid[7][b7mrp[1][0]]; B7(1,1) = grid[7][b7mrp[1][1]]; B7(1,2) = grid[7][b7mrp[1][2]];      B7(2,0) = grid[8][b7mrp[2][0]]; B7(2,1) = grid[8][b7mrp[2][1]]; B7(2,2) = grid[8][b7mrp[2][2]];      if ( b4mc[0] & ((1<<B7(0,0))+(1<<B7(1,0))+(1<<B7(2,0))) ) continue;      if ( b4mc[1] & ((1<<B7(0,1))+(1<<B7(1,1))+(1<<B7(2,1))) ) continue;      if ( b4mc[2] & ((1<<B7(0,2))+(1<<B7(1,2))+(1<<B7(2,2))) ) continue;#undef B7      if (!--ct) goto got_b7;    }  }got_b7:    /* complete B1 then, if successful, solve for the rest of the grid */    {    int b1opts[3][3];    int i, j;    for (i=0; i<3; i++) for (j=0; j<3; j++) b1opts[i][j] = 0x3fe;    for (i=0; i<3; i++) for (j=3; j<9; j++) {      b1opts[i][0] &= ~(1<<soln[i][j]);      b1opts[i][1] &= ~(1<<soln[i][j]);      b1opts[i][2] &= ~(1<<soln[i][j]);      b1opts[0][i] &= ~(1<<soln[j][i]);      b1opts[1][i] &= ~(1<<soln[j][i]);      b1opts[2][i] &= ~(1<<soln[j][i]);    }    for (i=0; i<9; i++) if (!b1opts[i/3][i%3]) break;    if (i == 9) {      char printme[82];      for (i=0; i<9; i++) {        for (j=1; j<=9; j++) if (b1opts[i/3][i%3]&(1<<j)) break;        soln[i/3][i%3] = j;      }      solve_b12347(soln,segnr,gridnr);      for (i=0; i<81; i++) printme[i] = ((int*)soln)[i] + '0';      printme[81] = '\n';      fwrite(printme,1,82,stdout);    } else {      goto restart;    }  }    return;}#undef PLOOP/* ---------------------------------------------------------------------------------- * read my file of 2865 representatives of B2347 * ---------------------------------------------------------------------------------- */ static void die_read_err(char *path, int nr) {  fprintf(stderr,"read error: line %d, file %s\n",nr+1,path);  /* could be more helpful ... */  exit(1);}static void read_b2347(char *path) {  FILE *fp;  int  nr;  int  i;    if (NULL == (fp=fopen(path,"r"))) {    perror(path);    exit(1);  }    for (nr=0; nr<280; nr++) {    char s[20];    int  d[4];    if (5!=fscanf(fp,"%d %s %d %d %d",d,s,d+1,d+2,d+3) || d[0]!=nr || strcmp(s,"minirows")) die_read_err(path,nr);    mrow[nr][0] = d[1];    mrow[nr][1] = d[2];    mrow[nr][2] = d[3];  }      for (nr=0; nr<2865; nr++) {    char   s[3][20];    int    d[18];    double f;    if (22!=fscanf(fp,"%d %s %d %s %lf %s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",d,s[0],d+1,s[1],&f,s[2],d+2,d+3,d+4,d+5,d+6,d+7,d+8,d+9,d+10,d+11,d+12,d+13,d+14,d+15,d+16,d+17) || d[0]!=nr || strcmp(s[0],"b347") || strcmp(s[1],"cdf") || strcmp(s[2],"ext")) die_read_err(path,nr+280);    rep[nr].b347 = d[1];    rep[nr].cdf  = f;    for (i=0; i<16; i++) rep[nr].ext[i] = d[2+i];  }  fclose(fp);}/* ---------------------------------------------------------------------------------- * kluge to get around the fact that rds_liuniform returns only 32-bit nrs! * ---------------------------------------------------------------------------------- */static u64 kluge_rds_liuniform(mt_state *state, u64 lower, u64 upper) {  u64 mask, rnd;    if (lower) {    fprintf(stderr,"ERROR: lower!=0 not implemented\n");    exit(1);  }    mask  = upper - 1;  mask |= mask >> 1;  mask |= mask >> 2;  mask |= mask >> 4;  mask |= mask >> 8;  mask |= mask >> 16;  mask |= mask >> 32;    while ( (rnd=mts_llrand(state)&mask) >= upper ) ;    return rnd;}/* ---------------------------------------------------------------------------------- * PROGRAM START * ---------------------------------------------------------------------------------- */ int main(void) {  double   tsum = 0.0;  int      i;    /* seed the Mersenne Twister from the time of day   */  mts_goodseed(&state);      /* read "gang of 2865" file   */  read_b2347("gen2865.out");    /* generate grids!   */  fprintf(stderr,"Go!\n");#if 1  for (;;) {#else  for (i=0; i<100000; i++) {#endif    u64 rnd;    int classnr, segnr, gridnr;    clock_t c = -clock();    rnd = kluge_rds_liuniform(&state,(u64)0,rep[2864].cdf);    for (classnr=0; classnr<2865 && rep[classnr].cdf<=rnd; classnr++) ;    gridnr = rds_iuniform(&state,0,rep[classnr].ext[15]);    for (segnr=0; segnr<16 && rep[classnr].ext[segnr]<=gridnr; segnr++) ;    if (segnr) gridnr -= rep[classnr].ext[segnr-1];    solve_rep(rep+classnr,segnr,gridnr);    c += clock(); tsum += c / (double)CLOCKS_PER_SEC;  }  fprintf(stderr,"Stop!  Rate = %f grids/sec\n",i/tsum);       return 0;}`

... which requires the presence, in the same directory, of gen2856.out, created by this code:
Hidden Text: Show
Code: Select all
`/* this program does the data prep for "urg2" */ #include <stdio.h>#include <stdlib.h>#include <mem.h>typedef unsigned char u8;typedef unsigned int  u32;typedef struct {  u32    b347;      /* max value is <280^3 */  double cdf;       /* max value is 510645067241472 */  short  ext[16];   /* max value is 11664 */} rep_t;static rep_t rep[2865];static int boxlookup[65536];static int colBits[84];static int cblookup[449];static int colcompat[84][20][2];static int boxcompat[280][56];static int ijk[280][3];static int ijklookup[1024];static u32 tmpl[10][16];static int ntmpl[10];static int extlist[16];#define CBCRUNCH(B) (B)#define CBLOOKUP(B) cblookup[CBCRUNCH(B)]#define ORDER3(A,B,C)                         \  if (A<B) {                                  \    if (B<C)      { int t=A; A=C; C=t;      } \    else if (A<C) { int t=A; A=B; B=C; C=t; } \    else          { int t=A; A=B; B=t;      } \  } else {                                    \    if (B>C)      { /* do nothing */        } \    else if (A>C) { int t=B; B=C; C=t;      } \    else          { int t=A; A=C; C=B; B=t; } \  }#define IJKCRUNCH(I,J,K) ( ( ((J)<<5)+(K) ) & 0x3ff )#define IJKLOOKUP(I,J,K) ijklookup[IJKCRUNCH(I,J,K)]#define PLOOP(P) for (P=&perm[0]; P<&perm[18]; P+=3)/* ----------------------------------------------------------------------- * Set up various useful lookup tables: * *   colBits[0..83] = bitwise representation of each of the 84 *                    possible sets of minicolumn contents  * *   colcompat[C][0..19][0..1] : 84 lists, one per minicolumn C. *                  Each list is 20 pairs of numbers; each pair *                  indexes the two minicolumns compatible (sharing *                  no values) with the first minicolumn  * *   ijk[0..279][0..2] = descriptions of the 280 miniboxes with *                  ordered columns; values are numbers 0..83, *                  ordered s.t. ijk[B][0] > ijk[B][1] > ijk[B][2] * *   boxcompat[B][0..55] = 280 lists, one per minibox B.  Each list *                  is a number 0..279 of a box that can go below *                  this one.  The lists are ordered so that *                  boxcompat[B][i] and boxcompat[B][i^1] are mutually *                  compatible with box B (you can stack all three, *                  with suitable minicol ordering, on top of each *                  other), with even numbered indices always being *                  the lesser index of the pair * * We also set up a couple of reverse maps that let us go from one * representation of a column or box to another: * *   cblookup[0..1023] : maps a 9-bit "colBits" value back to an *                  index 0..83 * *   ijklookup[0..1023] : maps an ordered triple i>j>k of column *                  indices 0..83 to a box index 0..279; but first *                  you have to use IJKCRUNCH() to reduce i,j,k * ----------------------------------------------------------------------- */static void pretabulate(void) {  int i, j, k;  int n = 0;  /* colBits[]   */  for (i=0; i<7; i++)    for (j=i+1; j<8; j++)      for (k=j+1; k<9; k++)        colBits[n++] = (1<<i) + (1<<j) + (1<<k);  /* colcompat[][][]   */  for (i=0; i<84; i++) {    cblookup[CBCRUNCH(colBits[i])] = i;    for (j=n=0; j<84; j++) {      if (!(colBits[i]&colBits[j])) {        for (k=0; k<84; k++)          if (!((colBits[i]|colBits[j])&colBits[k]))            break;        colcompat[i][n][0]   = j;        colcompat[i][n++][1] = k;      }    }  }  /* ijk[][]   */  for (i=2,k=0; i<84; i++) {    for (j=0; j<20; j++) {      if (i>colcompat[i][j][0] && colcompat[i][j][0]>colcompat[i][j][1]) {        ijk[k][0] = i;        ijk[k][1] = colcompat[i][j][0];        ijk[k][2] = colcompat[i][j][1];        ijklookup[IJKCRUNCH(i,ijk[k][1],ijk[k][2])] = k;#define BUP(I,J) boxlookup[0xffff&(colBits[ijk[k][I]]+256*colBits[ijk[k][J]])] = k;        BUP(0,1) BUP(0,2) BUP(1,0) BUP(1,2) BUP(2,0) BUP(2,1)        k++;      }    }  }  /* boxcompat[][]   */  for (k=0; k<280; k++) {    int k0    =         ijk[k][0];    int bits1 = colBits[ijk[k][1]];    int bits2 = colBits[ijk[k][2]];    int n;    for (i=n=0; i<20; i++) {      int colsA0, colsA[3];      int colsB0, colsB[3];      if ( (colsA0=colcompat[k0][i][0]) < (colsB0=colcompat[k0][i][1]) ) {        for (j=0; j<20; j++) {          if (  !(bits1&colBits[colsA[1]=colcompat[colsA0][j][0]]) \             && !(bits2&colBits[colsA[2]=colcompat[colsA0][j][1]]) ) {            colsA[0] = colsA0;            colsB[0] = colsB0;            colsB[1] = CBLOOKUP(0x1ff^bits1^colBits[colsA[1]]);            colsB[2] = CBLOOKUP(0x1ff^bits2^colBits[colsA[2]]);            ORDER3(colsA[0],colsA[1],colsA[2]);            boxcompat[k][n++] = IJKLOOKUP(colsA[0],colsA[1],colsA[2]);            ORDER3(colsB[0],colsB[1],colsB[2]);                        boxcompat[k][n++] = IJKLOOKUP(colsB[0],colsB[1],colsB[2]);          }        }      }    }  }  return;}/* ----------------------------------------------------------------------- * permute(b2,b3,b4,b7,seen) loops over renumberings of b2 from which * renumberings of b3,b4,b7 are deduced, and marks seen[] accordingly * ----------------------------------------------------------------------- */static int permute(int b2, int b3, int b4, int b7, u8 *seen) {  int perm[18] = { 0,1,2, 0,2,1, 1,0,2, 1,2,0, 2,0,1, 2,1,0 };  int *rp[3], *cp;  int box2[3][3], box3[3][3], box4[3][3], box7[3][3];  int ct, i, j, k;  i = colBits[ijk[b2][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][0] = j;  i = colBits[ijk[b2][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][1] = j;  i = colBits[ijk[b2][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][2] = j;  i = colBits[ijk[b3][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][0] = j;  i = colBits[ijk[b3][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][1] = j;  i = colBits[ijk[b3][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][2] = j;  i = colBits[ijk[b4][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[0][k++] = j;  i = colBits[ijk[b4][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[1][k++] = j;  i = colBits[ijk[b4][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[2][k++] = j;  i = colBits[ijk[b7][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[0][k++] = j;  i = colBits[ijk[b7][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[1][k++] = j;  i = colBits[ijk[b7][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[2][k++] = j;  ct = 0;  PLOOP(cp) PLOOP(rp[0]) PLOOP(rp[1]) PLOOP(rp[2]) {    int bits[10];    for (i=0; i<3; i++) for (j=0; j<3; j++) bits[ box2[rp[j][i]][cp[j]] ] = 1 << (i+3*j);    i = cblookup[ bits[box3[0][0]] + bits[box3[1][0]] + bits[box3[2][0]] ];    j = cblookup[ bits[box3[0][1]] + bits[box3[1][1]] + bits[box3[2][1]] ];    k = cblookup[ bits[box3[0][2]] + bits[box3[1][2]] + bits[box3[2][2]] ];    ORDER3(i,j,k) b3 = ijklookup[IJKCRUNCH(i,j,k)];    i = cblookup[ bits[box4[0][0]] + bits[box4[0][1]] + bits[box4[0][2]] ];    j = cblookup[ bits[box4[1][0]] + bits[box4[1][1]] + bits[box4[1][2]] ];    k = cblookup[ bits[box4[2][0]] + bits[box4[2][1]] + bits[box4[2][2]] ];    ORDER3(i,j,k) b4 = ijklookup[IJKCRUNCH(i,j,k)];    i = cblookup[ bits[box7[0][0]] + bits[box7[0][1]] + bits[box7[0][2]] ];    j = cblookup[ bits[box7[1][0]] + bits[box7[1][1]] + bits[box7[1][2]] ];    k = cblookup[ bits[box7[2][0]] + bits[box7[2][1]] + bits[box7[2][2]] ];    ORDER3(i,j,k) b7 = ijklookup[IJKCRUNCH(i,j,k)];    if (b4 > b7) { i = b4; b4 = b7; b7 = i; }    i = b3*39340 + (b4*(559-b4))/2 + b7;    if (!(seen[i>>3] & (1<<(i&7)))) {      seen[i>>3] |= 1 << (i&7);      ct += 1 + (b4 != b7);    }  }  return ct;}/* ----------------------------------------------------------------------- * count nr completions from B1/2/3/4/7 to a full grid * ----------------------------------------------------------------------- */static void solve_b12347(int (*soln)[9]) {  int xok[6], yok[6], next;  int goodx[2][2], goody[2][2];  int ii, jj;  int i, j;    for (i=3; i<9; i++) yok[i-3] = 0x3fe ^ (1<<soln[i][0]) ^ (1<<soln[i][1]) ^ (1<<soln[i][2]);  for (j=3; j<9; j++) xok[j-3] = 0x3fe ^ (1<<soln[0][j]) ^ (1<<soln[1][j]) ^ (1<<soln[2][j]);        for (i=3,j=0; i<9; i++) if (yok[i-3]&2) { goody[j/2][j%2] = i; j++; }  for (j=3,i=0; j<9; j++) if (xok[j-3]&2) { goodx[i/2][i%2] = j; i++; }    for (next=0; next<16; next++) {      /* part 0: set 1s     */    soln[ goody[0][  (next>>0)&1] ][ goodx[0][  (next>>2)&1] ] = 1;    soln[ goody[0][1^(next>>0)&1] ][ goodx[1][  (next>>3)&1] ] = 1;    soln[ goody[1][  (next>>1)&1] ][ goodx[0][1^(next>>2)&1] ] = 1;    soln[ goody[1][1^(next>>1)&1] ][ goodx[1][1^(next>>3)&1] ] = 1;    /* part 1: get u32 templates for digits 2...9     * (needs to be fast, so we didn't use recursion)     */    {      int pos[6][6];      int digit;      for (ii=i=0; ii<6; ii++)        for (jj=0; jj<6; jj++)          if (soln[3+ii][3+jj] == 0)            pos[ii][jj] = i++;      for (digit=2; digit<=8; digit++) {        int iposs[4], jposs[4];        int i1, i2, j1, j2;        for (ii=i=0; ii<6; ii++) if (yok[ii]&(1<<digit)) iposs[i++] = ii;        for (jj=i=0; jj<6; jj++) if (xok[jj]&(1<<digit)) jposs[i++] = jj;        ntmpl[digit] = 0;        for (i1=0; i1<2; i1++) {          for (j1=0; j1<2; j1++) {            if (soln[3+iposs[i1]][3+jposs[j1]] == 0) {              for (i2=2; i2<4; i2++) {                if (soln[3+iposs[i2]][3+jposs[j1^1]] == 0) {                  for (j2=2; j2<4; j2++) {                    if (  soln[3+iposs[i1^1]][3+jposs[j2  ]]==0 \                       && soln[3+iposs[i2^1]][3+jposs[j2^1]]==0 ) {                      tmpl[digit][ntmpl[digit]++] =                   \                            (((u32)1)<<pos[iposs[i1  ]][jposs[j1  ]]) \                          | (((u32)1)<<pos[iposs[i1^1]][jposs[j2  ]]) \                          | (((u32)1)<<pos[iposs[i2  ]][jposs[j1^1]]) \                          | (((u32)1)<<pos[iposs[i2^1]][jposs[j2^1]]) ;                    }                  } /* j2 */                }              } /* i2 */            }          } /* j1 */        }  /* i1 */      } /* digit */    }  /* end of template-packing block */    /* part 2: enumerate solutions     * (quicker not to use compression for this tiny problem)     */    {      u32 mask;      int i2, i3, i4, i5, i6, i7, i8;      int ct = 0;      for (i2=0; i2<ntmpl[2]; i2++) {        mask = tmpl[2][i2];        for (i3=0; i3<ntmpl[3]; i3++) {          if (0 == (mask&tmpl[3][i3])) {            mask |= tmpl[3][i3];            for (i4=0; i4<ntmpl[4]; i4++) {              if (0 == (mask&tmpl[4][i4])) {                mask |= tmpl[4][i4];                for (i5=0; i5<ntmpl[5]; i5++) {                  if (0 == (mask&tmpl[5][i5])) {                    mask |= tmpl[5][i5];                    for (i6=0; i6<ntmpl[6]; i6++) {                      if (0 == (mask&tmpl[6][i6])) {                        mask |= tmpl[6][i6];                        for (i7=0; i7<ntmpl[7]; i7++) {                          if (0 == (mask&tmpl[7][i7])) {                            mask |= tmpl[7][i7];                            for (i8=0; i8<ntmpl[8]; i8++) {                              if (0 == (mask&tmpl[8][i8])) {                                ct++;                              }                            } /* i8 */                            mask ^= tmpl[7][i7];                          }                        } /* i7 */                        mask ^= tmpl[6][i6];                      }                    } /* i6 */                    mask ^= tmpl[5][i5];                  }                }  /* i5 */                mask ^= tmpl[4][i4];              }            }  /* i4 */            mask ^= tmpl[3][i3];          }        } /* i3 */      } /* i2 */            extlist[next] = ct + (next ? extlist[next-1] : 0);    }  /* end of part 2 */    /* part 3: reset 1s     */    soln[ goody[0][  (next>>0)&1] ][ goodx[0][  (next>>2)&1] ] = 0;    soln[ goody[0][1^(next>>0)&1] ][ goodx[1][  (next>>3)&1] ] = 0;    soln[ goody[1][  (next>>1)&1] ][ goodx[0][1^(next>>2)&1] ] = 0;    soln[ goody[1][1^(next>>1)&1] ][ goodx[1][1^(next>>3)&1] ] = 0;  } /* next */  return;}static void count_b12347(int b2, int b3, int b4, int b7) {  int grid[9][9] = {{0}};  int i, j, k;  /* unpack b*index into actual boxes */    i = colBits[ijk[b2][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][3] = j;  i = colBits[ijk[b2][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][4] = j;  i = colBits[ijk[b2][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][5] = j;  i = colBits[ijk[b3][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][6] = j;  i = colBits[ijk[b3][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][7] = j;  i = colBits[ijk[b3][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[k++][8] = j;  i = colBits[ijk[b4][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[3][k++] = j;  i = colBits[ijk[b4][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[4][k++] = j;  i = colBits[ijk[b4][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[5][k++] = j;  i = colBits[ijk[b7][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[6][k++] = j;  i = colBits[ijk[b7][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[7][k++] = j;  i = colBits[ijk[b7][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) grid[8][k++] = j;  /* count */    solve_b12347(grid);  return;}/* ----------------------------------------------------------------------- * count number of actual B2/3/4/7 corresponding to this representative * ----------------------------------------------------------------------- */static int b1mc[288][3], b1mr[288][3];  /* too big for stack, probably */static int resolutions(int b2, int b3, int b4, int b7) {  int perm[18] = { 0,1,2, 0,2,1, 1,0,2, 1,2,0, 2,0,1, 2,1,0 };  int box2[3][3], box3[3][3], box4[3][3], box7[3][3];  int b3mc[3], *b2mcp[2], b7mr[3], *b4mrp[2];  int nb23, nb47;  int ct;  int i, j, k;  /* unpack b*index into actual boxes */    i = colBits[ijk[b2][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][0] = j;  i = colBits[ijk[b2][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][1] = j;  i = colBits[ijk[b2][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box2[k++][2] = j;  i = colBits[ijk[b3][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][0] = j;  i = colBits[ijk[b3][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][1] = j;  i = colBits[ijk[b3][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box3[k++][2] = j;  i = colBits[ijk[b4][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[0][k++] = j;  i = colBits[ijk[b4][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[1][k++] = j;  i = colBits[ijk[b4][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box4[2][k++] = j;  i = colBits[ijk[b7][0]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[0][k++] = j;  i = colBits[ijk[b7][1]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[1][k++] = j;  i = colBits[ijk[b7][2]]; for (j=1,k=0; k<3; j++,i>>=1) if (i&1) box7[2][k++] = j;    /* find all valid B2/B3 combos when B2 minicol 2 is fixed (at most 36*8 = 288 ways) */    b3mc[0] = (1<<box3[0][0]) + (1<<box3[1][0]) + (1<<box3[2][0]);  b3mc[1] = (1<<box3[0][1]) + (1<<box3[1][1]) + (1<<box3[2][1]);  b3mc[2] = (1<<box3[0][2]) + (1<<box3[1][2]) + (1<<box3[2][2]);  nb23=0; PLOOP(b2mcp[0]) PLOOP(b2mcp[1]) {    int b2mr[3], b3mr[3], *b3mcp[3];    b2mr[0] = (1<<box2[b2mcp[0][0]][0]) + (1<<box2[b2mcp[1][0]][1]) + (1<<box2[0][2]);    b2mr[1] = (1<<box2[b2mcp[0][1]][0]) + (1<<box2[b2mcp[1][1]][1]) + (1<<box2[1][2]);    b2mr[2] = (1<<box2[b2mcp[0][2]][0]) + (1<<box2[b2mcp[1][2]][1]) + (1<<box2[2][2]);    PLOOP(b3mcp[0]) PLOOP(b3mcp[1]) PLOOP(b3mcp[2]) {      b3mr[0] = (1<<box3[b3mcp[0][0]][0]) + (1<<box3[b3mcp[1][0]][1]) + (1<<box3[b3mcp[2][0]][2]);      b3mr[1] = (1<<box3[b3mcp[0][1]][0]) + (1<<box3[b3mcp[1][1]][1]) + (1<<box3[b3mcp[2][1]][2]);      b3mr[2] = (1<<box3[b3mcp[0][2]][0]) + (1<<box3[b3mcp[1][2]][1]) + (1<<box3[b3mcp[2][2]][2]);      if ( (b2mr[0]&b3mr[0]) || (b2mr[1]&b3mr[1]) || (b2mr[2]&b3mr[2]) ) continue;      b1mr[nb23][0] = 0x3fe ^ b2mr[0] ^ b3mr[0];      b1mr[nb23][1] = 0x3fe ^ b2mr[1] ^ b3mr[1];      b1mr[nb23][2] = 0x3fe ^ b2mr[2] ^ b3mr[2];      nb23++;    }  }    /* find all valid B4/B7 combos when B4 minirow 2 is fixed (at most 36*8 = 288 ways) */    b7mr[0] = (1<<box7[0][0]) + (1<<box7[0][1]) + (1<<box7[0][2]);  b7mr[1] = (1<<box7[1][0]) + (1<<box7[1][1]) + (1<<box7[1][2]);  b7mr[2] = (1<<box7[2][0]) + (1<<box7[2][1]) + (1<<box7[2][2]);  nb47=0; PLOOP(b4mrp[0]) PLOOP(b4mrp[1]) {    int b4mc[3], b7mc[3], *b7mrp[3];    b4mc[0] = (1<<box4[0][b4mrp[0][0]]) + (1<<box4[1][b4mrp[1][0]]) + (1<<box4[2][0]);    b4mc[1] = (1<<box4[0][b4mrp[0][1]]) + (1<<box4[1][b4mrp[1][1]]) + (1<<box4[2][1]);    b4mc[2] = (1<<box4[0][b4mrp[0][2]]) + (1<<box4[1][b4mrp[1][2]]) + (1<<box4[2][2]);    PLOOP(b7mrp[0]) PLOOP(b7mrp[1]) PLOOP(b7mrp[2]) {      b7mc[0] = (1<<box7[0][b7mrp[0][0]]) + (1<<box7[1][b7mrp[1][0]]) + (1<<box7[2][b7mrp[2][0]]);      b7mc[1] = (1<<box7[0][b7mrp[0][1]]) + (1<<box7[1][b7mrp[1][1]]) + (1<<box7[2][b7mrp[2][1]]);      b7mc[2] = (1<<box7[0][b7mrp[0][2]]) + (1<<box7[1][b7mrp[1][2]]) + (1<<box7[2][b7mrp[2][2]]);      if ( (b4mc[0]&b7mc[0]) || (b4mc[1]&b7mc[1]) || (b4mc[2]&b7mc[2]) ) continue;      b1mc[nb47][0] = 0x3fe ^ b4mc[0] ^ b7mc[0];      b1mc[nb47][1] = 0x3fe ^ b4mc[1] ^ b7mc[1];      b1mc[nb47][2] = 0x3fe ^ b4mc[2] ^ b7mc[2];      nb47++;    }  }  /* ------------------------------ */  /* main counting loop starts here */  /* ------------------------------ */       ct = 0;  for (i=0; i<nb23; i++) for (j=0; j<nb47; j++) {    if (!(b1mr[i][0] & b1mc[j][0])) continue;    if (!(b1mr[i][0] & b1mc[j][1])) continue;    if (!(b1mr[i][0] & b1mc[j][2])) continue;    if (!(b1mr[i][1] & b1mc[j][0])) continue;    if (!(b1mr[i][1] & b1mc[j][1])) continue;    if (!(b1mr[i][1] & b1mc[j][2])) continue;    if (!(b1mr[i][2] & b1mc[j][0])) continue;    if (!(b1mr[i][2] & b1mc[j][1])) continue;    if (!(b1mr[i][2] & b1mc[j][2])) continue;    ct++;  }    return ct;}/* ----------------------------------------------------------------------- * PROGRAM START * ----------------------------------------------------------------------- */int main(void) {  u8  *seen;  int b3, b4, b7;  int classnr;    pretabulate();  if (NULL == (seen=calloc(1376900,1))) {    perror("calloc(seen)");    exit(1);  }    classnr = 0;  for (b3=0; b3<280; b3++) for (b4=0; b4<280; b4++) for (b7=b4; b7<280; b7++) {    int x = b3*39340 + (b4*(559-b4))/2 + b7;    if (!(seen[x>>3] & (1<<(x&7)))) {      count_b12347(279,b3,b4,b7);      rep[classnr].b347 = 280*280*b3 + 280*b4 + b7;      rep[classnr].cdf  = permute(279,b3,b4,b7,seen) + permute(b3,279,b4,b7,seen) \                        + permute(b4,b7,279,b3,seen) + permute(b7,b4,279,b3,seen) ;      rep[classnr].cdf *= resolutions(279,b3,b4,b7) * (double)extlist[15];      if (classnr) rep[classnr].cdf += rep[classnr-1].cdf;      for (x=0; x<16; x++) rep[classnr].ext[x] = extlist[x];      classnr++;    }  }    {    int i, j, k;    for (i=0; i<280; i++) {      printf("%d minirows ",i);      for (j=colBits[ijk[i][0]],k=0; k<9; k++) if (j&(1<<k)) printf("%d",k+1); printf(" ");      for (j=colBits[ijk[i][1]],k=0; k<9; k++) if (j&(1<<k)) printf("%d",k+1); printf(" ");      for (j=colBits[ijk[i][2]],k=0; k<9; k++) if (j&(1<<k)) printf("%d",k+1); printf("\n");    }    for (i=0; i<2865; i++) {      printf("%d b347 %d cdf %.0f ext",i,rep[i].b347,rep[i].cdf);      for (j=0; j<16; j++) printf(" %d",rep[i].ext[j]); printf("\n");    }      }    return 0;}`

If you spot bugs (I last compiled this years ago), let me know. Thanks.
Red Ed

Posts: 633
Joined: 06 June 2005

Thanks! The generator works great since I finally get proper results for the subset method.

Just a small note: I had to add the "-lm" option for the compiler to work on Linux.
Afmob

Posts: 130
Joined: 28 June 2011

### Re: Unbiased grid generation: B2347 method

I must be getting old! (well I am, but anyway ...)

I've been away from C for some time (and don't miss it at all). What on earth does indexing some var with 1 ^ (next >> 0) & 1 do?

BTW: I get a bunch of "warnings: suggest parentheses around arithmetic" from gcc wrt these

Mathimagics
2017 Supporter

Posts: 545
Joined: 27 May 2015

### Re: Unbiased grid generation: B2347 method

(next >> 0) is for readability and does nothing. 0 means which bit (which power of 2) we are interested in.
x & 1 takes the less significant bit and clears others. It is equivalent of (x modulo 2).
1 ^ x inverts the less significant bit.
1 ^ (next >> 0) & 1 is equivalent of 1 - (next modulo 2). It returns 0 for odd next and 1 for even next.
dobrichev
2016 Supporter

Posts: 1369
Joined: 24 May 2010

### Re: Unbiased grid generation: B2347 method

Ok, I forgot "^" is XOR operator.

So he means 1 - (next & 1), it seems!

Mathimagics
2017 Supporter

Posts: 545
Joined: 27 May 2015