Latin squares are very interesting, thank you for sharing this!

For comparison, here is a Prolog solution that uses CLP(ℤ), constraint logic programming over integers, to express the task:

    :- 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).

To generate random solutions, you only need to implement a suitable labeling strategy, or use one of the built-in strategies if available.

For example, here are two different solutions for N = 5:

    ?- 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].

and:

    ?- 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].

Using the Prolog formulations, you can also generate solutions for much larger squares. For example, here is a solution for N = 12:

    ?- 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)

For best results, try a Prolog system with a fast constraint solver.