\\**************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]]]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[], 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]=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],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]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]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(s0, 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(s0,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]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); ) }