# Latin squares

A question came up on StackExchange Proof Assistants and noted Prolog and Latin squares.

OEIS A274806 states that for a square of size 4 with a diagonal and anti-diagonal there should be 48 Latin squares.

Since I could not find a listing of all 48 Latin squares of size 4 with a diagonal and anti-diagonal, created some Prolog code to generate them.

This is not optimized code and could probably be done using constraints. This was done to quickly get a listing of the 48 squares.

Click triangle to expand section.

Prolog source code

Note: The directory can be changed but the user must have write access to the directory.

Name: `example.pl`
Directory: `C:/Users/Groot`

``````:- module(example,
[
latin_square/2,
is_normal/1
]).

% Sxy   x - column, y - row from top to bottom (think over and down)
%   11 21 31 41
%   12 22 32 42
%   13 23 33 43
%   14 24 34 44

latin_square(N,Latin_square) :-
N == 4,
!,
sequence(N,Values),
main_diagonal(Values,MD),

sequence(N,L) :-
numlist(1,N,L).

main_diagonal(L,P) :-
permutation(L,P).

% Given main diagonal, generate valid antidiagonal
anti_diagonal(L,P) :-
permutation(L,P),
anti_diagonal_row_diff(L,P),
anti_diagonal_col_diff(L,P).

anti_diagonal_row_diff(A,B) :-
seq_diff(A,B).
anti_diagonal_col_diff(A,B) :-
reverse(B,C),
seq_diff(A,C).

seq_diff([],[]) :- !.
seq_diff([H|_],[H|_]) :-
!,
fail.
seq_diff([_|T0],[_|T1]) :-
seq_diff(T0,T1).

% Given diagonal and anti-diagonal finish Latin square.
finish(S,[S11,S22,S33,S44],[S41,S32,S23,S14],[[S11,S21,S31,S41],[S12,S22,S32,S42],[S13,S23,S33,S43],[S14,S24,S34,S44]]) :-
member(S21,S),
S21 \= S11,       true,   var(S31), S21 \= S41,
S21 \= S22,
S21 \= S23,
var(S24),

member(S31,S),
S31 \= S11, S31 \= S21,       true, S31 \= S41,
S31 \= S32,
S31 \= S33,
var(S34),

member(S12,S),
S12 \= S11,
true, S12 \= S22, S12 \= S32,   var(S42),
var(S13),
S12 \= S14,

member(S42,S),
S42 \= S41,
S42 \= S12, S42 \= S22, S42 \= S32,       true,
var(S43),
S42 \= S44,

member(S13,S),
S13 \= S11,
S13 \= S12,
true, S13 \= S23, S13 \= S33,   var(S43),
S13 \= S14,

member(S43,S),
S43 \= S41,
S43 \= S42,
S43 \= S13, S43 \= S23, S43 \= S33,       true,
S43 \= S44,

member(S24,S),
S24 \= S21,
S24 \= S22,
S24 \= S23,
S24 \= S14,       true,   var(S34), S24 \= S44,

member(S34,S),
S34 \= S31,
S34 \= S32,
S34 \= S33,
S34 \= S14, S34 \= S24,       true, S34 \= S44,
true.

:- multifile
user:portray/1.

user:portray([[S11,S21,S31,S41],[S12,S22,S32,S42],[S13,S23,S33,S43],[S14,S24,S34,S44]]) :-
nl,
format('~d ~d ~d ~d~n',[S11,S21,S31,S41]),
format('~d ~d ~d ~d~n',[S12,S22,S32,S42]),
format('~d ~d ~d ~d~n',[S13,S23,S33,S43]),
format('~d ~d ~d ~d~n',[S14,S24,S34,S44]).

check_uniqueness :-
findall(LS,latin_square(4,LS),List),
length(List,LL),
list_to_set(List,Set),
length(Set,LL).

is_normal([H|_]) :-
length(H,N),
numlist(1,N,H).
``````
Example run
``````Welcome to SWI-Prolog (threaded, 64 bits, version 8.5.5)

