## 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
pos_map = defaultdict(list)
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):
model = cp_model.CpModel()
x = [[model.NewIntVar(1, 2*n, 'x[%i][%i]' % (i,j)) for j in range(2)] for i in range(n)]

model.AddAllDifferent([x[i][j] for i in range(n) for j in range(2)])
for i in range(n):

solver = cp_model.CpSolver()
status = solver.Solve(model)
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
out = [0]*(2*n+1)
for i in range(n):
for j in range(2):
out[solver.Value(x[i][j])] = i+1
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)$:

-123456
111
211
311
411
522
622
722
833
933

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:

1. We choose only $1$ row among all rows containing the number $k$ in the matrix where $1 \leq k \leq n$.
2. 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):
solver = pywraplp.Solver.CreateSolver('SCIP')

n_rows, n_cols = sum(range(n-1, 2*n-1)), 2*n
matrix = [[0 for j in range(n_cols)] for i in range(n_rows)]
out = [0 for i in range(n_cols)]

# setting up the covering matrix
j = 0
for i in range(n):
for k in range(2*n-i-2):
matrix[j][k] = i + 1
matrix[j][k+i+2] = i + 1
j += 1

x = [solver.IntVar(0, 1, 'x[%i]' % j) for j in range(n_rows)]

# row constraints
j = 0
for i in range(n):
solver.Add(sum([x[k] for k in range(j, j + 2*n-i-2)])==1)
j += 2*n-i-2

# column constraints
for i in range(n_cols):
inds = []
for j in range(n_rows):
if matrix[j][i]:
inds.append(j)

solver.Minimize(sum([x[i] for i in  range(n_rows)]))

status = solver.Solve()
if status == pywraplp.Solver.OPTIMAL:
for i in range(n_rows):
if x[i].solution_value():
for j in range(n_cols):
out[j] += matrix[i][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):
x = ceil(n/4)
a, b, c, d = 2*x-1, 4*x-2, 4*x-1, 4*x
p = [i for i in range(1, a) if i % 2==1]
q = [i for i in range(2, a) if i % 2==0]
r = [i for i in range(a+2, b) if i % 2==1]
s = [i for i in range(a+1, b) if i % 2==0]
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