Langford problem

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):
        model.Add(x[i][1] - x[i][0] == i+2)

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

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.Add(sum([x[k] for k in inds])==1)

    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