function [ X, m ] = A264041( n, varargin ) %A264041 Use cplex to get optimal configuration on n x n grid % Returns optimal configuration X and number of occupied cells m % If optional second argument is present, produces graphical output % To produce the sequence A264041(1..N), use % A = zeros(1,N); % for n = 1:N % [~,A(n)] = A264041(n); % end % A % nvars = 2*n^2; % first n^2 constraints: x(ro,col,1) + x(ro,col,2) <= 1 Ai = []; Aj = []; conscount = 0; for ro=1:n for col=1:n conscount = conscount+1; Ai(end+[1,2]) = conscount; Aj(end+[1,2]) = [varnum(ro,col,1,n),varnum(ro,col,2,n)]; end end % next constraints: horizontal for ro = 1:n for col=1:n-1 conscount = conscount+2; Ai(end+[1:4]) = [conscount-1,conscount-1,conscount,conscount]; Aj(end+[1:4]) = [varnum(ro,col,1,n),varnum(ro,col+1,2,n),varnum(ro,col,2,n),varnum(ro,col+1,1,n)]; end end % next constraints: vertical for ro = 1:n-1 for col=1:n conscount = conscount+2; Ai(end+[1:4]) = [conscount-1,conscount-1,conscount,conscount]; Aj(end+[1:4]) = [varnum(ro,col,1,n),varnum(ro+1,col,2,n),varnum(ro,col,2,n),varnum(ro+1,col,1,n)]; end end % next constraints: nw to se for ro = 1:n-1 for col = 1:n-1 conscount = conscount+1; Ai(end+[1:2]) = conscount; Aj(end+[1:2]) = [varnum(ro,col,2,n),varnum(ro+1,col+1,2,n)]; end end % last constraints: ne to sw for ro = 1:n-1 for col = 2:n conscount = conscount+1; Ai(end+[1:2]) = conscount; Aj(end+[1:2]) = [varnum(ro,col,1,n),varnum(ro+1,col-1,1,n)]; end end As = ones(1,2*conscount); Aineq = sparse(Ai,Aj,As); bineq = ones(conscount,1); obj = -ones(nvars,1); [x,m, ~, output] = cplexbilp(obj,Aineq,bineq); if output.cplexstatus ~= 101 error(output.cplexstatusstring); end X = zeros(n,n); m = -m; for i=1:nvars if x(i) > 0.5 [ro,col,p] = unpackvar(i,n); X(ro,col) = 1 - 2*(p-1); end end if nargin == 2 fprintf('a(%d) = %d\n',n,m); disp(showconf(X)); end end function xn = varnum(ro,col,p,n) % variable number for x(ro,col,p) xn = (ro-1)*n + col + (p-1)*n^2; end function [ro,col,p] = unpackvar(xn,n) % row, column, p for variable #n if xn > n^2 p = 2; xp = xn - n^2; else p = 1; xp = xn; end col = mod(xp-1,n)+1; ro = (xp - col)/n+1; end function S = showconf(X) n = size(X,1); S = char(' '*ones(n,n)); S(X==1) = '/'; S(X==-1) = '\'; end