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],
[1,3,7,8],[1,4,5,9],[1,4,6,8],[1,4,7,10],[1,5,6,7],[1,8,9,10],
[2,3,5,6],[2,3,7,10],[2,3,8,9],[2,4,5,7],[2,4,6,9],[2,4,8,10],
[2,5,9,10],[2,6,7,8],[3,4,5,8],[3,4,6,7],[3,4,9,10],[3,5,7,9],
[3,6,8,10],[4,5,6,10],[4,7,8,9],[5,6,8,9],[5,7,8,10],[6,7,9,10]]):-!.
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],
[1,3,5,7],[1,3,6,9],[1,3,8,11],[1,3,10,13],[1,3,12,14],[1,4,5,10],
[1,5,9,11],[1,5,8,14],[1,5,12,13],[1,6,8,10],[1,7,10,12],[1,4,6,12],
[1,8,9,12],[1,4,8,13],[1,7,9,13],[1,4,9,14],[1,4,7,11],[1,6,7,14],
[1,6,11,13],[1,10,11,14],[2,3,5,9],[2,3,6,8],[2,3,7,12],[2,3,11,13],
[2,3,10,14],[2,5,12,14],[2,6,9,12],[2,6,11,14],[2,4,8,14],[2,7,9,14],
[2,4,9,11],[2,8,9,13],[2,8,10,12],[2,4,12,13],[2,4,5,7],[2,4,6,10],
[2,5,8,11],[2,5,10,13],[2,6,7,13],[2,7,10,11],[3,6,7,10],[3,7,9,11],
[3,8,9,14],[3,7,8,13],[3,4,7,14],[3,5,11,14],[3,6,13,14],[3,4,5,13],
[3,4,6,11],[3,5,6,12],[3,5,8,10],[3,4,8,12],[3,4,9,10],[3,9,12,13],
[3,10,11,12],[4,5,11,12],[4,5,6,14],[4,5,8,9],[4,6,7,8],[4,6,9,13],
[4,7,9,12],[4,7,10,13],[4,8,10,11],[4,10,12,14],[4,11,13,14],
[5,6,8,13],[5,7,8,12],[5,9,10,12],[5,6,7,9],[5,6,10,11],[5,7,10,14],
[5,7,11,13],[5,9,13,14],[6,7,11,12],[6,8,9,11],[6,8,12,14],
[6,9,10,14],[6,10,12,13],[7,8,9,10],[7,8,11,14],[7,12,13,14],
[8,10,13,14],[8,11,12,13],[9,10,11,13],[9,11,12,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).