Quelle  blockdesign.g   Sprache: unbekannt

############################################################################
##
##    blockdesign.g             Design 1.8.2 Package         Leonard Soicher
##
##
# Functions to study, construct and classify (sub)block designs 
# with given properties, and to partition block designs into 
# block designs with given properties.
# At present there is limited support for non-binary block designs.
# 
# Copyright (C) 2003-2024 Leonard H. Soicher
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#

BindGlobal("DESIGN_VERSION",InstalledPackageVersion("design"));
BindGlobal("DESIGN_INDETERMINATE_RATIONALS_1",Indeterminate(Rationals,1));

BindGlobal("IsSylowTowerGroupOfSomeComplexion",function(G)
#
# Suppose  G  is a finite group, and  p_1,...,p_k  are the distinct 
# (positive) primes dividing |G|, in some order. We say that 
# G is a *Sylow tower group of complexion* (p_1,...,p_k) if 
# G has a normal series 1=T_0<T_1<...<T_n=G such that, for 
# i=1,...,n, T_i/T_{i-1} is isomorphic to a Sylow p_i-subgroup
# of G.  
#
# This function returns  true  if G is a Sylow tower group of some complexion; 
# else  false  is returned. 
#
# See: L.H. Soicher, On classifying objects with specified groups of
# automorphisms, friendly subgroups, and Sylow tower groups, Port. Math. 74
# (2017), 233-242.
# 
local P,T,U,indices,usedindices,found,i;
if not IsGroup(G) or not IsFinite(G) then
   Error("<G> must be a finite group");
fi;
if IsNilpotentGroup(G) then
   return true;
elif not IsSolvableGroup(G) then 
   return false;
fi;
P:=List(Reversed(Set(FactorsInt(Size(G)))),p->SylowSubgroup(G,p));
T:=TrivialSubgroup(G); 
indices:=[1..Length(P)];
usedindices:=[];
while Length(usedindices)<Length(indices)-1 do
   #
   # Loop invariant: T is a normal subgroup of G,
   # and T is a Sylow tower group of complexion 
   # (p_{j_1},...,p_{j_m}), where  usedindices=[j_1,...,j_m],
   # and P[j] is a Sylow p_j-subgroup of G for j=1,...,Length(indices).
   #
   found:=false;
   for i in Difference(indices,usedindices) do
      U:=ClosureGroup(T,P[i]);
      if IsNormal(G,U) then
         found:=true;
         T:=U;
         Add(usedindices,i);
         break;
      fi;
   od;
   if not found then
      return false;
   fi;
od;
return true;
end);

BindGlobal("OnMultisetsRecursive",function(x,g)
#
# Action on multisets of multisets of ... of multisets of points.
# A multiset is given as a sorted list (of multisets or points).
#
local onmultisetsrecursive;

onmultisetsrecursive := function(x,g) 
local C;
if not IsList(x) then
   # x is a point
   return x^g;
fi;
C:=List(x,y->onmultisetsrecursive(y,g));
Sort(C);
return C;
end;

return onmultisetsrecursive(x,g);
end);

BindGlobal("HasLargerEntry",function(v,w)
#
# This boolean function assumes that v and w are equal-length integer vectors,
# and returns true iff v[i]>w[i] for some 1<=i<=Length(v).
#
local i;
for i in [1..Length(v)] do
   if v[i]>w[i] then
      return true;
   fi;
od;
return false;
end);

BindGlobal("BlockIntersectionPolynomial",function(x,m,lambdavec)
#
# Suppose  x  is a rational number or an indeterminate over the rationals, 
# and  m  and  lambdavec  are non-empty lists whose elements are rational 
# numbers  (and/or indeterminates over the rationals), such that 
# Length(lambdavec)<=Length(m). 
# 
# Then this function returns the block intersection polynomial 
# B(x) = B(x,m,lambdavec),  defined in:
# P.J. Cameron and L.H. Soicher, Block intersection polynomials, 
# Bull. London Math. Soc. 39 (2007), 559-564.
# 
local P,s,t,s1,s2,i,j;

P := function(x,k)
local prod,i;
prod:=1*x^0;
for i in [0..k-1] do
   prod:=prod*(x-i);
od;
return prod;
end;

if not ((IsRat(x) or IsPolynomialFunction(x)) and IsList(m) and IsList(lambdavec)) then
   Error("usage: BlockIntersectionPolynomial( <Rat> or <PolynomialFunction>, <List>, <List> )");
fi;
if m=[] or lambdavec=[] or Length(lambdavec)>Length(m) then
   Error("must have <m> and <lambdavec> non-empty, and Length(<lambdavec>)<=Length(<m>)");
fi;
if not ForAll(m,x->IsRat(x) or IsPolynomialFunction(x)) then
   Error("elements of <m> must be rationals or polynomials over the rationals");
fi;
if not ForAll(lambdavec,x->IsRat(x) or IsPolynomialFunction(x)) then
   Error("elements of <lambdavec> must be rationals or polynomials over the rationals");
fi;
s:=Length(m)-1;
t:=Length(lambdavec)-1;
s1:=0*x^0;
for j in [0..t] do
   s2:=0*x^0;
   for i in [j..s] do
      s2:=s2+P(i,j)*m[i+1];
   od;
   s1:=s1+Binomial(t,j)*P(-x,t-j)*(P(s,j)*lambdavec[j+1]-s2);
od;
return s1;
end);

BindGlobal("BlockIntersectionPolynomialCheck",function(m,lambdavec)
#
# Suppose  m  is a list of non-negative integers,
# and  lambdavec  is a list of non-negative rational numbers,
# with Length(lambdavec) odd, and Length(lambdavec)<=Length(m). 
# 
# Then this function checks whether the block intersection polynomial 
# B(x,m,lambdavec) is a polynomial over the integers and 
# has a non-negative value whenever x is an integer. 
# Returns true if all this is so; else false is returned.
#
local c,x,B,D,i,xmin,closestintxmin;
if not IsList(m) or not IsList(lambdavec) then
   Error("usage: BlockIntersectionPolynomialCheck( <List>, <List> )");
fi;
if Length(lambdavec) mod 2 = 0 or Length(lambdavec)>Length(m) then
   Error("Length(<lambdavec>) must be odd and at most Length(<m>)");
fi;
if not ForAll(m,x->IsInt(x) and x>=0) then
   Error("elements of <m> must be non-negative integers");
fi;
if not ForAll(lambdavec,x->IsRat(x) and x>=0) then
   Error("elements of <lambdavec> must be non-negative rationals");
fi;
x:=DESIGN_INDETERMINATE_RATIONALS_1;
B:=BlockIntersectionPolynomial(x,m,lambdavec);
c:=ShallowCopy(CoefficientsOfUnivariatePolynomial(B));
if Length(c)=0 then
   # B is the zero polynomial
   return true;
elif not ForAll(c,IsInt) then
   return false;
elif c[Length(c)]<0 or c[1]<0 or Length(c) mod 2 = 0 then
   # B has negative leading coef or B(0)<0 or B has odd degree.
   return false;
elif Length(c)=1 then
   # B is a constant positive polynomial.
   return true;
elif Length(c)=3 then
   # B is a quadratic with positive leading term.
   # Handle this case as suggested by Rhys Evans.
   xmin:=-c[2]/(2*c[3]); #  B has minimum value at xmin. 
   closestintxmin:=BestQuoInt(NumeratorRat(xmin),DenominatorRat(xmin));
   return Value(B,closestintxmin)>=0;
fi;
# Now apply the bound in [Soicher, MCompSci Thesis] to the absolute
# values of the zeros of B.  See also:  
# P.J. Cameron and L.H. Soicher, Block intersection polynomials, 
# Bull. London Math. Soc. 39 (2007), 559-564.
for i in [1..Length(c)-1] do
   c[i]:=-AbsoluteValue(c[i]);
od;
D:=UnivariatePolynomial(Rationals,c,1);
i:=1;
while Value(D,i)<=0 do
   if Value(B,i)<0 or Value(B,-i)<0 then
      return false;
   fi;
   i:=i+1;
od;
return true;
end);

BindGlobal("PartialDesignCheck",function(designinfo,block,blockbags)
#
# This boolean function returns  false  only if  blockbags  
# (and the points involved in blockbags) cannot be a subdesign 
# of a design with the given  designinfo  and  block.  
# Note: block must be in blockbags.  
#
local s,lambdavec,m,c,bag,intsize;
if not ForAny(blockbags,x->x.comb=block) then
   Error("<block> must be a block in a bag in <blockbags>");
fi;
s:=Length(block); # s>0
lambdavec:=ShallowCopy(designinfo.lambdavec);
if Length(lambdavec)<3 then
   if Length(lambdavec)=0 or not IsBound(designinfo.lambdamat) then
      Error("this should not happen");
   fi;
   if Length(lambdavec)=1 then
      lambdavec[2]:=Sum(List(block,i->designinfo.lambdamat[i][i]))/s;
   fi;
   # at this point, Length(lambdavec)=2
   if s>1 then
      lambdavec[3]:=0;
      for c in Combinations(block,2) do
         lambdavec[3]:=lambdavec[3]+designinfo.lambdamat[c[1]][c[2]];
      od;
      lambdavec[3]:=lambdavec[3]/Binomial(s,2);
   fi;
fi;
if Length(lambdavec)>s+1 then
   lambdavec:=lambdavec{[1..s+1]};
fi;
if Length(lambdavec) mod 2 = 0 then
   lambdavec:=lambdavec{[1..Length(lambdavec)-1]};
fi;
m:=ListWithIdenticalEntries(s+1,0);
for bag in blockbags do    
   intsize:=Size(Intersection(bag.comb,block));
   m[intsize+1]:=m[intsize+1]+bag.mult;
od;
if IsBound(designinfo.blockmmat) and IsBound(designinfo.blockmmat[s]) then
   # designinfo.blockmmat[s] is an integer vector of length s+1 and 
   # designinfo.blockmmat[s][i] is a lower bound on the number
   # of blocks intersecting a block of size s in exactly i-1 points.
   m:=List([1..s+1],i->Maximum(m[i],designinfo.blockmmat[s][i]));
fi;
return BlockIntersectionPolynomialCheck(m,lambdavec);
end);

