Matlab script

The script given below, written in Matlab, is intended to carry it out. An input, DVec, is a list of 27-bit numbers. Each number represents one case of the forbidden 2*** points. The output is ParetoCell, which lists the Pareto statistics for each DVec value. A more compact form of the output is ParetoInd.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DATE : 6 June, 2009
% AUTHOR : Michael Peake
% PROJECT : Polymath1
% Lookup table
% There is assumed a vector of 27-bit DISALLOWED numbers in DVEC
% Its output is PARETOCELL which lists the Pareto statistics
% for each DISALLOWED number
if ~exist('eee');

% List of Disallowed 27-bit numbers
if ~exist('DVec')
DVec = 0;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Find the 230 line-free subsets of ^2
[a,b,c,d,e,f,g,h,k] = ndgrid(0:1);
abcd = [a(:) b(:) c(:) d(:) e(:) f(:) g(:) h(:) k(:)];
All512 = abcd;
% Delete those with complete lines in them
for n=512:-1:1,
v=abcd(n,:);
if any(all(v([1 2 3;4 5 6;7 8 9;1 4 7;2 5 8;3 6 9;1 5 9;3 5 7]),2)),
abcd(n,:)=[];
end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some numbers needed repeatedly, so precalculate them
TwoEight = 2.^[0:8];
Eight = abcd*TwoEight';
TwoFourEight = 512*Eight;
% Student Matlab won't allow 230x230 array, but will allow cell arrays
Eights = cell(1,230);
% 27-bit numbers built from first layer and third layer
for n=1:230,Eights{n}=Eight(n)+512*512*Eight;end;

% find symmetry group of each one
%for n=1:230,
%v=reshape(abcd(n,:),3,3);
%w=v';
%x=flipud(v);
%y=fliplr(v);
%z=flipud(fliplr(v));
%a=flipud(w);
%b=fliplr(w);
%c=flipud(fliplr(w));
%e=[v(:) w(:) x(:) y(:) z(:) a(:) b(:) c(:)];
%same(n) = sum(all(e==(v(:)*ones(1,8))));
%end;

% calculate the statistics for each line-free square
for n=1:230,
v=reshape(abcd(n,:),3,3);
Stats2(n,:) = [sum(v([1 3 7 9])) sum(v([2 4 6 8])) v(5)];
end;

% Convert the stats to a unique index
Stats2Ind = Stats2*[1;9;9*13];
% The centre-slice has a different statistic
Stats2Shift = Stats2*[9;9*13;9*13*7];

% count the restricted patterns
for n=1:512,
WithinCell{n} = find(all(abcd <= All512(n*ones(1,230),:),2));
Within(n)     = length(WithinCell{n});
% Sort them so that Stats2Shift is increasing;
% will make it easier to ignore repetitions later
[a,b] = sort(Stats2Shift(WithinCell{n}));
WithinCell{n} = WithinCell{n}(b);
end;
Total = 0;Count = zeros(1,48);

% Compare every Statistics for the cube with every other Statistics,
% to see if one dominates the other - to keep track of Pareto maxima
Beaters = cell(9,13,7,2);
for a=1:9,for b=1:13,for c=1:7,for d=1:2,
Beaters{a,b,c,d} = zeros(9,13,7,2);
for e=1:9,
for f=1:13,
for g=1:7,
for h=1:2,
if all([a b c d]>=[e f g h]),
Beaters{a,b,c,d}(e,f,g,h)=1;
elseif all([a b c d]<=[e f g h])
Beaters{a,b,c,d}(e,f,g,h)=-1;
end;
end;end;end;end;
Beaters{a,b,c,d}(a,b,c,d)=-1;
end;end;end;end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For each Layer1 and Layer3, precalculate which cells
% are forbidden in Layer2
for a=1:230,
% Convert square of 01s to 9-bit number
Layer1 = abcd(a,:);
Num1 = TwoEight*Layer1';
% Some bits are needed for sloping lines
X123 = bitand(Num1,7)*8;
X789 = bitand(Num1,8*56)/8;
X147 = bitand(Num1,73)*2;
X369 = bitand(Num1,4*73)/2;

for b=1:230,

Layer3 = abcd(b,:);
Num2 = TwoEight*Layer3';

