# Generate a Graeco-Latin square

disclaimer: I'm not aware of any non-bruteforce solutions

A Graeco-Latin square is, for two sets of same length $$\n\$$, a $$\n \times n\$$ arrangement of cells, each containing a unique (across the entire square) pair of a element of the first set and a element of the second set, such that all first elements and all second elements of the pairs are unique in their row and column. The most common sets used are, as one could guess, the first $$\n\$$ letters of the Greek and the Latin alphabets.

Here is a picture of a 4x4 Graeco-Latin square:

Graeco-Latin squares are as useful as they sound (the Wikipedia article mentions "design of experiments, tournament scheduling and constructing magic squares"). Your task is, given a positive integer $$\n\$$, to generate a $$\n\times n\$$ Graeco-Latin square.

## Input

A positive integer $$\n > 2\$$; it is guaranteed that a $$\n\times n\$$ Graeco-Latin square exists (that is, $$\n \ne 6\$$).

## Output

A Graeco-Latin square with side length n as a two-dimensional array, a array of arrays, a flattened array or outputted directly.

## Notes

• You do not have to use the Greek and Latin alphabets specifically; for example, outputting pairs of positive integers is allowed as well.
• If you choose to use a alphabet that can't be extended arbitrarily, you have to (theoretically; your code doesn't have to finish before the heat death of the universe) support a maximal side length of at least 20.

This is , so the shortest code wins!

• Sandbox link Jun 4 '19 at 1:52
• Do we have to output a single square, or is it ok to output all possible squares as a list? Jun 4 '19 at 22:40

# 05AB1E, 2623 22 bytes

-3 bytes thanks to Emigna

-1 byte thanks to Kevin Cruijssen

Lãœ.ΔIôDζ«D€í«ε€нÙgQ}P


Try it online!

• n<ÝI‰ can be <Ýã Jun 4 '19 at 9:47
• ...and <Ý can be L. Thanks! Jun 4 '19 at 9:53
• ê}DIùQ can be ÙgQ}P to save a byte. Jun 4 '19 at 12:21
• @KevinCruijssen thanks! I edited that in. Jun 4 '19 at 12:26

# R, 164 148 bytes

-many bytes thanks to Giuseppe.

n=scan()
!=function(x)sd(colSums(2^x))
m=function()matrix(sample(n,n^2,1),n)
while(T)T=!(l=m())|!(g=m())|!t(l)|!t(g)|1-all(1:n^2%in%(n*l+g-n))
l
g


Try it online!

Dramatically inefficient - I think it's even worse than other brute force approaches. Even for n=3, it will probably time out on TIO. Here is an alternate version (155 bytes) which works for n=3 in about 1 second.

Works by rejection. The function m draws a random matrix of integers between $$\1\$$ and $$\n\$$ (without forcing each integer to appear exactly $$\n\$$ times, or in different columns - this explains the slowness of the code, and is the only thing changed in the faster version). Call this function twice, to get the "latin" and "greek" squares l and g. Then we check:

1. all(1:n^2%in%(n*l+g-n)): are there $$\n^2\$$ different pairs of values in l $$\\times\$$ g?
2. are l and g latin squares?

To check that a matrix is a latin square, use the function !. Since point 1. above has been validated, we know that each integer appears $$\n\$$ times in each of l and g. Compute the column-wise sums of 2^l: these sums are all equal if and only if each integer appears once in each column (and the value of the sum is then $$\2^{n+1}-2\$$). If this is true for both l and its transpose t(l), then l is a latin square; same for g. The function sd, which computes the standard deviation, is an easy way to check whether all values of a vector are equal. Note that it doesn't work for the trivial cases $$\n=0\$$ and $$\n=1\$$ , which is OK according to the OP.

A final note: as often in R code golf, I used the variable T, which is initialized as TRUE, to gain a few bytes. But this means that when I needed the actual value TRUE in the definition of m (parameter replace in sample), I had to use 1 instead of T. Similarly, since I am redefining ! as a function different from negation, I had to use 1-all(...) instead of !all(...).