BindGlobal("CC_PermutationGroupProperties",function(G,n)
#
# Determines the properties rank, isGenerouslyTransitive, isMultiplicityFree,
# and isStratifiable of the permutation group  G  acting on [1..n].
#
local prop,A,P,PP,i,j;
if not IsPermGroup(G) or not IsInt(n) then
   Error("usage: CC_PermutationGroupProperties( <PermGroup>, <Int> )");
fi;   
prop:=rec();
if not IsTransitive(G,[1..n]) then
   prop.rank:=Length(Orbits(G,Cartesian([1..n],[1..n]),OnTuples));
   prop.isGenerouslyTransitive:=false;
   prop.isMultiplicityFree:=false;
   prop.isStratifiable:=false;
   return prop;
fi;
A:=OrbitalGraphColadjMats(G);
prop.rank:=Length(A);
prop.isGenerouslyTransitive:=ForAll([1..Length(A)],i->A[i][i][1]=1);
if prop.isGenerouslyTransitive then
   prop.isMultiplicityFree:=true;
   prop.isStratifiable:=true;
   return prop;
fi;
P:=List(A,TransposedMat); 
# P  is a list of the intersection matrices in the order determined by  A.
prop.isMultiplicityFree:=ForAll(Combinations(P,2),x->x[1]*x[2]=x[2]*x[1]);
if prop.isMultiplicityFree then
   prop.isStratifiable:=true;
   return prop;
fi;
PP:=[];
for i in [1..Length(A)] do
   j:=First([i..Length(A)],j->A[i][j][1]=1); 
   # The orbital corresponding to A[i] is paired with that corresp. to A[j].
   if j<>fail then
      if j=i then 
         # A[i] correpsonds to a self-paired orbital.
         Add(PP,P[i]);
      else
         Add(PP,P[i]+P[j]);
      fi; 
   fi;   
od;   
prop.isStratifiable:=ForAll(Combinations(PP,2),x->x[1]*x[2]=x[2]*x[1]);
return prop;
end);

BindGlobal("IsBlockDesign",function(D)
if not (IsRecord(D) and IsBound(D.isBlockDesign) and D.isBlockDesign=true) then 
   return false;
fi;
# Now perform some easy checks.
# More complete checking is done by the function  BlockDesign.
if not (IsBound(D.v) and IsInt(D.v) and D.v>0) then
   Error("<D>.v must be a positive integer");
fi;
if not (IsBound(D.blocks) and IsList(D.blocks) and D.blocks<>[]) then
   Error("<D>.blocks must be a non-empty (sorted) list");
fi;
if not ForAll(D.blocks,x->Length(x)>0) then
   # We do not allow empty blocks.
   Error("Each block must be a non-empty (sorted) list");
fi;
return true;
end);

BindGlobal("BlockDesign",function(arg)
#
# Let  v=arg[1],  blocks=arg[2],  and  G=arg[3]  (default: Group(())). 
# Then this function returns the block design with point-set {1,...,v},
# and whose block-list is the sorted concatenation of the G-orbits of the
# the elements of  blocks. 
#
# Note: It is no longer required that  blocks  be sorted.
#
local v,block,blocks,G,newblocks;
v:=arg[1];
blocks:=arg[2];
if IsBound(arg[3]) then 
   G:=arg[3];
else
   G:=Group(());
fi;
if not (IsInt(v) and IsList(blocks) and IsPermGroup(G)) then
   Error("usage: BlockDesign( <Int>, <List> [, <PermGroup> ] )");
fi;
if v<1 then
   Error("<v> must be positive");
fi;
if v<LargestMovedPoint(GeneratorsOfGroup(G)) then 
   Error("<G> must be a permutation group on [1..<v>]");
fi;
if blocks=[] then
   Error("<blocks> must be non-empty");
fi;
if not ForAll(blocks,x->x<>[] and IsSortedList(x)) then
   Error("each element of <blocks> must be a non-empty sorted list");
fi;
if not IsSubset([1..v],Set(Flat(blocks))) then
   Error("not all elements of <blocks> are in [1..<v>]");
fi;
if IsTrivial(G) then
   return rec(isBlockDesign:=true, v:=v, blocks:=AsSortedList(blocks));
fi;
# Now handle the case where G is nontrivial.
newblocks:=[];
for block in blocks do 
   if IsSet(block) then
      Append(newblocks,Orbit(G,block,OnSets));
   else 
      Append(newblocks,Orbit(G,block,OnMultisetsRecursive));
   fi;
od;
return rec(isBlockDesign:=true, v:=v, blocks:=AsSortedList(newblocks),
   autSubgroup:=G);
end);

BindGlobal("BlockDesignPoints",function(D)
if not IsBlockDesign(D) then
   Error("usage: BlockDesignPoints( <BlockDesign> )");
fi;
return Immutable([1..D.v]);
end); 

BindGlobal("BlockDesignBlocks",function(D)
if not IsBlockDesign(D) then
   Error("usage: BlockDesignBlocks( <BlockDesign> )");
fi;
return Immutable(D.blocks);
end); 

BindGlobal("NrBlockDesignPoints",function(D)
if not IsBlockDesign(D) then
   Error("usage: NrBlockDesignPoints( <BlockDesign> )");
fi;
return D.v;
end); 

BindGlobal("NrBlockDesignBlocks",function(D)
if not IsBlockDesign(D) then
   Error("usage: NrBlockDesignBlocks( <BlockDesign> )");
fi;
return Length(D.blocks);
end); 

BindGlobal("IsBinaryBlockDesign",function(D)
if not IsBlockDesign(D) then
   Error("usage: IsBinaryBlockDesign( <BlockDesign> )");
fi;
if not IsBound(D.isBinary) then
   D.isBinary:=ForAll(D.blocks,IsSet);
fi;
return D.isBinary;
end); 

BindGlobal("IsSimpleBlockDesign",function(D)
if not IsBlockDesign(D) then
   Error("usage: IsSimpleBlockDesign( <BlockDesign> )");
fi;
if not IsBound(D.isSimple) then
   D.isSimple:=IsSet(D.blocks);
fi;
return D.isSimple;
end); 

BindGlobal("BlockSizes",function(D)
if not IsBlockDesign(D) then
   Error("usage: BlockSizes( <BlockDesign> )");
fi;
if not IsBound(D.blockSizes) then
   D.blockSizes:=AsSet(List(D.blocks,Length));
fi;
return Immutable(D.blockSizes);
end); 

BindGlobal("BlockNumbers",function(D)
if not IsBlockDesign(D) then
   Error("usage: BlockNumbers( <BlockDesign> )");
fi;
if not IsBound(D.blockNumbers) then
   D.blockNumbers:=
      Immutable(List(BlockSizes(D),x->Number(D.blocks,y->Length(y)=x)));
fi;
return Immutable(D.blockNumbers);
end); 

BindGlobal("TSubsetLambdasVector",function(D,t)
local c,i,j,ii,xx,tsubsets,tvector,block;
if not IsBlockDesign(D) or not IsInt(t) then
   Error("usage: TSubsetLambdasVector( <BlockDesign>, <Int> )");
fi;
if t<0 then
   Error("<t> must be non-negative");
elif t>D.v then
   return [];
fi;
tsubsets:=Combinations([1..D.v],t);
tvector:=ListWithIdenticalEntries(Length(tsubsets),0); # initialization
for block in D.blocks do
   if not IsSet(block) then
      xx:=ListWithIdenticalEntries(D.v,0);
      for i in block do 
         xx[i]:=xx[i]+1;
      od;
   fi;
   for c in Combinations(Set(block),t) do
      if IsSet(block) then
         ii:=1;
      else
         ii:=Product(xx{c});
      fi;
      j:=PositionSorted(tsubsets,c);
      tvector[j]:=tvector[j]+ii; 
   od;
od;
return tvector;
end);

BindGlobal("TSubsetStructure",function(D,t)
local tsubsets,tvector,lambdas,partition,tsubsetstructure,i;
if not IsBlockDesign(D) or not IsInt(t) then
   Error("usage: TSubsetStructure( <BlockDesign>, <Int> )");
fi;
if t<0 then
   Error("<t> must be non-negative");
fi;
if IsBound(D.tSubsetStructure) and t=D.tSubsetStructure.t then 
   return Immutable(D.tSubsetStructure);
fi;
if IsBound(D.allTDesignLambdas) and Length(D.allTDesignLambdas)>=t+1 then
   return Immutable(rec(t:=t,lambdas:=[D.allTDesignLambdas[t+1]]));
fi;
tsubsets:=Combinations([1..D.v],t);
tvector:=TSubsetLambdasVector(D,t);
lambdas:=Set(tvector);
tsubsetstructure:=rec(t:=t,lambdas:=lambdas);
if Length(lambdas)>1 then
   partition:=List(lambdas,x->[]);  # initialization
   for i in [1..Length(tvector)] do
      Add(partition[PositionSorted(lambdas,tvector[i])],tsubsets[i]); 
   od;
   tsubsetstructure.partition:=partition;
   return Immutable(tsubsetstructure);
elif not IsBound(D.tSubsetStructure) and
   t>1 and Length(lambdas)=1 and lambdas[1]>0 then
   # this is probably info worth recording 
   D.tSubsetStructure:=Immutable(tsubsetstructure);
   return D.tSubsetStructure;
else
   return Immutable(tsubsetstructure);
fi;
end);

