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.
