Steiner Quadruple Systems

This is based on a description given in "A Short Course in Combinatorial Designs" by Ian Anderson and Iiro Honkala.

 * Steiner Quadruple Systems
 * A Steiner quadruple system of order v, denoted SQS(v), is a pair (X,B)
 * where X is a set of cardinality v, and B is a set of 4-subsets of X
 * (called blocks), with the property that any 3-subset of X is contained in
 * a unique block. SQS(v) exists if and only if v=2 or 4 (mod 6).

% sqs(V, SQS) constructs the blocks of SQS(4), SQS(10), SQS(14) and,
%   recursively, SQS(2U) and SQS(3U-2) provided that SQS(U) can be constructed.
%   Valid values of V<100 are:
%   4,8,10,14,16,20,22,28,32,40,44,46,56,58,64,80,82,88,92,94.
sqs( 4, [[1,2,3,4]]):-!.
sqs(10, [[1,2,3,4],[1,2,5,8],[1,2,6,10],[1,2,7,9],[1,3,5,10],[1,3,6,9],
sqs(14, [[1,2,3,4],[1,2,5,6],[1,2,7,8],[1,2,9,10],[1,2,11,12],[1,2,13,14],
sqs(V, SQSv):-
  V > 0,
  0 =:= V mod 2,         % V=2U
  U is V // 2,
  sqs(U, SQSu), !,
  findall(Block, sqs_1(SQSu, U, Block), SQSv0),
  sort(SQSv0, SQSv).
sqs(V, SQSv):-
  V > 4,
  0 =:= (V + 2) mod 3,   % V=3U-2
  U is (V + 2) // 3,
  sqs(U, SQSu),
  sqs(10, SQS10),
  sqs3u2(SQSu, SQS10, [], SQSv1),
  findall(I, for(2, U, I), Is),
  findall(J, for(0, 2, J), Js),
  findall(P, point(Is, Js, P), Ps),
  renumber(SQSv1, [p(0,0)|Ps], SQSv2),
  sort(SQSv2, SQSv).     % Sorts, and removes duplicate blocks

% On backtracking, returns all blocks of SQS(2U)
sqs_1(SQSu, _, Block):-
  member(Block, SQSu).
sqs_1(SQSu, U, Block):-
  add_1(SQSu, U, Ds),
  member(Block, Ds).
sqs_1(_, U, Block):-
  design(U, Ps),
  add(Ps, U, Qs),
  member_2(Ps, Qs, P, Q),
  member(A, P),
  member(B, Q),
  append(A, B, Block).

% e.g. add([[[1,2],[3,4]],[[5,6],[7,8]]], 4, [[[5,6],[7,8]],[[9,10],[11,12]]]).
add([], _, []).
add([Xss|Xsss], Y, [Zss|Zsss]):-
  add_1(Xss, Y, Zss),
  add(Xsss, Y, Zsss).

% e.g. add_1([[1,2],[3,4]], 4, [[5,6],[7,8]]).
add_1([], _, []).
add_1([Xs|Xss], Y, [Zs|Zss]):-
  add_2(Xs, Y, Zs),
  add_1(Xss, Y, Zss).

% e.g. add_2([1,2], 4, [5,6]).
add_2([], _, []).
add_2([X|Xs], Y, [Z|Zs]):-
  Z is X + Y,
  add_2(Xs, Y, Zs).

% e.g. member_2([1,2,3], [4,5,6], X, Y)
member_2([X|_], [Y|_], X, Y).
member_2([_|Xs], [_|Ys], X, Y):-
  member_2(Xs, Ys, X, Y).

% Constructs an sqs(3U-2) containing "pairs" from an sqs(U) and an sqs(10)
sqs3u2([], _, SQSv, SQSv).
sqs3u2([[1|As]|SQSu], SQS10, SQSv0, SQSv):-!,
  % Handle a block of SQSu which contains a 1
  findall(I, for(0, 2, I), Is),
  findall(P, point(As, Is, P), Ps),
  new_blocksA(SQS10, [p(0,0)|Ps], SQSv0, SQSv1),
  sqs3u2(SQSu, SQS10, SQSv1, SQSv).
sqs3u2([As|SQSu], SQS10, SQSv0, SQSv):-
  % Handle a block of SQSu which does not contain a 1
  findall(Zs, ijkl([0,1,2], Zs), Zss),
  new_blocksB(Zss, As, SQSv0, SQSv1),
  sqs3u2(SQSu, SQS10, SQSv1, SQSv).

% e.g. point([1,5,8], [1,2], X) gives on backtracking
%   X=p(A,I) for each element of As and each element of Is
point(As, Is, p(A1,I)):-
  member(A, As),
  A1 is A - 1,
  member(I, Is).

% e.g. new_blocksA([[1,3,5],[1,4,6],[2,3,6],[2,4,5]], 
%                 [p(1,1),p(1,2),p(5,1),p(5,2),p(8,1),p(8,2)], [], X) gives
%   X=[[p(1,1),p(5,1),p(8,1)],[p(1,1),p(5,2),p(8,2)],
%      [p(1,2),p(5,1),p(8,2)],[p(1,2),p(5,2),p(8,1)]]
new_blocksA([], _, SQS, SQS).
new_blocksA([As|SQS10], Ps, SQS0, [Bs|SQS]):-
  renumber_1(Bs, Ps, As),
  new_blocksA(SQS10, Ps, SQS0, SQS).

% e.g. new_blocksB([[0,0,0,0],[0,0,1,2]], [2,4,5,6], [], X) gives
%   X=[[p(2,0),p(4,0),p(5,1),p(6,2)],[p(2,0),p(4,0),p(5,0),p(6,0)]]
new_blocksB([], _, SQS, SQS).
new_blocksB([Zs|Zss], As, SQS0, SQS):-
  new_blocksB_1(As, Zs, Bs),
  new_blocksB(Zss, As, [Bs|SQS0], SQS).

% e.g. new_blocksB_1([2,4,5,6], [0,0,0,0], X) gives
%   X=[p(2,0),p(4,0),p(5,0),p(6,0)]
new_blocksB_1([], [], []). 
new_blocksB_1([A|As], [Z|Zs], [p(A1,Z)|Bs]):-
  A1 is A - 1,
  new_blocksB_1(As, Zs, Bs).

% e.g. ijkl([0,1,2], [0,0,0,0])
ijkl(As, [I,J,K,L]):-
  member(I, As),
  member(J, As),
  member(K, As),
  member(L, As),
  0 =:= (I+J+K+L) mod 3.

% Replaces each pair in Gss by its index in Ps.  Or, replaces each index in
%   Iss by the I-th element of Ps.
renumber([], _, []).
renumber([Gs|Gss], Ps, [Is|Iss]):-
  renumber_1(Gs, Ps, Is),
  renumber(Gss, Ps, Iss).

% e.g. renumber_1([p(1,1),p(5,1),p(8,1)],
%                 [p(1,1),p(1,2),p(5,1),p(5,2),p(8,1),p(8,2)], X)
%   gives X=[1,3,5]
% e.g. renumber_1(X, [p(1,1),p(1,2),p(5,1),p(5,2),p(8,1),p(8,2)], [1,3,5])
%   gives X=[p(1,1),p(5,1),p(8,1)]
renumber_1([], _, []):-!.
renumber_1([G|Gs], Ps, [I|Is]):-
  nth_member(Ps, I, 1, G), !,
  renumber_1(Gs, Ps, Is).

 * Resolvable Designs
 * A (v,k,lambda) design is a pair (X,B) where X is a set of cardinality v,
 * and B is a set of k-subsets of X (called blocks), with the property that
 * any 2-subset of X is contained in exactly lambda blocks. It is required
 * that 2 <= k < v and lambda > 0.
 * The number of blocks is denoted by b. The number of blocks in which each
 * element of X occurs is denoted by r.
 * A (v,k,1) design is resolvable if its blocks can be partitioned into
 * r collections Bi each consisting of b/r = v/k of the blocks, and every
 * element of X appears in exactly one block in Bi for all i. The subsets
 * Bi are called parallel classes.

% Creates the partitions of a (N,2,1) resolvable design for even N
% e.g. design(4, [[[1,4],[2,3]],[[1,3],[2,4]],[[1,2],[3,4]]]).
design(N, Partitions):-
  N mod 2 =:= 0,
  integers(1, N, Ints),
  findall(Comb, comb(2, Ints, Comb), Combs),
  two_i_mod_two_n_minus_1(1, N, Is),
  add_class_index(Combs, N, Is, Yss),
  sort(Yss, Zss),
  resolve(Zss, Partitions).

% Creates 2I mod (2N - 1) for I = 1 to 2N - 1
% e.g. two_i_mod_two_n_minus_1(1, 6, [2,4,1,3,0]).
two_i_mod_two_n_minus_1(N, N, []):-!.
two_i_mod_two_n_minus_1(I, N, [J|Js]):-
  J is (2 * I) mod (N - 1),
  I1 is I + 1,
  two_i_mod_two_n_minus_1(I1, N, Js).

% Puts the index of the parallel class in front of each A,B
add_class_index([], _, _, []).
add_class_index([[A,N]|Xss], N, Is, [[A,A,N]|Yss]):-!,
  % The index is a if b = 2n
  add_class_index(Xss, N, Is, Yss).
add_class_index([[A,B]|Xss], N, Is, [[D,A,B]|Yss]):-
  % Otherwise the index is the value of i for which
  %   2i mod (2n-1) = (a+b) mod (2n-1)
  C is (A + B) mod (N - 1),
  nth_member(Is, D, 1, C), !,
  add_class_index(Xss, N, Is, Yss).

% e.g. resolve([[1,2,3],[1,1,4],[2,2,4],[2,1,3],[3,3,4],[3,1,2]],
%              [[[2,3],[1,4]], [[2,4],[1,3]], [[3,4],[1,2]]]).
resolve([[_|As]], [[As]]):-!.
resolve([[A|As],[A|Bs]|Bss], [[As|Dss]|Dsss]):-!,
  resolve([[A|Bs]|Bss], [Dss|Dsss]).
resolve([[_|As],Bs|Bss], [[As],Dss|Dsss]):-
  resolve([Bs|Bss], [Dss|Dsss]).

 * Library Procedures

/* for(I, J, K) is true if K is an integer between I and J inclusive.        */
for(I, J, I):-I =< J.
for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K).

/* comb(N, Xs, Ys) is true if Ys is a subset of cardinality N resulting      */
/*   from picking N things from the set Xs.                                  */
comb(0, _, []):-!.
comb(I, Xs, [Y|Ys]):-
  I1 is I - 1,
  rest(Y, Xs, Xs1),
  comb(I1, Xs1, Ys).

/* integers(M, N, Is) is true if Is is the list of integers from M to N      */
/*   inclusive.                                                              */
integers(N, N, [N]):-!.
integers(I, N, [I|Is]):-I < N, I1 is I + 1, integers(I1, N, Is).

/* nth_member(+Xs, ?N, ?B, ?X) is true if X is the N-th (base B) element of  */
/*   the list Xs.                                                            */
nth_member(Xs, N, B, X):-nth_member_1(Xs, X, B, N).

nth_member_1([X|_], X, I, I).
nth_member_1([_|Xs], X, I0, I):-
  I1 is I0 + 1,
  nth_member_1(Xs, X, I1, I).

/* rest(X, Ys, Zs) is true if X is a member of the list Ys, and the          */
/*   list Zs is the rest of the list following X.                            */
rest(X, [X|Ys], Ys).
rest(X, [_|Ys], Zs):-rest(X, Ys, Zs).

/* append(Xs, Ys, Zs) is true if Zs is the result of appending the list Xs   */
/*   to the list Ys.                                                         */
%append([], Ys, Ys).
%append([X|Xs], Ys, [X|Zs]):-append(Xs, Ys, Zs).

/* length(Xs, L) is true if L is the number of elements in the list Xs.      */
%length(Xs, L):-length_1(Xs, 0, L).

/* length_1(Xs, L0, L) is true if L is equal to L0 plus the number of        */
/*   elements in the list Xs.                                                */
%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).

/* member(X, Xs) is true if the element X is contained in the list Xs.       */
%member(X, [X|_]).
%member(X, [_|Xs]):-member(X, Xs).

LPA Index     Home Page

Valid HTML 4.0!