/* 

  Square sum problem in Picat.

  http://community.wolfram.com/groups/-/m/t/1264240
  """
  The question is: if you have the integers 1 through n, can you arrange that list in
  such a way that every two adjacent ones sum to a square number. 
  """

  Also see 
  * Numberphile: The Square-Sum Problem: https://www.youtube.com/watch?v=G1m7goLCJDY
  * Numberphile: The Square-Sum Problem (extra footage) https://www.youtube.com/watch?v=7_ph5djCCnM

  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 ?=>
  N = 15,
  square_sum(N,X),
  println(X),
  
  fail,
  
  nl.

go => true.


go2 => 
  nolog,
  foreach(N in 1..25) 
    All=find_all(X,square_sum(N,X)),
    println(N=All.len=All)
  end,
  nl.



go3 => 
  nolog,
  foreach(N in 1..125) 
    if square_sum(N,X) then
        println(N=X)
    end
  end,
  nl.

% [8,1,15,10,6,3,13,12,4,5,11,14,2,7,9]
go4 =>
  N = 15,
  Squares = [I*I : I in 1..N+1],
  Dom = new_list(N),
  foreach(I in 1..N) 
    Dom[I] := [J : J in 1..N, member(I+J,Squares)]
  end,

  X = new_list(N),
  X :: 1..N,
  Path = new_list(N),
  Path :: 1..N,

  foreach(I in 1..N) 
    println(i=I=Dom[I]),
    X[I] :: Dom[I]
  end,
  % Path[1] #< Path[N],
  all_different(X),
  all_different(Path),

  hamiltonian_path(N,X,Path),

  Vars = X ++ Path,
  solve($[ff,split],Vars),

  println(x=X),   
  println(path=Path),
  % fail,

  nl.

square_sum(N, X) =>
  X = new_list(N),
  X :: 1..N,
  all_distinct(X),
  X[1] #< X[N],
  foreach(I in 1..N-1)
     Y :: 0..N**2,
     X[I]+X[I+1] #= Y*Y
  end,

  solve([degree,split],X).


% This is based on my circuit_path/2 decomposition
% See http://hakank.org/picat/circuit.pi
%
hamiltonian_path(N,X,Z) =>
   println($hamiltonian_path(N,X,Z)),
   % Z = new_list(N),
   % Z :: 1..N,

   %
   % The main constraint is that Z[I] must not be 1 
   % until I = N, and for I = N it must be 1.
   %
   all_different(X),
   all_different(Z),
   % all_distinct(X), % slower
   % all_distinct(Z), % slower

   % put the orbit of x[1] in in z[1..n]
   % X[1] #= Z[1],
   
   % when I = N it must be 1
   % Note: This is not needed for a hamiltonial path
   % Z[N] #= 1,

   %
   % Get the orbit of Z.
   %
   foreach(I in 2..N)
     println(i=I),
     element(Z[I-1],X,Z[I])
   end.