[Question]
  [
 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:[![enter image description here](https://i.stack.imgur.com/jVIJv.png)](https://i.stack.imgur.com/jVIJv.png)
Graeco-Latin squares are as useful as they sound (the [Wikipedia article](https://en.wikipedia.org/wiki/Graeco-Latin_square) 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 [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'"), so the shortest code wins!
  
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), ~~26~~ ~~23~~ 22 bytes
*-3 bytes thanks to Emigna*
*-1 byte thanks to Kevin Cruijssen*
```
Lãœ.ΔIôDζ«D€í«ε€нÙgQ}P
```
[Try it online!](https://tio.run/##yy9OTMpM/f/f5/Dio5P1zk3xPLzF5dy2Q6tdHjWtObz20OpzW4GMC3sPz0wPrA34XwuUPrT7vzEA "05AB1E – Try It Online")
[Answer]
# [R](https://www.r-project.org/), ~~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!](https://tio.run/##TY3BCoMwEETv@Yp6EHbbeNDeCvkKPYvB1jSwWcUoVWq/PdW2YC/D8GaG6UOjmpHrwbYMszz4WvMOcMZnYLVBQFFF1R5N6K9Qt5SPzkNWTojC/Q2dHno7gdeuoxuw5DKTKUpG8bjblRRYqAhIOUBcIjA/MwB91OCSJpoI0su6jC3HwEc6mYTXHxImdL3lrSy@xqB4iQbOGN4 "R – 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](https://tio.run/##TYzBDoIwDIbvewo5kLQKB/BmwlPImbBMmEu6QjaIEPHZJ6gJXpqv39/@LrRFO7IaTMcwJwevJO8CZ3wGLjYJKOqo3qMJ/Q1UR9fResirCVHYv0crB2cmcE1PRsmhAU68tD2tgJgwisfdrEuJZREBFRYQlwj0Dwagz9S4ZKkkguzCVR4bjoGPdNLpWiJI6NA7w9ux@IJG8RItnDG8AQ) 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(...)`.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~21~~ 20 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
-1 thanks to Nick Kennedy (flat output option allows a byte save of `ż"þ`ẎẎQƑ$Ƈ` \$\rightarrow\$ `F€p`Z€QƑƇ`)
```
Œ!ṗ⁸Z€Q€ƑƇF€p`Z€QƑƇḢ
```
**[Try it online!](https://tio.run/##y0rNyan8///oJMWHO6c/atwR9ahpTSAQH5t4rN0NSBckgEVA3Ic7Fv3//98YAA "Jelly – Try It Online")** (Too slow for `4` in 60s on TIO, but if we replace the Cartesian power, `ṗ`, with Combinations, `œc`, [it will complete](https://tio.run/##y0rNyan8///oJMWjk5MfNe6IetS0JhCIj0081u4GpAsSwCIg7sMdi/7//28CAA "Jelly – Try It Online") - although 5 certainly will not!)
### How?
```
Œ!ṗ⁸Z€Q€ƑƇF€p`Z€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)
```
[Answer]
# 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!](https://tio.run/##LY3BboMwEETv/YrNgWDXBlK1p8CS3zCiVIoSQLbCuoK0ApV@O10rPaxGu292xp2/z9NltJ/3hPy13TrcCEvRoxi0w4M2SM8ksXQF62lnkmRdh3TyQytEPeul0ZapMIizVKJCXGQpXEYfNqP93kW0QxuR1H3I00bqCk1GKwdHSPLkj72o0zQddG101TTaqRd59DhI4bFu5Nb5URAgvOZAUCC8sSol4ecJYGwnJp0gmfMWnP7rzpc41mBZF55DzlIA5RCYYvhOjJeQET7mh2d@eOZw/jdyem2VakBBDHFouHia/K1Nb74PTdz6u/0B "JavaScript (Node.js) – 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 = []
```
[Answer]
# [Haskell](https://www.haskell.org/), ~~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(`notElem`r)p$r:m]
```
[Try it online!](https://tio.run/##Vc0xb4QgFAfwuX4KiN4FUuR06GKkW7d26khIj7vDiAIienHhu1tMpy4veS///@/1chmVMfuOPJkxRJLcMPMXJs/n@cJumWKdNqsKWQdcNG3Ja0qdYL2Sj6LKObqTEcd7WxoypiEKLrZeBZW9IKTJgE8e2zgw5hjnQsT6vWJ8bkKc29ITaQyaIeYWQgnhEGXiK0p1WQtBQlsewmuNT@pI@cKKpOr8EPV/MTRT/MtX6WNhyXQsqZordHXT@mGUvYZEhMaK3UrtAANW@q8fgPxz/V7Dp6NPt03hsdB0B0s/bRgUoANv@y8 "Haskell – 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
[Answer]
# [Haskell](https://www.haskell.org/), 187 bytes
```
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!](https://tio.run/##XcuxbgMhDAbgnafwcEOiwElVVujUMVWHjMSD76AC5SAInO3e/UIqRZG6@P9tfQ7Urn5Zti2mcqsMX8Q0nmJj8f8wnssSWbDhSrmVW/NiBgIPwdi2Fq1INq3mHV0utuB@5w@HZwZJ2dlqTL5PUFenVdDDJ0OTSSvrJINDWbVKiP3Bo/iFbIInN8zWkpxwpe4@xjGjnF4NLQ5zuOdr@@l8Y98YMhhIVL6h1JgZBniDvvQpEsUn@tPH7QE "Haskell – 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
```
[Answer]
## C#, 520 506 494 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](https://codegolf.stackexchange.com/users/52210/kevin-cruijssen), used "a" instead of "args".
Update: thanks to [gastropner](https://codegolf.stackexchange.com/users/75886/gastropner) and [someone](https://codegolf.stackexchange.com/users/36445/someone).
[Answer]
# [Octave](https://www.gnu.org/software/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!](https://tio.run/##ZY4xDoMwEAR7fyR3BEWOUiLzEeTCmMNYsg@F4ISCvzumCUWKLWZWWu1sV/OmnMfEdvUzd1arERiFVZ2spRafyQcCTpECWGwDsVsnSOyfiYrAfV8MO4LOR@PgleIh6xNul4ILmfDrTjg6ja1EYdToNyhLA3DNWDE2/b@yylzvvuobQTwcyVs5@8D8BQ "Octave – Try It Online")
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 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!](https://tio.run/##VY1BS8QwEIXv@RUDhR4kUlx7Wo3MgiiKYHRzGwYJy1ALbbomKSwsu3@9Rvbk6b0HH98bff6W0ed@55fFGitxnHOZU0jKGRd9SPspierWpqoRtzLILpOjl5AlptILSQdjyeKnD51gpY/VibVr8MBnN@8HSeeVXt1fP7SsH8suT1meosgHvU59QKy45nqxsQ@ZOmxZKdWn57ftD/kvhrWBPwwQwUPTwH/FnXoPuElJYlaXoM0wuDgLHW/0rW5PGi4ywCvomJdf "Wolfram Language (Mathematica) – 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];
```
[Answer]
# [Python 3](https://docs.python.org/3/), ~~271~~ ~~267~~ 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](https://codegolf.stackexchange.com/users/85686/alexz02) for golfing 26 bytes and to [ceilingcat](https://codegolf.stackexchange.com/users/52904/ceilingcat) for golfing 4 bytes.
[Try it online!](https://tio.run/##fU/LasMwELzrK/aoRws16YO6@EuMDiaRGwV5JVYrSCntr7uS0wbaQy/L7szOMJPe@Bhxt64zxQU8O@IYQwa/pEisxcHNMEtUvYA80ISvrh4vYQgOBcyRgMAjJEdL4Yl9xCwTxUPZs8w3WTUd@BmmEOQY5Pt5vLNNdW4qGr3G/vazTvuh9EZ3/9E0njQab79NTu0r/6G6X9QwoNb3DfEbYlVPjgtV91UgDLATIhZOheveeorr70/Z1iCRR5Z63B9Jtgrm@UmZy9FZ8/ig4Jr64rZFl950qoa3Sqxf "Python 3 – 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
```
]