# JavaScript (ES6),  159 147  140 bytes

Returns a flat array of $$\n\times n\$$ pairs of non-negative integers.

This is a simple brute force search, and therefore very slow.

n=>(g=(m,j=0,X=n*n)=>j<n*n?!X--||m.some(([x,y],i)=>(X==x)+(Y==y)>(j/n^i/n&&j%n!=i%n),g(m,j,X),Y=X/n|0,X%=n)?o:g([...m,[X,Y]],j+1):o=m)(o=[])


Try it online! (with prettified output)

### Commented

n => (                      // n = input
g = (                     // g is the recursive search function taking:
m,                      //   m[] = flattened matrix
j = 0,                  //   j   = current position in m[]
X = n * n               //   X   = counter used to compute the current pair
) =>                      //
j < n * n ?             // if j is less than n²:
!X-- ||               //   abort right away if X is equal to 0; decrement X
m.some(([x, y], i) => //   for each pair [x, y] at position i in m[]:
(X == x) +          //     yield 1 if X is equal to x OR Y is equal to y
(Y == y)            //     yield 2 if both values are equal
//     or yield 0 otherwise
>                   //     test whether the above result is greater than:
( j / n ^ i / n &&  //       - 1 if i and j are neither on the same row
j % n != i % n    //         nor the same column
),                  //       - 0 otherwise
//     initialization of some():
g(m, j, X),         //       do a recursive call with all parameters unchanged
Y = X / n | 0,      //       start with Y = floor(X / n)
X %= n              //       and X = X % n
) ?                   //   end of some(); if it's falsy (or X was equal to 0):
o                   //     just return o[]
:                     //   else:
g(                  //     do a recursive call:
[...m, [X, Y]],   //       append [X, Y] to m[]
j + 1             //       increment j
)                   //     end of recursive call
:                       // else:
o = m                 //   success: update o[] to m[]
)(o = [])                   // initial call to g with m = o = []

• 144? (On my phone, so not entirely sure it works) Jun 4 '19 at 11:16
• I don't think you need o, either; you can just return m at the end for 141 Jun 4 '19 at 11:17
• @Shaggy Both versions would fail for $n=5$ (probably because it's the first size where we actually need to backtrack -- but I didn't really check). Jun 4 '19 at 12:26

# Jelly,  21  20 bytes

