4
votes

I am new to Prolog programming and I wrote a code for solving 4 by 4 magic squares but when I run the program, the program does not give any output; it just keeps on executing(busy) and eventually I have to quit SWI-Prolog. Kindly guide me regarding this issue.

Code:

:- use_module(library(clpfd)).

magic_square(Puzzle,Sum) :-
    Puzzle = [S11,S12,S13,S14,
              S21,S22,S23,S24,
              S31,S32,S33,S34,
              S41,S42,S43,S44],

    Puzzle ins 1..16,
    all_different(Puzzle),
    labeling([],Puzzle),

    R1 = [S11,S12,S13,S14],           % rows
    R2 = [S21,S22,S23,S24],
    R3 = [S31,S32,S33,S34],
    R4 = [S41,S42,S43,S44],

    C1 = [S11,S21,S31,S41],           % columns
    C2 = [S12,S22,S32,S42],
    C3 = [S13,S23,S33,S43],
    C4 = [S14,S24,S34,S44],

    Diag1 = [S11,S22,S33,S44],        % diagonals
    Diag2 = [S14,S23,S32,S41],

    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.

I tried changing sum variable to 64, but still nothing happens. I am using SWI-Prolog for this.

1
for basic usage, label(Vars) (or labeling(Vars,[])) should be put after the constraintsCapelliC
but the same code works when edited for 3 by 3 magic square!Pooja Gupta
i got it to work but it is not printing full puzzle variable values,it is printing P = [1, 2, 15, 16, 12, 14, 3, 5, 13|...], as solution instead of 4 by 4 valuesPooja Gupta
As long as there are choice points left, you can also simply type "w" to print the solution.mat

1 Answers

4
votes

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:

  1. Declare the initial domain of all variables used
  2. State all relations that should hold via constraints
  3. 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!