BindGlobal("AllTDesignLambdas",function(D)
local b,k,m,alltdesignlambdas;

alltdesignlambdas:=function(D)
# This function does all the real work.
# It is assumed that D is a binary block design with 
# constant block size.
local all,t,v,k,lambda,b,T,G,DD;
v:=D.v;
k:=Length(D.blocks[1]);
b:=Length(D.blocks);
if IsBound(D.autGroup) then
   G:=D.autGroup;
elif IsBound(D.autSubgroup) then
   G:=D.autSubgroup;
else
   G:=Group(());
fi;
if IsTransitive(G,[1..v]) then
   if k=1 then 
      D.allTDesignLambdas:=Immutable([b,b/v]); 
      return D.allTDesignLambdas;
   fi;
   DD:=BlockDesign(v-1,List(Filtered(D.blocks,x->x[k]=v),y->y{[1..k-1]}));
   # DD is the derived design of D wrt the point v.
   DD.autSubgroup:=Stabilizer(G,v);
   D.allTDesignLambdas:=Immutable(Concatenation([b],alltdesignlambdas(DD)));
   return D.allTDesignLambdas;
fi;
# Now handle the case where G is not point-transitive.
all:=[b];
for t in [1..k] do
   lambda:=b*Binomial(k,t)/Binomial(v,t);
   if not IsInt(lambda) then
      D.allTDesignLambdas:=Immutable(all);
      return D.allTDesignLambdas;
   fi;
   T:=TSubsetLambdasVector(D,t);
   if not ForAll(T,x->x=lambda) then
      D.allTDesignLambdas:=Immutable(all);
      return D.allTDesignLambdas;
   fi;      
   # If we get here then we know that D is a t-(v,k,lambda) design.
   Add(all,lambda);
od;    
D.allTDesignLambdas:=Immutable(all);
return D.allTDesignLambdas;
end;

if not IsBlockDesign(D) then
   Error("usage: AllTDesignLambdas( <BlockDesign> )");
fi;
if IsBound(D.allTDesignLambdas) then
   return Immutable(D.allTDesignLambdas);
fi;
if Length(BlockSizes(D))>1 or not IsBinaryBlockDesign(D) then
   # D is not a t-design for any t.
   D.allTDesignLambdas:=Immutable([]); 
   return D.allTDesignLambdas;
fi;
b:=Length(D.blocks);
k:=BlockSizes(D)[1];
m:=b/Binomial(D.v,k);
if IsInt(m) and ForAll(Collected(D.blocks),x->x[2]=m) then
   # D is m times the trivial design.
   D.allTDesignLambdas:=Immutable(List([0..k],
      i->b*Binomial(k,i)/Binomial(D.v,i))); 
   return D.allTDesignLambdas;
else 
   return alltdesignlambdas(D);
fi;
end); 
   