-1 thanks to Nick Kennedy (flat output option allows a byte save of ż"þẎẎQƑ$Ƈ $$\\rightarrow\$$ F€pZ€QƑƇ) Œ!ṗ⁸Z€Q€ƑƇF€pZ€QƑƇḢ  Try it online! (Too slow for 4 in 60s on TIO, but if we replace the Cartesian power, ṗ, with Combinations, œc, it will complete - although 5 certainly will not!) ### How? Œ!ṗ⁸Z€Q€ƑƇF€pZ€QƑƇḢ - Link: integer, n Œ! - all permutations of [1..n] ⁸ - chain's left argument, n ṗ - Cartesian power (that is, all ways to pick n of those permutations, with replacement, not ignoring order) Z€ - transpose each Ƈ - filter, keeping those for which: Ƒ - invariant under: Q€ - de-duplicate each F€ - flatten each  - use this as both arguments of: p - Cartesian product Z€ - transpose each Ƈ - filter, keeping those for which: Ƒ - invariant under: Q - de-duplicate (i.e. contains all the possible pairs) Ḣ - head (just one of the Latin-Greaco squares we've found)  • Here's a 20. I originally wrote this independently of yours, but ended up with something pretty similar, and then took some inspiration from your use of Cartesian power in place of a permutation dyad, so it's probably best to use it to improve yours. Note you've misspelled Graeco in your explanation. Jun 4 '19 at 23:42 • Thanks Nick, I didn't notice we were allowed to output a flattened version. Jun 5 '19 at 10:11 # Haskell, 207 143 233 bytes (p,q)!(a,b)=p/=a&&q/=b e=filter f n|l<-[1..n]=head$0#[(c,k)|c<-l,k<-l]$[]where ((i,j)%p)m|j==n=[[]]|1>0=[q:r|q<-p,all(q!)[m!!a!!j|a<-[0..i-1]],r<-(i,j+1)%e(q!)p$m]
(i#p)m|i==n=[[]]|1>0=[r:o|r<-(i,0)%p$m,o<-(i+1)#e(notElemr)p$r:m]


Try it online!

OK, I think I finally got it this time. It works fine for n=5, n=6 times out on TIO but I think that might just be because this new algorithm is INCREDIBLY inefficient and basically checks all possibilities until it finds one that works. I'm running n=6 on my laptop now to see if it terminates with some more time.

Thanks again to @someone for pointing out the bugs in my previous versions

• I do not know Haskell, but this seems to error for me when I change the "4" in the footer to 5. Am I invoking this correctly? Jun 5 '19 at 8:13
• @someone Good catch, I should've tested that.I'm actually not sure what's going wrong here, this might take a while to debug Jun 5 '19 at 13:49
• I think this still has a bug; when run for n=5, the tuple (1,1) appears twice. Jun 5 '19 at 14:58
• @someone Man, this problem is a lot harder than I thought. I just can't find a reliable way to lock down all the constraints at once. As soon as I focus on one another one slips out of my grasp. I'm going to mark as non-competing for now until I can find some more time to work on this. Sorry for not testing as thoroughly as I should have Jun 5 '19 at 15:56

import Data.List
import Data.List.Split
t=transpose
c a e h=[s|p<-a,s<-c(a\\[p])(e++[p])h,and[r==nub r|d<-h<$>t s,m<-[d,t d],r<-m]]++[e] f n=head$c[[a,b]|a<-[1..n],b<-[1..n]][]$chunksOf n  Try it online! This is equivalent to a DFS Brute-Force Approach. Times out for n>3 on TIO, so only n=1 and n=3 work correctly. ## Explanation import Data.List -- used for nub, transpose, and difference (\\) import Data.List.Split -- used for chunksOf, shorter than a custom definition using splitAt or take/drop -- we store each pair as [a,b]::[Int] instead of (a,b)::(Int,Int) to save bytes t=transpose -- c :: [[Int]] -> [[Int]] -> Int -> [[[Int]]] -- recursive algorithm that yields a (lazy) list of all possible GL squares using -- a = available pairs left (avoids duplicates) -- e = existing square state, flattened (first length e pairs of the square in LTR reading order) -- h = chunksOf n c a e h=[ s| -- where each possible square comes from p<-a, -- some available pair to add s<-c -- each square comes from recursion with (a\\[p]) -- that pair taken from the list of available pairs (e++[p]) -- that pair added to the existing square state h, and[r==nub r| -- such that every "row" contains no duplicates d<-h<$>t s, -- d<-(list of [the rows of the first element of each pair, the rows of the second element of each pair])
m<-[d,t d], -- m<-(rows of the square, columns of the square)
r<-m -- row r<-m
]
]
++[e] -- base case: only reached if a==[] or if no valid possibilities exist (n==2 || n==6, not required to be handled)
f :: Int -> [[Int]]
f n=head$c -- the first valid GL square found from the initial conditions: [[a,b]|a<-[1..n],b<-[1..n]] -- all possible pairs [] -- empty existing GL square$chunksOf n -- given chunksOf function


## C#, 520506494 484 bytes

