View source for Second lookup table C code

// building a new lookup table (storing the Pareto-optimal sets directly) from the old one

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#define NUM_LINES_3D 25   // 25 non-horizontal lines in [3]^3
// the bitmasks of the 25 lines in [3]^3
const long line_bitmasks[NUM_LINES_3D] =
{262657l,263172l,266304l,270592l,525314l,532608l,1049601l,1050628l,1056832l,1065216l,2101256l,2105376l,4202512l,8396808l,8405024l,16781313l,16785412l,16810048l,16843008l,33562626l,33620096l,67117057l,67125252l,67174464l,67240192l};


/* stat[i] is 0,1,2,3 depending on whether i (represented base 3) is of type a,b,c,d. */

const int stat[27] = {0,1,0, 1,2,1, 0,1,0,
                      1,2,1, 2,3,2, 1,2,1,
					  0,1,0, 1,2,1, 0,1,0};


long *data[9][13][7][2];  // data[a][b][c][d] will point to all the (a,b,c,d) Moser sets

long numstat[9][13][7][2];  /* numstat[a][b][c][d] is the number of Moser sets with statistics (a,b,c,d) */
int statlabel[9][13][7][2];  /* statlabel[a][b][c][d] is the index of (a,b,c,d) in the sorted list of 499 non-zero statistics*/


