\\**************Instructions**************
\\
\\Copy and paste this file into PARI/GP. 
\\This will define a function CoveringMatrix(q,n,m).
\\This function either returns a minimal n x m covering 
\\matrix, or proves that no such matrix exists.
\\For example, enter CoveringMatrix(3,5,12) to find the
\\5 x 12 matrix for q = 3.

OneSearch(line)=  
{
gettime;

if(CheckKupon(line),,onetime+=[1,gettime];return(0));
if(kind==#Kinder,Kinder=Double(Kinder));kind++;Kinder[kind]=line;
onetime+=[1,gettime];
return(1)
}

CompleteSearch(permline,sym,Sym)= 
{
local(BS,BL,C,OTT,antalperm,www,A,SP);
BS=BlokStart(sym);BL=BlokLength(sym);OTT=MakeOTT(permline,sym);SP=SamePriority(OTT);
www=vector(#BS,j,if(j==1,RestrictedMulti(OTT[j],vector(BL[j],i,Sym[BS[j]-1+i])),Multi(OTT[j],vector(BL[j],i,Sym[BS[j]-1+i]))));
antalperm=prod(j=1,#BS,#www[j]);
if(#www>1,
A=matrix(#SP,#www,i,blok,
  vector(#www[blok],j,ComparePriority(www[blok][j],vector(BL[blok],ii,Sym[BS[blok]-1+ii]),SP[i][1],SP[i][2])));
);
  forvec(vvv=vector(#www,j,[1,#www[j]]),
  taeller++;

if(#www>1,
  for(i=1,#SP,
   for(blok=1,#vvv,
    if(A[i,blok][vvv[blok]]==0,for(j=blok+1,#www,vvv[j]=#www[j]);oktime+=[1,gettime];next(3));
    if(A[i,blok][vvv[blok]]==2,next(2));
   );
  );
);

  C=connect(vector(#www,j,www[j][vvv[j]]));
  if(CheckKupon(C),,for(blok=sym[fail]+1,#www,vvv[blok]=#www[blok]);next); 
  if(kind==#Kinder,Kinder=Double(Kinder));kind++;Kinder[kind]=Canonise(C,Sym);
  );
return(1);
}

CheckKupon(Values)=   
{
local(Howmany,t);
gettime;
t=X;for(i=1,X,if(Values[i]==Husk[1][i],,t=i-1;break)); 
if(t>=Husk[2],fail=Husk[2];checktime+=[1,gettime];return(0));
for(i=SubsetsOrder[l-1][t+1],#Subsets[l-1],
Howmany=vector(b);  
  forstep(j=#Subsets[l-1][i],1,-1,
  Howmany[Values[Subsets[l-1][i][j]]]++;
  if(#(Subsets[l-1][i])-Howmany[Values[Subsets[l-1][i][j]]]<mustremain,fail=Subsets[l-1][i][j];Husk=[Values,fail];checktime+=[1,gettime];return(0));
  );
);
Husk=[Values,X+1];checktime+=[1,gettime];return(1)
}


FindLineSet(q)=  
{
print("Last row: "LineSet[l-1][Counter[l-1]]);
    for(j=2,if(l>1,#ExtraSymmetry(l-2),0), 
    if(Lexicographic(LineSet[l-1][Counter[l-1]],
        Canonise(Permutation(LineSet[l-1][Counter[l-1]],ExtraSymmetry(l-2)[j]),Symmetry(l-2))),,
         LineSet[l]=[];print("Last row rejected");return(0))
    );
print("Last row accepted");

Kandidater=vector(100);  
StammKandidater=vector(100); 
antalkandidater=0;
mustremain=if(l==N,1,Omega[N-l]);
symmetri=Symmetry(l-2);
gotoonesearch=(symmetri==vector(X,i,i));
Husk=[vector(X),1];  

for(last=Counter[l-1],#LineSet[l-1],
Kinder=vector(100);kind=0;
if(gotoonesearch, OneSearch(LineSet[l-1][last]), CompleteSearch(LineSet[l-1][last],symmetri,Symmetri[l-1]));

Kinder=vecsort(RemoveZeros(Kinder),,2);
for(i=1,#Kinder,    for(j=1,l-1,for(t=1,#ExtraSymmetry(j-1),if(Lexicographic(Kupon[j],
                    Canonise(Permutation(Kinder[i],ExtraSymmetry(j-1)[t]),Symmetry(j-1))),,next(3))));
                if(antalkandidater==#Kandidater,Kandidater=Double(Kandidater);StammKandidater=Double(StammKandidater));
                antalkandidater++;Kandidater[antalkandidater]=Kinder[i];
                StammKandidater[antalkandidater]=vector(l,j,if(j<l,Stammbaum[l-1][last][j],i)));
);

LineSet[l]=vector(antalkandidater);Stammbaum[l]=vector(antalkandidater);
for(i=1,antalkandidater,LineSet[l][i]=Kandidater[i];Stammbaum[l][i]=StammKandidater[i];
    print(LineSet[l][i]" "Stammbaum[l][i]));
print("Found "#LineSet[l]" candidates for row "l);
return(1)
}

kill(CoveringMatrix);
CoveringMatrix(b,N,X)=   
{
Start=[1,1,1];
Omega=[0,[2,4,8,16,32,64,128],[2,3,5,8,12,18,29,44],[2,3,4,7,10],[2,3,4,5,8,11],[2,3,4,5,6,10],[2,3,4,5,6,7,11]][b];
kuponer=0; 
Kandidater=vector(100);StammKandidater=vector(100);
gettime;totalline=0;taeller=0;linesets=0;
checktime=[0,0];maketime=[0,0];multitime=[0,0];canonicaltime=[0,0];removetime=[0,0];oktime=[0,0];matrixtime=[0,0];
restrtime=[0,0];extratime=[0,0,0];makeextratime=[0,0];Ttime=[0,0];TTtime=[0,0];onetime=[0,0];
Kupon=vector(N);
Subsets=vector(N,l,vector(b^l)); 
SubsetsOrder=vector(N);Poolsets=vector(N);PoolsetsOrder=vector(N);
Counter=vector(N);LineSet=vector(N);Stammbaum=vector(N);Symmetri=vector(N);ExtraSymmetri=vector(N);
maxantalkandidater=0;maxkinder=0;
l=1;

n=0;LineSet[1]=vector(100);
maks=X-Omega[N-1];mini=X-(b-1)*maks;if(mini<0,mini=0);
forvec(v=vector(b,i,[mini,maks]),if(sum(i=1,b,v[i])==X,,next);
n++;LineSet[1][n]=ReverseOTT([vector(b,i,v[b-i+1])])
,1);
LineSet[1]=vecsort(RemoveZeros(LineSet[1]),,2);
totalline=#LineSet[1];linesets=1;
Stammbaum[1]=vector(#LineSet[1],i,[i]);
print("Searching for row l=1:");for(i=1,#LineSet[1],print(LineSet[1][i]" "Stammbaum[1][i]));

print("Found "#LineSet[1]" candidates for row 1");
if(#LineSet[1]==0,print("NO COVERING MATRIX FOUND");break(5));
Counter[1]=Start[1];
Kupon[1]=LineSet[1][Counter[1]];
MakeSym(1);
MakeSubsets(1);
MakeExtraSym(1);
l=2;

while(1,
print1("Searching for row "l" given ");for(j=1,l-1,print1([Counter[j],#LineSet[j]]));print();

FindLineSet(9);totalline+=#LineSet[l];linesets++;

if(l==N & LineSet[l]<>[], 
 for(i=1,#LineSet[l], 
    for(j=2,#ExtraSymmetry(l-1), 
        if(Lexicographic(LineSet[l][i],Canonise(Permutation(LineSet[l][i],ExtraSymmetry(l-1)[j]),Symmetry(l-1))),,
        next(2))
    );
 Kupon[l]=LineSet[l][i];Counter[l]=i;
 if(Dobbelttjek(Kupon),,next);
 kuponer++;
 print("MINIMAL COVERING MATRIX FOUND:");for(j=1,N,print(Kupon[j]));break(5);
 );
LineSet[l]=[]);

if(LineSet[l]<>[],Counter[l]=0;l++); 
while(Counter[l-1]==#LineSet[l-1],l--;if(l==1,break(2)));
Counter[l-1]++;
if(l==3&Counter[1]==Start[1]&Counter[2]<Start[2],Counter[2]=Start[2]); 
if(l==4&Counter[1]==Start[1]&Counter[2]==Start[2]&Counter[3]<Start[3],Counter[3]=Start[3]);
Kupon[l-1]=LineSet[l-1][Counter[l-1]];MakeSubsets(l-1);MakeSym(l-1);MakeExtraSym(l-1);
);

print("NO COVERING MATRICES FOUND")
}

MakeSym(l)=
{
if(l==1,Symmetri[1]=Kupon[1];return(1));
Symmetri[l]=vector(X);
Symmetri[l][1]=1;
for(i=2,X,
 if(Symmetri[l-1][i]==Symmetri[l-1][i-1] & Kupon[l][i]==Kupon[l][i-1],
  Symmetri[l][i]=Symmetri[l][i-1],Symmetri[l][i]=Symmetri[l][i-1]+1);
)
}

Symmetry(l)= 
{
return(if(l==0,vector(X,i,1),Symmetri[l]))
}

BlokStart(Sym)=  
{
local(v);
v=vector(Sym[#Sym]);
forstep(i=#Sym,1,-1,v[Sym[i]]=i);
return(v)
}

BlokSlut(Sym)=  
{
local(v);
v=vector(Sym[#Sym]);
for(i=1,#Sym,v[Sym[i]]=i);
return(v)
}

BlokLength(Sym)= 
{
local(v);
v=vector(Sym[#Sym]);
for(i=1,#Sym,v[Sym[i]]++);
return(v);
}

MakeSubsets(ll)=
{
gettime;
if(ll==1,Subsets[1]=vector(b,i,vector(X,j,X+1-j));
         for(j=1,X,Subsets[1][Kupon[1][j]][X+1-j]=0);
         for(j=1,b,Subsets[1][j]=RemoveZeros(Subsets[1][j]));
         Subsets[1]=vecsort(Subsets[1],,2);
         SubsetsOrder[ll]=vector(X,i,1);
         for(i=2,X,for(j=SubsetsOrder[ll][i-1],#Subsets[ll],if(Subsets[ll][j][1]>=i,SubsetsOrder[ll][i]=j;next(2))));
         return(1));
for(i=1,#Subsets[ll-1],
for(j=1,b,Subsets[ll][b*(i-1)+j]=Subsets[ll-1][i]);
for(j=1,#(Subsets[ll-1][i]),Subsets[ll][b*(i-1)+Kupon[ll][Subsets[ll-1][i][j]]][j]=0);
for(j=1,b,Subsets[ll][b*(i-1)+j]=RemoveZeros(Subsets[ll][b*(i-1)+j]));
);
Subsets[ll]=vecsort(Subsets[ll],,2);
SubsetsOrder[ll]=vector(X+1,i,1);
for(i=2,X,for(j=SubsetsOrder[ll][i-1],#Subsets[ll],if(Subsets[ll][j][1]>=i,SubsetsOrder[ll][i]=j;next(2))));
SubsetsOrder[ll][X+1]=#Subsets[ll]+1;
maketime+=[1,gettime];
}

RemoveZeros(vektor)=
{
return(vecextract(vektor,sum(i=1,#vektor,if(vektor[i],shiftmul(1,i-1),0))))
}

Double(A)=  
{
return(vector(2*#A,i,if(i<=#A,A[i],0)))
}

Transpose(A)= 
{
return(vector(#A[1],i,vector(#A,j,A[j][i])))
}

Bagfra(v)=
{
return(vector(#v,i,v[#v+1-i]))
}

kill(connect)     
connect(uu)=
{
lenn=sum(i=1,#uu,#(uu[i]));
wek=vector(lenn);
pri=0;
for(i=1,#(uu),for(j=1,#(uu[i]),pri++;wek[pri]=uu[i][j]));
return(wek)
}

Lexicographic(v,w)= 
{
for(i=1,#v,if(v[i]<w[i],return(1));if(v[i]>w[i],fail=i;return(0)));
return(1)
}

MakeOTT(v,sym)= 
{
local(OTT);
OTT=vector(sym[#sym],i,vector(b));
for(i=1,#sym,OTT[sym[i]][v[i]]++);
return(OTT);
}

ReverseOTT(OTT)=
{
local(n,v);
n=0;
v=vector(sum(i=1,#OTT,sum(j=1,#OTT[i],OTT[i][j])));
for(i=1,#OTT,for(j=1,#OTT[i],for(k=1,OTT[i][j],n++;v[n]=j)));
return(v)
}

Dobbelttjek(Kupon)=
{
for(i=1,N,for(j=i,N, 
 for(t=1,#ExtraSymmetry(i-1),
 if(Lexicographic(Kupon[i],Canonise(Permutation(Kupon[j],ExtraSymmetry(i-1)[t]),Symmetry(i-1))),,
 return(1/0));
 );
));
return(1)
}

IsIncreasing(v,sym,T=1000000)=  
{
for(i=2,min(#v,T),if(sym[i]==sym[i-1] & v[i]<v[i-1],return(0)));
return(1)
}

MakeIncreasing(v,sym)=
{
local(w,OTT,n);
w=vector(#v);OTT=MakeOTT(v,sym);n=0;
for(blok=1,sym[#sym],for(i=1,b,for(j=1,OTT[blok][i],n++;w[n]=i)));
return(w)
}

IsCanonical(v,Sym,T=X)=
{
local(OTT,n);gettime;
OTT=MakeOTT(v,Sym);n=Sym[T];if(T<X & Sym[T]==Sym[T+1],n--);fail=n+1;
 for(i=1,b-1,for(blok=1,fail-1,
 if(OTT[blok][i]<OTT[blok][i+1],fail=blok;break);
 if(OTT[blok][i]>OTT[blok][i+1],break);
 ));
canonicaltime+=[1,gettime];
if(fail==n+1,return(1),fail=BlokStart(Sym)[fail];return(0))
}

Priority(v,sym,a)= 
{
local(P);
P=vector(sym[#sym]);
for(j=1,#sym,if(v[j]==a,P[sym[j]]++));
return(P)
}

ComparePriority(v,sym,i,j)= 
{
local(P,Q);
P=Priority(v,sym,i);Q=Priority(v,sym,j);
for(k=1,#P,  if(P[k]>Q[k],return(2));  if(P[k]<Q[k],return(0))  );
return(1)
}

SamePriority(OTT)= 
{
local(T,L);
T=Transpose(OTT);
L=listcreate(b^2);
for(i=1,b,for(j=i+1,b,if(T[i]==T[j],listput(L,[i,j]))));
return(L)
}

Canonise(v,sym)=  
{
return(ReverseOTT(Transpose(vecsort(Transpose(MakeOTT(v,sym)),,6))))
}

Permutation(line,si)=  
{
local(v);
v=vector(#line);
for(i=1,#line,v[si[i]]=line[i]);
return(v)
}

InducedPermutation(line,si)=
{
local(perm,A,B);
perm=vector(#line,iii,iii);
 for(i=1,b,
 if(si[i]==i,next);
 A=SymbolPladser(line,i);B=SymbolPladser(line,si[i]); 
 for(j=1,#A,perm[B[j]]=A[j]);
 );
return(perm)
}

SymbolPladser(line,i)=   
{
local(A);
A=listcreate(#line);
for(j=1,#line,if(line[j]==i,listput(A,j)));
return(A)
}

ExtraSymmetry(l)=  
{
return(if(l==0,[vector(X,i,i)],ExtraSymmetri[l]))
}

MakeExtraSym(l)=  
{
local(A,B,sym,T,K,Values);
gettime;
A=vector(l);

 for(i=1,l,
 T=Transpose(MakeOTT(Kupon[i],Symmetry(i-1)));
 B=listcreate(b!);  
  for(j=1,b!,
  si=numtoperm(b,j);
  for(k=1,b,if(T[k]==T[si[k]],,next(2)));
  listput(B,si);
  );
 A[i]=B;  
 );

B=listcreate(prod(i=1,l,#A[i]));  
 forvec(index=vector(l,i,[1,#A[i]]),                    
 K=vector(l,i,vector(X,j,A[i][index[i]][Kupon[i][j]])); 
 Values=vector(X,k,k);                                  
  for(i=1,l,    
  T=InversePermutation(vecsort(-Transpose(MakeOTT(K[i],Symmetry(i-1))),,3));
  for(j=1,X,K[i][j]=T[K[i][j]]);  
  T=vecsort(vector(X,j,[Symmetry(i-1)[j],K[i][j]]),,3);  
  for(j=i,l,K[j]=vector(X,k,K[j][T[k]]));             
  Values=vector(X,k,Values[T[k]]);                    
  if(K[i]<>Kupon[i],next(2));
  );
 if(K<>vector(l,i,Kupon[i]),1/0);
 listput(B,InversePermutation(Values));
 );  

ExtraSymmetri[l]=listsort(B,1);
makeextratime+=[1,gettime];
return(1)
}

PermutationProduct(si,tau)=  
{
local(rho);
rho=vector(#si);
for(j=1,#rho,rho[j]=tau[si[j]]);
return(rho);
}

InversePermutation(perm)=
{
local(v);
v=vector(#perm);
for(i=1,#perm,v[perm[i]]=i);
return(v)
}


Extra(si)=  
{
local(sym,OTT,A,B);
sym=if(i==1,vector(X,j,1),Symmetri[i-1]);
OTT=MakeOTT(Kupon[i],sym);
for(j=1,#si, for(ii=1,#OTT, if(OTT[ii][j]<>OTT[ii][si[j]],return(0))));
for(j=i+1,l, if(Kupon[j]<>Canonise(vector(X,k,si[Kupon[j][k]]),Symmetri[j-1]), return(0)));
return(1);
}

kill(multinom);multinom(a,b,c)=(a+b+c)!/a!/b!/c!;

Multi(OTT,sym)= 
{
local(abc,Values,pointer,ott,n,C,r,s,CCCC);gettime;
CCCC=vector(100);
abc=sum(i=1,#OTT,OTT[i]);
pointer=abc;
ott=OTT;    
Values=ReverseOTT([OTT]);
CCCC[1]=Values;n=1;

while(pointer>0,
if(Values[pointer]==b,ott[b]--;pointer--;next);
ott[Values[pointer]]--;Values[pointer]++;ott[Values[pointer]]++;
for(i=2,b,if(ott[i]>OTT[i],next(2)));

r=0;while(pointer+r+1<=abc & sym[pointer+r+1]==sym[pointer],r++);  
s=sum(i=Values[pointer],b,OTT[i]-ott[i]); 
if(s<r,next);
for(i=pointer+1,pointer+r,for(symbol=Values[i-1],b,if(ott[symbol]<OTT[symbol],Values[i]=symbol;ott[symbol]++;next(2))));
for(i=pointer+r+1,abc,    for(symbol=1          ,b,if(ott[symbol]<OTT[symbol],Values[i]=symbol;ott[symbol]++;next(2))));
pointer=abc;ott=OTT;n++;CCCC[n]=Values;if(n==#CCCC,CCCC=Double(CCCC));
);

C=vector(n);for(i=1,n,C[i]=CCCC[i]);multitime+=[1,gettime];return(C)
}

RestrictedMulti(OTT,sym)= 
{
local(abc,Values,pointer,ott,n,C,r,s,t,CCCC);gettime;
CCCC=vector(100);
abc=sum(i=1,#OTT,OTT[i]);
pointer=abc;
ott=OTT;    
Values=ReverseOTT([OTT]);
CCCC[1]=Values;n=1;

while(pointer>0,
if(Values[pointer]==b,ott[b]--;pointer--;next);
ott[Values[pointer]]--;Values[pointer]++;ott[Values[pointer]]++;
for(i=2,b,if(ott[i]>OTT[i],next(2)));

r=0;while(pointer+r+1<=abc & sym[pointer+r+1]==sym[pointer],r++);  
s=sum(i=Values[pointer],b,OTT[i]-ott[i]); 
if(s<r,next);

t=pointer;while(t>0,if(t==#sym,break);if(sym[t+1]>sym[t],break);t--); 
if(t>0,for(i=1,#SP,if(ComparePriority(vector(t,j,Values[j]),vector(t,j,sym[j]),SP[i][1],SP[i][2]),,next(2))));

for(i=pointer+1,pointer+r,for(symbol=Values[i-1],b,if(ott[symbol]<OTT[symbol],Values[i]=symbol;ott[symbol]++;next(2))));
for(i=pointer+r+1,abc,    for(symbol=1          ,b,if(ott[symbol]<OTT[symbol],Values[i]=symbol;ott[symbol]++;next(2))));
pointer=abc;ott=OTT;

t=pointer;while(t>0,if(t==#sym,break);if(sym[t+1]>sym[t],break);t--); 
if(t>0,for(i=1,#SP,if(ComparePriority(vector(t,j,Values[j]),vector(t,j,sym[j]),SP[i][1],SP[i][2]),,next(2))));

n++;CCCC[n]=Values;if(n==#CCCC,CCCC=Double(CCCC));
);

C=vector(n);for(i=1,n,C[i]=CCCC[i]);restrtime+=[1,gettime];return(C)
}

MagicSquare(OTT,BL)=
{
local(L,v,colgulv,colloft,rowgulv,rowloft);
L=listcreate(100000);
v=vector(#BL,i,vector(#OTT));
rowgulv=vector(#BL+1,i,vector(#OTT));
rowloft=vector(#BL+1,i,vector(#OTT));
for(i=1,#OTT,rowgulv[1][i]=OTT[i]-sum(j=2,#BL,BL[j]));
for(i=1,#OTT,rowloft[1][i]=OTT[i]);
colgulv=vector(#BL,i,vector(#OTT+1));
colloft=vector(#BL,i,vector(#OTT+1));
for(i=1,#BL,colgulv[i][1]=BL[i]-sum(j=2,#OTT,OTT[j]));
for(i=1,#BL,colloft[i][1]=BL[i]);
loop(1,1);
return(L)
} 

loop(row,col)=
{
for(i=max(0,max(rowgulv[col][row],colgulv[col][row])),min(rowloft[col][row],colloft[col][row]),
v[col][row]=i;
if(col<#BL,rowgulv[col+1][row]=rowgulv[col][row]-i+BL[col+1];
           rowloft[col+1][row]=rowloft[col][row]-i);
if(row<#OTT,colgulv[col][row+1]=colgulv[col][row]-i+OTT[row+1];
            colloft[col][row+1]=colloft[col][row]-i);
if(row<#OTT,            loop(row+1,col  );next);
if(row==#OTT & col<#BL, loop(1,    col+1);next);
listput(L,v);
)
}