Construction of some (v,5,1)-BIBDs

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).

LPA Index     Home Page

Valid HTML 4.0!