/* 

  Samuraï puzzle in Picat.

  From Constraint 2011 exam (Polytechnique)
  via Nicolas Beldiceanu.
  
  For every sqrt(n) by sqrt(n) sub square there should be 
  exacly one 1, also there should be just one 1 on 
  each line/column.

  Here are all 16 possible solutions for n=4
  [0,2,1,3],
  [0,2,3,1],
  [0,3,1,2],
  [0,3,2,1],
  [1,2,0,3],
  [1,2,3,0],
  [1,3,0,2],
  [1,3,2,0],
  [2,0,1,3],
  [2,0,3,1],
  [2,1,0,3],
  [2,1,3,0],
  [3,0,1,2],
  [3,0,2,1],
  [3,1,0,2],
  [3,1,2,0]]

  Here is one solution for N=4:
   x = [1,3,2,4]
   1 _  _ _  
   _ _  3 _  

   _ 2  _ _  
   _ _  _ 4  


   1 0  0 0  
   0 0  1 0  

   0 1  0 0  
   0 0  0 1  




  Here is one (of many many) solutions for N=25:
    x = [1,6,11,16,21,2,7,12,17,22,3,8,13,18,23,4,9,14,19,24,5,10,15,20,25]
    1 _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  6 _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  11 _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  16 _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  21 _ _ _ _  

    _ 2 _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ 7 _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ 12 _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ 17 _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ 22 _ _ _  

    _ _ 3 _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ 8 _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ 13 _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ 18 _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ 23 _ _  

    _ _ _ 4 _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ 9 _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ 14 _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ 19 _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ 24 _  

    _ _ _ _ 5  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ 10  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ 15  _ _ _ _ _  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ 20  _ _ _ _ _  
    _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ _  _ _ _ _ 25  


  Using symmetry breaking with alldifferent_interval reduces
  the number of solutions. See go2/0 for more on this.

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

*/

import cp.

main => go.

go =>
  N = 4, % must be a square
  Symmetry = false,
  samurai_puzzle(N,Symmetry, X,Y),
  print_solution(N,X,Y),
  nl,
  fail,  
  nl.


go2 =>
  member(Symmetry,[false,alldifferent_interval,alldifferent_interval2]),
  println(symmetry=Symmetry),
  N = 4, % must be a square
  samurai_puzzle(N,Symmetry, X,Y),
  print_solution(N,X,Y),
  nl,
  fail,  
  nl.

/*

  Number of solutions for n
   n  #solutions
   -------------
   1  1
   4  16
   9  46656
   25 (too many solutions...)


  It seems to be http://oeis.org/A185141
  1, 16, 46656, 110075314176, 619173642240000000000,...
   (n!)^(2n)
  a(n) is the number of "templates", or ways of placing a single digit within an 
  n^2*n^2 Sudoku puzzle so that all rows, columns, and n*n blocks have exactly one 
  copy of the digit.
  Reference: G. Dahl, "Permutation matrices related to Sudoku, Linear Algebra and its 
  Applications", 430 (2001), 2457-2463.

  Number of solutions with the alldifferent_interval(X, M):
   n  #solutions
   -------------
   1  1
   4  8
   9  1296
   16 = 7962624

  This seems to be https://oeis.org/A091868
  "a(n) = (n!)^(n+1). 
  Let f(x) be a monic polynomial of degree n. Let u be any number and let m be the matrix 
  of values f(u+i-j) for i,j=1..n. Then the determinant of m is a(n)."

  Number of solutions with alldifferent_interval2(X, M)
   n  #solutions
   -------------
   1  1
   4  2
   9  2000
   25 
  
   Not in OEIS.


*/
go3 =>
  member(Symmetry,[alldifferent_interval,alldifferent_interval2,false]),
  println(symmetry=Symmetry),
  foreach(N in [I*I : I in 1..4])
    Count = count_all(samurai_puzzle(N,Symmetry, _X,_Y)),
    println(N=Count)
  end,
  fail,
  nl.


print_solution(N,X,Y) =>
  M = round(sqrt(N)),
  println(x=X),
  % println(y=Y),
  foreach(I in 1..N)
    foreach(J in 1..N)
      if X[I] == J then
        printf("%d ",J)
      else
        print("_ ")
      end,
      if J mod M == 0 then
        print(" ")
      end
    end,
    nl,
    if I mod M == 0 then
      nl
    end
  end,
  nl,
  foreach(I in 1..N)
    foreach(J in 1..N)
      printf("%d ",Y[I,J]),
      if J mod M == 0 then
        print(" ")
      end
    end,
    nl,
    if I mod M == 0 then
      nl
    end    
  end,
  
  nl.

samurai_puzzle(N,Symmetry, X,Y) =>

  M = ceiling(sqrt(N)),
  
  % decision variables
  X = new_list(N),
  X :: 1..N,

  % 0/1 version of the grid
  Y = new_array(N,N),
  Y :: 0..1,

  all_different(X),

  % channeling between x and y
  foreach(I in 1..N, J in 1..N) 
    Y[I,J] #= 1 #<=> X[I] #= J
  end,
  
  % each sub square has just one 1 (the rest is 0's)
  foreach(I in 0..M-1, J in 0..M-1) 
    sum([Y[R,C] : R in I*M+1..I*M+M, C in J*M+1..J*M+M]) #= 1
  end,

  if Symmetry == alldifferent_interval then
    alldifferent_interval(X, M)
  end,
  
  if Symmetry == alldifferent_interval2 then
    alldifferent_interval2(X, M)
  end,

  Vars = X ++ Y.vars,
  solve(Vars).



%
% alldifferent_interval(x, m)
%
% which requires that the differences of all elements in the mxm 
% subsquares should be >= m
%
% For n=4 it gives these 8 solutions
% x: [0, 2, 1, 3]
% x: [0, 2, 3, 1]
% x: [1, 3, 0, 2]
% x: [1, 3, 2, 0]
% x: [2, 0, 1, 3]
% x: [2, 0, 3, 1]
% x: [3, 1, 0, 2]
% x: [3, 1, 2, 0]
% 
alldifferent_interval(X, M) =>
  foreach(I in 0..M-1)
    foreach(J in I*M+1..I*M+M, K in J+1..I*M+M)
      abs(X[J]-X[K]) #>= M
    end
  end.


%
% alldifferent_interval2(x, m)
%
% A variant where the difference of all _neighboring_ elements 
% must be >= m.
%
% For n=4 it gives these two solutions:
% x: [1, 3, 0, 2]
% x: [2, 0, 3, 1]
%
alldifferent_interval2(X, M) =>
  foreach(I in 1..X.len-1)
    abs(X[I]-X[I+1]) #>= M
  end.