Solution: Magic Squares with initial constraints

The code...

declare
fun {Nuss N Data}
   NN=N*N
   L1N={List.number 1 N 1} 
in
   proc { $ Square}
      fun {Field I J}
         Square.((I-1)*N + J)
      end
      proc {Assert F}
         {FD.sum {Map L1N F} '=:' Sum}
      end
      Sum = {FD.decl}
in
      {FD.tuple square NN 1#NN Square}
      {FD.distinct Square}
      % Diagonals
      {Assert fun {$ I} {Field I I} end}
      {Assert fun {$ I} {Field I N+1-I} end}
      % Columns
      {For 1 N 1
       proc {$ I} {Assert fun {$ J} {Field I J} end} end}
      % Rows
      {For 1 N 1
       proc {$ J} {Assert fun {$ I} {Field I J} end} end}
      % Redundant Constraint
      NN*(NN+1) div 2 =: N*Sum
      %
      {ForAll Data
       proc {$ P}
          case P
          of p(X Y Val) then Square.((X-1)*N + Y) =: Val
          [] l(X Y) then Square.((X-1)*N + Y) <: 9
          end
       end
      }
      {FD.distribute split Square}
   end
end
   

Data1=[p(1 4 5) p(2 1 13) p(2 4 8) p(3 1 3) p(4 3 7) l(1 2) l(4 1)]
Data2=[p(2 2 6) p(2 4 4) p(3 1 10) p(4 1 2) p(4 3 3) l(1 2) l(3 4)]
Data3=[p(1 1 7) p(1 3 2) p(2 2 4) p(4 1 1) p(4 2 11)]
Data4=[p(1 2 3) p(2 1 15) p(2 2 6) p(3 3 5) p(4 1 4) l(3 1)]
Data5=[p(2 1 15) p(2 2 1) p(3 4 3) p(4 1 6) p(4 3 7) l(1 4) l(2 3)]
Data6=[p(1 1 7) p(2 4 8) p(3 3 3) p(4 2 5) l(1 2) l(4 1)]


{ExploreAll {Nuss 4 Data1}}
          

square(16  4  9  5
       13  1 12  8
        3 15  6 10
        2 14  7 11)