Unbiased grid generation: B2347 method

Programs which generate, solve, and analyze Sudoku puzzles

Unbiased grid generation: B2347 method

Postby Red Ed » Mon Jul 01, 2013 5:20 pm

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

Postby Afmob » Tue Jul 02, 2013 10:58 am

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

Postby Mathimagics » Sat Jan 09, 2016 7:01 am

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? :roll:

BTW: I get a bunch of "warnings: suggest parentheses around arithmetic" from gcc wrt these
User avatar
Mathimagics
2017 Supporter
 
Posts: 335
Joined: 27 May 2015

Re: Unbiased grid generation: B2347 method

Postby dobrichev » Sat Jan 09, 2016 7:50 am

(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: 1311
Joined: 24 May 2010

Re: Unbiased grid generation: B2347 method

Postby Mathimagics » Sat Jan 09, 2016 11:40 am

Ok, I forgot "^" is XOR operator.

So he means 1 - (next & 1), it seems!
User avatar
Mathimagics
2017 Supporter
 
Posts: 335
Joined: 27 May 2015


Return to Software