Section 3.3 of [Generalised 2-designs with Block Size 3(Andy L. Drizen)] describes a method of generating Latin Squares of order n attributed to Jacobson and Matthews. The purpose of this task is to produce a function which given a valid Latin Square transforms it to another using this method.

part 1
Use one of the 4 Latin Squares in reduced form of order 4 as X0 to generate 10000 Latin Squares using X(n-1) to generate X(n). Convert the resulting Latin Squares to their reduced form, display them and the number of times each is produced.

part 2
As above for order 5, but do not display the squares. Generate the 56 Latin Squares in reduced form of order 5, confirm that all 56 are produced by the Jacobson and Matthews technique and display the number of each produced.

part 3
Generate 750 Latin Squares of order 42 and display the 750th.

part 4
Generate 1000 Latin Squares of order 256. Don't display anything but confirm the approximate time taken and anything else you may find interesting

function shuffleCube(sequence c)
    integer n = length(c), rx, ry, rz
    bool bProper = true
    while true do
        rx = rand(n)
        ry = rand(n)
        rz = rand(n)
        if c[rx][ry][rz] == 0 then exit end if
    end while
    while true do
        integer ox, oy, oz
        for i=1 to n do
            ox = i
            if c[ox][ry][rz] == 1 then exit end if
        end for
        if not bProper and rand(2)==2 then
            for i=ox+1 to n do
                ox = i
                if c[ox][ry][rz] == 1 then exit end if
            end for
        end if
        for i=1 to n do
            oy = i
            if c[rx][oy][rz] == 1 then exit end if
        end for
        if not bProper and rand(2)==2 then
            for i=oy+1 to n do
                oy = i
                if c[rx][oy][rz] == 1 then exit end if 
            end for
        end if
        for i=1 to n do
            oz = i
            if c[rx][ry][oz] == 1 then exit end if
        end for
        if not bProper and rand(2)==2 then
            for i=oz+1 to n do
                oz = i
                if c[rx][ry][oz] == 1 then exit end if
            end for
        end if
 
        c[rx][ry][rz] += 1
        c[rx][oy][oz] += 1
        c[ox][ry][oz] += 1
        c[ox][oy][rz] += 1
 
        c[rx][ry][oz] -= 1
        c[rx][oy][rz] -= 1
        c[ox][ry][rz] -= 1
        c[ox][oy][oz] -= 1
 
        if c[ox][oy][oz] < 0 then
            {rx, ry, rz} = {ox, oy, oz}
            bProper = false
        else
            bProper = true
            exit
        end if
    end while
    return c
end function
 
function toMatrix(sequence c)
    integer n = length(c)
    sequence m = repeat(repeat(0,n),n)
    for i=1 to n do
        for j=1 to n do
            for k=1 to n do
                if c[i][j][k] != 0 then
                    m[i][j] = k
                    exit
                end if
            end for
        end for
    end for
    return m
end function
 
function toReduced(sequence m)
    integer n := length(m)
    m = deep_copy(m)
    for j=1 to n-1 do
        if m[1][j]!=j then
            for k=j+1 to n do
                if m[1][k]==j then
                    for i=1 to n do
                        atom mij = m[i][j],
                             mik = m[i][k]
                        m[i][j] = mik
                        m[i][k] = mij

                    end for
                    exit
                end if
            end for
        end if
    end for
    for i=2 to n-1 do
        if m[i][1]!=i then
            for k=i+1 to n do
                if m[k][1]==i then
                    for j=1 to n do
                        atom mij = m[i][j],
                             mkj = m[k][j]
                        m[i][j] = mkj
                        m[k][j] = mij
                    end for
                    exit
                end if
            end for
        end if
    end for
    return m
end function
 
function makeCube(object orig, integer n)
    sequence c = repeat(repeat(repeat(0,n),n),n)
    for i=1 to n do
        for j=1 to n do
            integer k = iff(orig==NULL?mod(i+j,n)+1:orig[i][j])
            c[i][j][k] = 1
        end for
    end for
    return c
end function
 
procedure main()
 
    printf(1,"Part 1: 10,000 latin Squares of order 4 in reduced form:\n\n")
    sequence orig = {{1, 2, 3, 4}, {2, 1, 4, 3}, {3, 4, 1, 2}, {4, 3, 2, 1}},
             c := makeCube(orig, 4), m, rm, fk
    integer freq = new_dict()
    for i=1 to 10000 do
        c = shuffleCube(c)
        m = toMatrix(c)
        rm = toReduced(m)
        setd(rm,getd(rm,freq)+1,freq)
    end for
    fk = getd_all_keys(freq)
    for i=1 to length(fk) do
        printf(1,"%v occurs %d times\n", {fk[i],getd(fk[i],freq)})
    end for
 
    printf(1,"\nPart 2: 10,000 latin squares of order 5 in reduced form:\n\n")
    orig = {{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}}
    c = makeCube(orig, 5)
    destroy_dict(freq, justclear:=true)
    for i=1 to 10000 do
        c = shuffleCube(c)
        m = toMatrix(c)
        rm = toReduced(m)
        setd(rm,getd(rm,freq)+1,freq)
    end for
    fk = getd_all_keys(freq)
    for i=1 to length(fk) do
        fk[i] = sprintf("%2d(%3d)", {i,getd(fk[i],freq)})
    end for
    puts(1,join_by(fk,8,7," ","\n"))
    destroy_dict(freq)
 
    -- part 3
    printf(1,"\nPart 3: 750 latin squares of order 42, showing the last one:\n\n")
    c = makeCube(NULL, 42)
    for i=1 to 750 do
        c = shuffleCube(c)
    end for
    m = toMatrix(c)
    integer n := length(m)
    for i=1 to n do
        for j=1 to n do
            m[i,j] = sprintf("%2d",m[i,j])
        end for
        m[i] = join(m[i]," ")
    end for
    printf(1,"%s\n",join(m,"\n"))
 
    -- part 4
    printf(1,"\nPART 4: 1000 latin squares of order 256:\n\n")
    atom t0 = time()
    c = makeCube(NULL, 256)
    for i=1 to 1000 do
        c = shuffleCube(c)
    end for
    printf(1,"Generated in %s\n", elapsed(time()-t0))
end procedure
main()