void init_numstat(void)  // This was precomputed
 {
 	int a,b,c,d,s[4];

	/* Initialise numstat */
	for (a=0; a<9; a++)
		for (b=0; b<13; b++)
			for (c=0; c<7; c++)
				for (d=0; d<2; d++)
					numstat[a][b][c][d]=0;

    numstat[0][0][0][0] = 1;
    numstat[0][0][0][1] = 1;
    numstat[0][0][6][0] = 1;
    numstat[0][12][0][0] = 1;
    numstat[8][0][0][0] = 1;
    numstat[4][12][0][0] = 2;
    numstat[0][0][1][0] = 6;
    numstat[0][0][1][1] = 6;
    numstat[0][0][5][0] = 6;
    numstat[4][0][5][0] = 6;
    numstat[0][0][3][1] = 8;
    numstat[0][6][6][0] = 8;
    numstat[0][9][3][0] = 8;
    numstat[1][0][0][0] = 8;
    numstat[1][0][0][1] = 8;
    numstat[1][0][6][0] = 8;
    numstat[1][12][0][0] = 8;
    numstat[3][12][0][0] = 8;
    numstat[5][9][0][0] = 8;
    numstat[7][0][0][0] = 8;
    numstat[7][3][0][0] = 8;
    numstat[0][0][2][1] = 12;
    numstat[0][1][0][0] = 12;
    numstat[0][1][0][1] = 12;
    numstat[0][1][6][0] = 12;
    numstat[0][11][0][0] = 12;
    numstat[6][0][2][0] = 12;
    numstat[0][0][2][0] = 15;
    numstat[0][0][4][0] = 15;
    numstat[2][0][6][0] = 16;
    numstat[2][12][0][0] = 16;
    numstat[4][0][0][1] = 16;
    numstat[6][6][0][0] = 16;
    numstat[0][0][3][0] = 20;
    numstat[0][10][1][0] = 24;
    numstat[2][0][0][1] = 24;
    numstat[4][11][0][0] = 24;
    numstat[4][9][2][0] = 24;
    numstat[5][0][3][0] = 24;
    numstat[6][0][1][0] = 24;
    numstat[7][1][0][0] = 24;
    numstat[7][2][0][0] = 24;
    numstat[2][0][0][0] = 28;
    numstat[6][0][0][0] = 28;
    numstat[3][0][0][1] = 32;
    numstat[3][9][3][0] = 32;
    numstat[4][0][3][1] = 32;
    numstat[0][5][6][0] = 48;
    numstat[0][9][2][0] = 48;
    numstat[1][0][1][0] = 48;
    numstat[1][0][1][1] = 48;
    numstat[1][0][5][0] = 48;
    numstat[3][0][5][0] = 48;
    numstat[4][1][5][0] = 48;
    numstat[4][4][5][0] = 48;
    numstat[4][6][2][1] = 48;
    numstat[6][3][2][0] = 48;
    numstat[6][4][1][0] = 48;
    numstat[0][2][6][0] = 54;
    numstat[4][0][1][1] = 54;
    numstat[3][0][0][0] = 56;
    numstat[5][0][0][0] = 56;
    numstat[0][2][0][1] = 60;
    numstat[4][0][4][0] = 60;
    numstat[6][1][2][0] = 60;
    numstat[0][6][0][1] = 64;
    numstat[0][6][3][1] = 64;
    numstat[1][0][3][1] = 64;
    numstat[1][6][6][0] = 64;
    numstat[1][9][3][0] = 64;
    numstat[0][10][0][0] = 66;
    numstat[0][2][0][0] = 66;
    numstat[4][0][0][0] = 70;
    numstat[0][1][1][0] = 72;
    numstat[0][1][1][1] = 72;
    numstat[0][1][5][0] = 72;
    numstat[4][0][2][1] = 72;
    numstat[2][6][6][0] = 80;
    numstat[0][1][3][1] = 96;
    numstat[0][8][3][0] = 96;
    numstat[1][0][2][1] = 96;
    numstat[1][1][0][0] = 96;
    numstat[1][1][0][1] = 96;
    numstat[1][1][6][0] = 96;
    numstat[1][11][0][0] = 96;
    numstat[5][0][2][0] = 96;
    numstat[5][4][3][0] = 96;
    numstat[5][7][1][0] = 96;
    numstat[5][8][0][0] = 96;
    numstat[6][2][2][0] = 96;
    numstat[2][9][3][0] = 104;
    numstat[3][0][3][1] = 104;
    numstat[0][4][6][0] = 108;
    numstat[0][6][5][0] = 108;
    numstat[0][7][4][0] = 108;
    numstat[2][0][5][0] = 108;
    numstat[6][5][0][0] = 108;
    numstat[0][3][6][0] = 112;
    numstat[3][6][3][1] = 112;
    numstat[1][0][2][0] = 120;
    numstat[1][0][4][0] = 120;
    numstat[3][11][0][0] = 120;
    numstat[4][6][1][1] = 120;
    numstat[5][0][1][0] = 120;
    numstat[5][6][2][0] = 120;
    numstat[6][1][1][0] = 120;
    numstat[2][0][1][1] = 132;
    numstat[4][2][5][0] = 132;
    numstat[4][6][4][0] = 132;
    numstat[0][1][2][1] = 144;
    numstat[2][0][3][1] = 144;
    numstat[3][0][1][1] = 144;
    numstat[4][1][0][1] = 144;
    numstat[4][3][5][0] = 144;
    numstat[4][9][1][0] = 144;
    numstat[2][0][1][0] = 156;
    numstat[6][1][0][0] = 156;
    numstat[0][3][0][1] = 160;
    numstat[1][0][3][0] = 160;
    numstat[4][10][0][0] = 162;
    numstat[5][1][3][0] = 168;
    numstat[6][3][1][0] = 168;
    numstat[0][1][2][0] = 180;
    numstat[0][1][4][0] = 180;
    numstat[2][1][6][0] = 180;
    numstat[0][5][0][1] = 192;
    numstat[1][10][1][0] = 192;
    numstat[3][10][1][0] = 192;
    numstat[4][4][3][1] = 192;
    numstat[2][11][0][0] = 204;
    numstat[4][0][3][0] = 212;
    numstat[0][6][1][1] = 216;
    numstat[0][6][2][1] = 216;
    numstat[0][9][1][0] = 216;
    numstat[3][0][2][1] = 216;
    numstat[4][6][0][1] = 216;
    numstat[6][2][1][0] = 216;
    numstat[0][3][0][0] = 220;
    numstat[0][9][0][0] = 220;
    numstat[0][1][3][0] = 240;
    numstat[0][4][0][1] = 240;
    numstat[2][0][2][1] = 240;
    numstat[3][0][4][0] = 240;
    numstat[4][0][1][0] = 246;
    numstat[3][0][1][0] = 264;
    numstat[3][6][5][0] = 264;
    numstat[4][1][3][1] = 264;
    numstat[2][1][0][1] = 276;
    numstat[4][7][3][0] = 288;
    numstat[2][0][4][0] = 300;
    numstat[6][4][0][0] = 300;
    numstat[3][9][2][0] = 312;
    numstat[2][1][0][0] = 324;
    numstat[0][2][5][0] = 336;
    numstat[3][1][0][1] = 336;
    numstat[4][0][2][0] = 336;
    numstat[0][2][1][1] = 348;
    numstat[2][0][2][0] = 360;
    numstat[4][5][2][1] = 360;
    numstat[6][2][0][0] = 360;
    numstat[0][2][1][0] = 384;
    numstat[0][5][3][1] = 384;
    numstat[1][5][6][0] = 384;
    numstat[1][9][2][0] = 384;
    numstat[2][10][1][0] = 384;
    numstat[4][8][2][0] = 384;
    numstat[5][3][3][0] = 384;
    numstat[5][2][3][0] = 408;
    numstat[0][2][3][1] = 432;
    numstat[1][2][6][0] = 432;
    numstat[3][7][4][0] = 432;
    numstat[5][1][0][0] = 432;
    numstat[2][0][3][0] = 440;
    numstat[6][3][0][0] = 440;
    numstat[4][1][1][1] = 456;
    numstat[0][8][2][0] = 480;
    numstat[1][2][0][1] = 480;
    numstat[3][8][3][0] = 480;
    numstat[3][0][3][0] = 488;
    numstat[0][4][0][0] = 495;
    numstat[0][8][0][0] = 495;
    numstat[0][5][5][0] = 504;
    numstat[3][0][2][0] = 504;
    numstat[3][1][5][0] = 504;
    numstat[5][7][0][0] = 504;
    numstat[1][6][0][1] = 512;
    numstat[1][6][3][1] = 512;
    numstat[4][1][4][0] = 516;
    numstat[1][10][0][0] = 528;
    numstat[1][2][0][0] = 528;
    numstat[2][5][6][0] = 528;
    numstat[4][2][0][1] = 540;
    numstat[1][1][1][0] = 576;
    numstat[1][1][1][1] = 576;
    numstat[1][1][5][0] = 576;
    numstat[3][1][0][0] = 600;
    numstat[4][1][2][1] = 600;
    numstat[5][5][2][0] = 648;
    numstat[4][1][0][0] = 660;
    numstat[0][2][2][1] = 672;
    numstat[4][5][4][0] = 672;
    numstat[5][1][2][0] = 696;
    numstat[2][9][2][0] = 720;
    numstat[5][6][1][0] = 720;
    numstat[4][2][3][1] = 744;
    numstat[4][5][1][1] = 744;
    numstat[2][2][6][0] = 756;
    numstat[0][3][5][0] = 768;
    numstat[0][7][3][0] = 768;
    numstat[1][1][3][1] = 768;
    numstat[1][8][3][0] = 768;
    numstat[2][6][3][1] = 768;
    numstat[4][9][0][0] = 772;
    numstat[0][5][0][0] = 792;
    numstat[0][5][1][1] = 792;
    numstat[0][7][0][0] = 792;
    numstat[4][3][3][1] = 792;
    numstat[4][5][0][1] = 792;
    numstat[0][6][4][0] = 816;
    numstat[3][10][0][0] = 816;
    numstat[0][3][1][1] = 864;
    numstat[0][4][3][1] = 864;
    numstat[1][4][6][0] = 864;
    numstat[1][6][5][0] = 864;
    numstat[1][7][4][0] = 864;
    numstat[0][2][4][0] = 870;
    numstat[0][8][1][0] = 870;
    numstat[5][1][1][0] = 888;
    numstat[0][4][5][0] = 894;
    numstat[0][3][3][1] = 896;
    numstat[1][3][6][0] = 896;
    numstat[3][6][0][1] = 896;
    numstat[0][6][0][0] = 924;
    numstat[0][2][2][0] = 930;
    numstat[0][5][2][1] = 1008;
    numstat[3][6][2][1] = 1008;
    numstat[3][1][3][1] = 1056;
    numstat[3][5][3][1] = 1056;
    numstat[4][3][0][1] = 1088;
    numstat[1][1][2][1] = 1152;
    numstat[2][6][0][1] = 1152;
    numstat[0][4][1][1] = 1158;
    numstat[2][10][0][0] = 1188;
    numstat[0][2][3][0] = 1200;
    numstat[0][3][1][0] = 1200;
    numstat[2][1][5][0] = 1224;
    numstat[4][4][0][1] = 1254;
    numstat[1][3][0][1] = 1280;
    numstat[2][4][6][0] = 1296;
    numstat[2][6][5][0] = 1296;
    numstat[2][2][0][1] = 1320;
    numstat[2][8][3][0] = 1344;
    numstat[2][7][4][0] = 1404;
    numstat[4][8][1][0] = 1416;
    numstat[1][1][2][0] = 1440;
    numstat[1][1][4][0] = 1440;
    numstat[4][4][2][1] = 1440;
    numstat[2][3][6][0] = 1456;
    numstat[3][2][0][1] = 1464;
    numstat[5][2][0][0] = 1464;
    numstat[3][1][1][1] = 1488;
    numstat[4][2][1][1] = 1488;
    numstat[2][1][1][1] = 1512;
    numstat[5][6][0][0] = 1512;
    numstat[0][3][2][1] = 1536;
    numstat[1][5][0][1] = 1536;
    numstat[3][5][5][0] = 1608;
    numstat[2][1][3][1] = 1632;
    numstat[4][2][4][0] = 1662;
    numstat[2][2][0][0] = 1716;
    numstat[1][6][1][1] = 1728;
    numstat[1][6][2][1] = 1728;
    numstat[1][9][1][0] = 1728;
    numstat[1][3][0][0] = 1760;
    numstat[1][9][0][0] = 1760;
    numstat[0][4][2][1] = 1788;
    numstat[2][1][1][0] = 1800;
    numstat[4][4][4][0] = 1806;
    numstat[3][6][1][1] = 1824;
    numstat[4][2][2][1] = 1824;
    numstat[5][4][2][0] = 1824;
    numstat[4][1][3][0] = 1896;
    numstat[1][1][3][0] = 1920;
    numstat[1][4][0][1] = 1920;
    numstat[4][4][1][1] = 1920;
    numstat[0][7][2][0] = 1932;
    numstat[5][2][2][0] = 1944;
    numstat[3][2][5][0] = 2016;
    numstat[0][7][1][0] = 2064;
    numstat[4][6][3][0] = 2096;
    numstat[3][9][1][0] = 2160;
    numstat[0][3][4][0] = 2172;
    numstat[3][1][2][1] = 2208;
    numstat[0][5][4][0] = 2232;
    numstat[4][1][1][0] = 2280;
    numstat[5][5][1][0] = 2352;
    numstat[4][3][1][1] = 2376;
    numstat[0][4][1][0] = 2436;
    numstat[4][3][2][1] = 2472;
    numstat[4][3][4][0] = 2496;
    numstat[3][1][4][0] = 2520;
    numstat[4][8][0][0] = 2634;
    numstat[4][7][2][0] = 2640;
    numstat[5][3][2][0] = 2640;
    numstat[1][2][5][0] = 2688;
    numstat[0][6][3][0] = 2696;
    numstat[0][3][2][0] = 2712;
    numstat[5][2][1][0] = 2712;
    numstat[2][1][2][1] = 2736;
    numstat[1][2][1][1] = 2784;
    numstat[4][2][0][0] = 2802;
    numstat[3][1][1][0] = 2808;
    numstat[5][3][0][0] = 2856;
    numstat[5][5][0][0] = 2856;
    numstat[3][2][0][0] = 2928;
    numstat[0][4][4][0] = 2988;
    numstat[2][6][2][1] = 3024;
    numstat[1][2][1][0] = 3072;
    numstat[1][5][3][1] = 3072;
    numstat[3][5][0][1] = 3072;
    numstat[4][1][2][0] = 3072;
    numstat[0][6][1][0] = 3192;
    numstat[0][3][3][0] = 3248;
    numstat[3][9][0][0] = 3320;
    numstat[0][5][1][0] = 3360;
    numstat[2][3][0][1] = 3360;
    numstat[3][3][0][1] = 3392;
    numstat[2][1][4][0] = 3420;
    numstat[1][2][3][1] = 3456;
    numstat[2][6][1][1] = 3456;
    numstat[5][4][0][0] = 3528;
    numstat[3][4][5][0] = 3600;
    numstat[2][5][0][1] = 3648;
    numstat[2][9][1][0] = 3672;
    numstat[1][8][2][0] = 3840;
    numstat[3][3][5][0] = 3840;
    numstat[3][2][3][1] = 3888;
    numstat[3][8][2][0] = 3936;
    numstat[1][4][0][0] = 3960;
    numstat[1][8][0][0] = 3960;
    numstat[3][4][3][1] = 3960;
    numstat[1][5][5][0] = 4032;
    numstat[2][1][2][0] = 4140;
    numstat[3][6][4][0] = 4176;
    numstat[2][9][0][0] = 4180;
    numstat[5][4][1][0] = 4224;
    numstat[0][6][2][0] = 4248;
    numstat[3][4][0][1] = 4416;
    numstat[5][3][1][0] = 4440;
    numstat[2][4][0][1] = 4800;
    numstat[0][4][2][0] = 4911;
    numstat[0][5][3][0] = 4920;
    numstat[2][5][3][1] = 4992;
    numstat[2][1][3][0] = 5040;
    numstat[3][7][3][0] = 5112;
    numstat[3][1][3][0] = 5136;
    numstat[0][4][3][0] = 5172;
    numstat[3][1][2][0] = 5328;
    numstat[1][2][2][1] = 5376;
    numstat[2][2][5][0] = 5376;
    numstat[2][3][0][0] = 5500;
    numstat[0][5][2][0] = 5712;
    numstat[3][5][2][1] = 6096;
    numstat[3][2][1][1] = 6120;
    numstat[1][3][5][0] = 6144;
    numstat[1][7][3][0] = 6144;
    numstat[3][3][3][1] = 6176;
    numstat[1][5][0][0] = 6336;
    numstat[1][5][1][1] = 6336;
    numstat[1][7][0][0] = 6336;
    numstat[4][7][1][0] = 6336;
    numstat[4][7][0][0] = 6384;
    numstat[1][6][4][0] = 6528;
    numstat[2][5][5][0] = 6552;
    numstat[4][2][3][0] = 6720;
    numstat[4][5][3][0] = 6720;
    numstat[1][3][1][1] = 6912;
    numstat[1][4][3][1] = 6912;
    numstat[2][2][3][1] = 6912;
    numstat[1][2][4][0] = 6960;
    numstat[1][8][1][0] = 6960;
    numstat[2][2][1][1] = 6960;
    numstat[4][3][0][0] = 7064;
    numstat[1][4][5][0] = 7152;
    numstat[1][3][3][1] = 7168;
    numstat[1][6][0][0] = 7392;
    numstat[1][2][2][0] = 7440;
    numstat[2][8][2][0] = 7680;
    numstat[3][5][1][1] = 8016;
    numstat[1][5][2][1] = 8064;
    numstat[3][2][2][1] = 8592;
    numstat[3][3][0][0] = 8600;
    numstat[3][8][0][0] = 9000;
    numstat[4][2][1][0] = 9132;
    numstat[2][2][1][0] = 9216;
    numstat[1][4][1][1] = 9264;
    numstat[1][2][3][0] = 9600;
    numstat[1][3][1][0] = 9600;
    numstat[4][6][2][0] = 9804;
    numstat[2][8][0][0] = 9900;
    numstat[3][2][4][0] = 10464;
    numstat[3][8][1][0] = 10776;
    numstat[4][6][0][0] = 11004;
    numstat[2][6][4][0] = 11424;
    numstat[2][3][5][0] = 11520;
    numstat[2][7][3][0] = 11520;
    numstat[4][2][2][0] = 11652;
    numstat[4][4][0][0] = 11742;
    numstat[2][4][0][0] = 11880;
    numstat[4][4][3][0] = 12000;
    numstat[2][2][2][1] = 12096;
    numstat[2][4][3][1] = 12096;
    numstat[4][3][3][0] = 12120;
    numstat[1][3][2][1] = 12288;
    numstat[2][4][5][0] = 12516;
    numstat[3][3][1][1] = 12768;
    numstat[3][2][1][0] = 13152;
    numstat[2][3][3][1] = 13440;
    numstat[2][5][1][1] = 13464;
    numstat[4][5][0][0] = 13512;
    numstat[3][4][1][1] = 14208;
    numstat[1][4][2][1] = 14304;
    numstat[3][4][2][1] = 14352;
    numstat[3][5][4][0] = 14688;
    numstat[2][5][2][1] = 15120;
    numstat[1][7][2][0] = 15456;
    numstat[2][2][4][0] = 15660;
    numstat[2][8][1][0] = 15660;
    numstat[3][3][2][1] = 15888;
    numstat[2][3][1][1] = 16416;
    numstat[1][7][1][0] = 16512;
    numstat[4][6][1][0] = 16536;
    numstat[2][7][0][0] = 16632;
    numstat[3][4][0][0] = 16920;
    numstat[3][7][0][0] = 17136;
    numstat[1][3][4][0] = 17376;
    numstat[1][5][4][0] = 17856;
    numstat[2][5][0][0] = 18216;
    numstat[1][4][1][0] = 19488;
    numstat[3][7][2][0] = 19896;
    numstat[2][6][0][0] = 20328;
    numstat[2][2][2][0] = 20460;
    numstat[4][3][1][0] = 20760;
    numstat[2][4][1][1] = 20844;
    numstat[1][6][3][0] = 21568;
    numstat[4][5][2][0] = 21576;
    numstat[1][3][2][0] = 21696;
    numstat[3][3][4][0] = 21936;
    numstat[3][2][3][0] = 22176;
    numstat[3][6][3][0] = 22904;
    numstat[3][5][0][0] = 23472;
    numstat[3][6][0][0] = 23520;
    numstat[1][4][4][0] = 23904;
    numstat[4][3][2][0] = 23904;
    numstat[3][2][2][0] = 23952;
    numstat[2][2][3][0] = 24000;
    numstat[3][4][4][0] = 24672;
    numstat[1][6][1][0] = 25536;
    numstat[1][3][3][0] = 25984;
    numstat[2][3][2][1] = 26112;
    numstat[1][5][1][0] = 26880;
    numstat[4][5][1][0] = 27360;
    numstat[2][3][1][0] = 27600;
    numstat[2][4][2][1] = 28608;
    numstat[4][4][2][0] = 29040;
    numstat[4][4][1][0] = 29550;
    numstat[3][7][1][0] = 31272;
    numstat[2][7][2][0] = 32844;
    numstat[2][5][4][0] = 33480;
    numstat[1][6][2][0] = 33984;
    numstat[3][3][1][0] = 35736;
    numstat[2][3][4][0] = 36924;
    numstat[2][7][1][0] = 39216;
    numstat[1][4][2][0] = 39288;
    numstat[1][5][3][0] = 39360;
    numstat[1][4][3][0] = 41376;
    numstat[2][6][3][0] = 43136;
    numstat[1][5][2][0] = 45696;
    numstat[2][4][4][0] = 47808;
    numstat[3][3][3][0] = 50912;
    numstat[3][5][3][0] = 52272;
    numstat[2][4][1][0] = 53592;
    numstat[3][6][2][0] = 54216;
    numstat[2][3][2][0] = 56952;
    numstat[3][6][1][0] = 58368;
    numstat[3][3][2][0] = 59952;
    numstat[2][3][3][0] = 61712;
    numstat[3][4][1][0] = 62400;
    numstat[2][6][1][0] = 63840;
    numstat[3][4][3][0] = 67416;
    numstat[2][5][1][0] = 70560;
    numstat[3][5][1][0] = 73176;
    numstat[2][6][2][0] = 76464;
    numstat[2][5][3][0] = 83640;
    numstat[3][5][2][0] = 88944;
    numstat[3][4][2][0] = 91824;
    numstat[2][4][3][0] = 93096;
    numstat[2][4][2][0] = 98220;
    numstat[2][5][2][0] = 108528;


// compute statlabel;

    int count=0;
	for (a=0; a<9; a++)
		for (b=0; b<13; b++)
			for (c=0; c<7; c++)
				for (d=0; d<2; d++)
					if (numstat[a][b][c][d])
						{
							statlabel[a][b][c][d] = count;
//							printf("%d:(%d %d %d %d) ", count,a,b,c,d);
							count++;
						}
					else
 						statlabel[a][b][c][d]=0;

 }


