No exit in prologue Program: clp (fd) too slow

First of all, the first!

With clp (fd), always defer enumeration targets such as labeling/2

until all constraints are specified. In general, good practice works like this:

• Declare the source domain of all used variables
• List all relationships to be done with constraints
• Launch a search for a solution using labeling/2
  
  or other target enumeration target domains.

In the comment, you noted the same as with the 3x3 magic squares; Sure, but it didn’t work because of the reverse ordering of labels and message-constraints — it worked despite this .

The reason you got the solutions in a reasonable amount of time is because finding 3x3 magic squares is an order of magnitude easier than finding 4x4.

So, change the order of the targets in your predicate magicSquare4x4_withSum/2

like this:

:- use_module(library(clpfd)).

magicSquare4x4_withSum(Zs,Sum) :-
    Zs = [S11,S12,S13,S14,
          S21,S22,S23,S24,
          S31,S32,S33,S34,
          S41,S42,S43,S44], 
    Zs ins 1..16,                   % state the initial finite domain
    S11 + S12 + S13 + S14 #= Sum,   % rows
    S21 + S22 + S23 + S24 #= Sum,
    S31 + S32 + S33 + S34 #= Sum,
    S41 + S42 + S43 + S44 #= Sum,    
    S11 + S21 + S31 + S41 #= Sum,   % columns
    S12 + S22 + S32 + S42 #= Sum,
    S13 + S23 + S33 + S43 #= Sum,
    S14 + S24 + S34 + S44 #= Sum,
    S11 + S22 + S33 + S44 #= Sum,   % diagonals
    S14 + S23 + S32 + S41 #= Sum,   
    all_different(Zs).              % no two variables shall have the same value

      

        
        
        
      

    

The most common query magicSquare4x4_withSum/2

deforms deterministically in the blink of an eye; however, the answer provided does not provide specific meanings for any 4x4 magic square; instead, it presents pending constraints that a solution must adhere to in order to exist.

?- time(magicSquare4x4_withSum(Zs,Sum)).
% 7,780 inferences, 0.001 CPU in 0.001 seconds (99% CPU, 8396967 Lips)
Zs = [_G9481, _G9484, _G9487, _G9490, _G9493, _G9496, _G9499, _G9502|...],
_G9481 in 1..16,
all_different([_G9481, _G9484, _G9487, _G9490, _G9493, _G9496, _G9499|...]),
_G9481+_G9496+_G9511+_G9526#=Sum,
_G9481+_G9493+_G9505+_G9517#=Sum,
_G9481+_G9484+_G9487+_G9490#=Sum,
_G9484 in 1..16,
_G9484+_G9496+_G9508+_G9520#=Sum,
_G9487 in 1..16,
_G9487+_G9499+_G9511+_G9523#=Sum,
_G9490 in 1..16,
_G9490+_G9499+_G9508+_G9517#=Sum,
_G9490+_G9502+_G9514+_G9526#=Sum,
_G9493 in 1..16,
_G9493+_G9496+_G9499+_G9502#=Sum,
_G9496 in 1..16,
_G9499 in 1..16,
_G9502 in 1..16,
_G9505 in 1..16,
_G9505+_G9508+_G9511+_G9514#=Sum,
_G9508 in 1..16,
_G9511 in 1..16,
_G9514 in 1..16,
_G9517 in 1..16,
_G9517+_G9520+_G9523+_G9526#=Sum,
_G9520 in 1..16,
_G9523 in 1..16,
_G9526 in 1..16,
Sum in 4..64.

      

        
        
        
      

    

Now let's get started!

?- time((magicSquare4x4_withSum(Zs,Sum), labeling([],Zs))).
% 8,936,459 inferences, 1.025 CPU in 1.025 seconds (100% CPU, 8720473 Lips)
Zs = [1, 2, 15, 16, 12, 14, 3, 5, 13|...],
Sum = 34 ;
% 37,098 inferences, 0.006 CPU in 0.006 seconds (100% CPU, 6073990 Lips)
Zs = [1, 2, 15, 16, 13, 14, 3, 4, 12|...],
Sum = 34                                    % and the search goes on...

      

        
        
        
      

    

What can we do to speed up the search? First, we are trying to use different heuristics / emulation options.

?- time((magicSquare4x4_withSum(Zs,Sum), labeling([ff],Zs))).
% 5,056,298 inferences, 0.578 CPU in 0.578 seconds (100% CPU, 8749040 Lips)
Zs = [1, 2, 15, 16, 12, 14, 3, 5, 13|...],
Sum = 34 ;
% 36,783 inferences, 0.006 CPU in 0.006 seconds (100% CPU, 5733914 Lips)
Zs = [1, 2, 15, 16, 13, 14, 3, 4, 12|...],
Sum = 34                                    % and the search goes on...

      

        
        
        
      

    

What else can we do? Help me find by putting the amount as a constant.

?- time((magicSquare4x4_withSum(Zs,34), labeling([],Zs))).
% 106,296 inferences, 0.017 CPU in 0.017 seconds (100% CPU, 6242045 Lips)
Zs = [1, 2, 15, 16, 12, 14, 3, 5, 13|...] ;
% 36,858 inferences, 0.009 CPU in 0.009 seconds (100% CPU, 4076658 Lips)
Zs = [1, 2, 15, 16, 13, 14, 3, 4, 12|...] ;
% 209,206 inferences, 0.028 CPU in 0.028 seconds (100% CPU, 7430044 Lips)
Zs = [1, 2, 16, 15, 13, 14, 4, 3, 12|...]   % and the search goes on...

      

        
        
        
      

    

Performance, a priori limiting the sum to 34 has a huge effect!

Edit 2015-04-24

What else can we do?

We can add additional constraints that limit the search and solution space by excluding solutions that are simply "flips" or "rotations" of each other.

Let's look at 3x3 magic squares: There are 8 different ones, but in fact they all "flip" / "rotate" the same magic square. Of these 8, we can safely get rid of 7 by knowing how to build them if we leave the rest.

Symmetry breaking limits the size of both the search space and the solution space, we can express it using order constraints between the four corners of the magic square:

magicSquare4x4withRestrictedSymmetries(Zs) :-
    Zs = [S11,_,_,S14,
            _,_,_,_,
            _,_,_,_,
          S41,_,_,S44], 
    S11 #< S14,
    S11 #< S44,
    S11 #< S41,
    S14 #< S41.

      

        
        
        
      

    

These additional constraints may or may not help us find the first solution faster, but they definitely help us in finding all solutions.

First, let's do it without additional symmetry-breaking constraints:

?- time((findall(t,(magicSquare4x4_withSum(Zs,34),
                    labeling([ff],Zs)),Ts),
         length(Ts,N_Sols))).
% 1,526,766,108 inferences, 152.1 CPU in 152.1 seconds (100% CPU, 10033768 Lips)
Ts = [t, t, t, t, t, t, t, t, t|...],
N_Sols = 7040.

      

        
        
        
      

    

Now with additional restrictions:

?- time((findall(t,(magicSquare4x4_withSum(Zs,34),
                    magicSquare4x4withRestrictedSymmetries(Zs),
                    labeling([ff],Zs)),Ts),
         length(Ts,N_Sols))).
% 129,527,384 inferences, 12.580 CPU in 12.578 seconds (100% CPU, 10295893 Lips)
Ts = [t, t, t, t, t, t, t, t, t|...],
N_Sols = 880.

      

        
        
        
      

    

Great victory! Finding all solutions more than 10 times faster. As expected, the number of solutions is divisible by exactly 8 when symmetry is broken (880 * 8 = 7040). NTN!