Restrict search in Prolog - Magic Sqare

I want to solve the Most Perfect Magic Square with a Prolog program.

Wiki page: https://en.wikipedia.org/wiki/Most-perfect_magic_square

When I input the query "magic_square(4, [[7, 12, 1, 14], [2, 13, 8, 11], [16, 3, 10, 5], [9, 6, 15, 4]])." (which is a valid Magic Square) my program returns true. So I assume my rule basis is correct.

Unfortunately, if more than 9 Values are unknown, it takes forever to find a solution.

I need help to restrict my search, so the Program finds a solution in a reasonable time. Ideally, it should also work with a 12 x 12 grid(and also every other multiple of 4) and no Values given: magic_square(12, Matrix).

Thanks a lot!

Here my Code:

:- use_module(library(clpfd)).

diag2_sum(0, _, _, _).
diag2_sum(I0, C1, Row1, Row3) :-
    I0 > 0,
    nth1(I0,Row1,A),
    (I0 =:= 2 -> I2 = 4 ; I2 is mod(I0 + 2,4)),
    nth1(I2,Row3,B),
    C1 =:= A + B,
    I1 is I0 - 1,
    diag2_sum(I1, C1, Row1, Row3).

diag_sum([_,_], _, _).
diag_sum([Row1|Tail], C1, N) :-
    nth1(2,Tail,Row3),
    diag2_sum(N, C1, Row1,Row3),
    diag_sum(Tail, C1, N).

square_sum_x(_, _, _, 0).
square_sum_x(Row1, Row2, C2, I0) :-
    (I0 =:= 3 -> I2 = 4 ; I2 is mod(I0 + 1,4)),
    nth1(I0,Row1,Elem1),
    nth1(I2,Row1,Elem2),
    nth1(I0,Row2,Elem3),
    nth1(I2,Row2,Elem4),
    C2 =:= Elem1 + Elem2 + Elem3 + Elem4,
    I1 is I0 - 1,
    square_sum_x(Row1, Row2, C2, I1).


square_sum_y(_, _, 0, _).
square_sum_y(Matrix, C2, I0, N) :-
    (I0 =:= 3 -> I2 = 4 ; I2 is mod(I0 + 1,4)),
    nth1(I0,Matrix,Row1),
    nth1(I2,Matrix,Row2),
    
    square_sum_x(Row1,Row2, C2, N),
    I1 is I0 - 1,
    square_sum_y(Matrix, C2, I1, N).

magic_square(N, Matrix) :-
    Nmax is N * N,
    C1 is Nmax + 1,
    C2 is C1 * 2,
    write(C1),nl,write(C2),nl,
    length(Matrix, N),
    maplist(same_length(Matrix), Matrix),
    append(Matrix, Vs),
    Vs ins 1..Nmax, all_different(Vs),
    maplist(label, Matrix),
    %write(Matrix),nl,
    diag_sum(Matrix, C1, N),
    square_sum_y(Matrix, C2, N, N).

% some queries i tryed out:
% magic_square(4, [[7, 12, 1, 14], [2, 13, 8, 11], [16, 3, 10, 5], [9, 6, 15, 4]]).    
% magic_square(4, [[7, 12, 1, 14], [2, 13, 8, B4], [C1, C2, C3, C4], [D1, D2, D3, D4]]).
% magic_square(4, [[7, A2, A3, 14], [2, B2, 8, B4], [C1, C2, 10, C4], [D1, D2, D3, 4]]).

Solution 1:

I confirm that a query with 10 holes seems to take very long, and even with 9 holes it's not very fast in SWI-Prolog:

?- time(magic_square(4, [[_, _, _, _], [_, _, 8, 11], [_, _, 10, 5], [_, 6, 15, 4]])).
17
34
% 88,086,790 inferences, 3.449 CPU in 3.449 seconds (100% CPU, 25543070 Lips)
true .

So where does that time go? Let's use SWI-Prolog's profiler:

?- profile(magic_square(4, [[_, _, _, _], [_, _, 8, 11], [_, _, 10, 5], [_, 6, 15, 4]])).
17
34
true.

This opens a GUI window. If we click on magic_square/2 in the list of predicates, we see:

Note the line for __aux_maplist/2_label+0/1: 94.7% of your execution time is spent there! This is your maplist(label, Matrix) line, where you get stuck before ever getting to the interesting part with the diag_sum and square_sum_y predicates that actually define what a magic square is.

Basically, by posting all_different(Vs) followed immediately by labeling it, you are asking Prolog to enumerate all permutations of different numbers that could fill the holes in the input term. The number of such permutations, and therefore your execution time, grows exponentially.

You should not use labeling in this way, before all constraints are known. Labeling should be the last thing you do in your predicate, or even better, it should be done by a different predicate from the one where you define your constaints.

So you could set up things like this:

% "Core relation" defining the shape of the solution and constraints.
magic_square_(N, Matrix) :-
    Nmax is N * N,
    C1 is Nmax + 1,
    C2 is C1 * 2,
    write(C1),nl,write(C2),nl,
    length(Matrix, N),
    maplist(same_length(Matrix), Matrix),
    append(Matrix, Vs),
    Vs ins 1..Nmax, all_different(Vs),
    diag_sum(Matrix, C1, N),
    square_sum_y(Matrix, C2, N, N).

magic_square(N, Matrix) :-
    magic_square_(N, Matrix),
    maplist(label, Matrix).

Now diag_sum and square_sum_x will be called with unbound (but constrained) variables C1 and C2, respectively. This is fine! This is what you want when you use CLP(FD)! But you need to change the numerical equality goals for these variables from =:= (which expects ground values) to #= (which actually works with constraint variables).

And now you're fast with 9 variables:

?- Matrix = [[_, _, _, _], [_, _, 8, 11], [_, _, 10, 5], [_, 6, 15, 4]], time(magic_square(4, Matrix)).
17
34
% 8,473 inferences, 0.006 CPU in 0.006 seconds (100% CPU, 1506518 Lips)
Matrix = [[7, 12, 1, 14], [2, 13, 8, 11], [16, 3, 10, 5], [9, 6, 15, 4]] ;
% 33 inferences, 0.000 CPU in 0.000 seconds (91% CPU, 749778 Lips)
false.

And with 10:

?- Matrix = [[_, _, _, _], [_, _, 8, 11], [_, _, 10, 5], [_, _, 15, 4]], time(magic_square(4, Matrix)).
17
34
% 8,933 inferences, 0.002 CPU in 0.002 seconds (100% CPU, 4338676 Lips)
Matrix = [[7, 12, 1, 14], [2, 13, 8, 11], [16, 3, 10, 5], [9, 6, 15, 4]] .

And with even more:

?- Matrix = [[_, _, _, _], [_, _, 8, _], [_, _, 10, _], [_, _, _, 4]], time(magic_square(4, Matrix)).
17
34
% 24,710 inferences, 0.009 CPU in 0.009 seconds (99% CPU, 2658996 Lips)
Matrix = [[7, 2, 11, 14], [12, 13, 8, 1], [6, 3, 10, 15], [9, 16, 5, 4]] .