# Langford problem

Three different solutions.
integer programming
puzzles
Published

November 1, 2021

## 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)$$:

- 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:

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