An NxN Latin square consists of N sets of the numbers 1 to N arranged in such a way that no row or column contains the same number twice. Two Latin squares are orthogonal if no pair of corresponding elements occurs more than once. A set of NxN Latin squares is mutually orthogonal if every pair of Latin squares from the set is orthogonal.
Some of the procedures below require the Affine Planes of Prime Power Order and Resolvability of Affine Planes packages.
/*
* Constructs Q-1 MOLS of prime power order Q
*/
% n_1_mols(Q, MOLS) is true if Q is prime power and MOLS contains Q-1
% mutually orthogonal Latin squares of order Q derived from an affine
% plane of order Q.
% This works for primes, and prime powers less than or equal to 25.
% e.g. n_1_mols(3, X) gives
% X=[[1,2,3,
% 3,1,2,
% 2,3,1],[1,2,3,
% 2,3,1,
% 3,1,2]]
n_1_mols(Q, MOLS):-
affine_plane(Q, Blocks),
parallel_classes(Blocks, [_,_|Classes]),
Q2 is Q * Q,
findall(Point, for(1, Q2, Point), Points),
n_1_mols_1(Classes, Points, MOLS).
% e.g. n_1_mols_1([[[1,5,9],[2,6,7],[3,4,8]],[[1,6,8],[2,4,9],[3,5,7]]],
% [1,2,3,4,5,6,7,8,9], X) gives
% X=[[1,2,3,3,1,2,2,3,1],[1,2,3,2,3,1,3,1,2]]
n_1_mols_1([], _, []).
n_1_mols_1([Class|Classes], Points, [LS|LSs]):-
n_1_mols_2(Points, Class, LS),
n_1_mols_1(Classes, Points, LSs).
% e.g. n_1_mols_2([1,2,3,4,5,6,7,8,9],
% [[1,5,9],[2,6,7],[3,4,8]],
% [1,2,3,3,1,2,2,3,1]).
n_1_mols_2([], _, []).
n_1_mols_2([Point|Points], Class, [L|LS]):-
row_containing_point(Class, Point, 1, L),
n_1_mols_2(Points, Class, LS).
% e.g. row_containing_point([[1,2,3],[4,5,6],[7,8,9]], 9, 1, L) gives
% L=3
row_containing_point([Row|_], Point, L, L):-
member(Point, Row), !.
row_containing_point([_|Class], Point, L0, L):-
L1 is L0 + 1,
row_containing_point(Class, Point, L1, L).
% member(X, Ys) is true if the element X is contained in the list Ys.
% member(X, [X|_]).
% member(X, [_|Ys]):-member(X, Ys).
/*
* Constructs 3 MOLS of order N
*/
% three_mols(N, MOLS) is true if N is not divisible by 2 or 3, and MOLS
% contains 3 mutually orthogonal Latin squares of order N.
% e.g. three_mols(5, X) gives
% X=[[1,2,3,4,5,2,3,4,5,1,3,4,5,1,2,4,5,1,2,3,5,1,2,3,4],
% [1,2,3,4,5,3,4,5,1,2,5,1,2,3,4,2,3,4,5,1,4,5,1,2,3],
% [1,2,3,4,5,5,1,2,3,4,4,5,1,2,3,3,4,5,1,2,2,3,4,5,1]]
three_mols(N, [Ls,Ms,Ns]):-
N > 1,
N mod 2 > 0,
N mod 3 > 0,
findall(Column, column(N, Column), Columns0),
develop(N, Columns0, Columns),
transpose(Columns, [Ls,Ms,Ns]).
% Construct a column
column(N, [X,Y,Z]):-
N1 is N - 1,
for(0, N1, I),
X is I mod N,
Y is (2 * I) mod N,
Z is (4 * I) mod N.
% Develop each column
% e.g. develop(3, [[0,0,0],[1,2,3],[3,1,2]], X) gives
% X=[[1,1,1],[2,2,2],[3,3,3],[2,3,1],[3,1,2],[1,2,3],[1,2,3],[2,3,1],[3,1,2]]
develop(N, Xss, Yss):-
develop_1(Xss, N, [], Yss).
develop_1([], _, Yss, Yss).
develop_1([Xs|Xss], N, Yss0, Yss):-
develop_1(Xss, N, Yss0, Yss1),
develop_2(0, N, Xs, Yss1, Yss).
% e.g. develop_2(0, 3, [1,2,3], [[1,2,3],[2,3,1],[3,1,2]], X) gives
% X=[[2,3,1],[3,1,2],[1,2,3],[1,2,3],[2,3,1],[3,1,2]]
develop_2(N, N, _, Yss, Yss):-!.
develop_2(I, N, Xs, Yss0, [Ys|Yss]):-
I1 is I + 1,
develop_3(Xs, I, N, [], Ys),
develop_2(I1, N, Xs, Yss0, Yss).
% e.g. develop_3([1,2,3], 1, 3, [], X) gives X=[3,1,2]
develop_3([], _, _, Ys, Ys).
develop_3([X|Xs], I, N, Ys0, [Y|Ys]):-
Y is (X + I) mod N + 1,
develop_3(Xs, I, N, Ys0, Ys).
transpose([], [[],[],[]]).
transpose([[A,B,C]|Columns], [[A|Row1],[B|Row2],[C|Row3]]):-
transpose(Columns, [Row1,Row2,Row3]).
/*
* Constructs 3 MOLS of order L*M
*/
% three_mols(L, M, MOLS) is true if L and M are greater than 3, and either
% are not divisible by 2 or 3, or are primes, or are prime powers less
% than or equal to 25.
% e.g. three_mols(4, 5, X)
three_mols(L, M, MOLS):-
three_mols_1(L, Ls),
three_mols_1(M, Ms),
minus_one(Ls, [A,B,C|_]),
minus_one(Ms, [D,E,F|_]),
direct_product(A, D, P),
direct_product(B, E, Q),
direct_product(C, F, R),
are_orthogonal(P, Q),
are_orthogonal(P, R),
are_orthogonal(Q, R),
plus_one([P,Q,R], MOLS).
three_mols_1(L, Ls):-
three_mols(L, Ls), !.
three_mols_1(L, Ls):-
n_1_mols(L, Ls).
% Subtracts 1 from each element of a list of lists
minus_one([], []).
minus_one([Xs|Xss], [Ys|Yss]):-
minus_one_1(Xs, Ys),
minus_one(Xss, Yss).
minus_one_1([], []).
minus_one_1([X|Xs], [Y|Ys]):-
Y is X - 1,
minus_one_1(Xs, Ys).
/*
* Constructs 2 MOLS of order N, where N is equal to 4 mod 6
*/
% two_mols_a(N, MOLS) is true if N is equal to 4 mod 6, and MOLS contains 2
% mutually orthogonal Latin squares of order N.
% e.g. two_mols_a(4, X) gives
% X=[[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]]
two_mols_a(N, [Ls,Ms]):-
N > 0,
4 =:= N mod 6, !,
Q is (N - 1) // 3,
P is N - Q,
retractall(ls1(_,_,_)),
retractall(ls2(_,_,_)),
top_left(P), % Assert PxP top-left identical Latin squares
bottom_right(Q, P), % Assert QxQ bottom-right orthogonal Latin squares
construct_mols(0, Q, P), % Assert NxN mutually orthogonal Latin squares
setof(L, als1(L), Ls0),
setof(M, als2(M), Ms0),
extract(Ls0, Ls),
extract(Ms0, Ms).
% Constructs two copies of a Latin square of order P, formed using integers
% 1,...,P and placed at the top-left of NxN squares
top_left(P):-
P > 1,
for(1, P, I),
for(1, P, J),
L is (I + J - 2) mod P + 1,
assert(ls1(I,J,L)),
assert(ls2(I,J,L)),
fail, !.
top_left(_).
% Constructs two orthogonal Latin squares of order Q, formed using integers
% P+1,...,N and placed at the bottom-right of NxN squares
bottom_right(Q, P):-
R is P + 1,
Q > 0,
Q_1 is Q - 1,
for(0, Q_1, I),
for(0, Q_1, J),
L is (I + J) mod Q + R,
M is (Q + I - J) mod Q + R,
I1 is I + R,
J1 is J + R,
assert(ls1(I1,J1,L)),
assert(ls2(I1,J1,M)),
fail, !.
bottom_right(_, _).
% Used by findall/3
als1([I,J,K]):-ls1(I, J, K).
% Used by findall/3
als2([I,J,K]):-ls2(I, J, K).
% Constructs the Latin squares using Q transversals in each square
construct_mols(Q, Q, _):-!.
construct_mols(J, Q, P):-
J1 is J + 1,
construct_mols_1(0, P, J1),
construct_mols(J1, Q, P).
% In each square, replaces a transversal by a repeated new element, and adds a
% new row and a new column containing the transversal
construct_mols_1(P, P, _):-!.
construct_mols_1(C, P, J):-
R1 is (J + C) mod P + 1,
R2 is (P - J + C) mod P + 1,
X is J + P,
C1 is C + 1,
retract(ls1(R1, C1, T1)), % A transversal element
assert(ls1(R1, C1, X)), % Replace the transversal element
assert(ls1(X, C1, T1)), % New row element
assert(ls1(R1, X, T1)), % New column element
retract(ls2(R2, C1, T2)), !, % A transversal element
assert(ls2(R2, C1, X)), % Replace the transversal element
assert(ls2(X, C1, T2)), % New row element
assert(ls2(R2, X, T2)), % New column element
construct_mols_1(C1, P, J).
% Extracts the elements of the Latin square
extract([], []).
extract([[_,_,L]|Xs], [L|Ls]):-
extract(Xs, Ls).
/*
* Constructs 2 MOLS of order N, where N is not equal to 2 mod 4
*/
% two_mols_b(N, MOLS) is true if N is greater than 2 and is not equal to
% 2 mod 4, and MOLS contains 2 mutually orthogonal Latin squares of
% order N.
% e.g. two_mols_b(3, X) gives X=[[1,2,3,2,3,1,3,1,2], [1,3,2,2,1,3,3,2,1]]
two_mols_b(N, MOLS):-
two_mols_1(N, MOLS0),
plus_one(MOLS0, MOLS).
two_mols_1(4, [Ls,Ms]):-!,
Ls=[0,2,3,1,3,1,0,2,1,3,2,0,2,0,1,3],
Ms=[0,3,1,2,2,1,3,0,3,0,2,1,1,2,0,3].
two_mols_1(8, [Ls,Ms]):-!,
Ls=[0,1,2,3,4,5,6,7,1,0,3,2,5,4,7,6,2,3,0,1,6,7,4,5,3,2,1,0,7,6,5,4,
4,5,6,7,0,1,2,3,5,4,7,6,1,0,3,2,6,7,4,5,2,3,0,1,7,6,5,4,3,2,1,0],
Ms=[0,1,2,3,4,5,6,7,2,3,0,1,6,7,4,5,4,5,6,7,0,1,2,3,6,7,4,5,2,3,0,1,
5,4,7,6,1,0,3,2,7,6,5,4,3,2,1,0,1,0,3,2,5,4,7,6,3,2,1,0,7,6,5,4].
two_mols_1(N, [Ls,Ms]):-
N > 1,
1 =:= N mod 2, !,
two_mols_odd(N, [Ls,Ms]).
two_mols_1(N, [Ls,Ms]):-
N > 8,
is_power_of_2(N), !,
N1 is N // 4,
two_mols_1(4, [Ps,Qs]),
two_mols_1(N1, [Rs,Ss]),
direct_product(Ps, Rs, Ls),
direct_product(Qs, Ss, Ms).
two_mols_1(N, [Ls,Ms]):-
N > 1,
N =\= 2 mod 4, !,
factors(N, 1, N1, N2),
two_mols_1(N1, [Ps,Qs]),
two_mols_1(N2, [Rs,Ss]),
direct_product(Ps, Rs, Ls),
direct_product(Qs, Ss, Ms).
% Constructs two orthogonal Latin squares of order N, where N > 1 and odd.
% e.g. two_mols_odd(3, X) gives X=[[0,1,2,1,2,0,2,0,1],[0,2,1,1,0,2,2,1,0]]
two_mols_odd(N, Lss):-
N > 1,
1 =:= N mod 2,
two_mols_odd_1(0, N, Lss).
two_mols_odd_1(I, N, [[],[]]):-
I =:= N * N, !.
two_mols_odd_1(I, N, [[L|Ls],[M|Ms]]):-
X is I // N,
Y is I mod N,
L is (X + Y) mod N,
M is (N + X - Y) mod N,
I1 is I + 1,
two_mols_odd_1(I1, N, [Ls,Ms]).
% Is true if N is a power of 2
is_power_of_2(1):-!.
is_power_of_2(N):-
N > 1,
0 =:= N mod 2,
N1 is N // 2,
is_power_of_2(N1).
% factors(N, 1, F, G) is true if N=F*G, F is a power of 2, and G is odd.
% e.g. factors(32384, 1, F, G) gives F=128, G=253
factors(N, F, F, N):-
N mod 2 =:= 1, !.
factors(N, F0, F, G):-
F1 is 2 * F0,
N1 is N // 2,
factors(N1, F1, F, G).
% e.g. direct_product([2,0,1,1,2,0,0,1,2],[0,2,3,1,3,1,0,2,1,3,2,0,2,0,1,3],X).
direct_product(Ls, Ms, Ps):-
length(Ls, M0),
length(Ms, N0),
M is sqrt(M0),
N is sqrt(N0),
mn_indices(M, N, Xs),
direct_product_1(Xs, Xs, Ls, Ms, M, N, [], Qs),
indices(Qs, Xs, Ps).
direct_product_1([], _, _, _, _, _, Ps, Ps).
direct_product_1([X|Xs], Ys, Ls, Ms, M, N, Ps0, Ps):-
direct_product_2(Ys, X, Ls, Ms, M, N, Ps0, Ps1),
direct_product_1(Xs, Ys, Ls, Ms, M, N, Ps1, Ps).
direct_product_2([], _, _, _, _, _, Ps, Ps).
direct_product_2([p(J1,J2)|Ys], p(I1,I2), Ls, Ms, M, N, Ps0, Ps):-
I1J1 is I1 * M + J1,
I2J2 is I2 * N + J2,
nth_member(Ls, I1J1, A),
nth_member(Ms, I2J2, B), !,
direct_product_2(Ys, p(I1,I2), Ls, Ms, M, N, [p(A,B)|Ps0], Ps).
% e.g. mn_indices(2, 3, [p(0,0),p(0,1),p(0,2),p(1,0),p(1,1),p(1,2)])
mn_indices(M, N, Ps):-
M1 is M - 1,
N1 is N - 1,
findall(X, for(0, M1, X), Xs),
findall(Y, for(0, N1, Y), Ys),
mn_indices_1(Xs, Ys, Ps, []).
mn_indices_1([], _, Ps, Ps).
mn_indices_1([X|Xs], Ys, Ps0, Ps):-
mn_indices_1(Xs, Ys, Ps1, Ps),
mn_indices_2(Ys, X, Ps0, Ps1).
mn_indices_2([], _, Ps, Ps).
mn_indices_2([Y|Ys], X, [p(X,Y)|Ps0], Ps):-
mn_indices_2(Ys, X, Ps0, Ps).
indices([], _, []).
indices([P|Ps], Is, [Q|Qs]):-
nth_member(Is, Q, P), !,
indices(Ps, Is, Qs).
% Adds 1 to each element of a list of lists
plus_one([], []).
plus_one([Xs|Xss], [Ys|Yss]):-
plus_one_1(Xs, Ys),
plus_one(Xss, Yss).
plus_one_1([], []).
plus_one_1([X|Xs], [Y|Ys]):-
Y is X + 1,
plus_one_1(Xs, Ys).
% for(I, J, K) is true if K is an integer between I and J inclusive.
for(I, J, I):-I =< J.
for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K).
% length(Xs, L) is true if L is the number of elements in the list Xs.
% length(Xs, L):-length_1(Xs, 0, L).
% length_1(Xs, L0, L) is true if L is equal to L0 plus the number of
% elements in the list Xs.
% length_1([], L, L).
% length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).