?- working_directory(_,'C:/Users/Groot').
true.

?- [example].
true.

?- tell('Latin Squares.txt'),findall(_,(latin_square(4,LS),print(LS)),S),length(S,N),told.
S = [_, _, _, _, _, _, _, _, _|...],
N = 48.

?- example:check_uniqueness.
true.

?- latin_square(4,LS),is_normal(LS).
LS =
1 2 3 4
4 3 2 1
2 1 4 3
3 4 1 2
;
LS =
1 2 3 4
3 4 1 2
4 3 2 1
2 1 4 3
;
false.
``````
File: Latin Squares.txt
``````
1 3 4 2
4 2 1 3
2 4 3 1
3 1 2 4

1 4 2 3
3 2 4 1
4 1 3 2
2 3 1 4

1 4 3 2
3 2 1 4
2 3 4 1
4 1 2 3

1 3 2 4
4 2 3 1
3 1 4 2
2 4 1 3

1 2 4 3
4 3 1 2
3 4 2 1
2 1 3 4

1 4 3 2
2 3 4 1
4 1 2 3
3 2 1 4

1 4 2 3
2 3 1 4
3 2 4 1
4 1 3 2

1 2 3 4
4 3 2 1
2 1 4 3
3 4 1 2

1 2 3 4
3 4 1 2
4 3 2 1
2 1 4 3

1 3 4 2
2 4 3 1
3 1 2 4
4 2 1 3

1 3 2 4
2 4 1 3
4 2 3 1
3 1 4 2

1 2 4 3
3 4 2 1
2 1 3 4
4 3 1 2

2 3 4 1
4 1 2 3
1 4 3 2
3 2 1 4

2 4 1 3
3 1 4 2
4 2 3 1
1 3 2 4

2 4 3 1
3 1 2 4
1 3 4 2
4 2 1 3

2 3 1 4
4 1 3 2
3 2 4 1
1 4 2 3

2 1 4 3
4 3 2 1
3 4 1 2
1 2 3 4

2 4 3 1
1 3 4 2
4 2 1 3
3 1 2 4

2 4 1 3
1 3 2 4
3 1 4 2
4 2 3 1

2 1 3 4
4 3 1 2
1 2 4 3
3 4 2 1

2 1 3 4
3 4 2 1
4 3 1 2
1 2 4 3

2 3 4 1
1 4 3 2
3 2 1 4
4 1 2 3

2 3 1 4
1 4 2 3
4 1 3 2
3 2 4 1

2 1 4 3
3 4 1 2
1 2 3 4
4 3 2 1

3 2 4 1
4 1 3 2
1 4 2 3
2 3 1 4

3 4 1 2
2 1 4 3
4 3 2 1
1 2 3 4

3 4 2 1
2 1 3 4
1 2 4 3
4 3 1 2

3 2 1 4
4 1 2 3
2 3 4 1
1 4 3 2

3 1 4 2
4 2 3 1
2 4 1 3
1 3 2 4

3 4 2 1
1 2 4 3
4 3 1 2
2 1 3 4

3 4 1 2
1 2 3 4
2 1 4 3
4 3 2 1

3 1 2 4
4 2 1 3
1 3 4 2
2 4 3 1

3 1 2 4
2 4 3 1
4 2 1 3
1 3 4 2

3 2 4 1
1 4 2 3
2 3 1 4
4 1 3 2

3 2 1 4
1 4 3 2
4 1 2 3
2 3 4 1

3 1 4 2
2 4 1 3
1 3 2 4
4 2 3 1

4 2 3 1
3 1 4 2
1 3 2 4
2 4 1 3

4 3 1 2
2 1 3 4
3 4 2 1
1 2 4 3

4 3 2 1
2 1 4 3
1 2 3 4
3 4 1 2

4 2 1 3
3 1 2 4
2 4 3 1
1 3 4 2

4 1 3 2
3 2 4 1
2 3 1 4
1 4 2 3

4 3 2 1
1 2 3 4
3 4 1 2
2 1 4 3

