/*

Magic squares in B-Prolog.

Model created by Hakan Kjellerstrand, hakank@gmail.com

*/

% Licenced under CC-BY-4.0 : http://creativecommons.org/licenses/by/4.0/

%
% Reporting both time and backtracks
%
time2(Goal):-
cputime(Start),
statistics(backtracks, Backtracks1),
call(Goal),
statistics(backtracks, Backtracks2),
cputime(End),
T is (End-Start)/1000,
Backtracks is Backtracks2 - Backtracks1,
format('CPU time ~w seconds. Backtracks: ~d\n', [T, Backtracks]).

%
% Simple run
%
go :-
time2(magic(5,_Square)).

%
% Run for different sizes
%
% ff/updown solves N=6 in 0.06s and 6856 backtracks
%                   10 in 2.99s and 289939 backtracks
% ffd/updown solves N=10 in 2.37s and 306336 backtracks
go2 :-
foreach(N in 3..11, [Square], time2(magic(N,Square))).

%
% All solutions.
%
go3 :-
N = 4,
findall(Square,magic(N,Square),L),
length(L,Len),
writeln('Len':Len).

magic(N,Square) :-
format('\n\nN: ~d\n', [N]),
NN is N*N,
Sum is N*(NN+1)//2,% magical sum
format('Sum = ~d\n', [Sum]),

new_array(Square,[N,N]),
array_to_list(Square, Vars),
Vars :: 1..NN,

alldifferent(Vars),
foreach(I in 1..N,
(
Sum #= sum([ T : J in 1..N, [T], T @= Square[I,J]]),% rows
Sum #= sum([ T : J in 1..N, [T], T @= Square[J,I]]) % column
)
),

Sum #= sum(Square^diagonal1),% diagonal sums
Sum #= sum(Square^diagonal2),% diagonal sums

% Symmetry breaking
Square[1,1] #< Square[1,N],
Square[1,1] #< Square[N,N],
Square[1,1] #< Square[N,1],
Square[1,N] #< Square[N,1],

labeling([ffd,updown],Vars),
print_square(Square).

print_square(Square) :-
N is Square^length,
foreach(I in 1..N,
(
foreach(J in 1..N,
[S],
(
S @= Square[I,J],
format('~3d', [S])
)
),
nl
)
),
nl.