#define NUM_ACD 49 // number of interesting possible acd statistics for 3D Moser sets

int a_base[NUM_ACD], c_base[NUM_ACD], d_base[NUM_ACD];

void init_pareto(void) // this was precomputed
 {
  int i=0;

  a_base[i] = 2; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 5; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 6; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 7; c_base[i] = 0; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 5; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 6; c_base[i] = 1; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 5; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 6; c_base[i] = 2; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 3; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 3; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 3; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 3; d_base[i] = 0; i++;
  a_base[i] = 5; c_base[i] = 3; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 4; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 4; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 4; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 4; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 5; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 5; d_base[i] = 0; i++;
  a_base[i] = 3; c_base[i] = 5; d_base[i] = 0; i++;
  a_base[i] = 4; c_base[i] = 5; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 6; d_base[i] = 0; i++;
  a_base[i] = 2; c_base[i] = 6; d_base[i] = 0; i++;
  a_base[i] = 0; c_base[i] = 0; d_base[i] = 1; i++;
  a_base[i] = 2; c_base[i] = 0; d_base[i] = 1; i++;
  a_base[i] = 3; c_base[i] = 0; d_base[i] = 1; i++;
  a_base[i] = 4; c_base[i] = 0; d_base[i] = 1; i++;
  a_base[i] = 0; c_base[i] = 1; d_base[i] = 1; i++;
  a_base[i] = 2; c_base[i] = 1; d_base[i] = 1; i++;
  a_base[i] = 3; c_base[i] = 1; d_base[i] = 1; i++;
  a_base[i] = 4; c_base[i] = 1; d_base[i] = 1; i++;
  a_base[i] = 0; c_base[i] = 2; d_base[i] = 1; i++;
  a_base[i] = 2; c_base[i] = 2; d_base[i] = 1; i++;
  a_base[i] = 3; c_base[i] = 2; d_base[i] = 1; i++;
  a_base[i] = 4; c_base[i] = 2; d_base[i] = 1; i++;
  a_base[i] = 0; c_base[i] = 3; d_base[i] = 1; i++;
  a_base[i] = 2; c_base[i] = 3; d_base[i] = 1; i++;
  a_base[i] = 3; c_base[i] = 3; d_base[i] = 1; i++;
  a_base[i] = 4; c_base[i] = 3; d_base[i] = 1; i++;

 }


