First things first!
With clp(fd), always postpone enumeration goals like labeling/2
until all constraints have been stated. In general, it is good practise to proceed as follows:
- Declare the initial domain of all variables used
- State all relations that should hold via constraints
- Trigger search for a solution by
labeling/2
or other finite-domain enumeration goals
In a comment you noted the same worked with 3x3 magic squares; Granted, but it worked not because of the reversed ordering of labeling and constraint-posting---it worked in spite of it.
The reason you got solutions in reasonable time is that finding magic squares of size 3x3 is orders of magnitude easier than finding the 4x4 ones.
So reorder the goals 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 general query of magicSquare4x4_withSum/2
succeeds deterministically in no time; however, the answer given does not show concrete values of some 4x4 magic square; instead it presents the pending constraints that a solution must obey 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 it 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 search? First, we try different labeling heuristics / 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? Aid the search by supplying the sum 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-wise, a-priori constraining the sum to the value 34 has a massive effect!
Edit 2015-04-24
What else can we do?
We can add additional constraints that restrict the search and the solution space by excluding solutions that are mere "flips" or "rotations" of each other.
Let's consider magic squares of size 3x3: There are 8 different ones, but in fact they are all "flips"/"rotations" of the same magic square.
Of these 8, we can eliminate 7 safely in knowing how to construct them when given the remaining one.
Breaking symmetry restricts the size of
both the search space and the solution space, we can express it by using ordering 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 when searching for all solutions.
First, let's do it without the additional constraints breaking symmetries:
?- 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 the additional constraints:
?- 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.
Big win! The search for all solutions is more than 10 times faster. As expected, the number of solutions is divided exactly by 8 when breaking symmetries (880 * 8 = 7040). HTH!