magicSquare(N) :-
	makeSquare(N, SQ),
	solve(N, SQ),
	printSquare(SQ).

solve(N, SQ) :-
	% Let Sum be the sum that all rows and cols must have
	getRows(SQ, Rows),
	getCols(SQ, N, Cols),
	setSumConstraint(Rows, Sum),
	setSumConstraint(Cols, Sum),
	getVars(SQ, Vars),
	fd_all_different(Vars),
	Nsq is N * N,
	fd_domain(Vars, 1, Nsq),
	fd_labeling(Vars).

getRows(SQ, SQ).

getCols(SQ, N, Cols) :-
	getCols(SQ, 0, N, Cols).

getCols(_, N, N, []) :- !.
getCols(SQ, I, N, [Col|Cols]) :-
	I1 is I + 1,
	getCols(SQ, I1, N, Cols),
	getCol(SQ, I, Col).

getCol([], _, []) :- !.
getCol([Row|Rows], I, [X|Xs]) :-
	getCol(Rows, I, Xs),
	getEle(Row, I, X).

getEle([X|_], 0, X) :- !.
getEle([_|Xs], I, X) :-
	I1 is I - 1,
	getEle(Xs, I1, X).

setSumConstraint([], _) :- !.
setSumConstraint([L|Ls], Sum) :-
	getSum(L, S), 
	S #= Sum,
	setSumConstraint(Ls, Sum).

getSum([], 0) :- !.
getSum([X|Xs], S) :-
	getSum(Xs, Ss),
	S = X + Ss.

getVars([], []).
getVars([Row|Rows], Vars) :-
	getVars(Rows, Vars1),
	append(Row, Vars1, Vars).


makeSquare(N, SQ) :-
	makeSquare(N, N, SQ).

makeSquare(_, 0, []) :- !.
makeSquare(N, RowsLeft, [Row|Rows]) :-
	makeRow(N, Row),
	RowsLeft1 is RowsLeft - 1,
	makeSquare(N, RowsLeft1, Rows).

makeRow(0, []) :- !.
makeRow(N, [_|Rest]) :-
	N1 is N - 1,
	makeRow(N1, Rest).

printSquare([]) :- !.
printSquare([Row|Rows]) :-
	printRow(Row), printSquare(Rows).

printRow([]) :- nl, !.
printRow([X|Xs]) :- write(X), write('\t '), printRow(Xs).

