## The Langford Problem

In combinatorial mathematics, a \(\textbf{Langford pairing}\), also called a \(\textbf{Langford sequence}\), is a permutation of the sequence of \(2n\) numbers \(1, 1, 2, 2, ..., n, n\) in which the two \(1\)s are one unit apart, the two \(2\)s are two units apart, and more generally the two copies of each number \(k\) are \(k\) units apart. Langford pairings are named after C. Dudley Langford, who posed the problem of constructing them in 1958. \(\textbf{Langford's problem}\) is the task of finding Langford pairings \(L(n)\) for a given value of \(n\).

### Solution using Constraint Programming

Let \((x_{is}, x_{ie})\) be the tuple of decision variables indicating the starting and ending position of number \(i\) in the Langford sequence. The decision variables need to satisfy the following constraints:

\[ 1 \leq x_{is}, x_{ie} \leq 2n , 1 \leq i \leq n \\\\ x_{ie} = x_{is} + (i + 1) , 1 \leq i \leq n \\\\ \text{All members of the set $\\{x_{it} | 1 \leq i \leq n, t \in {s, e} \\}$ are different}. \]

Here is the Python code implementing the above model using the fantastic Google \(\textbf{OR-Tools}\) library:

```
from ortools.sat.python import cp_model
from collections import defaultdict
def langford_seq_checker(seq):
if not seq:
return False
= defaultdict(list)
pos_map for i, n in enumerate(seq):
pos_map[n].append(i)for n,p in pos_map.items():
if len(p) != 2:
return False
if (p[1] - p[0]) != n + 1:
return False
return True
def langford_cp(n):
= cp_model.CpModel()
model = [[model.NewIntVar(1, 2*n, 'x[%i][%i]' % (i,j)) for j in range(2)] for i in range(n)]
x
for i in range(n) for j in range(2)])
model.AddAllDifferent([x[i][j] for i in range(n):
1] - x[i][0] == i+2)
model.Add(x[i][
= cp_model.CpSolver()
solver = solver.Solve(model)
status if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
= [0]*(2*n+1)
out for i in range(n):
for j in range(2):
= i+1
out[solver.Value(x[i][j])] return out[1:]
else:
return None
```

Here is the sequence \(L(100)\) calculated using the above code:

`[51, 79, 80, 82, 1, 30, 1, 64, 70, 87, 95, 50, 21, 66, 33, 4, 29, 97, 20, 15, 4, 69, 22, 52, 59, 28, 100, 81, 46, 26, 36, 57, 49, 8, 21, 15, 30, 91, 31, 20, 83, 86, 8, 34, 23, 22, 29, 68, 33, 40, 53, 14, 51, 72, 28, 84, 26, 74, 89, 63, 32, 55, 50, 75, 98, 76, 14, 36, 23, 7, 31, 48, 64, 93, 96, 46, 52, 7, 34, 70, 66, 79, 49, 80, 59, 65, 82, 92, 71, 57, 40, 69, 94, 32, 99, 78, 61, 87, 67, 90, 6, 43, 62, 88, 53, 19, 95, 6, 60, 81, 35, 77, 24, 85, 73, 97, 68, 55, 56, 11, 48, 54, 58, 63, 83, 19, 72, 100, 86, 91, 25, 11, 74, 13, 27, 47, 17, 24, 45, 75, 84, 10, 76, 38, 41, 43, 35, 13, 89, 9, 44, 65, 10, 42, 17, 37, 25, 39, 61, 9, 71, 16, 27, 98, 12, 62, 67, 93, 3, 60, 2, 96, 3, 2, 78, 56, 54, 12, 16, 18, 92, 58, 38, 47, 45, 5, 41, 94, 73, 77, 90, 5, 88, 37, 99, 44, 42, 39, 18, 85]`

### Solution using Integer Programming

Another way to solve the Langford problem is to treat it as a set covering problem. To visualize this we make use of the following array for \(L(3)\):

- | 1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|

1 | 1 | 1 | ||||

2 | 1 | 1 | ||||

3 | 1 | 1 | ||||

4 | 1 | 1 | ||||

5 | 2 | 2 | ||||

6 | 2 | 2 | ||||

7 | 2 | 2 | ||||

8 | 3 | 3 | ||||

9 | 3 | 3 |

To solve the problem, we need to select one row for the \(1\)’s in the sequence, one row for the \(2\)’s and one row for the \(3\)’s, such that if we stack these rows on top of each other, no column contains more than one number.

In case of \(L(n)\) it is easy to see that the number of columns in the matrix will be \(2n\) and the number of rows will be \(r = 2n-2 + ... + n-1\). Let \(\{x_i|1 \leq i \leq r\}\) be the set of decision variables, one for each row in the matrix such that \(x_i \in \{0,1\}\). We have the following constraints:

- We choose only \(1\) row among all rows containing the number \(k\) in the matrix where \(1 \leq k \leq n\).
- For each column, we choose only \(1\) row among all rows containing non zero values.

The Python code implementing the above model using the Google \(\textbf{OR-Tools}\) library is given below:

```
from ortools.linear_solver import pywraplp
def langford_ip(n):
= pywraplp.Solver.CreateSolver('SCIP')
solver
= sum(range(n-1, 2*n-1)), 2*n
n_rows, n_cols = [[0 for j in range(n_cols)] for i in range(n_rows)]
matrix = [0 for i in range(n_cols)]
out
# setting up the covering matrix
= 0
j for i in range(n):
for k in range(2*n-i-2):
= i + 1
matrix[j][k] +i+2] = i + 1
matrix[j][k+= 1
j
= [solver.IntVar(0, 1, 'x[%i]' % j) for j in range(n_rows)]
x
# row constraints
= 0
j for i in range(n):
sum([x[k] for k in range(j, j + 2*n-i-2)])==1)
solver.Add(+= 2*n-i-2
j
# column constraints
for i in range(n_cols):
= []
inds for j in range(n_rows):
if matrix[j][i]:
inds.append(j)sum([x[k] for k in inds])==1)
solver.Add(
sum([x[i] for i in range(n_rows)]))
solver.Minimize(
= solver.Solve()
status if status == pywraplp.Solver.OPTIMAL:
for i in range(n_rows):
if x[i].solution_value():
for j in range(n_cols):
+= matrix[i][j]
out[j] return out
else:
return None
```

### Beautiful analytical solution

For the positive cases (\(n = 4k\) or \(4k+3\)) an algorithm for calculating the sequence can be found here. This beautiful algorithm was discovered by Roy Davies in \(1959\).

Here are the details, where \(R\) denotes the reversal of a sequence.

\[ x = Ceiling[n/4] \\\\ \{a, b, c, d\} = \{2x-1, 4x-2, 4x-1, 4x\} \\\\ p = \text{odds in $[1, a-1]$} \\\\ q = \text{evens in $[2, a-1]$} \\\\ r = \text{odds in $[a+2, b-1]$} \\\\ s = \text{evens in $[a+1, b-1]$} \]

If \(4\) divides \(n\), the sequence is \(\{R[s], R[p], b, p, c, s, d, R[r], R[q], b, a, q, c, r, a, d\}\).

If \(n \equiv 3 (\mod 4)\), it is \(\{R[s], R[p], b, p, c, s, a, R[r], R[q], b, a, q, c, r\}\).

The Python code implementing the above algorithm is given below

```
from math import ceil
def R(l):
return list(reversed(l))
def langford_davies(n):
= ceil(n/4)
x = 2*x-1, 4*x-2, 4*x-1, 4*x
a, b, c, d = [i for i in range(1, a) if i % 2==1]
p = [i for i in range(2, a) if i % 2==0]
q = [i for i in range(a+2, b) if i % 2==1]
r = [i for i in range(a+1, b) if i % 2==0]
s if n%4 == 0:
return R(s) + R(p) + [b] + p + [c] + s + [d] + R(r) + R(q) + [b,a] + q + [c] + r + [a, d]
if n%4 == 3:
return R(s) + R(p) + [b] + p + [c] + s + [a] + R(r) + R(q) + [b,a] + q + [c] + r
return None
```