void init_data(void)
{

	int i,j,k;

	long NUM_SETS_3D = 1;
    for (i=0; i<27; i++) NUM_SETS_3D *= 2;

	// allocate space

	init_numstat();       // load statistics
    init_pareto();        // load pareto counts

	int aa, bb, cc, dd;

	for (aa=0; aa<9; aa++)
		for (bb=0; bb<13; bb++)
			for (cc=0; cc<7; cc++)
				for (dd=0; dd<2; dd++)
					if (numstat[aa][bb][cc][dd] > 0)
						{
							data[aa][bb][cc][dd] = malloc(numstat[aa][bb][cc][dd]*sizeof(long));
					    }


    /* enumerate the 230 line-free sets in [3]^2 by brute force */

    int moser_2d[512];
	int lincount=0;

    for (i=0; i < 512; i++)
     {
      moser_2d[i] = 0;
	  if ((i & 7) == 7) continue;
	  if ((i & 56) == 56) continue;
	  if ((i & 73) == 73) continue;
	  if ((i & 84) == 84) continue;
	  if ((i & 146) == 146) continue;
	  if ((i & 273) == 273) continue;
	  if ((i & 292) == 292) continue;
	  if ((i & 448) == 448) continue;
	  moser_2d[i] = 1;
	  lincount++;
	}

	printf("Number of 2D Moser sets: %d\n", lincount);

    /* Now enumerate all line-free sets in [3]^3 by brute force. */
    long set, tmpset;
	long slice1, slice2, slice3;
	int s[4];
	long temp_numstat[9][13][7][2];
	long num_linefree=0;

	/* Initialise temp_numstat */
	for (aa=0; aa<9; aa++)
		for (bb=0; bb<13; bb++)
			for (cc=0; cc<7; cc++)
				for (dd=0; dd<2; dd++)
					temp_numstat[aa][bb][cc][dd]=0;

    printf("Computing 3D Moser sets...\n");

	for (slice1 = 0; slice1 < 512; slice1++) if (moser_2d[slice1])
	 for (slice2 = 0; slice2 < 512; slice2++) if (moser_2d[slice2])
	  for (slice3 = 0; slice3 < 512; slice3++) if (moser_2d[slice3])
		{
			set = slice1 | (slice2 << 9) | (slice3 << 18);

			if ((set & line_bitmasks[0]) == line_bitmasks[0]) continue;
			if ((set & line_bitmasks[1]) == line_bitmasks[1]) continue;
			if ((set & line_bitmasks[2]) == line_bitmasks[2]) continue;
			if ((set & line_bitmasks[3]) == line_bitmasks[3]) continue;
			if ((set & line_bitmasks[4]) == line_bitmasks[4]) continue;
			if ((set & line_bitmasks[5]) == line_bitmasks[5]) continue;
			if ((set & line_bitmasks[6]) == line_bitmasks[6]) continue;
			if ((set & line_bitmasks[7]) == line_bitmasks[7]) continue;
			if ((set & line_bitmasks[8]) == line_bitmasks[8]) continue;
			if ((set & line_bitmasks[9]) == line_bitmasks[9]) continue;
			if ((set & line_bitmasks[10]) == line_bitmasks[10]) continue;
			if ((set & line_bitmasks[11]) == line_bitmasks[11]) continue;
			if ((set & line_bitmasks[12]) == line_bitmasks[12]) continue;
			if ((set & line_bitmasks[13]) == line_bitmasks[13]) continue;
			if ((set & line_bitmasks[14]) == line_bitmasks[14]) continue;
			if ((set & line_bitmasks[15]) == line_bitmasks[15]) continue;
			if ((set & line_bitmasks[16]) == line_bitmasks[16]) continue;
			if ((set & line_bitmasks[17]) == line_bitmasks[17]) continue;
			if ((set & line_bitmasks[18]) == line_bitmasks[18]) continue;
			if ((set & line_bitmasks[19]) == line_bitmasks[19]) continue;
			if ((set & line_bitmasks[20]) == line_bitmasks[20]) continue;
			if ((set & line_bitmasks[21]) == line_bitmasks[21]) continue;
			if ((set & line_bitmasks[22]) == line_bitmasks[22]) continue;
			if ((set & line_bitmasks[23]) == line_bitmasks[23]) continue;
			if ((set & line_bitmasks[24]) == line_bitmasks[24]) continue;

			/* now compute the a,b,c,d stats */
			s[0]=s[1]=s[2]=s[3]=0;
		    tmpset = set;
			for (j=0; j < 27; j++)
				{
					if (tmpset & 1l) s[stat[j]]++;
					tmpset /= 2;
				}

       		data[s[0]][s[1]][s[2]][s[3]][temp_numstat[s[0]][s[1]][s[2]][s[3]]] = set;

			temp_numstat[s[0]][s[1]][s[2]][s[3]]++;
			num_linefree++;
			if (num_linefree % 10000 == 0)  printf("*");
		}

	printf("\nNumber of 3D Moser sets: %ld\n", num_linefree);
// checksum
	for (aa=0; aa<9; aa++)
		for (bb=0; bb<13; bb++)
			for (cc=0; cc<7; cc++)
				for (dd=0; dd<2; dd++)
					if (numstat[aa][bb][cc][dd] != temp_numstat[aa][bb][cc][dd])
						printf("Inconsistency at %d %d %d %d!  %ld != %ld\n", aa,bb,cc,dd,numstat[aa][bb][cc][dd],temp_numstat[aa][bb][cc][dd]);
}





