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 [3]^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;