/*

  Sicherman Dice in Picat.

  From http://en.wikipedia.org/wiki/Sicherman_dice
  """ 
  Sicherman dice are the only pair of 6-sided dice which are not normal dice, 
  bear only positive integers, and have the same probability distribution for 
  the sum as normal dice.
  
  The faces on the dice are numbered 1, 2, 2, 3, 3, 4 and 1, 3, 4, 5, 6, 8.
  """

  I read about this problem in a book/column by Martin Gardner long
  time ago, and got inspired to model it now by the WolframBlog post
  "Sicherman Dice": http://blog.wolfram.com/2010/07/13/sicherman-dice/

  This model gets the two different ways, first the standard way and
  then the Sicherman dice:
  
  x1 = array1d(1..6, [1, 2, 3, 4, 5, 6]);
  x2 = array1d(1..6, [1, 2, 3, 4, 5, 6]);
  ----------
  x1 = array1d(1..6, [1, 2, 2, 3, 3, 4]);
  x2 = array1d(1..6, [1, 3, 4, 5, 6, 8]);


  Extra: If we also allow 0 (zero) as a valid value then the 
  following two solutions are also valid:
  
  x1 = array1d(1..6, [0, 1, 1, 2, 2, 3]);
  x2 = array1d(1..6, [2, 4, 5, 6, 7, 9]);
  ----------
  x1 = array1d(1..6, [0, 1, 2, 3, 4, 5]);
  x2 = array1d(1..6, [2, 3, 4, 5, 6, 7]);
  
  These two extra cases are mentioned here:
  http://mathworld.wolfram.com/SichermanDice.html

  Also see Sicherman's page which includes some
  more history on this:
  https://userpages.monmouth.com/~colonel/sdice.html
  It includes a Perl program that is quite fast but is - IMHO -
  more complex using cyclotomic polynomials etc.

  
  go2/0 shows the solutions for dice with sides 2..10 with min values 1 and 0.

  Here are the solutions for N=2..10 (only for min value = 1):

  n = 2
  [1,2]
  [1,2]

  n = 3
  [1,2,3]
  [1,2,3]

  n = 4
  [1,2,3,4]
  [1,2,3,4]

  [1,2,2,3]
  [1,3,3,5]

  n = 5
  [1,2,3,4,5]
  [1,2,3,4,5]

  n = 6
  [1,2,3,4,5,6]
  [1,2,3,4,5,6]

  [1,2,2,3,3,4]
  [1,3,4,5,6,8]

  n = 7
  [1,2,3,4,5,6,7]
  [1,2,3,4,5,6,7]

  n = 8
  [1,2,2,3,3,4,4,5]
  [1,3,5,5,7,7,9,11]
        
  [1,2,2,3,5,6,6,7]
  [1,3,3,5,5,7,7,9]

  [1,2,3,4,5,6,7,8]
  [1,2,3,4,5,6,7,8]

  [1,2,3,3,4,4,5,6]
  [1,2,5,5,6,6,9,10]

  n = 9
  [1,2,2,3,3,3,4,4,5]
  [1,4,4,7,7,7,10,10,13]

  [1,2,3,4,5,6,7,8,9]
  [1,2,3,4,5,6,7,8,9]

  n = 10
  [1,2,2,3,3,4,4,5,5,6]
  [1,3,5,6,7,8,9,10,12,14]

  [1,2,3,4,5,6,7,8,9,10]
  [1,2,3,4,5,6,7,8,9,10]



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

*/

import sat.

main => go.

% The standard Sicherman's dice
go ?=>
   nolog,
   N = 6,
   show_all_solutions(N),
   nl.
go => true.

% For N=2..10 (number of sides per die)
go2 ?=>
   nolog,
   foreach(N in 2..10)
     show_all_solutions(N)
   end,
   nl.
go2 => true.


% N = 12
go3 ?=>
   nolog,
   N = 12,
   show_all_solutions(N),
   nl.
go3 => true.


%
% Show all solutions for N (number of sides per die).
%
show_all_solutions(N) =>
   println(n=N),
   StandardDist = distribution(N),   
   _ = findall(_, sicherman_dice(N,StandardDist,_X1a,_X2a,1)),
   _ = findall(_, sicherman_dice(N,StandardDist,_X1b,_X2b,0)),
   flush(stdout),
   nl.


%
% N: Number of side per dice
% M: Max values of each side
%
sicherman_dice(N,StandardDist, X1,X2, Min) =>
   printf("\nMin: %d\n", Min),
   DistLen = StandardDist.len,   
   println([standardDistribution=StandardDist, distLen=DistLen]),

   M = DistLen-1, % max value of each side

   % The two dice
   X1 = new_list(N),
   X1 :: Min..M,   
   X2 = new_list(N),
   X2 :: Min..M,

   % ensure standard distributions of the sums
   foreach(K in 1..DistLen)
      StandardDist[K] #= sum([(X1[I]+X2[J] #= K+1) : I in 1..N, J in 1..N]) 
   end,

   % Symmetry breaking
   increasing(X1),
   increasing(X2),

   % Symmetry breaking: x1 is less <= to x2
   foreach(I in 1..N) X1[I] #=< X2[I] end,

   Vars = X1 ++ X2,
   solve($[ffd,updown],Vars),

   println(X1),
   println(X2),
   % check_distribution(X1,X2),
   flush(stdout),
   nl.

%
% What is the distribution of 2 normal dice with
% the NumSides sides (with values 1..NumSides)?
%
distribution(N) = Dist =>
  Dist = 1..N ++ (N-1)..-1..1.

%
% Check the distribution of two dice and compare with the
% theoretical distribution.
%
check_distribution(D1,D2) =>
   Dist = distribution(D1.len),
   println('dist    '=Dist),
   DistMap = new_map(),
   foreach(I in D1, J in D2)
      S = I+J,
      DistMap.put(S,DistMap.get(S,0)+1)
   end,
   ThisDist = [V : _=V in DistMap.to_list.sort],
   println(thisDist=ThisDist),
   if Dist == ThisDist then
     println(ok)
   else
     println(not_ok)
   end.