Generate, Test and Intertwine

%% Generate, Test and Intertwine
%% -----------------------------

%% The aim here is to generate 4x4 magic squares. First of all, We start
%% with the 3x3 case, since that's easier.
%% This way any improvements can be reused, assuming they are general enough

%% I'll define a magic square as, an nxn table of numbers such that all the
%% rows, columns and diagonals sum to the same value, each number in the range
%% 1 .. n^2 appears once in the grid.

%% In Prolog, We'll represent squares as lists:

%: [4,8,2,
%:  3,5,7,
%:  8,1,6]

%% So to generate this, We'll select numbers from a list of digits and insert
%% them into the table and check that everything adds up.

saturn_digits([1,2,3,4,5,6,7,8,9]).

saturn([A,B,C,
        D,E,F,
        G,H,I]) :-
        saturn_digits(D1),
        select(A, D1, D2),
        select(B, D2, D3),
        select(C, D3, D4),
        select(D, D4, D5),
        select(E, D5, D6),
        select(F, D6, D7),
        select(G, D7, D8),
        select(H, D8, D9),
        select(I, D9, []),
        S is A + B + C,
        S is D + E + F,
        S is G + H + I,
        S is A + D + G,
        S is B + E + H,
        S is C + F + I,
        S is A + E + I,
        S is C + E + G.

%% That first attempt works fine:

%? ?- saturn(Square).
%? Square = [2, 7, 6, 9, 5, 1, 4, 3, 8] ;
%? Square = [2, 9, 4, 7, 5, 3, 6, 1, 8] ;
%? Square = [4, 3, 8, 9, 5, 1, 2, 7, 6] ;
%? Square = [4, 9, 2, 3, 5, 7, 8, 1, 6] ;
%? Square = [6, 1, 8, 7, 5, 3, 2, 9, 4] ;
%? Square = [6, 7, 2, 1, 5, 9, 8, 3, 4] ;
%? Square = [8, 1, 6, 3, 5, 7, 4, 9, 2] ;
%? Square = [8, 3, 4, 1, 5, 9, 6, 7, 2] ;

%% It's quite clear the program is correct,
%% The two important properties are
%%  * Soundness -- Generates only magic squares
%%  * Completeness -- Generates every magic square
%%
%% Since digits are selected from the list, It's possible for the square to
%% be filled in every possible way, and due to:
%:   select(I, D9, []),
%% It's ensured that every digit is used. Every row, column and diagonal is
%% summed, since S is used for all summations, this ensures that all the sums
%% are equal.


%% There is this pattern in the code though, which I don't very much like:

%: select(A, D1, D2),
%: select(B, D2, D3),
%: select(C, D3, D4),
%: select(D, D4, D5),
%: select(E, D5, D6),
%: select(F, D6, D7),
%: select(G, D7, D8),
%: select(H, D8, D9),
%: select(I, D9, []),

%% Prolog has a macro, DCG style, that lets us erase this pattern.
%% If we use --> instead of :- , then this is implicit, everything is sewn
%% together like we did manually, except automatically.
%% {} is used to escape it.
%% So now, It can be rewritten to the equivalent program:

saturn_2([A,B,C,
          D,E,F,
          G,H,I]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        select(E),
        select(F),
        select(G),
        select(H),
        select(I),
        { S is A + B + C },
        { S is D + E + F },
        { S is G + H + I },
        { S is A + D + G },
        { S is B + E + H },
        { S is C + F + I },
        { S is A + E + I },
        { S is C + E + G }.

%? ?- saturn_digits(D), saturn_2(Square, D, []).
%? D = [1, 2, 3, 4, 5, 6, 7, 8, 9],
%? Square = [2, 7, 6, 9, 5, 1, 4, 3, 8] ;
%? ..

%% It produces all the same answers in the same order. One of the fantastic
%% thing about Prolog is you could actually implement DCG in it, if it wasn't
%% already there, Or even your own macros.

%% Anyway, there is an optmisation which is quite blatant: select always
%% succeeds, It's the tests in {} which cause backtracking, So if we were to
%% intertwine the generating (done by select) with the testing, (as done by
%% summation) then backtracking which does occur, will happen sooner.

saturn_3([A,B,C,
          D,E,F,
          G,H,I]) -->
        select(A),
        select(B),
        select(C),
        { S is A + B + C },
        select(D),
        select(G),
        { S is A + D + G },
        select(E),
        { S is C + E + G },
        select(F),
        { S is D + E + F },
        select(H),
        { S is B + E + H },
        select(I),
        { S is G + H + I },
        { S is C + F + I },
        { S is A + E + I }.


%% So now, It's straightforward to use the same process for the 4x4 case.
%% We use DCG first in this case though.

jupiter_digits([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]).

jupiter([A,B,C,D,
         E,F,G,H,
         I,J,K,L,
         M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        select(E),
        select(F),
        select(G),
        select(H),
        select(I),
        select(J),
        select(K),
        select(L),
        select(M),
        select(N),
        select(O),
        select(P),
        
        { S is A + B + C + D },
        { S is E + F + G + H },
        { S is I + J + K + L },
        { S is M + N + O + P },
        
        { S is A + E + I + M },
        { S is B + F + J + N },
        { S is C + G + K + O },
        { S is D + H + L + P },
        
        { S is A + F + K + P },
        { S is D + G + J + M }.

%% And then rearrange

jupiter_2([A,B,C,D,
           E,F,G,H,
           I,J,K,L,
           M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        %% * * * *
        %% o o o o
        %% o o o o
        %% o o o o
        { S is A + B + C + D },
        
        select(E),
        select(I),
        select(M),
        %% * * * *
        %% * o o o
        %% * o o o
        %% * o o o
        { S is A + E + I + M },
        
        select(F),
        select(K),
        select(P),
        %% * * * *
        %% * * o o
        %% * o * o
        %% * o o *
        { S is A + F + K + P },
        
        select(G),
        select(H),
        %% * * * *
        %% * * * *
        %% * o * o
        %% * o o *
        { S is E + F + G + H },
        
        select(L),
        %% * * * *
        %% * * * *
        %% * o * *
        %% * o o *
        { S is D + H + L + P },
        
        select(J),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * o o *
        { S is I + J + K + L },
        { S is D + G + J + M },
        
        select(N),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * * o *
        { S is B + F + J + N },
        
        select(O),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * * * *
        { S is M + N + O + P },
        { S is C + G + K + O }.

%? ?- jupiter_digits(D), jupiter_2(Square, D, []).
%? 
%? D = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
%? Square = [1, 2, 15, 16, 12, 14, 3, 5, 13, 7, 10, 4, 8, 11, 6, 9] ;
%? ...

%% There is a huge speed difference in the two versions.
%% Here's the final program without comments:

jupiter_3([A,B,C,D,
           E,F,G,H,
           I,J,K,L,
           M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        { S is A + B + C + D },
        select(E),
        select(I),
        select(M),
        { S is A + E + I + M },
        select(F),
        select(K),
        select(P),
        { S is A + F + K + P },
        select(G),
        select(H),
        { S is E + F + G + H },
        select(L),
        { S is D + H + L + P },
        select(J),
        { S is I + J + K + L },
        { S is D + G + J + M },
        select(N),
        { S is B + F + J + N },
        select(O),
        { S is M + N + O + P },
        { S is C + G + K + O }.