4 3 1 2
1 2 4 3
2 1 3 4
3 4 2 1

4 1 2 3
3 2 1 4
1 4 3 2
2 3 4 1

4 1 2 3
2 3 4 1
3 2 1 4
1 4 3 2

4 2 3 1
1 3 2 4
2 4 1 3
3 1 4 2

4 2 1 3
1 3 4 2
3 1 2 4
2 4 3 1

4 1 3 2
2 3 1 4
1 4 2 3
3 2 4 1
``````

EDIT

In the original question on StackExchange Proof Assistants is

I’m interested in automating that part of some discrete problems which reduces the search space “without loss of generality” by exploiting various symmetries.

“Graph Rewriting for Natural Deduction and the Proper Treatment of Variables” by Willem Heijltjes (pdf)

is

There are three kinds of transformations: contractions, permutations and simplifications.

• Contractions are used to remove a consecutive introduction and elimination of the same connective.
• Simplifications remove parts of a proof that are unused, which occurs when a closed assumption class of an ∃E- or ∨E-application is empty.
• Permutations are used to shift ∨E- and ∃E applications down over elimination rules to expose other contractions.

… the formal descriptions of the normalization rules can be found in Appendix A.

While the use of contractions, permutations and simplifications is common in Prolog (think term rewriting among other predicates), formal descriptions of the (natural deduction) normalization rules found in Appendix A is not so easily found.

3 Likes

It is typical to normalize latin squares by fixing the first row to be the identity permutation `[1,2,3,4]` or similar. This has the effect of reducing the count by a factor of 4! = 24, which makes it easy to verify the outcomes by inspection. (Transposition of the square across the main diagonal accounts for another factor of two in solutions, but it then becomes attractive to fix the main diagonal as entries `1,2,3,4`.)

Here is a fairly general model using clp(fd). See http://hakank.org/swi_prolog/latin_squares_diagonals.pl for the complete program.

``````
:- use_module(library(clpfd)).
:- use_module(hakank_utils).

%
% Show all solutions for N=4.
%
go :-
N = 4,
new_matrix(N,N,1..N,X),
latin_square_diagonals(X),
flatten(X,Vars),
labeling([ffc,up,bisect],Vars),
print_matrix(X),
fail,
nl.
go.

%
% Find and count the number of solutions for N=6.
%
go2 :-
N = 6,
new_matrix(N,N,1..N,X),
time(findall(X,(latin_square_diagonals(X),
flatten(X,Vars),
labeling([ffc,up,step],Vars)),L)),
length(L,Len),
writeln(len=Len),

nl.
go2.

%
% latin_square_diagonal(X)
%
% Ensure that X is a Latin square as well as the
% constraints that the two diagonals should be distinct.
%
latin_square_diagonals(X) :-
maplist(all_different,X),

transpose(X,XT),
maplist(all_different,XT),

diagonal1_slice(X,Slice1),
all_different(Slice1),

diagonal2_slice(X,Slice2),
all_different(Slice2).
``````

Note that the program uses some of my utilities from http://hakank.org/swi_prolog/hakank_utils.pl

• `new_matrix/4`: define a matrix with a domain
• `diagonal1_slice/2`: extract the elements for main diagonal
• `diagonal2_slice/2`: extract the elements for second diagonal
• `print_matrix/1`: pretty print a matrix

Here are the result from N=1…6 (by running `go2/0`):

