Bluskov's Construction

Constructs covering designs using Bluskov's construction. This requires the Transversal Designs and Mutually Orthogonal Latin Squares packages.

```/*
* See "Infinite Classes of Covering Numbers" by Iliya Bluskov
*
* Theorem 2.8:
* Extend each group of a td(n-m,n) with the same additional r points. For each
*   extended group form two new blocks of size n-m and a third block of size
*   2(r+m). Combine sets of the third blocks into blocks of size at most n-m
*   containing all pairs from the corresponding third blocks.
*/

% If n is a prime power, m is an integer such that n>=3m, and r is an integer
%   such that 0<=r<=(n-3m)/2, then
%   c(n^2-mn+r,n-m,2,2) <= n^2+2(n-m)+ceil((n-m)/floor((n-m-r)/(2m+r)))
%   But see restrictions on n_1_mols/2.
% e.g. bluskov(3, 0, 1) outputs
%   c(10,3,2,2,17)
%   [[1,2,3],[1,4,7],[1,4,10],[1,5,8],[1,6,9],[1,7,10],[2,3,10],[2,4,9],[2,5,7],
%    [2,6,8],[3,4,8],[3,5,9],[3,6,7],[4,5,6],[5,6,10],[7,8,9],[8,9,10]]
bluskov(N, M, R):-
0 =< M, M =< N/3,
0 =< R, R =< (N-3*M)/2,
M+R > 0, % Avoids division by zero
N_M is N - M,
n_1_mols(N, MOLS),
T is N_M - 2,
mols_oa(MOLS, T, N, OA),
oa_td(OA, Groups, Blocks),
J is N*(N-M)+1,
L is J+R-1,
findall(X, for(J, L, X), Xs), % Xs=[] if J>L, that is, R=0
bluskov_1(Groups, Xs, M, R, Blocks2, Blocks, Blocks1),
S is (N-M-R) // (M+M+R),
unite(Blocks2, Xs, S, S, [], Blocks3),
fill_all(Blocks3, N_M, Blocks1, Cover0),
sort(Cover0, Cover),
length(Cover, B),
V is N*(N-M)+R,
schonheim(V, N_M, 2, B),             % Solutions = Schonheim's lower bound
%  schonheim(V, N_M, 2, Sch), B > Sch, % Solutions > Schonheim's lower bound
write(c(V,N_M,2,2,B)), nl,
write(Cover), nl.

bluskov_1([], _, _, _, [], Yss, Yss).
bluskov_1([As|Ass], Xs, M, R, [Zs|Zss], Yss0, Yss):-
M_R is M+R,
M_M_R is M+M+R,
split([M_R], As, [_,As1]),
append(As1, Xs, Ys1),             % Ys1=a(r+m+1),...,a(n),x(1),...,x(r)
split([M_R,M], As, [As2,_,As3]),
append(As2, As3, Ys2),            % Ys2=a(1),...,a(r+m),a(r+2m+1),...,a(n)
split([M_M_R], As, [Zs,_]),       % Zs=a(1),...,a(r+2m) without x(1),...,x(r)
bluskov_1(Ass, Xs, M, R, Zss, [Ys1,Ys2|Yss0], Yss).

% e.g. unite([[1,2,3],[9,10,11],[17,18,19]], [57], 2, 2, [], X) gives
%   X=[[1,2,3,9,10,11,57],[17,18,19,57]]
unite([], Xs, _, _, Zs0, [Zs]):-!,
append(Zs0, Xs, Zs).
unite(Yss, Xs, 0, N0, Zs0, [Zs|Zss]):-!,
append(Zs0, Xs, Zs),
unite(Yss, Xs, N0, N0, [], Zss).
unite([Ys|Yss], Xs, N, N0, Zs, Zss):-
append(Zs, Ys, Zs1),
N1 is N - 1,
unite(Yss, Xs, N1, N0, Zs1, Zss).

% Adds distinct unused integers starting at 1 to each ordered list so that
%   the length of each is I.
% e.g. fill_all([[2,4,5]], 6, [], [[1,2,3,4,5,6]])
fill_all([], _, Yss, Yss).
fill_all([Xs|Xss], I, Yss0, Yss):-
length(Xs, L),
J is I-L,
fill(J, 1, Xs, Ys),
fill_all(Xss, I, [Ys|Yss0], Yss).

% Adds J distinct unused integers starting at I to the ordered list.
% e.g. fill(3, 1, [2,4,5], [1,2,3,4,5,6])
fill(0, _, Xs, Xs):-!.
fill(J, I, [X|Xs], [X|Ys]):-
I >= X, !,
I1 is I + 1,
fill(J, I1, Xs, Ys).
fill(J, I, Xs, [I|Ys]):-
J1 is J - 1,
I1 is I + 1,
fill(J1, I1, Xs, Ys).

% e.g. split([2,2,1], [1,2,3,4,5,6,7], [[1,2],[3,4],[5],[6,7]]).
split([], Xs, [Xs]).
split([N|Ns], Xs, [Ys|Yss]):-
split_1(N, Xs, Ys, Xs1),
split(Ns, Xs1, Yss).

split_1(0, Xs, [], Xs):-!.
split_1(N, [X|Xs], [X|Ys], Xs1):-
N1 is N - 1,
split_1(N1, Xs, Ys, Xs1).

% The Schonheim lower bound
% L = ceil((v/k)ceil(((v-1)/(k-1))...ceil((v-t+1)/(k-t+1))))
schonheim(V, K, T, L):-
V >= K, K >= T,
V1 is V - T + 1,
K1 is K - T + 1,
schonheim_1(T, V1, K1, 1, L).

schonheim_1(0, _, _, L0, L):-!,
L=L0.
schonheim_1(T, V, K, L0, L):-
T1 is T - 1,
V1 is V + 1,
K1 is K + 1,
L1 is (L0 * V) / K,
ceil(L1, L2),
schonheim_1(T1, V1, K1, L2, L).

/* ceil(R, I) is true if I is the smallest integer not less than the         */
/*   positive number R.  Note that conversion from real to integer rounds.   */
ceil(X, Y):-X is int(X), !, Y=X.
ceil(X, Y):-Y is int(X) + 1.
```