int perms[48][27];  // The 48 symmetries of the cube
long pow2[27];      // The 27 powers of 2

char table[144][16384][NUM_ACD];

void buildperms(void)
{
   int i,j,k,p,x,y,z,u,v,w;
   int n=0;
   long power=1;

   for (i=0; i<27; i++) { pow2[i] = power; power *= 2; }

   for (x=0;x<2;x++)
	   for (y=0;y<2;y++)
		   for (z=0;z<2;z++)
			   for (p=0;p<6;p++)
				   {					   ;
					   for (i=0;i<3;i++)
						   for (j=0;j<3;j++)
							   for (k=0;k<3;k++)
								   {
								   switch (p)
							   		{
										case 0: u=i; v=j; w=k; break;
										case 1: u=i; v=k; w=j; break;
										case 2: u=j; v=i; w=k; break;
										case 3: u=j; v=k; w=i; break;
										case 4: u=k; v=j; w=i; break;
										case 5: u=k; v=i; w=j; break;
								 	 }
									if (x) u = 2-u;
									if (y) v = 2-v;
									if (z) w = 2-w;
									perms[n][i+3*j+9*k] = u+3*v+9*w;
								}
						n++;
					}
}


long perm(long set, int i) // permutes set by permutation i
{
	long in,out;
	int j;

	      in = set; out = 0;
		  for (j=0; j < 27; j++)
			  {
				 if (in & 1) out |= pow2[perms[i][j]];
				 in >>= 1;
			  }

	  return out;
	}


