bibd(21):-!,
develop(21, [[0,1,4,14,16]], Blocks),
write(c(21,5,2,2,21)), nl,
write(Blocks), nl.
bibd(41):-!,
develop(41, [[0,1,4,11,29],[0,2,8,17,22]], Blocks),
write(c(41,5,2,2,82)), nl,
write(Blocks), nl.
bibd(61):-!,
develop(61, [[0,1,3,13,34],[0,4,9,23,45],[0,6,17,24,32]], Blocks),
write(c(61,5,2,2,183)), nl,
write(Blocks), nl.
bibd(81):-!,
develop(81, [[0,1,3,7,33],[0,5,20,28,39],[0,9,21,52,65],[0,10,24,46,64]], Blocks),
write(c(81,5,2,2,324)), nl,
write(Blocks), nl.
% If Q is a prime power satisfying Q==1 (mod 4), then there exists a
% (V,5,1)-BIBD where V=5Q.
% See "On the Existence of (v,5,1)-BIBDs" by Clayton Smith.
% Note: Only Q=9 and prime Q are handled here.
% e.g. bibd(85) constructs an (85,5,1)-BIBD, i.e. c(85,5,2,2,357)
bibd(V):-
0 =:= V mod 5,
Q is V // 5,
gdd(Q, Gs, Bs),
append(Gs, Bs, Blocks0),
sort(Blocks0, Blocks),
V is 5 * Q,
B is Q + 5 * (Q - 1) // 4 * Q,
write(c(V,5,2,2,B)), nl,
write(Blocks), nl.
% If Q is a prime power satisfying Q==1 (mod 4), then there exists a
% group-divisible design with Q groups of size 5, and blocks of size 5.
% See "On the Existence of (v,5,1)-BIBDs" by Clayton Smith.
% Note: Only Q=9 and prime Q are handled here.
% e.g. gdd(17, Groups, Blocks)
gdd(9, Groups, Blocks):-!,
findall(X, point(9, X), Points),
findall(G, group(9, Points, G), Groups),
findall(B, block(9, 4, 2, Points, B), Blocks0),
sort(Blocks0, Blocks).
gdd(Q, Groups, Blocks):-
1 =:= Q mod 4,
is_prime(Q),
primitive_root(Q, A),
D is (Q - 1) // 4,
findall(X, point(Q, X), Points),
findall(G, group(Q, Points, G), Groups),
findall(B, block(Q, A, D, Points, B), Blocks0),
sort(Blocks0, Blocks).
% On backtracking, finds all points (as 2-tuples)
point(Q, [Z,C]):-
Q_1 is Q - 1,
for(0, Q_1, C),
for(0, 4, Z).
% On backtracking, finds all groups (of integers)
group(Q, Points, G):-
Q_1 is Q - 1,
for(0, Q_1, C),
G0=[[0,C],[1,C],[2,C],[3,C],[4,C]],
indices(G0, Points, G).
% On backtracking, finds all blocks (of integers)
block(Q, A, D, Points, Block):-
Q=9, A=4, D=2, !, % Arithmetic in GF(3^2)
Add=[0,1,2,3,4,5,6,7,8,
1,2,0,4,5,3,7,8,6,
2,0,1,5,3,4,8,6,7,
3,4,5,6,7,8,0,1,2,
4,5,3,7,8,6,1,2,0,
5,3,4,8,6,7,2,0,1,
6,7,8,0,1,2,3,4,5,
7,8,6,1,2,0,4,5,3,
8,6,7,2,0,1,5,3,4],
Sub=[0,2,1,6,8,7,3,5,4,
1,0,2,7,6,8,4,3,5,
2,1,0,8,7,6,5,4,3,
3,5,4,0,2,1,6,8,7,
4,3,5,1,0,2,7,6,8,
5,4,3,2,1,0,8,7,6,
6,8,7,3,5,4,0,2,1,
7,6,8,4,3,5,1,0,2,
8,7,6,5,4,3,2,1,0],
Pow=[1,0,0,0,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,
1,2,1,2,1,2,1,2,1,
1,3,2,6,1,3,2,6,1,
1,4,6,7,2,8,3,5,1,
1,5,3,8,2,7,6,4,1,
1,6,2,3,1,6,2,3,1,
1,7,3,4,2,5,6,8,1,
1,8,6,5,2,4,3,7,1],
D_1 is D - 1,
Q_1 is Q - 1,
for(0, 4, I),
Ip1 is (I + 1) mod 5,
Im1 is (5 + I - 1) mod 5,
for(0, D_1, J),
T1 is A*9+J+1, nth_member(Pow, T1, AJ), % A^J
JD is (J + D) mod Q,
T2 is A*9+JD+1, nth_member(Pow, T2, AJD), % A^(J+D)
for(0, Q_1, C),
T3 is C*9+AJ+1, nth_member(Add, T3, K), % C+A^J
T4 is C*9+AJ+1, nth_member(Sub, T4, L), % C-A^J
T5 is C*9+AJD+1, nth_member(Add, T5, M), % C+A^(J+D)
T6 is C*9+AJD+1, nth_member(Sub, T6, N), % C-A^(J+D)
Block0=[[I,C],[Ip1,K],[Ip1,L],[Im1,M],[Im1,N]],
indices(Block0, Points, Block1),
sort(Block1, Block).
block(Q, A, D, Points, Block):- % Arithmetic modulo Q
D_1 is D - 1,
Q_1 is Q - 1,
for(0, 4, I),
Ip1 is (I + 1) mod 5,
Im1 is (5 + I - 1) mod 5,
for(0, D_1, J),
powermod(A, J, Q, AJ), % (A^J) mod Q
JD is J + D,
powermod(A, JD, Q, AJD), % (A^(J+D)) mod Q
for(0, Q_1, C),
K is (AJ + C) mod Q, % (C+A^J) mod Q
L is (Q + C - AJ) mod Q, % (C-A^J) mod Q
M is (AJD + C) mod Q, % (C+A^(J+D)) mod Q
N is (Q + C - AJD) mod Q, % (C-A^(J+D)) mod Q
Block0=[[I,C],[Ip1,K],[Ip1,L],[Im1,M],[Im1,N]],
indices(Block0, Points, Block1),
sort(Block1, Block).
% Returns the indices in Yss of each list in Xss
indices([], _, []).
indices([Xs|Xss], Yss, [Z|Zs]):-
nth_member(Yss, Z, Xs), !,
indices(Xss, Yss, Zs).
/* primitive_root(N, Root) is true if Root is a primitive root of Z_N. It */
/* is assumed, but not checked, that N is prime. */
/* e.g. primitive_root(239, 7) */
primitive_root(N, Root):-
N1 is N - 1,
findall(X, for(1,N1,X), Ms),
primitive_root_1(Ms, Ms, N, Root).
primitive_root_1([Root|_], Ms0, N, Root):-
is_primitive_root(Ms0, Root, N), !.
primitive_root_1([_|Ms], Ms0, N, Root):-
primitive_root_1(Ms, Ms0, N, Root).
is_primitive_root(Ms0, M, N):-is_primitive_root_1(Ms0, M, M, N).
is_primitive_root_1([], _, _, _):-!.
is_primitive_root_1(Ms, M, M0, N):-
remove(M, Ms, Ms1), !,
M1 is (M * M0) mod N,
is_primitive_root_1(Ms1, M1, M0, N).
% Develops the blocks of the cover corresponding to the base blocks of the
% difference family in Zv.
% e.g. develop(13, [[1,3,9],[2,6,5]], X) gives X=[[1,2,11],...,[9,12,13]]
develop(V, Sets, Blocks):-
develop_1(Sets, V, [], Blocks0),
sort(Blocks0, Blocks).
develop_1([], _, Blocks, Blocks).
develop_1([Set|Sets], V, Blocks0, Blocks):-
develop_2(V, V, Set, Blocks0, Blocks1),
develop_1(Sets, V, Blocks1, Blocks).
develop_2(0, _, _, Blocks, Blocks):-!.
develop_2(I, V, Set, Blocks0, [Block|Blocks]):-
I1 is I - 1,
develop_3(Set, I1, V, Block0),
sort(Block0, Block),
develop_2(I1, V, Set, Blocks0, Blocks).
develop_3([], _, _, []).
develop_3([X|Set], Y, V, [Point|Block]):-
Point is (X + Y) mod V + 1,
develop_3(Set, Y, V, Block).
/* 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).
/* 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).
/* is_prime(P) is true if P is a prime number. */
is_prime(2):-!.
is_prime(P):-
P > 2,
P mod 2 =:= 1,
SqrtP is sqrt(P),
is_prime_1(3, P, SqrtP).
is_prime_1(I, _, SqrtP):-
I > SqrtP, !.
is_prime_1(I, P, SqrtP):-
P mod I > 0,
I1 is I + 2,
is_prime_1(I1, P, SqrtP).
/* nth_member(+Xs, ?N, ?X) is true if X is the N-th (base 1) element of */
/* the list Xs. */
nth_member(Xs, N, X):-nth_member_1(Xs, X, 1, 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).
/* powermod(X, N, M, Y) is true if Y is the result of raising X to the */
/* power N, mod M. */
powermod(X, N, M, Y):-N >= 0, M >= 0, powermod_1(N, X, M, 1, Y).
powermod_1(0, _, _, Y, Y):-!.
powermod_1(N, X, M, Y0, Y):-
1 =:= N mod 2, !,
N1 is N // 2,
X1 is X * X mod M,
Y1 is Y0 * X mod M,
powermod_1(N1, X1, M, Y1, Y).
powermod_1(N, X, M, Y0, Y):-
N1 is N // 2,
X1 is X * X mod M,
powermod_1(N1, X1, M, Y0, Y).
/* remove(X, Ys, Zs) is true if Zs is the result of removing one */
/* occurrence of the element X from the list Ys. */
%remove(X, [X|Ys], Ys).
%remove(X, [Y|Ys], [Y|Zs]):-remove(X, Ys, Zs).