/*
The University of York
CONSTRAINT SATISFACTION AND COMBINATORIAL OPTIMISATION
Open Examination

http://www-course.cs.york.ac.uk/cop/exam2000.pdf


An n-by-n magic square is an grid with n columns and n rows in which each element of
the grid contains a different integer and all rows, all columns and the main diagonals
sum to the same value. If S is a finite set of integers, a magic square for S is one in
which each element of the grid contains a member of S.

Write a program that given S determines whether or not there is a 3-by-3
magic square for S.
*/


/*
Solutions by George Rabanca, July, 2007
*/


magic_square1(S, Square):-
	Square = [A1, A2, A3, B1, B2, B3, C1, C2, C3], 
	
	Square in S, 
	all_different(Square), 
	
	A1 + A2 + A3 #= Sum, 
	B1 + B2 + B3 #= Sum, 
	C1 + C2 + C3 #= Sum, 
	A1 + B2 + C3 #= Sum, 
	C1 + B2 + A3 #= Sum, 
	A1 + B1 + C1 #= Sum, 
	A2 + B2 + C2 #= Sum, 
	A3 + B3 + C3 #= Sum, 
	
	labeling_ff(Square).
	
	
/*
We can make the program more efficient by breaking simmetry.
We will force the upper right corner to be smaller than all others to avoid rotation.
One more simmetry is left: rotate the square once and take the mirror immage.  Break this 
simmetry by choosing the upper right corner to be less than the lower left corner.
*/


magic_square2(S, Square):-
	Square = [A1, A2, A3, B1, B2, B3, C1, C2, C3], 
	
	Square in S, 
	all_different(Square), 
		
	min([A1, A3, C1, C3]) #= A1,
	A3 #< C1,
	
	A1 + A2 + A3 #= Sum, 
	B1 + B2 + B3 #= Sum, 
	C1 + C2 + C3 #= Sum, 
	A1 + B2 + C3 #= Sum, 
	C1 + B2 + A3 #= Sum, 
	A1 + B1 + C1 #= Sum, 
	A2 + B2 + C2 #= Sum, 
	A3 + B3 + C3 #= Sum, 
	
	labeling_ff(Square).
	
	
/*
A general magic square solution.
*/	

magic_square3(S, N, Square):-

	new_array(Square, [N, N]), 
	array_to_list(Square, L), 
	L in S,
	all_different(L),
		
	%calculate the sum of each row or column
	sum(L, SumTotal), 
	
	%Sum is N * (N*N + 1) // 2, 
	
	Sum #= SumTotal // N, 
	Sum * N #= SumTotal, 
	
	%all columns should sum to Sum
	Square^columns @= Columns, 
	all_sum_equal(Columns, Sum), 
	
	%all lines should sum to Sum
	Square^rows @= Rows, 
	all_sum_equal(Rows, Sum),
	
	%the two diagonals should also sum to Sum
	Square^diagonal1 @= Diag1, 
	Square^diagonal2 @= Diag2, 
	sum(Diag1, Sum), 
	sum(Diag2, Sum),
	
	%break simmetry
	Square^[1,1] @= A, 
	Square^[1,N] @= B, 
	Square^[N,N] @= C, 
	Square^[N,1] @= D, 
	min([A, B, C, D]) #= A,
	B #< D,
	
	labeling_ff(L),
	
	writeln(Square), fail.	
	

%tests that the sum of each list member of this list is Sum	
all_sum_equal([], _).
all_sum_equal([L|Ls], Sum):-
	sum(L, Sum), 
	all_sum_equal(Ls, Sum).
	
	
/*
A more efficient magic square general solution.
*/	

go:-
    statistics(runtime,[Start|_]),
    top,
    statistics(runtime,[End|_]),
    T is End-Start,
    write('execution time is '),write(T), write(milliseconds),nl.

top:-
    magic_square4([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25], 5, X).
top.


magic_square4(S, N, Square):-

	new_array(Square, [N, N]), 
	array_to_list(Square, L), 
	
	%get a subset of N*N elements to try to create the square with
	NN is N * N,
	eliminate_duplicate(S, S1),
	get_subset(TempElements, S1, NN), 
	L in TempElements,
	all_different(L),
		
	%calculate the sum of each row or column
	sum(TempElements, SumTotal), 
	
	Sum #= SumTotal // N, 
	Sum * N #= SumTotal, 
	
	%all columns should sum to Sum
	Square^columns @= Columns, 
	all_sum_equal(Columns, Sum), 
	
	%all lines should sum to Sum
	Square^rows @= Rows, 
	all_sum_equal(Rows, Sum),
	
	%the two diagonals should also sum to Sum
	Square^diagonal1 @= Diag1, 
	Square^diagonal2 @= Diag2, 
	sum(Diag1, Sum), 
	sum(Diag2, Sum),
	
	%break simmetry
	Square^[1,1] @= A, 
	Square^[1,N] @= B, 
	Square^[N,N] @= C, 
	Square^[N,1] @= D, 
	min([A, B, C, D]) #= A,
	B #< D,
	
	labeling_ff(L),
	
	writeln(Square).
	
	
%returns subsets of size N of Set; 
%pass the length of the remaining set for efficiency
get_subset(Subset, Set, N):-
	length(Set, Len),
	get_subset(Subset, Set, [], N, Len).
get_subset(Subset, [], Subset, 0, _).
get_subset(Subset, [X|Set], Acc, N, Len):-
	N1 is N - 1,
	Len1 is Len -1, 
	Acc1 = [X|Acc], 
	get_subset(Subset, Set, Acc1, N1, Len1).	
get_subset(Subset, [_|Set], Acc, N, Len):- 
	Len > N, 
	Len1 is Len -1,
	get_subset(Subset, Set, Acc,Len1).
	
	
%Sum is the sum of all elements in L	
sum(L, Sum):- sum(L, Sum, 0). 
sum([], Sum, Sum). 
sum([X|Xs], Sum, Acc):-
	Acc1 #= Acc + X, 
	sum(Xs, Sum, Acc1).