int index(long set)  // gives the permutation index of set of minimal size
 {
	 long best = set, tmp;
	 int i,j=0;

	 for (i=0; i < 48; i++)
		 {
			 tmp = perm(set, i);
			 if (tmp < best)
				 {
					 best = tmp;
					 j = i;
				}
			}
	  return j;
 }


#define NUM_CRUSH 144

const int expand[NUM_CRUSH] = {0, 1, 3, 6, 7, 15, 19, 20, 21, 22, 23, 28, 29, 31, 54, 55, 56, 57, 58, 59, 62, 63, 96, 97, 99, 101, 102, 103, 111,
112, 113, 115, 116, 117, 118, 119, 124, 125, 127, 240, 241, 243, 246, 247, 255, 270, 271, 284, 285, 286, 287, 318, 319, 324, 325, 327, 332,
333, 334, 335, 343, 348, 349, 350, 351, 360, 361, 362, 363, 364, 365, 366, 367, 377, 379, 380, 381, 382, 383, 449, 451, 454, 455, 457, 459,
462, 463, 468, 469, 470, 471, 473, 475, 476, 477, 478, 479, 497, 498, 499, 502, 503, 504, 505, 506, 507, 510, 511, 876, 877, 879, 895, 911,
924, 925, 927, 938, 939, 942, 943, 958, 959, 995, 998, 999, 1005, 1007, 1011, 1013, 1014, 1015, 1020, 1021, 1023, 1782, 1783, 1785, 1787,
1791, 2013, 2015, 2046, 2047, 4095};  // minimal b-labels