% Check for vertical lines
Out = bitand(Num1,Num2);

% Set up for sloping lines
Y123 = bitand(Num2,7);Y123=Y123*8;
Y789 = bitand(Num2,448);Y789=Y789/8;
Y147 = bitand(Num2,73);Y147=Y147*2;
Y369 = bitand(Num2,292);Y369=Y369/2;

% Check for sloping lines
Out = bitor(Out,bitand(X123,Y789));
Out = bitor(Out,bitand(X789,Y123));
Out = bitor(Out,bitand(X147,Y369));
Out = bitor(Out,bitand(X369,Y147));

% Check for major diagonals, and bit 5
v = 16*any(Layer1([1 3 7 9]) & Layer3([9 7 3 1]));
Out = bitor(Out,v);

NotOut = bitcmp(Out,9);
% Index
In=NotOut+1;

% 9-bit number shows which bits are allowed
Allowed{a}(b) = In;

end;end;%for a,for b

%%%%%%%%%%%%%%%%%%%%%%%%%
eee = 'done';
end;% if ~exist('eee');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up enough memory to keep track of Pareto sets
% These are the pareto-maxima of the second slice,
% will be subject to the DISALLOWED restriction

% Count of Disallowed 27-bit numbers
nd = length(DVec);

% List of Pareto-max Layer2s
ParetoInd = cell(1,nd);

% Provide enough space
[ParetoInd{:}] = deal(ones(1,230));

% Initialize
[ParetoInd{:}] = deal(ones(1));
NPareto = ones(1,nd);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DATE : 6 June, 2009
% AUTHOR : Michael Peake
% PROJECT : Polymath1
% Lookup table
% There is assumed a vector of 27-bit DISALLOWED numbers in DVEC
% Its output is PARETOCELL which lists the Pareto statistics
% for each DISALLOWED number
for a=1:230,

InFirst = Allowed{a};
NN0Vec = Eights{a};

for b=1:230,

% Statistics of Layer1 and Layer3 put together
Stats0Ind = Stats2Ind(a)+Stats2Ind(b)+1;

% List of possible Layer2s
Rvec0 = WithinCell{InFirst(b)};

% Stats of all three Layers put together
Stats3Ind0 = Stats0Ind+Stats2Shift(Rvec0);

% 27-bit number of Layer1 and Layer3
NN0 = NN0Vec(b);

% Pick which Disallowed are dead already
DOk = find(~bitand( NN0, DVec));

% 27-bit number including various possible Layer2s
NN = NN0+TwoFourEight(Rvec0);

% Only the sequel changes for different DISALLOWEDS
for di = DOk

% Pick a 27-bit number of DISALLOWED bits
Disallowed = DVec(di);

% Which Layer2s give a cube entirely within allowed bits
NOk = find(~bitand(NN,Disallowed));

% Don't bother if none apply
if length(NOk)

Stats3Ind = Stats3Ind0(NOk);

% Remove repeats (values are already sorted)
for c=[1; 1+find(diff(Stats3Ind))]',

NewStat = Stats3Ind(c);
% To compare this new Statistic with current Pareto set, look up the
% pre-calculated comparisons
Mat = Beaters{NewStat};
% Fetch the list of Pareto sets, to prevent repeated indexing
PInd = ParetoInd{di};
Comparison = Mat(PInd);

if (all(Comparison>=0)),
% Then this is a new Pareto set

% remove the dominated Pareto sets
Old = find(Comparison);
PInd(Old) = [];
% include the new Pareto set
PInd(end+1) = NewStat;
NPareto(di) = length(PInd);
% Put the updated Pareto list back in the cell array
ParetoInd{di} = PInd;
end;%if Comparison
end;%for c[1; 1+find(diff(Stats3Ind))]
end;%if length(NOk)
end;%for di
%disp(num2str([floor(toc) a b Total Count(find(Count))]));
end;%for b
%disp(num2str([a NPareto]));
end;% for a

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert from Pareto indexes to lists of Pareto statistics

for di = 1:nd,
[a,b,c,d] = ind2sub([9 13 7 2],ParetoInd{di}(1:NPareto(di)));
ParetoCell{di} = [a;b;c;d]'-1;
end;