class P{static void Main(string[]a){int n=int.Parse(a[0]);int[,,]m=new int[n,n,2];int i=n,j,k,p,I,J;R:for(;i-->0;)for(j=n;j-->0;)for(k=2;k-->0;)if((m[i,j,k]=(m[i,j,k]+ 1) % n)!=0)goto Q;Q:for(i=n;i-->0;)for(j=n;j-->0;){for(k=2;k-->0;)for(p=n;p-->0;)if(p!=i&&m[i,j,k]==m[p,j,k]||p!=j&&m[i,j,k]==m[i,p,k])goto R;for(I=i;I<n;I++)for(J=0;J<n;J++)if(I!=i&&J!=j&&m[i,j,0]==m[I,J,0]&&m[i,j,1]==m[I,J,1])goto R;}for(i=n;i-->0;)for(j=n;j-->0;)System.Console.Write(m[i,j,0]+"-"+m[i,j,1]+" ");}}


The algorithm of findinf a square is very simple. It is... bruteforce. Yeah, it's stupid, but code golf is not about speed of a program, right?

The code before making it shorter:

using System;

public class Program
{
static int[,,] Next(int[,,] m, int n){
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
for (int k = 0; k < 2; k++)
{
if ((m[i, j, k] = (m[i, j, k] + 1) % n) != 0)
{
return m;
}
}
}
}
return m;
}
static bool Check(int[,,] m, int n)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
for (int k = 0; k < 2; k++)
{
for (int p = 0; p < n; p++)
{
if (p != i)
if (m[i, j, k] == m[p, j, k])
return false;
}
for (int p = 0; p < n; p++)
{
if (p != j)
if (m[i, j, k] == m[i, p, k])
return false;
}
}
}
}

for (int i_1 = 0; i_1 < n; i_1++)
{
for (int j_1 = 0; j_1 < n; j_1++)
{
int i_2 = i_1;
for (int j_2 = j_1 + 1; j_2 < n; j_2++)
{
if (m[i_1, j_1, 0] == m[i_2, j_2, 0] && m[i_1, j_1, 1] == m[i_2, j_2, 1])
return false;
}
for (i_2 = i_1 + 1; i_2 < n; i_2++)
{
for (int j_2 = 0; j_2 < n; j_2++)
{
if (m[i_1, j_1, 0] == m[i_2, j_2, 0] && m[i_1, j_1, 1] == m[i_2, j_2, 1])
return false;
}
}
}
}
return true;
}
public static void Main()
{
int n = 3;
Console.WriteLine(n);
int maxi = (int)System.Math.Pow((double)n, (double)n*n*2);
int[,,] m = new int[n, n, 2];
Debug(m, n);
do
{
m = Next(m, n);
if (m == null)
{
Console.WriteLine("!");
return;
}
Console.WriteLine(maxi--);
} while (!Check(m, n));

Debug(m, n);
}

static void Debug(int[,,] m, int n)
{
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
Console.Write(m[i, j, 0] + "-" + m[i, j, 1] + " ");
}
Console.WriteLine();
}
Console.WriteLine();
}
}


Now, if you want to test it with n=3 you wil have to wait like an hour, so here is another version:

public static void Main()
{
int n = 3;
Console.WriteLine(n);
int maxi = (int)System.Math.Pow((double)n, (double)n*n*2);
int[,,] result = new int[n, n, 2];
Parallel.For(0, n, (I) =>
{
int[,,] m = new int[n, n, 2];
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
{
m[i, j, 0] = I;
m[i, j, 1] = I;
}
while (true)
{
m = Next(m, n);
if (Equals(m, n, I + 1))
{
break;
}
if (Check(m, n))
{
Debug(m, n);
}
}
});
}


Update: forgot to remove "public".

Update: used "System." instead of "using System;"; Also, thanks to Kevin Cruijssen, used "a" instead of "args".

Update: thanks to gastropner and someone.