int crush[4096];   // maps the label of a b-word to the label of its minimal permutation
int crushperm[4096];  // maps the label of a b-word to the permutation used to minimise it

const int bc[12] = { 1, 3, 5, 7, 9, 11, 15, 17, 19, 21, 23, 25 };  // the b points in [3]^3
const int ac[14] = { 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26 }; // the ac points in [3]^3

void buildhash(void)
{
   int i,j,k,count=0;
   long l,m;

   printf("Building hash...\n");
   for (i = 0; i<4096;i++)
	   {
		   j=i; l=0;
		   for (k=0;k<12;k++)
			   {
				   if (j & 1) l |= pow2[bc[k]];
				   j >>=1;
				}
			// l is now the b-word associated to i
		    crushperm[i] = index(l);
//			printf("%o (%lo) crushed by %d ", i, l, crushperm[i]);
		    m = perm(l, crushperm[i]);
//			printf("and permuted to %o ", m);

			j=0;
			for (k=0;k<12;k++)
				if (m & pow2[bc[k]]) j |= pow2[k];
			// j is now the label of the reduction of l

//			printf(" (i.e. %d)", j);
			for (k=0; k<NUM_CRUSH; k++)
				if (expand[k] == j) crush[i] = k;

			if (i == expand[crush[i]])
				{
					printf("%d ", i);
					count++;
				}

		}
  printf("Total mins: %d\n", count);
}

#define MAX_PARETO 27   // any forbidden set can generate at most 27 Pareto-optimals

FILE * newtable;

void read_lookup(void)
 {
	  FILE *input = fopen("lookup.dat", "r");

	  printf("Loading lookup table...\n");
	  int i, j, k;
	  for (i=0; i<144; i++)
		{
		  printf(".");
		  for (j=0; j<16384; j++)
			  fread(table[i][j],sizeof(char),NUM_ACD,input);
//			  for (k=0; k<NUM_ACD; k++)
//				  table[i][j][k] = (char) fgetc(input);
	    }
	  printf("\nTable loaded!\n");
 }

