024: Langford's number problem

/*

  Langford's number problem in Picat.

  Langford's number problem (CSP lib problem 24)
  http://www.csplib.org/prob/prob024/
  """
  Arrange 2 sets of positive integers 1..k to a sequence,
  such that, following the first occurence of an integer i, 
  each subsequent occurrence of i, appears i+1 indices later
  than the last. 
  For example, for k=4, a solution would be 41312432
  """
  
  * John E. Miller: Langford's Problem
    http://www.lclark.edu/~miller/langford.html
  
  * Encyclopedia of Integer Sequences for the number of solutions for each k
    http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=014552


  Note: For k=4 there are two different solutions:
     solution:[4,1,3,1,2,4,3,2]
     position:[2,5,3,1,4,8,7,6]
  and
     solution:[2,3,4,2,1,3,1,4]
     position:[5,1,2,3,7,4,6,8]

  Which this symmetry breaking

     Solution[1] #< Solution[K2],

  then just the second solution is shown.

  Note: There are only solutions when K mod 4 == 0 or K mod 4 == 3.


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

*/

% Licenced under CC-BY-4.0 : http://creativecommons.org/licenses/by/4.0/

import util.
import cp.

main => go.

go =>
    K = 12,
    writef("K: %d\n", K),
    time2($langford(K, Solution, Position)),
    writef("Solution: %w\n", Solution),
    writef("Position: %w\n", Position).

%
% Get the first (if any) solutions for K in 2..40
%
go2 => 
    foreach(K in 2..40) 
       writef("K: %d\n", K),
       if time2($langford(K, Solution, Position)) then
           if nonvar(Solution) then 
             writef("Solution: %w\n", Solution),
             writef("Position: %w\n", Position)
           else 
             writef("No solution\n")
           end,
          writef("\n")
       end
    end.

 
langford(K, Solution, Position) =>

        if not ((K mod 4 == 0; K mod 4 == 3)) then
           println("There is no solution for K unless K mod 4 == 0 or K mod 4 == 3"),
           fail
        end,

        K2 = 2*K,
        Position = new_list(K2),
        Position :: 1..K2,

        Solution = new_list(K2),
        Solution :: 1..K,

        all_different(Position),

        % symmetry breaking:
        Solution[1] #< Solution[K2],
        % Solution[1] #> Solution[K2],

        foreach(I in 1..K)
            % This "advanced" usage of element don't work
            /*
              Position[K+I] #= Position[I] + I+1,
              Solution[Position[I]] = I,
              Solution[Position[K+I]] = I
            */

            IK = I+K,
            element(IK, Position, PositionIK),
            element(I, Position, PositionI),
            PositionIK #= PositionI + I+1,
            element(PositionI,Solution,SolutionPositionI),
            SolutionPositionI #= I,
            element(PositionIK,Solution,SolutionPositionIK),
            SolutionPositionIK #= I
        end,

        Vars = Solution ++ Position,
        solve([ff,updown], Vars).