1function r = delsarte(n, d)
 2    % https://www.win.tue.nl/~aeb/codes/binary-1.html
 3
 4    if d > n
 5        r = 1;
 6        return
 7    end
 8
 9    if 1.5 * d > n
10        r = 2;
11        return
12    end
13    
14    m = n - d + 1;
15
16    % populate the objective function
17    c = -ones(1, m);
18    
19    A = zeros(n, m);
20    b = zeros(1, n);
21    for t = 1:n
22        % populate the A matrix
23        for j = 1:m
24            i = n - m + j;
25            A(t, j) = - k(t, n, i);
26        end
27        % populate the b vector
28        b(1, t) = k(t, n, 0);
29    end
30
31    % populate non-negativity constraints
32    lb = zeros(1, m);
33    
34    % fill other parameters
35    Aeq = [];
36    beq = [];
37
38    % find an optimal solution
39    s = sym(linprog(c, A, b, Aeq, beq, lb));
40
41    % add the value of x_0 and round the result
42    r = fix(sum(s) + 1);
43end