• args can be a :) Jun 4 '19 at 8:02
• Each for loop could be transformed from for(X = 0; X < Y; X++) to for(X = Y; X-->0; ), which should save a byte per loop. Jun 4 '19 at 8:49
• Have you tried the Visual C# Interactive Compiler? It can save bytes. You can also submit a anonymous function. You can also assign i = 0 in the defininion of i and save a byte. Jun 4 '19 at 8:52
• 405 bytes based on @someone's suggestion. Of course it times out after 60 sec on TIO, but it does save bytes by using a lambda and the Interactive Compiler with implicit System. Also, if((m[i,j,k]=(m[i,j,k]+ 1) % n)!=0) can be if((m[i,j,k]=-~m[i,j,k]%n)>0). Jun 4 '19 at 14:29
• @Kevin I don't really feel like reading through that code trying to golf it. Are you sure the printing part works right? It looks like it either should use Write or could save bytes by adding \n to the string inside the call or is otherwise broken. I think you can also return a array directly. Jun 4 '19 at 14:49

# Octave, 182 bytes

Brute force method, TIO keeps timing out and I had to run it a bunch of times to get output for n=3, but theoretically this should be fine. Instead of pairs like (1,2) it outputs a matrix of complex conjugates like 1+2i. This might be stretching the rule a bit, but in my opinion it stll fits the output requirements. There must be a better way to do the two lines under the functino declaration though, but I'm not sure at the moment.

function[c]=f(n)
c=[0,0]
while(numel(c)>length(unique(c))||range([imag(sum(c)),imag(sum(c.')),real(sum(c)),real(sum(c.'))])>0)
a=fix(rand(n,n)*n);b=fix(rand(n,n)*n);c=a+1i*b;
end
end


Try it online!

# Wolfram Language (Mathematica), 123 bytes

P=Permutations
T=Transpose
g:=#&@@Select[T[Intersection[x=P[P@Range@#,{#}],T/@x]~Tuples~2,2<->4],DuplicateFreeQ[Join@@#]&]&


Try it online!

I use TwoWayRule notation Transpose[...,2<->4] to swap the 2nd and 4th dimensions of an array; otherwise this is fairly straightforward.

Ungolfed:

(* get all n-tuples of permutations *)
semiLSqs[n_] := Permutations@Range@n // Permutations[#, {n}] &;

(* Keep only the Latin squares *)
LSqs[n_] := semiLSqs[n] // Intersection[#, Transpose /@ #] &;

isGLSq[a_] := Join @@ a // DeleteDuplicates@# == # &;

(* Generate Graeco-Latin Squares from all pairs of Latin squares *)
GLSqs[n_] :=
Tuples[LSqs[n], 2] // Transpose[#, 2 <-> 4] & // Select[isGLSq];


# Python 3, 271267 241 bytes

Brute-force approach: Generate all permutations of the pairs until a Graeco-Latin square is found. Too slow to generate anything larger than n=3 on TIO.

Thanks to alexz02 for golfing 26 bytes and to ceilingcat for golfing 4 bytes.

Try it online!

from itertools import*
def f(n):
s=range(n);l=len
for r in permutations(product(s,s)):
if all([l({x[0]for x in r[i*n:-~i*n]})*l({x[1]for x in r[i*n:-~i*n]})*l({r[j*n+i][0]for j in s})*l({r[j*n+i][1]for j in s})==n**4for i in s]):return r


Explanation:

from itertools import *  # We will be using itertools.permutations and itertools.product
def f(n):  # Function taking the side length as a parameter
s = range(n)  # Generate all the numbers from 0 to n-1
l = len  # Shortcut to compute size of sets
for r in permutations(product(s, s)):  # Generate all permutations of all pairs (Cartesian product) of those numbers, for each permutation:
if all([l({x[0] for x in r[i * n : (- ~ i) * n]})  # If the first number is unique in row i ...
* l({x[1] for x in r[i * n:(- ~ i) * n]})  # ... and the second number is unique in row i ...
* l({r[j * n + i][0] for j in s})  # ... and the first number is unique in column i ...
* l({r[j * n + i][1] for j in s})  # ... and the second number is unique in column i ...
== n ** 4 for i in s]):  # ... in all columns i:
return r  # Return the square
`
• -26 bytes
– user85686
Jun 5 '19 at 20:14