/*

Langford's number problem in Picat.

Langford's number problem (CSP lib problem 24)
http://www.csplib.org/Problems/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

*/

% 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).