``````N   #sols   time
--------------------
1      1           1,972 inferences,  0.000 CPU in  0.000 seconds (100% CPU, 7287805 Lips)
2      0           8,094 inferences,  0.001 CPU in  0.001 seconds (100% CPU, 11839322 Lips
3      0          22,259 inferences,  0.001 CPU in  0.001 seconds (100% CPU, 15460076 Lips)
4     48         205,097 inferences,  0.010 CPU in  0.010 seconds (100% CPU, 21095143 Lips)
5    960       4,728,770 inferences,  0.204 CPU in  0.204 seconds (100% CPU, 23172521 Lips)
6  92160   1,339,542,123 inferences, 59.874 CPU in 59.874 seconds (100% CPU, 22372698 Lips)
``````

The number of solutions for N=7 is 862848000, which would take a little too long time.

3 Likes

Yesterday I didn’t posted this snippet, because is not so interesting, but now, after seen your code, I wonder how I could profile it…
Indeed seems it’s functionally equivalent to your one, but I get 10 times slower execution (at best).
Maybe there is something evident I’m missing ?

``````:- module(latin_square,
[latin_square/2
,latin_square_double_diagonals/2
]).
:- use_module(library(clpfd)).

latin_square(N,LatinSquare) :-
length(LatinSquare,N),
maplist({N}/[Row]>>(length(Row,N),Row ins 1..N),LatinSquare),
maplist(alldiff,LatinSquare),
transpose(LatinSquare,Transposed),
maplist(alldiff,Transposed).

latin_square_double_diagonals(N,LatinSquareWithDiags) :-
latin_square(N,LatinSquareWithDiags),
diagonals(LatinSquareWithDiags,A,B),
alldiff(A),
alldiff(B).

diagonals(Square,A,B) :-
length(Square,N),
bagof(V,N^X^Row^(
between(1,N,X),
matrix_row_col_val(Square,X,X,V)
),A),
bagof(V,N^X^Y^Row^(
between(1,N,X),
Y is N+1-X,
matrix_row_col_val(Square,X,Y,V)
),B).

% evaluate which perform better (all_different vs all_distinct)
% ?- time(latin_square:count_lsd(4,N)).
alldiff(L) :- all_different(L).  % 2,443,915 inferences, 0.125 CPU in 0.125 seconds (100% CPU, 19551320 Lips)
%alldiff(L) :- all_distinct(L).   % 17,511,788 inferences, 1.516 CPU in 1.522 seconds (100% CPU, 11554169 Lips)

matrix_row_col_val(M,R,C,V) :-
nth1(R,M,U),
nth1(C,U,V).

% utility, shortcut a long CLI command...
count_lsd(N,M) :-
latin_square_double_diagonals(N,S),!,
flatten(S,F),
aggregate_all(count,label(F),M).
%why Hakank get so faster execution ?
%aggregate(count,F^labeling([ffc,up,bisect],F),M).

``````

@CapelliC My guess is that the culprit is the calculation of diagonals.

Here’s the time/1 for N=5 using your code:

``````88,982,449 inferences, 3.084 CPU in 3.084 seconds (100% CPU, 28856005 Lips)
``````

After replacing your `diagonals/3` with the following (using the predicates from my utils module: http://hakank.org/swi_prolog/hakank_utils.pl)

``````diagonals(Square,A,B) :-
diagonal1_slice(Square,A),
diagonal2_slice(Square,B).
``````

I got this time/1:

``````4,730,500 inferences, 0.196 CPU in 0.196 seconds (100% CPU, 24172257 Lips)
``````

Compare with the statistics for my original model, i.e. almost the same time and inferences.

`````` 4,728,770 inferences,  0.204 CPU in  0.204 seconds (100% CPU, 23172521 Lips)
``````

Here’s the test predicate:

``````go :-
N = 5,
time(findall(M,(latin_square_double_diagonals(N,M),flatten(M,Vars),labeling([ffc,up,step],Vars)),L)),
length(L,Len),
writeln(len=Len),
nl.
go.
``````
1 Like

Thanks, will try to understand how is possible that a small deterministic predicate like diagonals/3 does so much damage (note the cut in count_lsd). Maybe it introduces unconstrained vars. I checked, it should not, but evidently I missed some details…

Here’s the time when I tested N=5 using your `count_lsd/2` and my diagonal predicates:

`````` 8,040,082 inferences, 0.307 CPU in 0.307 seconds (100% CPU, 26154419 Lips)
``````

I.e. almost twice as many inferences and 0.1s slower than using `findall/2` and `length/2`.

Well, thanks again, seems that nth1/3 (I think) creates constraints.

``````?- time((latin_square(3,S),latin_square:diagonals(S,A,B))).
% 14,299 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
S = [[_A, _B, _C], [_D, _E, _F], [_G, _H, _I]],
A = [_A, _E, _I],
B = [_C, _E, _G],
clpfd:in(_A, ..(1, 3)),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:in(_B, ..(1, 3)),
clpfd:all_different([_B, _E, _H]),
clpfd:in(_C, ..(1, 3)),
clpfd:all_different([_C, _F, _I]),
clpfd:in(_D, ..(1, 3)),
clpfd:all_different([_D, _E, _F]),
clpfd:in(_E, ..(1, 3)),
clpfd:in(_F, ..(1, 3)),
clpfd:in(_G, ..(1, 3)),
clpfd:all_different([_G, _H, _I]),
clpfd:in(_H, ..(1, 3)),
clpfd:in(_I, ..(1, 3)).

?- time((latin_square(3,S),latin_square:diagonals_(S,A,B))).
% 29,077 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
S = [[_A, _B, _C], [_D, _E, _F], [_G, _H, _I]],
A = [_A, _E, _I],
B = [_C, _E, _G],
clpfd:in(_A, ..(1, 3)),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
clpfd:all_different([_A, _B, _C]),
clpfd:all_different([_A, _D, _G]),
...
clpfd:all_different([_G, _H, _I]),
clpfd:all_different([_G, _H, _I]),
clpfd:all_different([_G, _H, _I]),
clpfd:all_different([_G, _H, _I]),
clpfd:in(_H, ..(1, 3)),
clpfd:in(_I, ..(1, 3)).

``````

diagonals_/3 it’s my bagof+nth1 buggy code…

Ah, the infamous matrix_element issue! I didn’t noticed that.

Ideally one would want to `element/3` instead of `nth1/3` in clpfd models, but that is often not possible when working with matrices, since the first index is selecting the (entire) row.

Here are some of the implementations of `matrix_element/4`, that I tend to experiment with. However, they don’t give any performance boost in your model. (The name of them are of historical reason.)

``````matrix_element2(X, I, J, Val) :-
nth1(I, X, Row),
element(J, Row, Val).

matrix_element3(X, I, J, Val) :-
freeze(I, (nth1(I, X, Row),freeze(J,nth1(J,Row,Val)))).

matrix_element5(X, I, J, Val) :-
nth1(I, X, Row),
nth1(J, Row, Val).
``````

Indeed, a pretty subtle issue… I tried to rewrite completely the matrix’ element access by index, to avoid nth1/3 (I hadn’t still inspected your code, as I was more interested in understanding from where such strange behaviour is coming, and I use experimenting with code - where possible - to understand it)

``````/*
matrix_row_col_val(M,R,C,V) :-
nth1(R,M,U),
nth1(C,U,V).
*/
matrix_row_col_val(M,1,1,V) =>
M=[[V|_]|_].
matrix_row_col_val(M,1,C,V), C>1 =>
C1 is C-1,
M=[[_|T]|_],
matrix_row_col_val([T|_],1,C1,V).
matrix_row_col_val(M,R,C,V), R>1 =>
R1 is R-1,
M=[_|T],
matrix_row_col_val(T,R1,C,V).

``````

but it has absolutely no effect… now I’m going to try to postpone the constraint posting…

Seems that postponing the alldiff calls solved the issue: I reverted to the nth1, but changed the ‘driver’

``````latin_square_double_diagonals(N,LatinSquareWithDiags) :-
alloc_matrix(N,LatinSquareWithDiags),
domain_cells(LatinSquareWithDiags,N),
diagonals_(LatinSquareWithDiags,A,B),
mat_diff_rows_cols(LatinSquareWithDiags),
alldiff(A),
alldiff(B).
alloc_matrix(N,Mat) :-
length(Mat,N),
maplist({N}/[Row]>>length(Row,N),Mat).
domain_cells(Mat,N) :-
maplist({N}/[Row]>>(Row ins 1..N),Mat).
mat_diff_rows_cols(Mat) :-
maplist(alldiff,Mat),
transpose(Mat,Transposed),
maplist(alldiff,Transposed).
....

?- time(latin_square:count_lsd(4,N)).
% 152,576 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
N = 48.

?- time(latin_square:count_lsd(5,N)).
% 6,678,914 inferences, 0.453 CPU in 0.453 seconds (100% CPU, 14739672 Lips)
N = 960.

``````

And now I’m left with why your code is working

1 Like

What is your intuition regarding the performance between `bagof/3` and `findall/3`? I’m using `findall/3` deep in the diagonal logic in about the same way you use `bagof/3`

I’m using bagof/3 because - I think - it doesn’t copy_term/2 the Template, while findall/3 should do. Maybe I’m wrong… in my old interpreter, I had only findall/3, and of course I was copying Template every time…

Now I will try if findall/3 instead of bagof does change the behaviour… touching constrained vars can be a pain…

As suspected, findall/3 seems to copy variables:

``````diagonals__(Square,A,B) :-
length(Square,N),
findall(V,(
between(1,N,X),
matrix_row_col_val(Square,X,X,V)
),A),
findall(V,(
between(1,N,X),
Y is N+1-X,
matrix_row_col_val(Square,X,Y,V)
),B).

``````

and the result isn’t valid anymore

``````?- time(latin_square:count_lsd(4,N)).
% 776,688 inferences, 0.047 CPU in 0.047 seconds (100% CPU, 16569344 Lips)
N = 576.
``````

Also I have noticed such difference between findall/3 and bagof/3, which was sensitive to apply “pure literal rule” to set of clauses with prolog variables as object variables. I took time for debugging. In fact I posted it as “nice to know” category.

``````?- bagof(X, member(X, [a(A), b(B)]), S).
S = [a(A), b(B)].
?- findall(X, member(X, [a(A), b(B)]), S).
S = [a(_), b(_)].
``````

Both copy the template. Bagof extends the template with the variables shared with the environment and restores the variable binding by unifying the copied variable template with the original variables. AFAIK all bagof/3 implementations work that way. If you could figure out that a template instantiation is older than the bagof/3 call and not further instantiated by that call you could avoid the copying. That is rather hard. In simple cases you could derive that using static analysis and use a different bagof/3 implementation. I’m not aware of any Prolog system doing that. Maybe there are?

3 Likes

Thanks Jan for the clarification. Could also help me to implement bagof/3 in my old interpreter .

As I experimented with the SSU implementation of matrix_row_col_val/4 above, the core problem does not origin from nth1/3 (or bagof), but (I think) from the unification algorithm.

Modelling with CLP(FD), I often stumble in the problem of detecting the right ‘osmosis point’ at the interface between Prolog and library(clpfd) … my fault for sure, and I remember some discussions with Markus about it. Sorry for the noise, but since I’m going to dive deeper in s(CASP), I need to understand better about such details.

I don’t see it as noise but a side topic of practical value.

If you want I can split off the related replies into a new topic.

If you want that done I would need a reference to the first reply for the new topic, the associated replies, which could be all that follow or identify specific ones with references, and a title for the new topic.

The overall problem is that you create constraints and then use bagof/findall to copy them. The all solutions predicates (or anything that copies terms) do not play too well with constraints. We can fix this using this implementation for finding the diagonals:

``````diagonals(Square,A,B) :-
length(Square,N),
numlist(1,N,Xs),
maplist(matrix_diagonal_one(N, Square), Xs, A),
maplist(matrix_diagonal_two(N, Square), Xs, B).

matrix_diagonal_one(_, Square, X, V) :-
nth1(X, Square, Row),
nth1(X, Row, V).
matrix_diagonal_two(N, Square, X, V) :-
nth1(X, Square, Row),
Y is N+1-X,
nth1(Y, Row, V).
``````

As maplist/4 doesn’t copy, the problem vanishes, just like creating the constraints after the copying. Functional programmers avoid the two helpers using yall

6 Likes