/* 

  Magic square integer programmingin Picat.

  From GLPK:s example magic.mod
  """
  MAGIC, Magic Square

  Written in GNU MathProg by Andrew Makhorin <mao@mai2.rcnet.ru>

  In recreational mathematics, a magic square of order n is an
  arrangement of n^2 numbers, usually distinct integers, in a square,
  such that n numbers in all rows, all columns, and both diagonals sum
  to the same constant. A normal magic square contains the integers
  from 1 to n^2.

  (From Wikipedia, the free encyclopedia.)
  """

  This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import cp.


main => go.

go ?=>
  nolog,
  % square order
  N = 4,
  magic_square(N, Square,S),

  println(s=S),
  foreach(I in 1..N)
    println(Square[I])
  end,
  nl,
  % fail,
  nl.

go => true.


go2 =>
  nolog,
  foreach(N in 3..5)
    println(n=N),
    time2(magic_square(N, Square, S)),
    println(s=S),
    foreach(I in 1..N)
      println(Square[I])
    end,
    nl
    
  end,
  nl.


magic_square(N, Square,S) => 

  % integers to be placed
  N2 = N*N,

  % x[i,j,k] = 1 means that cell (i,j) contains integer k
  X = new_array(N,N,N2),
  X :: 0..1,

  Square = new_array(N,N),
  Square :: 1..N2,


  S :: 1..N*N*N, % the magic sum

  % each cell must be assigned exactly one integer
  foreach(I in 1..N, J in 1..N)
    sum([X[I,J,K] : K in 1..N2]) #= 1
  end,

  % each integer must be assigned exactly to one cell
  foreach(K in 1..N2) 
    sum([X[I,J,K] : I in 1..N, J in 1..N]) #= 1
  end,

  % the sum in each row must be the magic sum 
  foreach(I in 1..N) 
    sum([K * X[I,J,K] : J in 1..N, K in 1..N2]) #= S
  end,

  % the sum in each column must be the magic sum
  foreach(J in 1..N) 
    sum([K * X[I,J,K] : I in 1..N, K in 1..N2]) #= S
  end,

  % the sum in the diagonal must be the magic sum
  sum([K * X[I,I,K] : I in 1..N, K in 1..N2]) #= S,

  % the sum in the co-diagonal must be the magic sum
  sum([K * X[I,N-I+1,K] : I in 1..N, K in 1..N2]) #= S,

  % for output
  foreach(I in 1..N, J in 1..N)
    Square[I,J] #= sum([K * X[I,J,K] : K in 1..N2]) 
  end,


  Vars = X.vars ++ Square.vars, % ++ [S],
  solve([ff,split],Vars).