Latin squares are very interesting, thank you for sharing this!<p>For comparison, here is a Prolog solution that uses CLP(ℤ), constraint logic programming over integers, to express the task:<p><pre><code> :- use_module(library(clpz)).
latin_square(N, Rows) :-
length(Rows, N),
maplist(same_length(Rows), Rows),
append(Rows, Vs),
Vs ins 1..N,
maplist(all_distinct, Rows),
transpose(Rows, Columns),
maplist(all_distinct, Columns).
</code></pre>
To generate <i>random</i> solutions, you only need to implement a suitable labeling strategy, or use one of the built-in strategies if available.<p>For example, here are two different solutions for N = 5:<p><pre><code> ?- latin_square(5, Rows),
append(Rows, Vs),
labeling([ff, random_value(0)], Vs),
maplist(portray_clause, Rows).
[4, 2, 1, 3, 5].
[5, 1, 3, 4, 2].
[1, 4, 2, 5, 3].
[2, 3, 5, 1, 4].
[3, 5, 4, 2, 1].
</code></pre>
and:<p><pre><code> ?- latin_square(5, Rows),
append(Rows, Vs),
labeling([ff, random_value(1)], Vs),
maplist(portray_clause, Rows).
[4, 3, 5, 2, 1].
[5, 1, 3, 4, 2].
[2, 5, 1, 3, 4].
[3, 4, 2, 1, 5].
[1, 2, 4, 5, 3].
</code></pre>
Using the Prolog formulations, you can also generate solutions for much larger squares. For example, here is a solution for N = 12:<p><pre><code> ?- time((latin_square(12, Rows),
append(Rows, Vs),
labeling([ff], Vs),
maplist(portray_clause, Rows))).
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12].
[2, 12, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11].
[3, 1, 11, 12, 2, 4, 5, 6, 7, 8, 9, 10].
[4, 3, 9, 10, 11, 12, 1, 2, 5, 6, 7, 8].
[5, 4, 2, 7, 8, 11, 12, 9, 10, 1, 3, 6].
[6, 5, 4, 2, 12, 8, 9, 10, 3, 11, 1, 7].
[7, 6, 5, 1, 3, 2, 10, 11, 12, 4, 8, 9].
[8, 7, 6, 5, 9, 10, 3, 12, 11, 2, 4, 1].
[9, 8, 7, 6, 10, 1, 11, 3, 2, 12, 5, 4].
[10, 9, 8, 11, 1, 3, 2, 4, 6, 7, 12, 5].
[11, 10, 12, 8, 7, 9, 4, 5, 1, 3, 6, 2].
[12, 11, 10, 9, 6, 7, 8, 1, 4, 5, 2, 3].
% 8,146,020 inferences, 2.212 CPU in 2.283 seconds (97% CPU, 3683380 Lips)
</code></pre>
For best results, try a Prolog system with a fast constraint solver.