BindGlobal("PossibleTDesignLambdasFromParameters",function(t,v,k,lambda)
local b;
if not (IsInt(t) and IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: PossibleTDesignLambdasFromParameters( <Int>, <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 or lambda<=0 then
   Error("must have <v>,<k>,<lambda> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k,lambda>0.
#
b:=lambda*Binomial(v,t)/Binomial(k,t);
return List([0..k],i->b*Binomial(k,i)/Binomial(v,i)); 
end);

BindGlobal("IsPossiblySBIBD",function(v,k,lambda)
#
# This boolean function returns false only if there is no
# (v,k,lambda)-SBIBD.
#
local n,p,p1,p2,num;
if not (IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: IsPossiblySBIBD( <Int>, <Int>, <Int> )");
fi;
if not (2<=k and k<v and lambda>0) then
   Error("must have 2 <= <k> < <v>, and  <lambda> > 0"); 
fi;
if lambda*(v-1)/(k-1) <> k then
   return false;
fi;
# Parameters are admissible, and  r=k, b=v.
# Now check the Bruck-Ryser-Chowla condition.
n:=k-lambda;
if v mod 2 = 0 then
   # return true iff  n  is the square of an integer.
   return n=RootInt(n,2)^2; 
fi;
# For v odd, we test BRC by applying Theorem 2.5 in 
# E.S. Lander, Symmetric Designs: An Algebraic Approach, CUP, 1983. 
p1:=List(Filtered(Collected(FactorsInt(lambda)),x->x[2] mod 2 <> 0),y->y[1]);
# now p1 is a list of the distinct primes dividing the 
# squarefree part of  lambda.
p2:=List(Filtered(Collected(FactorsInt(n)),x->x[2] mod 2 <> 0),y->y[1]);
# now p2 is a list of the distinct primes dividing the 
# squarefree part of  n.
for p in Difference(p1,Union([2],p2)) do 
   if not Legendre(n,p)=1 then
      # BRC condition does not hold
      return false;
   fi;
od;
num:=(-1)^((v-1)/2)*Product(p1);
for p in Difference(p2,Union([2],p1)) do 
   if not Legendre(num,p)=1 then
      # BRC condition does not hold
      return false;
   fi;
od;
num:=(-1)^((v+1)/2)*Product(p1)*Product(p2);
for p in Difference(Intersection(p1,p2),[2]) do 
   if not Legendre(num/p^2,p)=1 then
      # BRC condition does not hold
      return false;
   fi;
od;
# BRC condition holds.
return true;
end);

BindGlobal("IsPossiblyBIBD",function(v,k,lambda)
#
# This boolean function returns false only if there is no
# (v,k,lambda)-BIBD.
#
local b,r;
if not (IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: IsPossiblyBIBD( <Int>, <Int>, <Int> )");
fi;
if not (2<=k and k<v and lambda>0) then
   Error("must have 2 <= <k> < <v>, and  <lambda> > 0"); 
fi;
if lambda*(v-1) mod (k-1) <> 0 then
   return false;
fi;
r:=lambda*(v-1)/(k-1);
if v*r mod k <> 0 then
   return false;
fi;
b:=v*r/k;
if b<v then 
   # Fisher's inequality does not hold.
   return false;
fi; 
if b=v then
   return IsPossiblySBIBD(v,k,lambda);
fi;
if r=k+lambda and lambda<=2 then
   # We are looking for embeddable quasiresidual BIBDs.
   # If lambda=1 then we are looking for affine planes 
   # (embeddable in projective planes), and if lambda=2 then
   # embeddablity is by the Hall-Connor theorem. 
   if not IsPossiblySBIBD(v+r,r,lambda) then
      # no SBIBDs to embed the required designs.
      return false;
   fi;
fi;
return true;
end);

BindGlobal("TDesignLambdas",function(t,v,k,lambda)
local lambdas;
if not (IsInt(t) and IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: TDesignLambdas( <Int>, <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 or lambda<=0 then
   Error("must have <v>,<k>,<lambda> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k,lambda>0.
#
lambdas:=List([0..t],s->lambda*Binomial(v-s,t-s)/Binomial(k-s,t-s));
if not ForAll(lambdas,IsInt) then
   # parameters t-(v,k,lambda) are inadmissible
   return fail;
fi;
return lambdas;
end);

BindGlobal("TDesignIntersectionTriangle",function(t,v,k,lambda)
local lambdas,T,i,j;
if not (IsInt(t) and IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: TDesignIntersectionTriangle( <Int>, <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 or lambda<=0 then
   Error("must have <v>,<k>,<lambda> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k,lambda>0.
#
lambdas:=TDesignLambdas(t,v,k,lambda);
if lambdas=fail then
   return fail;
fi;
# Make first column of T.
T:=List(lambdas,x->[x]);
# Make remaining columns of T.
for j in [1..t] do
   for i in [0..t-j] do
      T[i+1][j+1]:=T[i+1][j]-T[i+2][j];
   od;
od;
return T;
end);

BindGlobal("SteinerSystemIntersectionTriangle",function(t,v,k)
local lambdas,T,i,j;
if not (IsInt(t) and IsInt(v) and IsInt(k)) then
   Error("usage: SteinerSystemIntersectionTriangle( <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 then
   Error("must have <v>,<k> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k>0.
#
lambdas:=TDesignLambdas(t,v,k,1);
if lambdas=fail then
   return fail;
fi;
# Make first column of T.
T:=List(lambdas,x->[x]);
for i in [t+1..k] do
   T[i+1]:=[1];
od;
# Make remaining columns of T.
for j in [1..k] do
   for i in [0..k-j] do
      T[i+1][j+1]:=T[i+1][j]-T[i+2][j];
   od;
od;
return T;
end);

BindGlobal("TDesignLambdaMin",function(t,v,k)
#
# Returns the minimum lambda such that TDesignLambdas(t,v,k,lambda)<>fail.
#
local s,numer,denom,lambdamin;
if not (IsInt(t) and IsInt(v) and IsInt(k)) then
   Error("usage: TDesignLambdaMin( <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 then
   Error("must have <v>,<k> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k>0.
#
lambdamin:=1;
for s in [0..t] do
   numer:=lambdamin*Binomial(v-s,t-s);
   denom:=Binomial(k-s,t-s);
   lambdamin:=lambdamin*denom/Gcd(denom,numer);
od;
return lambdamin;
end);

BindGlobal("TDesignBlockMultiplicityBound",function(t,v,k,lambda)

local blockmultiplicitybound;

blockmultiplicitybound:=function(t,v,k,lambda)
#
# This function does most of the work.
# It is assumed that t,v,k,lambda are integers, with
# v,k,lambda>0 and 0<=t<=k<=v, but this is not checked.
#
local newlambda,lambdavec,multbound,m,lower,upper,mid;
if t=0 or k=v then
   return lambda;
elif k>v/2 and v-k>=t then
   # consider the design with the blocks complemented
   newlambda:=lambda*Binomial(v-k,t)/Binomial(k,t);
   if not IsInt(newlambda) then
      # no design
      return 0;
   else
      return blockmultiplicitybound(t,v,v-k,newlambda);
   fi;
fi;
lambdavec:=TDesignLambdas(t,v,k,lambda);
if lambdavec=fail then
   # parameters t-(v,k,lambda) inadmissible
   return 0;
elif t=1 then 
   return lambda;
fi;
# Now 2<=t<=k<v.
if t=2 then
   if not IsPossiblyBIBD(v,k,lambda) then
      return 0;
   fi;
   multbound:=lambda;
else
   # consider the derived design
   multbound:=blockmultiplicitybound(t-1,v-1,k-1,lambda);
fi;
if t mod 2 <> 0 or multbound=0 then
   return multbound;
fi;
# Now t>=2 is even.
# Refine bound using block intersection polynomials and binary search.
m:=ListWithIdenticalEntries(k+1,0);
lower:=0;
upper:=multbound+1;
# The bound on the multiplicity of a block is  >= lower  and  < upper.
while upper-lower>1 do
   mid:=Int((lower+upper)/2); 
   m[k+1]:=mid;
   if BlockIntersectionPolynomialCheck(m,lambdavec) then
      lower:=mid;
   else
      upper:=mid;
   fi;
od;
return lower;
end;

if not (IsInt(t) and IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: TDesignBlockMultiplicityBound( <Int>, <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 or lambda<=0 then
   Error("must have <v>,<k>,<lambda> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
return blockmultiplicitybound(t,v,k,lambda);
end);

BindGlobal("OARunMultiplicityBound",function(N,k,s,t)
#
# Suppose N is a positive integer, and k and s are non-empty
# lists of the same length, w say, of positive integers, with 
# each element of s >1, and t is a non-negative integer <= Sum(k). 
#
# Then this function returns an upper bound on the multiplicity of 
# any run in a (mixed) orthogonal array OA(N,s[1]^k[1],...,s[w]^k[w],t). 
#
# If k is given as an integer, it is replaced by [k].
# If s is given as an integer, it is replaced by [s].
#
# The technique used is based on the results in: 
# P.J. Cameron and L.H. Soicher, Block intersection polynomials, 
# Bull. London Math. Soc. 39 (2007), 559-564.
#

local runmultiplicitybound,S,i,j;

runmultiplicitybound:=function(N,S,t)
#
# This function does most of the work.
# S is a nonempty list whose length is the number of factors, and 
# with S[i] being the size of the symbol set for the i-th factor.
#
local lambdavec,m,i,C,c,count,nrtuples,nrfactors,num,
   symbolsetsizes,maxsymbolsetsize;
if t=0 then
   return N;
fi; 
nrfactors:=Length(S);
symbolsetsizes:=Set(S); # the set of the sizes of the symbol sets
maxsymbolsetsize:=symbolsetsizes[Length(symbolsetsizes)];
if t=1 then
   if ForAny(symbolsetsizes,size->N mod size <> 0) then
      return 0;
   else
      return N/maxsymbolsetsize;
   fi;
fi;
# Now Length(S) >= t >= 2. 
if t mod 2 <> 0 then
   return Minimum(List(symbolsetsizes,size->
      runmultiplicitybound(N/size,S{Difference([1..nrfactors],[Position(S,size)])},t-1)));
fi;
# Now t>=2 is even, and we determine the bound using
# block intersection polynomials.
num:=ListWithIdenticalEntries(maxsymbolsetsize,0);
for i in [1..nrfactors] do
   num[S[i]]:=num[S[i]]+1;
od;
# So now, num[j] is the number of factors having a symbol set of size j.
lambdavec:=[N];
for i in [1..t] do
   #
   # determine lambdavec[i+1]
   #
   lambdavec[i+1]:=0;
   C:=Combinations(S,i); # the set of submultisets of S of size i 
   for c in C do
      nrtuples:=Product(c);
      # nrtuples is the number of i-tuples where the first element
      # is from a symbol set of size c[1], the second from a symbol
      # set of size c[2], and so on.
      if N mod nrtuples <> 0 then
         return 0;
      fi;
      count:=Product(Collected(c),x->Binomial(num[x[1]],x[2]));
      # count is the number of times the multiset  c  occurs over 
      # all choices of  i  elements from  S.
      #
      # Now add in their contribution.
      # 
      lambdavec[i+1]:=lambdavec[i+1]+count*(N/nrtuples);
   od;
   # Now take the average contribution.
   lambdavec[i+1]:=lambdavec[i+1]/Binomial(nrfactors,i); 
od;
m:=ListWithIdenticalEntries(nrfactors+1,0);
m[nrfactors+1]:=1;
while BlockIntersectionPolynomialCheck(m,lambdavec) do
   m[nrfactors+1]:=m[nrfactors+1]+1;
od;
return m[nrfactors+1]-1;
end;

if IsInt(k) then
   k:=[k];
fi;
if IsInt(s) then
   s:=[s];
fi;
if not (IsInt(N) and IsList(k) and IsList(s) and IsInt(t)) then
   Error("usage: OARunMultiplicityBound( <Int>, <Int> or <List>, <Int> or <List>, <Int> )");
fi;
if N<1 or t<0 then
   Error("<N> must be positive and <t> must be non-negative");
fi;
if Length(s)<>Length(k) or Length(s)=0 then
   Error("<s> and <k> must have equal positive length");
fi;
if ForAny(k,x->not IsPosInt(x)) or t>Sum(k) then
   Error("<k> must be a list of positive integers, with Sum(<k>) >= <t>");
fi;
if ForAny(s,x->not IsInt(x) or x<2) then
   Error("each element of <s> must be an integer > 1");
fi;
S:=[];
for i in [1..Length(s)] do
   for j in [1..k[i]] do
      Add(S,s[i]);
   od;
od;
return runmultiplicitybound(N,S,t);
end);

BindGlobal("ResolvableTDesignBlockMultiplicityBound",function(t,v,k,lambda)

local multbound,r,lambdavec,m,lower,upper,mid;
if not (IsInt(t) and IsInt(v) and IsInt(k) and IsInt(lambda)) then
   Error("usage: ResolvableTDesignBlockMultiplicityBound( <Int>, <Int>, <Int>, <Int> )");
fi;
if v<=0 or k<=0 or lambda<=0 then
   Error("must have <v>,<k>,<lambda> > 0");
fi;
if 0>t or t>k or k>v then
   Error("must have 0 <= <t> <= <k> <= <v>");
fi;
#
# Now 0<=t<=k<=v and v,k,lambda>0.
#
if (v mod k <> 0) then
   # resolvable t-(v,k,lambda) design cannot exist
   return 0;
elif k=v then
   # complete blocks
   return lambda;
fi;
lambdavec:=TDesignLambdas(t,v,k,lambda);
if lambdavec=fail then
   # parameters t-(v,k,lambda) inadmissible
   return 0;
fi;
if t=0 then 
   r:=lambdavec[1]*k/v;
   if not IsInt(r) then
      # design cannot be resolvable
      return 0;
   else
      return r;
   fi;
elif t=1 then
   return lambda;
fi;
# Now 2<=t<=k<v. 
if lambdavec[1]<v+lambdavec[2]-1 then
   # Bose's inequality does not hold
   return 0;
fi;
multbound:=TDesignBlockMultiplicityBound(t,v,k,lambda); # initialize
if t mod 2 <> 0 or multbound=0 then
   return multbound;
fi;
# Now t>=2 is even.
# Refine bound using block intersection polynomials and binary search.
m:=ListWithIdenticalEntries(k+1,0);
lower:=0;
upper:=multbound+1;
# The bound on the multiplicity of a block is  >= lower  and  < upper.
while upper-lower>1 do
   mid:=Int((lower+upper)/2); 
   m[k+1]:=mid;
   m[1]:=mid*(v/k-1);
   if BlockIntersectionPolynomialCheck(m,lambdavec) then
      lower:=mid;
   else
      upper:=mid;
   fi;
od;
return lower;
end);

BindGlobal("IncidenceGraphSupportBlockDesign",function(D)
local G,supportblocks,gamma,pointnames;
if not IsBlockDesign(D) then
   Error("usage: IncidenceGraphSupportBlockDesign( <BlockDesign> )");
fi;
if IsBound(D.autGroup) then
   G:=D.autGroup;
elif IsBound(D.autSubgroup) then
   G:=D.autSubgroup;
else
    G:=Group(());
fi;
if IsSimpleBlockDesign(D) then
   supportblocks:=D.blocks;
else
   supportblocks:=Set(D.blocks);
fi;
gamma:=Graph(G,Concatenation([1..D.v],supportblocks),OnMultisetsRecursive,
             function(x,y)
             if IsInt(x) then
                return IsList(y) and x in y;
             else
                return IsInt(y) and y in x;
             fi;
             end,true);
if IsBound(D.pointNames) then
   pointnames:=D.pointNames;
else
   pointnames:=[1..D.v];
fi;
gamma.names:=Immutable(Concatenation(pointnames,
                          List(supportblocks,x->pointnames{x})));
return gamma;
end);
   
BindGlobal("IsConnectedBlockDesign",function(D)
if not IsBlockDesign(D) then
   Error("usage: IsConnectedBlockDesign( <BlockDesign> )");
fi;
if not IsBound(D.isConnected) then 
   D.isConnected:=IsConnectedGraph(IncidenceGraphSupportBlockDesign(D));
fi;
return D.isConnected;
end); 
   
BindGlobal("ReplicationNumber",function(D)
#
# If the block design  D  is equireplicate, then this function
# returns its replication number  r;  otherwise  fail  is returned.
#
local count,blocks,block,point;
if not IsBlockDesign(D) then
   Error("usage: ReplicationNumber( <BlockDesign> )");
fi;
if IsBound(D.r) then
   return D.r;
fi;   
count:=ListWithIdenticalEntries(D.v,0); 
for block in D.blocks do
   for point in block do
      count[point]:=count[point]+1;
   od;
od;
if ForAll(count,x->x=count[1]) then 
   # D is equireplicate
   D.r:=count[1];
   return D.r;
else
   # D is not equireplicate 
   return fail;
fi;   
end);

BindGlobal("PairwiseBalancedLambda",function(D)
local lambdas;
if not IsBinaryBlockDesign(D) then
   Error("usage: PairwiseBalancedLambda( <BinaryBlockDesign> )");
fi;
if D.v<2 then
   return fail;
fi;
lambdas:=TSubsetStructure(D,2).lambdas;
if Length(lambdas)=1 and lambdas[1]>0 then 
   # D is pairwise balanced
   return lambdas[1]; 
else
   # D is not pairwise balanced
   return fail;
fi;
end);

BindGlobal("IsVarianceBalanced",function(D)
local vec,j,c,twosubsets,block;
if not IsBinaryBlockDesign(D) then
   Error("<D> must be a binary block design");
fi;
if D.v=1 or not IsConnectedBlockDesign(D) then
   # this function is not applicable
   return fail;
fi;   
twosubsets:=Combinations([1..D.v],2);
vec:=ListWithIdenticalEntries(Length(twosubsets),0); # initialization
for block in D.blocks do 
   for c in Combinations(block,2) do
      j:=PositionSorted(twosubsets,c);
      vec[j]:=vec[j]+1/Length(block);
   od;
od;
return ForAll(vec,x->x=vec[1]);
end);

BindGlobal("AffineResolvableMu",function(D)
#
# A block design is *affine resolvable* if the design is resolvable 
# and any two blocks not in the same parallel class of a resolution 
# meet in a constant number mu of points. 
#
# If the block design  <D>  is affine resolvable, then this function
# returns its  mu;  otherwise  `fail'  is returned.
# 
# The value 0 is returned iff  <D>  consists of a single parallel class.
#
local blocks_remaining,mu,A,B,AA,m;
if not IsBlockDesign(D) then
   Error("usage: AffineResolvableMu( <BlockDesign> )");
fi;
if not IsBinaryBlockDesign(D) then
   #  D  not resolvable
   return fail;
fi;
if not IsSimpleBlockDesign(D) then
   if ForAll(D.blocks,x->Length(x)=D.v) then
      return D.v;
   else
      return fail;
   fi;
fi;
blocks_remaining:=ShallowCopy(D.blocks);
mu:=0;
while blocks_remaining <> [] do
   A:=blocks_remaining[1];
   RemoveSet(blocks_remaining,A);
   AA:=[]; # build in AA the elements other than A in the parallel class of A
   for B in blocks_remaining do
      m:=Length(Intersection(A,B));
      if m=0 then
         if ForAny(AA,x->Intersection(B,x)<>[]) then
            return fail;
         else
            Add(AA,B);
         fi;
      elif mu=0 then
         mu:=m;
      elif m<>mu then
         # non-constant mu
         return fail;
      fi; 
   od;
   if Length(A)+Sum(List(AA,Length))<>D.v then
      # [A] union AA is not a parallel class
      return fail;
   else
      SubtractSet(blocks_remaining,AA);
      # at this point mu>0 or blocks_remaining=[] 
      if ForAny(AA,x->ForAny(blocks_remaining,
                         y->Length(Intersection(x,y))<>mu)) then
         return fail;
      fi;
   fi;
od;
return mu;
end);

BindGlobal("AutGroupBlockDesign",function(D)
local collectedblocks,blockmultiplicities,vertexpartition,i,gamma;
if not IsBlockDesign(D) then
   Error("usage: AutGroupBlockDesign( <BlockDesign> )");
fi;
if IsBound(D.autGroup) then
   return D.autGroup;
fi;
if not IsBinaryBlockDesign(D) then
   Error("aut group not yet implemented for non-binary block designs");
fi;
gamma:=IncidenceGraphSupportBlockDesign(D);
collectedblocks:=Collected(D.blocks);
blockmultiplicities:=Set(List(collectedblocks,x->x[2]));
vertexpartition:=List(blockmultiplicities,x->[]); # initialization
for i in [1..Length(collectedblocks)] do
   Add(vertexpartition[PositionSorted(blockmultiplicities,collectedblocks[i][2])],D.v+i);
od;
vertexpartition:=Concatenation([[1..D.v]],vertexpartition);
D.autGroup:=Action(AutGroupGraph(gamma,vertexpartition),[1..D.v]);
return D.autGroup;
end); 
   
InstallOtherMethod(AutomorphismGroup,"for block design",[IsRecord],10,
function(D)
  if not IsBlockDesign(D) then
     TryNextMethod();
  fi;
  return AutGroupBlockDesign(D);
end);

BindGlobal("BlockDesignIsomorphismClassRepresentatives",function(L)
#
# Given a list  L  of binary block designs, this function returns
# a list containing pairwise non-isomorphic elements of  L,  representing
# all the isomorphism classes of elements of  L. 
#
local graphs,i,designs,reps;
if not IsList(L) then
   Error("usage: BlockDesignIsomorphismClassRepresentatives( <List> )");
fi;
if not ForAll(L,x->IsBlockDesign(x) and IsBinaryBlockDesign(x)) then
   Error("each element of <L> must be a binary block design");
fi;
designs:=ShallowCopy(L);
if Length(designs)<=1 then
   return designs;
fi;
graphs:=List(designs,A->Graph(Group(()),[1..A.v+Length(A.blocks)],
                 function(x,g) return x; end,
                 function (x,y) 
                 if x<=A.v then
                    if y<=A.v then 
                       return x=y;  # put loops on point-vertices
                    fi;
                    return x in A.blocks[y-A.v];
                 else 
                    return y<=A.v and y in A.blocks[x-A.v];
                 fi;
                 end,true));
for i in [1..Length(graphs)] do
   SetAutGroupCanonicalLabelling(graphs[i]);
   if not IsBound(designs[i].autGroup) then
      designs[i].autGroup:=Action(AutGroupGraph(graphs[i]),[1..designs[i].v]);
   fi;
   graphs[i]:=DirectedEdges(
      GraphImage(graphs[i],graphs[i].canonicalLabelling^(-1)));
   # Now  graphs[i]  is the set of directed edges of the canonical 
   # representative of the incidence graph of  designs[i]
   # (these directed edges determine this canonical graph since 
   # we assume no block design has an empty block).
od;
SortParallel(graphs,designs);
reps:=[designs[1]];
for i in [2..Length(graphs)] do
   if graphs[i]<>graphs[i-1] then
      # new isomorphism class representative
      Add(reps,designs[i]);
   fi;
od;
return reps;
end); 
   
BindGlobal("IsEqualBlockDesign",function(A,B)
#
# This boolean function returns true iff block designs  A  and  B 
# are equal (i.e. A.v=B.v and A and B have equal block-multisets).  
#
if not IsBlockDesign(A) or not IsBlockDesign(B) then
   Error("usage: IsEqualBlockDesign( <BlockDesign>, <BlockDesign> )");
fi;
return A.v=B.v and A.blocks=B.blocks;
end); 

BindGlobal("IsIsomorphicBlockDesign",function(A,B)
#
# This boolean function returns true iff block designs  A  and  B 
# are isomorphic.  The case where both  A  and   B   are non-binary is 
# not yet fully implemented.
#
if not IsBlockDesign(A) or not IsBlockDesign(B) then
   Error("usage: IsIsomorphicBlockDesign( <BlockDesign>, <BlockDesign> )");
fi;
if IsEqualBlockDesign(A,B) then
   return true;
fi;
if A.v<>B.v or BlockSizes(A)<>BlockSizes(B) or BlockNumbers(A)<>BlockNumbers(B) then 
   return false;
fi;
if IsBinaryBlockDesign(A) and not IsBinaryBlockDesign(B) or 
   IsBinaryBlockDesign(B) and not IsBinaryBlockDesign(A) then
   return false;
fi;
if not IsBinaryBlockDesign(A) and not IsBinaryBlockDesign(B) then
   Error("isomorphism test not yet fully implemented for non-binary block designs");
fi;
return Length(BlockDesignIsomorphismClassRepresentatives([A,B]))=1;
end); 

BindGlobal("DualBlockDesign",function(D)
#
# Suppose  D  is a block design for which every point
# lies on at least one block.  Then this function returns
# the dual of  D,  the block design in which the 
# roles of points and blocks are interchanged, but incidence
# (including repeated incidence) stays the same.  Note
# that, since lists of blocks are always sorted, the dual of the
# dual of  D  may not equal  D. 
#
local i,j,dualblocks,dual;
if not IsBlockDesign(D) then
   Error("usage: DualBlockDesign( <BlockDesign> )");
fi;
dualblocks:=List([1..D.v],x->[]);
for i in [1..Length(D.blocks)] do
   for j in D.blocks[i] do
      Add(dualblocks[j],i);
   od;
od;
if ForAny(dualblocks,x->Length(x)=0) then
   Error("each point of <D> must lie on at least one block");
fi;
dual:=BlockDesign(Length(D.blocks),dualblocks);
if IsBound(D.pointNames) then
   dual.pointNames:=Immutable(List(D.blocks,x->D.pointNames{x}));   
else 
   dual.pointNames:=Immutable(D.blocks);   
fi;
return dual;
end); 

BindGlobal("ComplementBlocksBlockDesign",function(D)
#
# Suppose  D  is a binary incomplete-block design.
# Then this function returns the block design on the same
# point-set as  D,  whose blocks are the complements of
# those of  D.
#
local newblocks,newdesign;
if not IsBinaryBlockDesign(D) then
   Error("usage: ComplementBlocksBlockDesign( <BinaryBlockDesign> )");
fi;
if ForAny(D.blocks,x->Length(x)=D.v) then
   Error("<D> must not contain a complete block");
fi;
newblocks:=List(D.blocks,x->Difference([1..D.v],x));
newdesign:=BlockDesign(D.v,newblocks);
if IsBound(D.pointNames) then
   newdesign.pointNames:=Immutable(D.pointNames);   
fi;
if IsBound(D.autGroup) then
   newdesign.autGroup:=D.autGroup;
fi;
return newdesign;
end); 

BindGlobal("DeletedPointsBlockDesign",function(D,Y)
#
# Suppose  D  is a block design, and  Y  is a proper subset of the 
# point-set of  D.
# 
# Then this function returns a block design  F  obtained from  D 
# by deleting the points in  Y  from the point-set of  D,  and
# removing the points in  Y  from each block of  D.
# It is an error if  F  then contains an empty block. 
# The points of  F  are relabelled with   [1..D.v-Length(Y)], 
# preserving the order of the corresponding points in  D. 
# The corresponding point-names of D become the point-names of F.
#
local newblocks,newpoints,newlabels,i,newdesign;
if not (IsBlockDesign(D) and IsSet(Y)) then
   Error("usage: DeletedPointsBlockDesign( <BlockDesign>, <Set> )");
fi;
if Y=[] then
   # no points to delete
   return StructuralCopy(D);
fi;
if not IsSubset([1..D.v],Y) or Length(Y)=D.v then 
   Error("<Y> must be a proper subset of [1..<D>.v]");
fi;
newblocks:=List(D.blocks,x->Filtered(x,y->not (y in Y)));
if ForAny(newblocks,x->x=[]) then
   Error("cannot create a block design with empty blocks");
fi;
newpoints:=Difference([1..D.v],Y);
newlabels:=[];
for i in [1..Length(newpoints)] do
   newlabels[newpoints[i]]:=i;
od;
# now relabel the points in the new blocks
newblocks:=List(newblocks,x->newlabels{x});
newdesign:=BlockDesign(Length(newpoints),newblocks);
if not IsBound(D.pointNames) then 
   newdesign.pointNames:=Immutable(newpoints);
else
   newdesign.pointNames:=Immutable(D.pointNames{newpoints});
fi;
return newdesign;
end); 

BindGlobal("DeletedBlocksBlockDesign",function(D,Y)
#
# Suppose  D  is a block design, and  Y  is a proper sublist of the 
# block-list of  D.  Y  need not be sorted.
# 
# Then this function returns a block design  F  obtained from  D 
# by deleting the blocks in  Y  (counting repeats) from the block-list 
# of  D. 
#
local blocks,newblocks,i,j,jj,newdesign;
if not (IsBlockDesign(D) and IsList(Y)) then
   Error("usage: DeletedBlocksBlockDesign( <BlockDesign>, <List> )");
fi;
if Y=[] then 
   # no blocks to delete
   return StructuralCopy(D);
fi;
blocks:=D.blocks;
if Length(Y)=Length(blocks) then 
   Error("<Y> must be a proper sublist of <D>.blocks");
fi;
if not IsSortedList(Y) then
   Y:=SortedList(Y);
fi;
newblocks:=[];
i:=1;
j:=1;
while i<=Length(Y) and j<=Length(blocks) do
   if Y[i]<blocks[j] then
      Error("<Y> is not a sublist of <D>.blocks");
   elif Y[i]=blocks[j] then
      # deleted block
      i:=i+1;
      j:=j+1;
   else
      # blocks[j] is not deleted
      Add(newblocks,blocks[j]);
      j:=j+1;
   fi;
od;
if i<=Length(Y) then
   # not all elements of  Y  are in  blocks.
   Error("<Y> is not a sublist of <D>.blocks");
fi;
for jj in [j..Length(blocks)] do
   Add(newblocks,blocks[jj]);
od;
newdesign:=BlockDesign(D.v,newblocks);
if IsBound(D.pointNames) then 
   newdesign.pointNames:=Immutable(D.pointNames);
fi;
return newdesign;
end); 

BindGlobal("AddedPointBlockDesign",function(arg)
#
# Suppose  D=arg[1]  is a block design, and  Y=arg[2]  is a sublist 
# of the block-list of  D.  Y  need not be sorted.
# 
# Then this function returns the block design  F  obtained from  D 
# by adding the new point  D.v+1,  and adding this point (once) to each 
# block in  Y  (including repeats).  
#
# The optional parameter  arg[3]  is used to specify a point-name for 
# the new point.
#
local D,Y,newpoint,blocks,i,j,newdesign,pointname;
D:=arg[1];
Y:=arg[2];
if not (IsBlockDesign(D) and IsList(Y)) then
   Error("usage: AddedPointBlockDesign( <BlockDesign>, <List> [, <Obj> ] )");
fi;
if not IsSortedList(Y) then
   Y:=SortedList(Y);
fi;
newpoint:=D.v+1;
blocks:=ShallowCopy(D.blocks);
i:=1;
j:=1;
while i<=Length(Y) and j<=Length(blocks) do
   if Y[i]<blocks[j] then
      Error("<Y> is not a sublist of <D>.blocks");
   elif Y[i]=blocks[j] then
      # add the new point to  blocks[j]
      blocks[j]:=SortedList(Concatenation(blocks[j],[newpoint]));
      i:=i+1;
      j:=j+1;
   else
      # blocks[j] is not altered
      j:=j+1;
   fi;
od;
if i<=Length(Y) then
   # not all elements of  Y  are in  blocks.
   Error("<Y> is not a sublist of <D>.blocks");
fi;
newdesign:=BlockDesign(D.v+1,blocks);
if IsBound(arg[3]) then
   pointname:=arg[3];
else
   pointname:=newpoint;
fi;
if IsBound(D.pointNames) then 
   newdesign.pointNames:=Immutable(Concatenation(D.pointNames,[pointname]));
elif IsBound(arg[3]) then
   newdesign.pointNames:=Immutable(Concatenation([1..D.v],[pointname]));
fi;
return newdesign;
end); 

BindGlobal("AddedBlocksBlockDesign",function(D,Y)
#
# Suppose  Y  is a list of multisets of points of the block design  D.
# Then this function returns a new block design, whose point-set
# is that of  D,  and whose block list is that of  D  with
# the elements of  Y  (including repeats)  added.
#
local newdesign;
if not (IsBlockDesign(D) and IsList(Y)) then
   Error("usage: AddedBlocksBlockDesign( <BlockDesign>, <List> )");
fi;
if Y=[] then 
   # no blocks to add
   return StructuralCopy(D);
fi;
newdesign:=BlockDesign(D.v,Concatenation(D.blocks,Y));
if IsBound(D.pointNames) then 
   newdesign.pointNames:=Immutable(D.pointNames);
fi;
return newdesign;
end); 

BindGlobal("DerivedBlockDesign",function(D,x)
#
# Suppose  D  is a block design, and  x  is a point or block of  D.
# Then this function returns the derived block design  DD  of  D,
# with respect to  x.
#
# If  x  is a point then  DD  is the block design whose blocks 
# are those of  D  containing  x,  but with  x  deleted from these
# blocks, and the points of  DD  are those which occur in some block
# of  DD.
#
# If  x  is a block, then the points of  DD  are the points in  x,  and 
# the blocks of  DD  are the blocks of  D  other than  x  containing 
# at least one point of  x,  but with all points not in  x  deleted from 
# these blocks.  Note that any repeat of  x,  but not  x  itself, 
# is a block of  DD.
#
# It is an error if the resulting block design  DD  has 
# no blocks or an empty block.
#
# The points of  DD  are relabelled with  1,2,...,
# preserving the order of the corresponding points in  D. 
# The corresponding point-names of D become the point-names of DD.
#
local newpoints,newblocks,DD; 
if not (IsBlockDesign(D) and (IsInt(x) or IsList(x))) then
   Error("usage: DerivedBlockDesign( <BlockDesign>, <Int> or <List> )");
fi;
if not (x in BlockDesignPoints(D)) and not (x in BlockDesignBlocks(D)) then
   Error("<x> must be a point or block of <D>");
fi; 
if IsInt(x) then
   # x is a point
   newblocks:=List(Filtered(D.blocks,A->x in A),A->Filtered(A,y->y<>x));
   if newblocks=[] or ForAny(newblocks,x->x=[]) then
      Error("cannot create a block design with no blocks or with empty blocks");
   fi;
   newpoints:=Union(newblocks); 
else 
   # x is a block
   newpoints:=ShallowCopy(x);
   newblocks:=ShallowCopy(D.blocks);
   Remove(newblocks,PositionSorted(newblocks,x));
   newblocks:=Filtered(List(newblocks,A->Filtered(A,y->y in x)),A->A<>[]); 
   if newblocks=[] then
      Error("cannot create a block design with no blocks");
   fi;
fi;
DD:=BlockDesign(D.v,newblocks);
if IsBound(D.pointNames) then 
   DD.pointNames:=Immutable(D.pointNames);
fi;
# now strip out the points not in newpoints, and relabel the points
DD:=DeletedPointsBlockDesign(DD,Difference(BlockDesignPoints(DD),newpoints));
return DD;
end); 

BindGlobal("ResidualBlockDesign",function(D,x)
#
# Suppose  D  is a block design, and  x  is a point or block of  D.
# Then this function returns the residual block design  RD  of  D,
# with respect to  x.
#
# If  x  is a point then  RD  is the block design whose blocks 
# are those of  D  not containing  x,  and the points of  RD  are those 
# which occur in some block of  RD.
#
# If  x  is a block, then the points of  RD  are those of  D  not in  x,  
# and the blocks of  RD  are the blocks of  D  (including repeats) 
# containing at least one point not in  x,  but with all points in  x  
# deleted from these blocks.  
#
# It is an error if the resulting block design  RD  has no blocks.
#
# The points of  RD  are relabelled with  1,2,...,
# preserving the order of the corresponding points in  D. 
# The corresponding point-names of D become the point-names of RD.
#
local newpoints,newblocks,RD; 
if not (IsBlockDesign(D) and (IsInt(x) or IsList(x))) then
   Error("usage: ResidualBlockDesign( <BlockDesign>, <Int> or <List> )");
fi;
if not (x in BlockDesignPoints(D)) and not (x in BlockDesignBlocks(D)) then
   Error("<x> must be a point or block of <D>");
fi; 
if IsInt(x) then
   # x is a point
   newblocks:=Filtered(D.blocks,A->not(x in A));
   if newblocks=[] then
      Error("cannot create a block design with no blocks");
   fi;
   newpoints:=Union(newblocks); 
else 
   # x is a block
   newpoints:=Difference(BlockDesignPoints(D),x);
   newblocks:=Filtered(List(D.blocks,A->Filtered(A,y->not (y in x))),A->A<>[]);
   if newblocks=[] then
      Error("cannot create a block design with no blocks");
   fi;
fi;
RD:=BlockDesign(D.v,newblocks);
if IsBound(D.pointNames) then 
   RD.pointNames:=Immutable(D.pointNames);
fi;
# now strip out the points not in newpoints, and relabel the points
RD:=DeletedPointsBlockDesign(RD,Difference(BlockDesignPoints(RD),newpoints));
return RD;
end); 

BindGlobal("TDesignFromTBD",function(D,t,K)
#
# Let  D  be a  t-(v,K,lambda)  design, where  t  is a positive integer 
# and  K  is a nonempty set of integers such that 
# t <= K[1] < K[2] < ... < K[Length(K)] <= v, and let  k = K[1].  
# Then this function returns a certain   t-(v,k,n*lambda)  design,
# where  n = Lcm(Binomial(K[1]-t,k-t),...,Binomial(K[Length(K)]-t,k-t)). 
#
# If  K  is not a set, then it must be an integer,  k  say, such that 
# t <= k <= Minimum(BlockSizes(D)),  and then this function returns 
# TDesignFromTBD(D, t, Union([k],BlockSizes(D) ),  which is the
# design  D*(t,k),  described and studied in
# J.P. McSorley and L.H. Soicher, Constructing t-designs from t-wise 
# balanced designs, European J. Combinatorics, to appear. 
# 
local k,lambdas,lambda,newblocks,block,n,c,j,newdesign;
if not (IsBinaryBlockDesign(D) and IsInt(t) and (IsInt(K) or IsSet(K))) then
   Error("usage: TDesignFromTBD( <BinaryBlockDesign>, <Int>, <Int> or <Set> )");
fi;
if t<=0 then 
   Error("<t> must be positive");
fi;
if IsInt(K) then 
   k:=K;
   if k<t or k>Minimum(BlockSizes(D)) then 
      Error("must have <t> <= <k> <= Minimum(BlockSizes(<D>))");
   fi;
   K:=Union([k],BlockSizes(D));
else
   #  K  is a set.
   if K=[] or not ForAll(K,IsInt) or K[1]<t or K[Length(K)]>D.v 
      or not IsSubset(K,BlockSizes(D)) then
      Error("invalid <K>");
   fi;
   k:=K[1];
fi;
lambdas:=TSubsetStructure(D,t).lambdas;
if Length(lambdas) <> 1 then
   Error("<D> is not <t>-wise balanced");
fi;
lambda:=lambdas[1];
n:=Lcm(List(K,x->Binomial(x-t,k-t)));
newblocks:=[];
for block in D.blocks do
   for c in Combinations(block,k) do
      for j in [1..n/Binomial(Length(block)-t,k-t)] do 
         Add(newblocks,c);
      od;
   od;
od;
newdesign:=BlockDesign(D.v,newblocks);
if IsBound(D.pointNames) then
   newdesign.pointNames:=Immutable(D.pointNames);
fi;
if lambda=1 and t<k and IsBound(D.autGroup) then
   newdesign.autGroup:=D.autGroup;
elif IsBound(D.autGroup) then
   newdesign.autSubgroup:=D.autGroup;
elif IsBound(D.autSubgroup) then
   newdesign.autSubgroup:=D.autSubgroup;
fi;
return newdesign;
end);

BindGlobal("PGPointFlatBlockDesign",function(n,q,d)
#
# Returns the block design whose points are the (projective) points
# of  PG(n,q)  and whose blocks are the  d-flats  of  PG(n,q) 
# (i.e. the subspaces of projective dimension d). 
# 
local F,V,points,i,flat,block,normedpoints,gens,frobaut,zerovec,
      flatbasis,G,D;
if not IsInt(n) or not IsInt(q) or not IsInt(d) then
   Error("usage: PGPointFlatBlockDesign( <Int>, <Int>, <Int> )");
fi;
if n<0 then 
   Error("<n> must be non-negative");
elif d<0 or d>n then
   Error("<d> must be in [0..<n>]");
fi;
F:=GaloisField(q);
V:=F^(n+1);
points:=AsSet(Subspaces(V,1));
normedpoints:=List(points,x->NormedRowVectors(x)[1]);
# Now calculate a block-transitive group  G  of automorphisms of the design.
gens:=ShallowCopy(GeneratorsOfGroup(Action(GL(n+1,q),normedpoints,OnLines)));
if not IsPrimeInt(q) then
   # Add to  gens  the Frobenius automorphism of  F,  acting on the
   # (indices of the) points.
   frobaut:=FrobeniusAutomorphism(F);
   Add(gens,PermList(List(normedpoints,x->
      PositionSorted(points,SubspaceNC(V,[List(x,y->y^frobaut)],"basis"))))); 
fi;
G:=Group(gens);
# Now make one block of the design.
zerovec:=ListWithIdenticalEntries(n+1,Zero(F));
flatbasis:=[];
for i in [1..d+1] do
   flatbasis[i]:=ShallowCopy(zerovec);
   flatbasis[i][i]:=One(F);
od;
flat:=Subspace(V,flatbasis,"basis");
block:=Filtered([1..Length(points)],i->normedpoints[i] in flat); 
D:=BlockDesign(Length(points),[block],G); 
D.pointNames:=points;
return D;
end);

BindGlobal("AGPointFlatBlockDesign",function(n,q,d)
#
# Returns the block design whose points are the points
# of  AG(n,q)  and whose blocks are the  d-flats  of  AG(n,q) 
# (i.e. the cosets of the subspaces of dimension d). 
# 
local F,V,G,points,zerovec,flatbasis,flat,i,block,gens,translation,
      frobaut,vec,D;
if not IsInt(n) or not IsInt(q) or not IsInt(d) then
   Error("usage: AGPointFlatBlockDesign( <Int>, <Int>, <Int> )");
fi;
if n<1 then 
   Error("<n> must be positive");
elif d<0 or d>n then
   Error("<d> must be in [0..<n>]");
fi;
F:=GaloisField(q);
V:=F^n;
points:=AsSet(Elements(V)); 
# Now calculate a block-transitive group  G  of automorphisms of the design.
gens:=ShallowCopy(GeneratorsOfGroup(Action(GL(n,q),points)));
# Now add to  gens  the translation by  [1,0,0,...,0]. 
vec:=ListWithIdenticalEntries(n,Zero(F));
vec[1]:=One(F); 
translation:=PermList(List(points,x->PositionSorted(points,x+vec)));
Add(gens,translation);
if not IsPrimeInt(q) then
   # Add to  gens  the Frobenius automorphism of  F,  acting on the
   # (indices of the) points.
   frobaut:=FrobeniusAutomorphism(F);
   Add(gens,PermList(List(points,x->
               PositionSorted(points,List(x,y->y^frobaut))))); 
fi;
G:=Group(gens);
# Now make one block of the design.
zerovec:=ListWithIdenticalEntries(n,Zero(F));
flatbasis:=[];
for i in [1..d] do
   flatbasis[i]:=ShallowCopy(zerovec);
   flatbasis[i][i]:=One(F);
od;
flat:=Subspace(V,flatbasis,"basis");
block:=Filtered([1..Length(points)],i->points[i] in flat); 
D:=BlockDesign(Length(points),[block],G); 
D.pointNames:=points;
return D;
end);

BindGlobal("WittDesign",function(n)
local M12,M24,W12,W24,W,i;
if not IsInt(n) then
   Error("usage: WittDesign( <Int> )");
fi;
if not (n in [9,10,11,12,21,22,23,24]) then
   Error("must have <n> in [9,10,11,12,21,22,23,24]");
fi;
if n<=12 then
   M12:=Group( [ (1,2,3,4,5,6,7,8,9,10,11), 
      (3,7,11,8)(4,10,5,6), (1,12)(2,11)(3,6)(4,8)(5,9)(7,10) ] );
   W12:=BlockDesign(12,[[1,2,3,4,5,7]],M12);
   W12.autGroup:=M12;
   Unbind(W12.autSubgroup);
   W12.allTDesignLambdas:=Immutable(TDesignLambdas(5,12,6,1));
   W:=W12;
   for i in Reversed([n+1..12]) do
      W:=DerivedBlockDesign(W,i);
   od;
else
   M24:=Group([ (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23), 
      (3,17,10,7,9)(4,13,14,19,5)(8,18,11,12,23)(15,20,22,21,16), 
      (1,24)(2,23)(3,12)(4,16)(5,18)(6,10)(7,20)(8,14)(9,21)(11,17)(13,22)
      (15,19) ]);
   W24:=BlockDesign(24,[[1,2,3,4,5,8,11,13]],M24);
   W24.autGroup:=M24;
   Unbind(W24.autSubgroup);
   W24.allTDesignLambdas:=Immutable(TDesignLambdas(5,24,8,1));
   W:=W24;
   for i in Reversed([n+1..24]) do
      W:=DerivedBlockDesign(W,i);
   od;
fi;
return W;
end);

BindGlobal("BlockDesigns",function(param)
#
# Function to classify block designs with given properties. 
#
# These block designs need not be simple and block sizes 
# need not be constant.  These block designs need not even be
# binary, although this is the default.
#
local t,v,b,blocksizes,k,lambda,blocknumbers,blockmaxmultiplicities,
   r,blockintersectionnumbers,isolevel,C,G,B,N,ntflags,
   act,rel,weightvector,gamma,KK,L,S,s,hom,GG,CC,NN,leastreps,
   clique,ans,blockbags,c,d,i,j,jj,issimple,allbinary,tsubsetstructure,
   blocks,blockdesign,A,kk,AA,tsubsets,maxlambda,targetvector,weightvectors,
   designinfo,lambdavec,lambdamat,m,T,bb,justone,isnormal,issylowtowergroup,
   testfurther;

act := function(x,g) 
# The boolean variable  allbinary  is global, and  allbinary=true 
# iff all possible blocks are sets.
if allbinary or IsSet(x.comb) then
   return rec(comb:=OnSets(x.comb,g),mult:=x.mult); 
else
   return rec(comb:=OnMultisetsRecursive(x.comb,g),mult:=x.mult); 
fi;
end;
             
weightvector := function(x)
# param, v, t, blocksizes, and tsubsets are global
local wv,c,i,xx;
wv:=ListWithIdenticalEntries(Length(tsubsets),0);
xx:=ListWithIdenticalEntries(v,0);
for i in x.comb do 
   xx[i]:=xx[i]+1;
od;
for c in Combinations(Set(x.comb),t) do
   wv[PositionSorted(tsubsets,c)]:=x.mult*Product(xx{c});
od;
if IsBound(param.r) then 
   Append(wv,x.mult*xx);
fi;
if IsBound(param.blockNumbers) then 
   for i in [1..Length(blocksizes)] do
      if Length(x.comb)=blocksizes[i] then
         Add(wv,x.mult);
      else
         Add(wv,0);
      fi;
   od;
fi;
if IsBound(param.b) then
   Add(wv,x.mult);
fi;
return wv;
end;

rel := function(x,y) 
# v, blocksizes, targetvector, blockbags, weightvectors, 
# and blockintersectionnumbers are global.
# The parameters x and y are indices into blockbags (and weightvectors).
local i,xx,yy,s;
if blockbags[x].comb=blockbags[y].comb then
   return false;
fi;
s:=weightvectors[x]+weightvectors[y];
if HasLargerEntry(s,targetvector) then
   return false;
fi;
if allbinary or (IsSet(x) and IsSet(y)) then
   s:=Size(Intersection(blockbags[x].comb,blockbags[y].comb));
else 
   xx:=ListWithIdenticalEntries(v,0);
   yy:=ListWithIdenticalEntries(v,0);
   for i in blockbags[x].comb do 
      xx[i]:=xx[i]+1;
   od;
   for i in blockbags[y].comb do 
      yy[i]:=yy[i]+1;
   od;
   s:=0;
   for i in [1..v] do
      s:=s+Minimum(xx[i],yy[i]);
   od;
fi;
return 
   s in blockintersectionnumbers[Position(blocksizes,Length(blockbags[x].comb))][Position(blocksizes,Length(blockbags[y].comb))]; 
end;

#
# begin BlockDesigns
# 
if not IsRecord(param) then
   Error("usage: BlockDesigns( <Record> )");
fi;
param:=ShallowCopy(param);
if not IsSubset(["v","blockSizes","tSubsetStructure","blockDesign",
                 "blockMaxMultiplicities","blockIntersectionNumbers",
   "r","blockNumbers","b","isoLevel","isoGroup",
   "requiredAutSubgroup"], 
                RecNames(param)) then
   Error("<param> contains an invalid component-name");
fi;   
if not IsSubset(RecNames(param),["v","blockSizes","tSubsetStructure"]) then
   Error("<param> missing a required component");
fi;   
v:=param.v;
if not IsPosInt(v) then
   Error("<param>.v must be positive integer");
fi;
blocksizes:=ShallowCopy(param.blockSizes);
if not (IsSet(blocksizes) and blocksizes<>[] and ForAll(blocksizes,IsPosInt)) then
   Error("<param>.blockSizes must be a non-empty set of positive integers"); 
fi;
if Length(blocksizes)=1 then 
   k:=blocksizes[1];
fi;
if not IsBound(param.blockDesign) then
   allbinary:=true;
else
   allbinary:=IsBinaryBlockDesign(param.blockDesign);
fi;
# Note: allbinary=true iff all possible blocks are sets.
tsubsetstructure:=ShallowCopy(param.tSubsetStructure);
if not ForAll(tsubsetstructure.lambdas,x->IsInt(x) and x>=0) then
   Error("all <param>.tSubsetStructure.lambdas must be non-negative integers");
fi;
if not IsDuplicateFree(tsubsetstructure.lambdas) then
   Error("<param>.tSubsetStructure.lambdas must not contain duplicates");
fi;
if IsBound(tsubsetstructure.partition) then
   if Length(tsubsetstructure.partition)<>Length(tsubsetstructure.lambdas) then
      Error("<param>.tSubsetStructure.partition must have the same length\n",
            "as <param>.tSubsetStructure.lambdas");
   fi;
elif Length(tsubsetstructure.lambdas)<>1 then
   Error("must have Length(<param>.tSubsetStructure.lambdas)=1\n",
         "if <param>.tSubsetStructure.partition is unbound");
fi;
t:=tsubsetstructure.t;
if not (IsInt(t) and t>=0 and t<=v) then
   Error("<t> must be an integer with 0<=<t><=<v>");
fi;
if not ForAll(blocksizes,x->x>=t) then 
   Error("each element of <blocksizes> must be >= <t>");
fi;
if IsBound(tsubsetstructure.partition) then
   # check it
   if not ForAll(tsubsetstructure.partition,x->IsSet(x) and x<>[]) then 
      Error("the parts of the t-subset partition must be non-empty sets");
   fi;
   c:=Concatenation(tsubsetstructure.partition);
   if not ForAll(c,x->IsSet(x) and Size(x)=t) then
      Error("the parts of the t-subset partition must be sets of t-subsets");
   fi;
   if Length(c)<>Binomial(v,t) or Length(Set(c))<>Binomial(v,t) then
      Error("t-subset partition is not a partition of the t-subsets");
   fi;
fi;
maxlambda:=Maximum(tsubsetstructure.lambdas);
if maxlambda<1 then
   Error("at least one element of <param>.tSubsetStructure.lambdas must be positive");
fi;
if Length(tsubsetstructure.lambdas)=1 then
   # constant lambda 
   lambda:=tsubsetstructure.lambdas[1];
fi;
tsubsets:=Combinations([1..v],t);
if IsBound(param.blockMaxMultiplicities) then
   blockmaxmultiplicities:=ShallowCopy(param.blockMaxMultiplicities);
   if not (IsList(blockmaxmultiplicities) and ForAll(blockmaxmultiplicities,x->IsInt(x) and x>=0)) then
      Error("<param>.blockMaxMultiplicities must be a list of non-negative integers");
   fi;   
   if Length(blockmaxmultiplicities)<>Length(blocksizes) then 
      Error("must have Length(<param>.blockMaxMultiplicities)=Length(<param>.blockSizes)");
   fi;
   blockmaxmultiplicities:=List(blockmaxmultiplicities,x->Minimum(x,maxlambda));
   # since *every* possible block is required to contain at least 
   # t  distinct points.
else 
   blockmaxmultiplicities:=ListWithIdenticalEntries(Length(blocksizes),maxlambda);
fi;
if IsBound(param.blockIntersectionNumbers) then
   blockintersectionnumbers:=StructuralCopy(param.blockIntersectionNumbers);
   if Length(blockintersectionnumbers)<>Length(blocksizes) then 
      Error("must have Length(<param>.blockIntersectionNumbers>)=Length(<param>.blockSizes>)");
   fi;
   blockintersectionnumbers:=List(blockintersectionnumbers,x->List(x,Set));
   if blockintersectionnumbers<>TransposedMat(blockintersectionnumbers) then
      Error("<blockintersectionnumbers> must be a symmetric matrix");
   fi;
else 
   blockintersectionnumbers:=List([1..Length(blocksizes)],x->List([1..Length(blocksizes)],
                      y->[0..Minimum(blocksizes[x],blocksizes[y])]));
fi;
if allbinary and maxlambda<=1 then 
   blockintersectionnumbers:=List(blockintersectionnumbers,x->List(x,y->Intersection(y,[0..t-1])));
fi;
# Compute the number  ntflags  of (t-subset,block)-flags 
# (counting multiplicities). 
if IsBound(lambda) then
   ntflags:=lambda*Binomial(v,t); 
else
   ntflags:=tsubsetstructure.lambdas*List(tsubsetstructure.partition,Length);
fi;
if IsBound(param.blockNumbers) then
   blocknumbers:=ShallowCopy(param.blockNumbers);
   # We will require  blocknumbers[i]  blocks of size  blocksizes[i] 
   # for i=1,...,Length(blocknumbers).
   if not (IsList(blocknumbers) and ForAll(blocknumbers,x->IsInt(x) and x>=0)) then
      Error("<param>.blockNumbers must be a list of non-negative integers"); 
   fi;
   if Length(blocknumbers)<>Length(blocksizes) then
      Error("must have Length(<param>.blockNumbers)=Length(<param>.blockSizes)");
   fi;
   b:=Sum(blocknumbers);
   if b=0 then
      Error("At least one element of <param>.blockNumbers must be positive");  
   fi;
   if allbinary and IsBound(ntflags) then 
--> --------------------

--> maximum size reached

--> --------------------