# Matlab script

From Polymath1Wiki

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;