/*
     Magic Square

?- [-'magic.pl'].

?- magic3CLP(-Vars,-Milliseconds).     % 3x3 MagicSquare CLP Version
?- magic3(-Var,-Milliseconds).         % 3x3 MagicSquare Prolog Version
?- mc(4).       % 4x4 MagicSquare CLP Version.    Arg is 3,4,5...
?- mp(4).       % 4x4 MagicSquare Prolog Version. Arg is 3,4,5...
?- mc(4),fail.  % All solution
*/

:- dlib_require(clp).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
magic3(Vars,Time):-
  statistics(runtime,[Begin,_]),
    Vars   = [X1,X2,X3, X4,X5,X6, X7,X8,X9],    
    Domain = [1,2,3,4,5,6,7,8,9],
    alldifferent(Vars,Domain),
    15 =:= X1+X2+X3,
    15 =:= X4+X5+X6,
    15 =:= X7+X8+X9,
    15 =:= X1+X4+X7,
    15 =:= X2+X5+X8,
    15 =:= X3+X6+X9,
    15 =:= X1+X5+X9,
    15 =:= X3+X5+X7,
  statistics(runtime,[End,_]),
  Time is End-Begin.

%%
magic3CLP(Vars,Time):-
  statistics(runtime,[Begin,_]),
    Vars   = [X1,X2,X3, X4,X5,X6, X7,X8,X9],    
    Vars in 1..9,
    alldifferent(Vars),
    15 #= X1+X2+X3,
    15 #= X4+X5+X6,
    15 #= X7+X8+X9,
    15 #= X1+X4+X7,
    15 #= X2+X5+X8,
    15 #= X3+X6+X9,
    15 #= X1+X5+X9,
    15 #= X3+X5+X7,
  labeling(Vars),
  statistics(runtime,[End,_]),
  Time is End-Begin.

%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% NxN %%%%%%%%%%

magic(N,Vars,Time):-
  statistics(runtime,[Begin,_]),
    Sum is integer((N*N*N+N)/2),
    Max is N*N,
    gen_var_dom(1,Max,Vars,Domain),
    make_constraint('=:=',Vars,Vars,N,N,Sum,P,P,Goal),
    alldifferent(Vars,Domain),
   call(Goal),
  statistics(runtime,[End,_]),
  Time is End-Begin.

%%%
magicCLP(N,Vars,Time):-
  statistics(runtime,[Begin,_]),
    Sum is integer((N*N*N+N)/2),
    Max is N*N,
    gen_var_dom(1,Max,Vars,_),
    make_constraint('#=',Vars,Vars,N,N,Sum,P,P,Goal),
    Vars in 1..Max,
    alldifferent(Vars),
    call(Goal),
   labeling(Vars),
  statistics(runtime,[End,_]),
  Time is End-Begin.

%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Prolog Utility %%%%%%

my_select([A|L],A,L).
my_select([B|L],A,[B|R]):- my_select(L,A,R).

alldifferent([],_).
alldifferent([X|L],Dom):- my_select(Dom,X,Other),alldifferent(L,Other).
%%%
gen_var_dom(M,M,[X],[M]):-!.
gen_var_dom(S,M,[X|R],[S|L]):- SS is S+1,gen_var_dom(SS,M,R,L).

%% 
make_constraint(Func,[],[],_,_,_,_,_,true):-!.
make_constraint(Func,[],[A|Vars],_,N,Sum,_,_,Goal):-
	!,nth_list(N,N,[A|Vars],Ans,[A|Tail]),
	MM is N+1,nth(N,MM,MM,Vars,Tail,T),
	SS is N-1,skip(SS,Vars,T,TL,O),nth(N,SS,SS,O,TL,[]),
	make_constraint(Func,Ans,[],N,N,Sum,G,G,Goal).
make_constraint(Func,[V|L],Vars,1,N,Sum,V,Q,','(Goal,Else)):- 
	!,Goal=..[Func,Sum,Q],
	make_constraint(Func,L,Vars,N,N,Sum,E,E,Else).
make_constraint(Func,[V|L],Vars,M,N,Sum,+(P,V),Q,G):- 
	MM is M-1,make_constraint(Func,L,Vars,MM,N,Sum,P,Q,G). 
%% 
nth(1,_,_, L, T,T):-!.       % End of Line
nth(N,1,P,[A|L],[A|S],T):-!,NN is N-1,nth(NN,P,P,L,S,T).
nth(N,M,P,[_|L],A,T):-MM is M-1,nth(N,MM,P,L,A,T).

nth_list(0,_,_,L,L):-!.
nth_list(N,M,[A|Vars],[A|Else],Tail):-
	nth(M,M,M,Vars,Else,T),
	NN is N-1,nth_list(NN,M,Vars,T,Tail).

skip(1,[A|L],[A|T],T,L):-!.
skip(N,[_|L],S,T,O):- NN is N-1,skip(NN,L,S,T,O).

%%%%%
mc(N):- e,magicCLP(N,A,B),u,
	write_magic(N,N,A),nl,write([N,'x',N,' magic(CLP) ',B,msec]),nl,nl.
mp(N):- e,magic(N,A,B),u,
	write_magic(N,N,A),nl,write([N,'x',N,' magic(Prolog) ',B,msec]),nl,nl.

e:- abolish('$reg'/1),assert('$reg'(1)).
u:- retract('$reg'(N)),!,NN is N+1,assert('$reg'(NN)),write('No'=N),nl.

write_magic(_,_,[]):-!.
write_magic(0,N,L):- !,nl,write_magic(N,N,L).
write_magic(M,N,[A|L]):- 
	MM is M-1,write(A),((N>3,A<10) -> tab(2);tab(1)),write_magic(MM,N,L).

%%%%%%%%%%%% Omake %%%%%%%%%%%
%% combination of 3X3 is !(9,N)  ==>   N = 362880.0
%% combination of 4X4 is !(16,N) ==>   N = 20922789888000.0    20 trillion 

!(N,M):- kaijo(N,1.0,M).
kaijo(0,X,X):-!.
kaijo(N,S,X):- NN is N-1,SS is S*N,kaijo(NN,SS,X).
%%%%

