function magicsquare = evenmagic(dim)
% EVENMAGIC magic square of of even numbered dimensions
%
% a = evenmagic(dim) produces a magic square of dimensions dim.
% evenmagic returns an error if dim is an odd number. for odd numbered
% magic squares use evenmagic. the smallest dimension accepted is 4.
%
% Algorithm developed by Sri. Samavedula Sita Rama Sastry
% http://www.me.berkeley.edu/~srinath/magic

if ( mod(dim, 2) ~= 0 )
	error('Error: evenmagic only works for odd numbers');
end
if ( dim < 4 )
	error('Error: Smallest even sided magic square possible is 4 x 4');
end

if ( dim == 4 )
	magicsquare = magic(4);
	return;
end


prenum = 1;

row = 2;
col = 1;

magicsquare(2, 1) = prenum;

if (mod(dim, 4) == 0)
	div4 = 1;
else
	div4 = 0;
end

if (~div4)
	magicsquare(1,2) = 2;
	prenum = 3;
else
	prenum = 2;
end

col = dim;
row = 3;

for i=prenum:2:prenum + (dim - 5)
	magicsquare(row, col) = i;
	magicsquare(row + 1, col) = i+1;
	row = row + 2;
	col = dim + 1 - col;
	prenum = prenum + 2;
end

magicsquare(row, col) = prenum;

lastsmall = 2*dim-2;
magicsquare(1,1) = lastsmall;
magicsquare(1,dim) = lastsmall - 1;

prenum = prenum + 1;

remain = [prenum:1:lastsmall-2];

reqdlowsum = sum([prenum:1:lastsmall])/2;
if ( ~div4 )
	reqdlowsum = reqdlowsum + 1;
end

%printf('reqdlowsum = %g', reqdlowsum);

if ( div4 )
	nlow = length(prenum:lastsmall)/2;
else
	nlow = (length(prenum:lastsmall) + 1)/2;
end

[lownums, hinums] = distribute(remain, reqdlowsum, nlow);

%printf('reqdlowsum = %d, nlow = %d', reqdlowsum, nlow);
%printf('remain :');
%disp(remain);
%printf('lownums :');
%disp(lownums);
%printf('hinums :');
%disp(hinums);

if ( div4 )
	magicsquare(dim, 2:1+length(lownums)) = lownums;
	col = 1+length(lownums);
	magicsquare(1,col+1:col+length(hinums)) = hinums;
else
	magicsquare(dim, 3:2+length(lownums)) = lownums;
	col = 2+length(lownums);
	magicsquare(1,col+1:col+length(hinums)) = hinums;
end

outersum = 1+dim^2;
for i=2:dim-1
	emptyrow = dim;
	if (magicsquare(1,i) == 0) 
		emptyrow = 1;
	end
	magicsquare(emptyrow, i) = outersum - magicsquare(dim-emptyrow+1, i);
	emptycol = dim;
	if (magicsquare(i,1) == 0) 
		emptycol = 1;
	end
	magicsquare(i, emptycol) = outersum - magicsquare(i, dim-emptycol+1);
end
magicsquare(dim,1) = outersum - magicsquare(1,dim);
magicsquare(dim,dim) = outersum - magicsquare(1,1);

magicsquare(2:dim-1, 2:dim-1) = evenmagic(dim-2) + lastsmall;
	
if ( length(unique(magicsquare)) ~= dim^2 )
	printf('Magic Square Error: not all numbers from %g to %g uniquely used', 1, dim^2);
	return;
end