/* 

  Integer partition in Picat.

  From SICStus Prolog partit.pl
  """
   Benchmark (Finite Domain)            INRIA Rocquencourt - ChLoE Project 
                                                                           
   Name           : partit.pl                                              
   Title          : integer partitionning                                  
   Original Source: Daniel Diaz - INRIA France                             
   Adapted by     :                                                        
   Date           : September 1993                                         
                                                                           
   Partition numbers 1,2,...,N into two groups A and B such that:          
     a) A and B have the same length,                                      
     b) sum of numbers in A = sum of numbers in B,                         
     c) sum of squares of numbers in A = sum of squares of numbers in B.   
                                                                           
   This problem admits a solution if N is a multiple of 8.                 
                                                                           
   Note: finding a partition of 1,2...,N into 2 groups A and B such that:  
                                                                           
       Sum (k^p) = Sum l^p                                                 
     k in A      l in B                                                    
                                                                           
   admits a solution if N mod 2^(p+1) = 0 (N is a multilple of 2^(p+1)).   
   Condition a) is a special case where p=0, b) where p=1 and c) where p=2.
                                                                           
   Two redundant constraints are used:                                     
                                                                           
     - in order to avoid duplicate solutions (permutations) we impose      
       A1<A2<....<AN/2, B1<B2<...<BN/2 and A1<B1. This achieves much more  
       pruning than only alldifferents(A) and alldifferents(B).            
                                                                           
     - the half sums are known                                             
                                N                                          
          Sum k^1 = Sum l^1 = (Sum i) / 2 = N*(N+1) / 4                    
         k in A    l in B      i=1                                         
                                N                                          
          Sum k^2 = Sum l^2 = (Sum i^2)/2 = N*(N+1)*(2*N+1) / 12           
         k in A    l in B      i=1                                         
                                                                           
   Solution:                                                               
   N=8  A=[1,4,6,7]                                                        
        B=[2,3,5,8]                                                        
                                                                           
   N=16 A=[1,2,7,8,11,12,13,14]                                            
        B=[3,4,5,6, 9,10,15,16]                                            
                                                                           
   N=24 A=[1,2,3,8,13,14,15,16,18,19,20,21]                                
        B=[4,5,6,7, 9,10,11,12,17,22,23,24]                                
                                                                           
   N=32 A=[1,3,4, 7, 8, 9,16,20,21,22,23,24,25,26,27,28]                   
        B=[2,5,6,10,11,12,13,14,15,17,18,19,29,30,31,32]                   

  """

  It seems that there are solutions IFF
  * (N mod 2**2) == 0 AND
  * (N mod 2**3) == 0 OR (N mod 2**2) == 4  (not just == 0)
  Except for N = 4, which has no solution.

  This works for N <= 80.

  Also, there is an experimental symmetry breaking:
   - we can assign numbers 1..N div 8 to partition 1,
   - and numberts 1+N div 8.1 .. 2*(N div 8) to partition 2

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

*/


import cp.

main => go.

%
% Testing partit/2.
%
go ?=>
  foreach(N in 2..60)
    println(n=N),
    NMod4 = (N mod 2**2),
    NMod8 = (N mod 2**3),
    println([nmod4=NMod4,nmod8=NMod8]),

    if NMod4 == 0, (NMod8 == 0 ; NMod8 == 4) then
      garbage_collect(300_000_000),
      if partit(N, X) then 
        foreach(I in 1..2)
          T = [J : J in 1..N, X[J] == I],
          println(I=T=len(T)=sum(T)=sum([J**2 : J in T]))
        end,
        nl
      else
        if N != 4 then
          print("THIS SHOULD BE A SOLUTION!"),
          halt
        end
      end
    end,
    nl,
    flush(stdout)
  end,
  % fail,
  
  nl.
go => true.


%
% Testing partit2/2
%
go2 ?=>
  foreach(N in 2..60)
     println(n=N),
     NMod4 = (N mod 2**2),
     NMod8 = (N mod 2**3),
     if NMod4 == 0, (NMod8 == 0 ; NMod8 == 4) then
        if partit2(N,X) then
         foreach(I in 1..2)
            T = [X[I,J] : J in 1..N div 2],
            println(I=T=len(T)=sum(T)=sum([J**2 : J in T]))
         end,
         nl
       else
         if N != 4 then
          print("THIS SHOULD BE A SOLUTION!"),
          halt
         end
       end
    end
  end,
  nl.
go2 => true.

%
% Using one list and assign the partition number.
% Slower: N=2..50: 2.032s
% 
partit(N, X) =>

   X = new_list(N),
   X :: 1..2,
   
   Sum = N*(N+1) div 4,
   SquaredSum = N*(N+1)*(2*N+1) div 12,
   PartitionLen = N div 2,
   println([len=PartitionLen,sum=Sum,squaredSum=SquaredSum, nDiv8=(N div 8)]),

   foreach(I in 1..2)
     PartitionLen #= sum([X[J] #= I : J in 1..N]),
     Sum #= sum([J*(X[J] #= I) : J in 1..N]),
     SquaredSum #= sum([(J**2)*(X[J] #= I) : J in 1..N])     
   end,

   % X[1] #= 1, % symmetry breaking
   % Experimental symmetry breaking:
   % - we can assign numbers 1..N div 8 to partition 1,
   % - and numbers 1+N div 8.1 .. 2*(N div 8) to partition 2
   Ndiv8 = N div 8,
   foreach(I in 1..Ndiv8)
     X[I] #= 1,
     X[I+Ndiv8] #= 2
   end,
   println(solve),
   solve($[degree,updown],X).

%
% Another approach: Use 2 separate lists .
% Faster: N=2..50: 0.359s
% 
partit2(N, X) =>
   N2 = N div 2,

   X = new_array(2,N2),
   X :: 1..N,

   all_different(X.vars),

   Sum = N*(N+1) div 4,
   SquaredSum = N*(N+1)*(2*N+1) div 12,
   println([len=N2,sum=Sum,squaredSum=SquaredSum, nDiv8=(N div 8)]),
   
   foreach(I in 1..2)
     increasing(X[I,1..N2]),
     Sum #= sum(X[I,1..N2]),
     SquaredSum #= sum([X[I,J]**2 : J in 1..N2]),
   end,

   % Same symmetry breaking as in partit/2.
   Ndiv8 = N div 8,
   foreach(J in 1..Ndiv8)
     X[1,J] #= J,
     X[2,J] #= J+Ndiv8,
   end,

   solve($[degree,updown],X).