int test( long forbidden )  // see what statistics are compatible with a forbidden set.  Note that the 222 bit is ignored.
 {
  int i,j,k,l,a,b,x,y,z,w;
  long forb, recon;

//  printf("Testing %lo (= %ld in decimal)\n", forbidden, forbidden);

  a=0;
  for (k=0;k<12;k++)
	 if (forbidden & pow2[bc[k]]) a |= pow2[k];

  forb = perm(forbidden, crushperm[a]);
  a = crush[a];
//  forb = forbidden;  // can do this b/c we are only working with minimal sets right now

//  printf("Permuted version: %lo (= %ld in decimal), using permutation %d\n", forb, forb, crushperm[a]);


  b=0;

  for (k=0;k<14;k++)
	  if (forb & pow2[ac[k]]) b |= pow2[k];

//  printf("Hashes: %o %o \n", a, b );

//  recon = 0;
//  for (k=0;k<12;k++)
 //    if (expand[a] & pow2[k]) recon |= pow2[bc[k]];

//  for (k=0;k<14;k++)
//     if (b & pow2[k]) recon |= pow2[ac[k]];

//  printf("Reconstituted permutation: %lo (= %ld in decimal)\n", recon, recon);


  			/* now compute the a,b,c,d stats */

  int s[4];
  long tmpset = forbidden;

	s[0]=s[1]=s[2]=s[3]=0;

	for (j=0; j < 27; j++)
		{
			if (tmpset & 1l) s[stat[j]]++;
			tmpset /= 2;
		}


  int feasible[9][13][7][2];   // a partial list of which (a,b,c,d) are feasible

  for (i=0;i<9; i++)
	  for (j=0;j<13;j++)
		  for(k=0;k<7;k++)
			  for(l=0;l<2;l++)
				  feasible[i][j][k][l]=0;

  for (i=0; i<NUM_ACD;i++)
	  if (table[a][b][i] != -1)
		 {
//		  printf("[%d %d %d %d] ", a_base[i], (int) table[a][b][i], c_base[i], d_base[i]);
		  feasible[a_base[i]][(int) table[a][b][i]][c_base[i]][d_base[i]] = 1;
		 }

  // do the easy cases with only one non-zero

  feasible[8-s[0]][0][0][0] = 1;
  feasible[0][12-s[1]][0][0] = 1;
  feasible[0][0][6-s[2]][0] = 1;
  feasible[0][0][0][1-s[3]] = 1;

  // fill in the (1,b,c,d) guys

  if (s[0] < 8)
	  for (j=0; j<13;j++)
		  for (k=0;k<7;k++)
			  for (l=0;l<2;l++)
				  if (feasible[0][j][k][l]) feasible[1][j][k][l] = 1;

  // fill in the (0,1,c,d) guys

  if (s[1] < 12)
		  for (k=0;k<7;k++)
			  for (l=0;l<2;l++)
				  if (feasible[0][0][k][l]) feasible[0][1][k][l] = 1;

  // fill in the (0,0,1,1) guys

  if (s[2] < 6)
				  if (feasible[0][0][0][1]) feasible[0][0][1][l] = 1;


  // eliminate non-optimals


  for (i=0;i<9; i++)
	  for (j=0;j<13;j++)
		  for(k=0;k<7;k++)
			  for(l=0;l<2;l++)
				  if (feasible[i][j][k][l])
						  for (x=0; x<= i; x++)
							  for (y=0; y<= j; y++)
								  for (z=0; z<= k; z++)
									  for (w=0; w<=l; w++)
										  if (x+y+z+w < i+j+k+l) feasible[x][y][z][w] = 0;

  int count = 0;

  short int labels[MAX_PARETO];

  for (i=0;i<MAX_PARETO;i++) labels[i]=-1;

  for (i=0;i<9; i++)
	  for (j=0;j<13;j++)
		  for(k=0;k<7;k++)
			  for(l=0;l<2;l++)
				  if (feasible[i][j][k][l])
					 {
						 labels[count] = statlabel[i][j][k][l];
//						 printf("(%d,%d,%d,%d) ",i,j,k,l);
//						 printf("%d ", statlabel[i][j][k][l]);
						 count++;
					 }

  fwrite(labels, sizeof(short int), MAX_PARETO, newtable);

//  printf("\nTest complete.\n");
  return count;
 }


int main()
{
   int i,j,k,l,m;
   long forb, forb_base;

   printf("sizeof(short int)=%d ", sizeof(short int));

   srand ( time(NULL) );

   buildperms();
   buildhash();
   init_pareto();
   init_numstat();
//   init_data();
   read_lookup();

   newtable = fopen("newlookup.dat","wb");

	int maxcount=0, count;
	  for (i=0; i<144; i++)
		{
		  printf("%d ",maxcount);
		  for (j=0; j<16384; j++)
			  {
			   forb=0;
               for (k=0;k<12;k++)
                 if (expand[i] & pow2[k]) forb |= pow2[bc[k]];
               for (k=0;k<14;k++)
			     if (j & pow2[k]) forb |= pow2[ac[k]];
			   count = test(forb);
		       if (maxcount < count) maxcount=count;
		      }
		}

   fclose(newtable);
//   for (i=0; i<1000; i++)
//	   {
//		   forb = rand()%32768 + 32768l * (rand()%4096);
//		   printf("%lo -> %d ", forb, test(forb));
 //      }
//   test(110555621l);
}