with(linalg): with(LinearAlgebra): with(combinat): ### Bit #### bit:=[1,-1]; ### CMM cmm:=proc(XX,YY) evalm(XX&*YY-YY&*XX) end: "cmm(X,Y) computes the commutator of two matrices"; ### GoLie goLie:=proc(AA,BB) local i,RK,T,j,eta,idx,CC,MaxRnk,OldRnk,k,kk,new,Rnk;global xi,ddm; idx:=2; xi:='xi'; MaxRnk:=2; OldRnk:=0; ddm:=2; new:=true; xi[1]:=evalm(AA); xi[2]:=evalm(BB): CC[1]:=xi[1]; CC[2]:=xi[2]; ## print(xi[1],xi[2]); for kk from 2 while new or kk<=ddm do: new:=false; for k from 1 to ddm-1 do: idx:=idx+1; CC[idx]:=cmm(xi[kk],xi[k]); if(CC[idx]=0) then idx:=idx-1; next fi; Rnk:=goRank(seq(CC[t],t=1..idx)); #print("AT ",kk," TRYING ",k," RANKS ",Rnk,MaxRnk," DIM ",ddm); if(Rnk>MaxRnk) then ddm:=ddm+1; xi[ddm]:=CC[idx]; MaxRnk:=Rnk; new:=true; print(ddm,"new basis element",rank(xi[ddm])); fi; od:od: for i from 1 to ddm do RK[i]:=rank(xi[i]); eta[i]:=evalm(xi[i]); od; T:=sort([seq(RK[t],t=1..ddm)]);print("ranks ",seq(T[ddm-t],t=0..ddm-1)); ## for i from 1 to ddm do ## for j from 1 to ddm do ## if rank(eta[j]) = T[1+ddm-i] then xi[i]:=evalm(eta[j]); eta[j]:=0; break; fi; ## od; ## od; "Lie algebra has dimension ", ddm end:"goLie generates the Lie algebra generated by 2 matrices"; ### GoQuikLie goQuikLie:=proc(RR,BB) local URAN,REP,TR,TB,nr,nb,NR,NB; global NLBAS; TR:=transpose(RR);TB:=transpose(BB); nr:=rowdim(RR)-rank(RR);nb:=rowdim(BB)-rank(BB); NR:=nullspace(TR);NB:=nullspace(TB); URAN:=stackmatrix(seq(NR[t],t=1..nr),seq(NB[s],s=1..nb)); NLBAS:=NullSpace(Matrix(URAN)); if nops(NLBAS)>1 then REP:=goRep([RR,BB],NLBAS) else return "The representation is one-dimensional" fi; goLie(REP[1],REP[2]) end:"goQuikLie calculates Lie algebra on the common range of R and B"; ### ### ExpVec # This expands a matrix in terms of a given matrix basis expVec:=proc(BB,v) global DM; DM:=doDIM(BB); expVecDim(BB,v,DM) end:" expVec(Basis,vector) returns the coefficients and assigns DM, dimension of the basis"; expVecDim:=proc(BB,v,dm) local i,vv,BX,vw,BY; global DM; DM:=dm; for i from 1 to dm do; vv[i]:=convert(evalm(BB[i]),vector); od; BX:=evalm(concat(seq(vv[ii],ii=1..dm))); vw:=convert(evalm(v),vector); BY:=evalm(evalm((transpose(BX)&*BX)^(-1))&*transpose(BX)); evalm(BY&*vw) end:" expVecDim(Basis,vector,dim of Basis) returns the coefficients and assigns DM, dimension of the basis"; ### GoRank goRank:=proc() local i,vxv; for i from 1 to nargs do; vxv[i]:=convert(Matrix(args[i]),vector); od; evalm(concat(seq(vxv[ii],ii=1..nargs))); rank(%) end: "goRank(matrix arguments) will give the linear rank of the set of arguments"; ### Kirlv kirlv:=proc() local i,j,n; global K; n:=rowdim(adjt[1]); K:=array(1..n,1..n,antisymmetric); for i from 1 to n-1 do: for j from 2 to n do: K[i,j]:=sum(x['k']*adjt[i]['k',j],'k'=1..n); od;od; evalm(K) end:"kirlv() computes the Kirillov matrix given adjt/after running Killing, goLie"; ### GoKilling goKilling:=proc() local i,KR,cc,j,k,Y,f,v,vv,i1; global Killing,adjt; Y:=evalm(sum(f[i1]*xi[i1],i1=1..ddm)): for i from 1 to ddm do: KR[i]:=expVecDim(xi,cmm(xi[i],Y),ddm);od: Killing:=array(1..ddm,1..ddm,symmetric); for i from 1 to ddm do:adjt[i]:=matrix(ddm,ddm,0);cc[i]:=vector(ddm,0);od: for i from 1 to ddm do: for j from 1 to ddm do: cc[j]:=map(diff,KR[i],f[j]); for k from 1 to ddm do: adjt[i][k,j]:=cc[j][k]; od;od;od; for i from 1 to ddm do: v[i]:=convert(adjt[i],vector); vv[i]:=convert(transpose(adjt[i]),vector);od; for i from 1 to ddm do: for j from 1 to i do: Killing[i,j]:=dotprod(v[i],vv[j]); od;od; trace(Killing), kernel(Killing), charpoly(Killing,lambda) end:"Killing computes the Killing form, its trace, nullspace, and characteristic polynomial"; ### GenBasis genBasis:=proc(n) local k,l,i,counter; global x,DZ,P; for k from 1 to n-1 do: x[k]:=matrix(n,n,(i,j)-> if j=1 then 1 elif j=k+1 then -1 else 0 fi); od; x[n]:=matrix(n,n,(i,j)-> if j=1 and i>1 then -1 elif i>1 and i=j then 1 else 0 fi); counter:=n; for k from 1 to n-1 do: for l from 1 to n-1 do: if k=l then next fi;counter:=counter+1; x[counter]:=matrix(n,n,(i,j)-> if i=k+1 and j=l+1 then 1 elif i=k+1 and j=1 then -1 else 0 fi); od;od; for k from 1 to n-2 do: counter:=counter+1; x[counter]:=matrix(n,n,(i,j)->if (i=k+1 and j=k+1) or (i=n and j=1) then 1 elif i=k+1 and j=1 then -1 elif i=j and i=n then -1 else 0 fi); od; DZ:=n^2-n; # for k from DZ+1 to n^2 do: # x[k]:=matrix(n,n,(i,j)-> if i=k-DZ and j=i then 1 else 0 fi); # od; DZ end:"genBasis computes the standard basis x of the ZRS(n) algebra of dimension DZ"; ### GenBasisAtState genBasisAtState:=proc(mm) local i,yy; global Q,m; m:=mm; Q:=matrix(n,n,(i,j)-> if i=j and i<>1 and i<>m then 1 elif i=1 and j=m then 1 elif i=m and j=1 then 1 else 0 fi); for i from 1 to n do: x||m[i]:=evalm(Q&*x[i]&*Q) od; for i from n+1 to DZ do: x||m[i]:=evalm(Q&*x[i]&*Q); od; DZ end:"genBasisAtState computes the standard basis x of the ZRS(n) algebra of dimension DZ based at state m"; ### GenPBasis genPBasis:=proc(n) local i; global P,DD,xp,y; P:=matrix(n,n,(i,j)-> if i=j then (n-1)/n else -1/n fi); genBasis(n); for i from n to DZ do: y[i+1-n]:=submatrix(x[i],2..n,2..n); od; for i from n to DZ do: xp[i+1-n]:=evalm(P&*x[i]); od; DD:=DZ+1-n; DD end:"genPBasis computes the projected standard basis xp of the ZRS(n) algebra of dimension DD"; ### ProjVec projVec:=proc(vv) local ww; ww:=expVec(xp,P&*vv); evalm(sum(ww['k']*y['k'],'k'=2..DD)) end:"projVec projects from PZRS to sl(n-1) with basis y"; ### GoRepRight goREPRight:=proc(ER,EE) local i,j,v,erdim,eedim; global REP; REP:='REP'; erdim:=doDIM(ER); eedim:=doDIM(EE); print(erdim," acting from the right on ",eedim); for i from 1 to erdim do: for j from 1 to eedim do: v[j]:=expVec(EE,evalm(EE[j]&*ER[i])); od; REP[i]:=map(simplify,evalm(concat(seq(v[p],p=1..eedim)))); print(`REP`[i]=evalm(REP[i])); od; [seq(REP[t],t=1..erdim)] end:"goREPRight computes the representation of arg1 on arg2 acting on the right and returns the set of matrices REP"; ### GoRep goREP:=proc(ER,EE) local i,j,v,erdim,eedim; global REP; REP:='REP'; erdim:=doDIM(ER); eedim:=doDIM(EE); print(erdim," acting on ",eedim); for i from 1 to erdim do: for j from 1 to eedim do: v[j]:=expVec(EE,evalm(ER[i]&*EE[j])); od; REP[i]:=map(simplify,evalm(concat(seq(v[p],p=1..eedim)))); print(`REP`[i]=evalm(REP[i])); od; [seq(REP[t],t=1..erdim)] end:"goREP computes the representation of arg1 on arg2 and returns the set of matrices REP"; ### DoKilling doKilling:=proc(BBB) local i,KR,cc,j,k,Y,f,v,vv; global Killing,adjt,eedim; eedim:=doDIM(BBB); print("Dimension is ",eedim); Y:=evalm(sum(f[i1]*BBB[i1],i1=1..eedim)): for i from 1 to eedim do: KR[i]:=expVec(BBB,cmm(BBB[i],Y));od: Killing:=array(1..eedim,1..eedim,symmetric); for i from 1 to eedim do:adjt[i]:=matrix(eedim,eedim,0);cc[i]:=vector(eedim,0);od: for i from 1 to eedim do: for j from 1 to eedim do: cc[j]:=map(diff,KR[i],f[j]); for k from 1 to eedim do: adjt[i][k,j]:=cc[j][k]; od;od;od; for i from 1 to eedim do: v[i]:=convert(adjt[i],vector); vv[i]:=convert(transpose(adjt[i]),vector);od; for i from 1 to eedim do: for j from 1 to i do: Killing[i,j]:=dotprod(v[i],vv[j]); od;od; trace(Killing), kernel(Killing), charpoly(Killing,lambda) end:"doKilling computes the Killing form for the input basis, its trace, nullspace, and characteristic polynomial"; ### DeComp Decomp:=proc(vv) local va,pva,vb,pvb; va:=expVec(GBAS,vv); pva:=evalm(sum(va['k']*GBAS['k'],'k'=1..DM)); vb:=expVec(IBAS,vv); pvb:=evalm(sum(vb['k']*IBAS['k'],'k'=1..DM)); evalm(pva),evalm(pvb) end:"Decomp gives the decomposition of a vector with respect to two subspaces: GBAS, IBAS"; ### Decomp ZRS DecompZRS:=proc(vv,mm) local vvv,pva,pvb; global gbas,ibas; vvv:=expVec(x||mm,vv); pva:=evalm(sum(vvv['k']*x||mm['k'],'k'=n+1..DZ)); pvb:=evalm(sum(vvv['k']*x||mm['k'],'k'=1..n)); evalm(pva),evalm(pvb) end:"DecompZRS gives the decomposition of a vector with respect to G(n) and I(n) in ZRS(n) at state m"; ### CommonNullspace commonNullspace:=proc(Vseq) local i,M; M:=stackmatrix(seq(evalm(Vseq[i]),i=1..doDIM(Vseq))); NullSpace(Matrix(M)) end:"commonNullspace returns the common nullspace of a sequence of matrices"; ### DODIM doDIM:=proc(Mseq) local edim; if nops(Mseq)>1 then return nops(Mseq) fi; # edim:=1; # while type(evalm(Mseq[edim]),matrix) do: edim:=edim+1;od; # edim-1 nops(op(Mseq)) end:"doDIM returns the length of the input sequence, dimension of the span"; ### ADJREP adjREPDim:=proc(ER,EE,dm) local i,j,v,erdim,eedim;global ADJREP; ADJREP:='ADJREP'; erdim:=dm; eedim:=doDIM(EE); print(erdim," acting on ",eedim); for i from 1 to erdim do: for j from 1 to eedim do: v[j]:=expVec(EE,cmm(ER[i],EE[j])); od; ADJREP[i]:=map(simplify,evalm(concat(seq(v[p],p=1..eedim)))); print(`ADJREP`[i]=evalm(ADJREP[i])); od; end:"adjREP computes the adjoint representation of arg1 on arg2 and returns the set of matrices ADJREP"; adjREP:=proc(ER,EE) local erdim; erdim:=doDIM(ER); adjREPDim(ER,EE,erdim) end: ### GenZFCBasis genZFCBasis:=proc(n) local i; global U,UU,z; U:=matrix(n,n,(i,j)-> if j=1 then 1 elif j>1 and i=j then -1 else 0 fi); genBasis(n); for i from 1 to DZ do: z[i]:=evalm(U&*x[i]&*U) od; DZ end:"genZFCBasis computes the ZRS(n) algebra of dimension DZ in the zero first column form"; ### Nullity nullity:=proc(X) rowdim(X)-rank(X) end:"Nullity returns the nullity of a matrix as a linear transformation."; ### Nullstack nullstack := proc() nullspace(stackmatrix(seq(args[t], t = 1 .. nargs))) end:"nullstack gives the nullspace of a set of matrices."; ### StdREP stdREP:=proc(RR,BB) local RCSP,BCSP,NRSP,NBSP; global RRBAS,NRBAS,RNBAS,NNBAS,FOURBAS,STDREP; RCSP:=colspace(RR); BCSP:=colspace(BB); NRSP:=nullspace(transpose(RR)); NBSP:=nullspace(transpose(BB)); RRBAS:=nullstack(seq(NRSP[t],t=1..nullity(RR)),seq(NBSP[s],s=1..nullity(BB))); NRBAS:=nullstack(seq(BCSP[t],t=1..rank(BB)),seq(NRSP[s],s=1..nullity(RR))); RNBAS:=nullstack(seq(NBSP[t],t=1..nullity(BB)),seq(RCSP[s],s=1..rank(RR))); NNBAS:=nullstack(seq(RCSP[t],t=1..rank(RR)),seq(BCSP[s],s=1..rank(BB))); FOURBAS:=[op(RRBAS),op(NRBAS),op(RNBAS),op(NNBAS)]; STDREP:=goRep([RR,BB],FOURBAS); STDREP[1], STDREP[2] end:"stdREP gives the matrix representation of R and B in terms of the common ranges and null perpendicular spaces."; ######## CNN CommonNullSpace ################ CNN:=proc() "Lie(R,B): ",commonNullspace([seq(xi[t],t=3..ddm)]) ,"Transposed L:", commonNullspace([seq(transpose(xi[t]),t=3..ddm)]) end:"CNN gives the common nullspace of the Basis of L from xi-3 on"; #################### GoLieE ##################### goLieE:=proc(AA,BB) local Q,idx,CC,MaxRnk,OldRnk,k,kk,new,Rnk;global xi,ddm; idx:=2; xi:='xi'; MaxRnk:=2; OldRnk:=0; ddm:=2; new:=true; xi[1]:=evalm(AA); xi[2]:=evalm(BB): CC[1]:=xi[1]; CC[2]:=xi[2]; ## print(xi[1],xi[2]); for kk from 2 while new or kk<=ddm do: new:=false; #print("COUNTER ONE ",kk); for k from 1 to ddm-1 do: idx:=idx+1; #print("COUNTER TWO ",k); CC[idx]:=ecliffe(cmm(xi[kk],xi[k])); if(CC[idx]=0) then idx:=idx-1; next fi; Rnk:=goRank(seq(CC[t],t=1..idx)); #print(idx," ",CC[idx]," AT ",kk," TRYING ",k," RANKS ",Rnk,MaxRnk," DIM ",ddm); if(Rnk>MaxRnk) then ddm:=ddm+1; xi[ddm]:=CC[idx]; MaxRnk:=Rnk; new:=true; print(ddm,"new basis element",rank(xi[ddm])); else next fi; Q:=map(simplify,evalm(xi[ddm])); xi[ddm]:=ecliffe(Q); od:od: "Lie algebra has dimension ", ddm end:"goLieE generates the Lie algebra generated by 2 matrices with e-factors"; ############################################# ######### CLIFF ROUTINES ####################### ############################################# cliffit:=proc() local k,T; T:=simplify(expand(args[1])); for k from 1 to n do T:=expand(algsubs(p[k]^2=p[k],T)) od; for k from 1 to n do T:=expand(algsubs(p[k]^2=p[k],T)) od; for k from 1 to n do T:=algsubs(p[k]^2=p[k],T) od; for k from 1 to n do T:=expand(algsubs(q[k]^2=q[k],T)) od; for k from 1 to n do T:=expand(algsubs(q[k]^2=q[k],T)) od; for k from 1 to n do T:=algsubs(q[k]^2=q[k],T) od; for k from 1 to n do T:=algsubs(p[k]*q[k]=0,expand(T)) od; for k from 1 to n do T:=algsubs(q[k]*p[k]=0,expand(T)) od; simplify(T) end:"cliffit reduces p-q expressions of an expression"; cliffe:=proc(X) local T; T:=map(cliffit,X); map(simplify,evalm(T)) end:"cliffe this clears the squares of the p's in a matrix"; cliffsubsmat:=proc(X,ff) local T; T:=map(cliffsubs,X,ff); map(simplify,evalm(T)) end:"cliffsubsmat substitutes for the p's in a matrix"; cliffsubs:=proc(X,ff) local k,T; T:=expand(X); for k from 1 to n do T:=algsubs(p[k]=ff[k],T) od; for k from 1 to n do T:=algsubs(q[k]=1-ff[k],T) od; simplify(expand(T)) end:"cliffsubs substitutes for the p's"; icliffit:=proc() local k,T; T:=simplify(expand(args[1])); for k from 1 to n do T:=expand(algsubs(e[k]^2=e[k],T)) od; for k from 1 to n do T:=expand(algsubs(e[k]^2=e[k],T)) od; for k from 1 to n do T:=algsubs(e[k]^2=e[k],T) od; simplify(T) end:"icliffit reduces 0-1 e's of an expression"; icliffe:=proc(X) local T; T:=map(cliffit,X); map(simplify,evalm(T)) end:"icliffe this clears the squares of the e's in a matrix"; icliffsubsmat:=proc(X,ff) local T; T:=map(cliffsubs,X,ff); map(simplify,evalm(T)) end:"icliffsubsmat substitutes for the e's in a matrix"; icliffsubs:=proc(X,ff) local k,T; T:=expand(X); for k from 1 to n do T:=algsubs(e[k]=ff[k],T) od; simplify(expand(T)) end:"icliffsubs substitutes for the e's"; ### ECLIFFs ################### ecliffit:=proc() local k,T; T:=expand(args[1]); for k from 1 to n do T:=(algsubs(e[k]^2=1,T)) od; for k from 1 to n do T:=algsubs(e[k]^2=1,T) od; simplify(T) end:"ecliffit clears e-squares of an expression"; ecliffe:=proc() local i,j,k,T; X:=args[1]; T:=map(ecliffit,X); evalm(T) end:"ecliffe this clears the squares of the e's in a matrix"; ecliffsubsmat:=proc(X,ff) local i,j,k,T; T:=map(ecliffsubs,X,ff); map(simplify,evalm(T)) end:"ecliffsubsmat substitutes for the e's in a matrix"; ecliffsubs:=proc(X,ff) local k,T; T:=expand(X); for k from 1 to n do T:=algsubs(e[k]=ff[k],T) od; simplify(expand(T)) end:"ecliffsubs substitutes for the e's"; ### ZCLIFFs ################### zcliffit:=proc() local k,T,nn; T:=expand(args[1]); if(nargs>1) then nn:=args[2] else nn:=n fi; for k from 1 to nn do T:=simplify(algsubs(f[k]^2=1,T)) od; for k from 1 to nn do T:=algsubs(f[k]^2=1,T) od; simplify(T) end:"zcliffit clears f-squares of an expression"; zcliffe:=proc() local i,j,k,T; X:=args[1]; if(nargs>1) then nn:=args[2] else nn:=n fi; T:=map(zcliffit,X,nn); evalm(T) end:"zcliffe this clears the squares of the f's in a matrix"; zcliffsubsmat:=proc(X,ff) local i,j,k,T; T:=map(zcliffsubs,X,ff); map(simplify,evalm(T)) end:"zcliffsubsmat substitutes for the f's in a matrix"; zcliffsubs:=proc(X,ff) local k,T; T:=expand(X); for k from 1 to n do T:=algsubs(f[k]=ff[k],T) od; simplify(expand(T)) end:"zcliffsubs substitutes for the f's"; ### SYMMETRIC MULTIPLICATION ######### symm:=proc(X) local Y,Y0; Y:=map(simplify,evalm(transpose(X)&*X)); Y0:=cliffe(evalm(Y)); evalm(Y0) end:"symm symmetrizes with transpose on the left"; sym:=proc(X) local Y,Y0; Y:=map(simplify,evalm(X&*transpose(X))); #Y0:=cliffe(evalm(Y)); evalm(Y) end:"symmm symmetrizes with transpose on the right"; symmult:=proc(X,Y) local Y1,Y0; Y1:=map(simplify,evalm(transpose(X)&*Y&*X)); #Y0:=cliffe(evalm(Y1)); evalm(Y1) end:"symmult does symmetric multiplication of X by Y with transpose on the left"; symult:=proc(X,Y) local Y1,Y0; Y1:=map(simplify,evalm(X&*Y&*transpose(X))); #Y0:=cliffe(evalm(Y1)); evalm(Y1) end:"symult does symmetric multiplication of X by Y with transpose on the right"; ######## NORMS ########### nrrm:=proc(v) local z; cliffit(evalm(sum(v[z]^2,z=1..n))) end:"nrrm gives the p-reduced sum of squares"; enrrm:=proc(v) local z; ecliffit(evalm(sum(v[z]^2,z=1..n))) end:"enrrm gives the e-reduced sum of squares"; nm:=proc(v) local z; simplify(evalm(sum(v[z]^2,z=1..doDIM(v)))) end:"nm gives the raw sum of squares"; #### SymTrace ################ symtrace:=proc(v) local d,i,T,theta; d:=doDIM(v); setparams(d); T:=matrix(d,d,0): for i from 1 to d do: T[i,v[i]]:=1;od: symtracemat(T) end:"symtrace of a vector calculates the generating function of the symmetric traces of the corresponding transformation"; symtracemat:=proc(T) local d,i,theta; d:=rowdim(T); setparams(d); theta:=det(J-t*T)^(-1); theta,taylor(theta,t=0,10) end:"symtracemat of a matrix calculates the generating function of the symmetric traces"; #### Transmat ################ transmat:=proc(v) local i,d,X; d:=doDIM(v); X:=matrix(d,d,0): for i from 1 to d do: j:=convert(v[i], decimal, hex): X[i,j]:=1;od: setparams(d); evalm(X) end:"transmat produces the matrix corresponding to the transformation v"; setparams:=proc(dm) global n,J,u,h; n:=dm; J:=array(1..n,1..n,identity); u:=vector(n,1); end:"setparams sets n to the dimension and J to the n-by-n identity."; matrans:=proc(X) local i,j,r; global trf;r:=rowdim(X);trf:=vector(r); for i from 1 to r do for j from 1 to r do if X[i,j]=1 then trf[i]:=j fi;od;od; h:=map(convert,evalm(trf) ,hex); [seq(h[k],k=1..r)] end:"matrans returns the transformation corresponding to a given matrix."; ######### ColSpan ################## colSpan:=proc() local i,vxv; for i from 1 to nargs do; vxv[i]:=Matrix(args[i]); od; evalm(concat(seq(vxv[ii],ii=1..nargs))); colspace(%) end: "colSpan(matrix arguments) will give the column space of the set of arguments"; ######### RowSpan ################## rowSpan:=proc() local i,vxv; for i from 1 to nargs do; vxv[i]:= Matrix(args[i]); od; evalm(stackmatrix(seq(vxv[ii],ii=1..nargs))); rowspace(%) end: "rowSpan(matrix arguments) will give the row space of the set of arguments"; ################## Centralizer ###################### Big:=proc(Y) local JJ,nn; nn:=rowdim(Y); JJ:=IdentityMatrix(nn); evalm(kron(Y,JJ)-kron(JJ,transpose(Y))) end:"Big"; Centr:=proc() local k,c,i,j,ZA,vxv,Y,nn; global Z,h,dc; h:='h'; Z:='Z'; nn:=rowdim(args[1]); for i from 1 to nargs do; vxv[i]:= Matrix(args[i]); od; Y:=evalm(stackmatrix(seq(Matrix(Big(vxv[ii])),ii=1..nargs))); ZA:=NullSpace(Matrix(Y)); for k from 1 to doDIM(ZA) do c:=1; Z[k]:=matrix(nn,nn); for i from 1 to nn do for j from 1 to nn do Z[k][i,j]:=ZA[k][c];c:=c+1; od; od; od; dc:=doDIM(ZA); evalm(sum(h[t]*evalm(Z[t]),t=1..dc)) end:"Centr(matrix) gives a basis for the centralizer of the argument"; ### Kronecker Product kron:= proc(AA,BB) local i,j,CX,UX,ii,jj,n,m; n:=rowdim(AA);m:=coldim(AA); for i to n do: for j to m do: CX[j]:=scalarmul(BB,AA[i,j]) od; UX[i]:=concat(seq(CX[jj],jj=1..m)) od; stackmatrix(seq(UX[ii],ii=1..n)) end:"kron(matrix,matrix) returns the Kronecker product of the matrices"; proj:=proc() local BX; BX:=evalm(concat(seq(args[ii],ii=1..nargs))); map(simplify,evalm(evalm(BX&*(htranspose(BX)&*BX)^(-1))&*htranspose(BX))); end:" proj(seq of vectors) returns the orthogonal projection onto their span"; goProj:=proc() local i;global ASEV,AKEV,DSEV,DKEV,P,Q,S,T,alpha,beta,lambda,mu; ASEV:=[eigenvectors(AS)];AKEV:=[eigenvectors(AK)]; DSEV:=[eigenvectors(DS)];DKEV:=[eigenvectors(DK)]; for i from 1 to doDIM(ASEV) do P[i]:=proj(seq(ASEV[i][3][t],t=1..ASEV[i][2]));od; for i from 1 to doDIM(AKEV) do Q[i]:=proj(seq(AKEV[i][3][t],t=1..AKEV[i][2]));od; for i from 1 to doDIM(DSEV) do S[i]:=proj(seq(DSEV[i][3][t],t=1..DSEV[i][2]));od; for i from 1 to doDIM(DKEV) do T[i]:=proj(seq(DKEV[i][3][t],t=1..DKEV[i][2]));od; alpha:='alpha';beta:='beta';lambda:='lambda';mu:='mu'; map(simplify,evalm(add(lambda[t]*P[t],t=1..doDIM(ASEV)))), map(simplify,evalm(add(mu[t]*Q[t],t=1..doDIM(AKEV)))), map(simplify,evalm(add(alpha[t]*S[t],t=1..doDIM(DSEV)))), map(simplify,evalm(add(beta[t]*T[t],t=1..doDIM(DKEV)))) end:"goProj() calculates the projection operators for AS, AK, DS, DK, namely P,Q,S,T"; ### EXTERIOR POWERS asym2:=proc(X) local d,nn,X2,i,j,k,l,aa,bb; d:=rowdim(X); nn:=binomial(d,2): X2:=matrix(nn,nn); aa:=0;bb:=0; for i from 1 to d-1 do for j from i+1 to d do aa:=aa+1; bb:=0; for k from 1 to d-1 do for l from k+1 to d do bb:=bb+1; X2[aa,bb]:=abs(det(submatrix(X,[i,j],[k,l]))); od; od; od; od; map(simplify,evalm(X2)) end:"asym2 returns the 2nd exterior power of a matrix with positive entries"; extpow:=proc(X,p) local nn,d,X2,CC,i,j; nn:=rowdim(X); d:=binomial(nn,p); X2:=matrix(d,d); CC:=choose(nn,p); for i from 1 to d do for j from 1 to d do X2[i,j]:=det(submatrix(X,CC[i],CC[j])) od; od; map(simplify,evalm(X2)) end:"extpow(matrix,power) returns the p-th exterior power of the matrix"; extpowdiag:=proc(X,p) local nn,d,X2,CC,i; nn:=rowdim(X); d:=binomial(nn,p); X2:=vector(d); CC:=choose(nn,p); for i from 1 to d do X2[i]:=det(submatrix(X,CC[i],CC[i])) od; map(simplify,evalm(X2)) end:"extpow(matrix,power) returns the diagonal of the p-th exterior power of the matrix"; extpowf:=proc(X,p) local nn,d,X2,CC,i,j; nn:=rowdim(X); d:=binomial(nn,p); X2:=matrix(d,d); CC:=choose(nn,p); for i from 1 to d do for j from 1 to d do X2[i,j]:=abs(det(submatrix(X,CC[i],CC[j]))) od; od; map(simplify,evalm(X2)) end:"extpow(matrix,power) returns the p-th exterior power of the matrix absolute val of entries"; extbas:=proc(m) local d,k,CC; global II,eta; eta:='eta'; d:=binomial(m,2); II:=IdentityMatrix(m); CC:=choose(m,2); for k from 1 to d do eta[k]:=wej(col(II,CC[k][1]),col(II,CC[k][2])) od; "basis of size", k-1, " generated" end:"extbas(int) generates eta, the second exterior power of the standard basis"; extpowshow:=proc(X,nn) local n,p,EX; n:=rowdim(X): for p to nn do EX:=extpowdiag(X,p): print(p," ",seq(EX[i],i=1..binomial(n,p))); od: end:"extpowshow prints the diagonal elements of the exterior powers of a matrix"; extlie:=proc(X) local X2,i,j,nn,d,k,ytemp,y; nn:=rowdim(X); d:=binomial(nn,2); X2:=matrix(d,d): k:=0; for i from 1 to nn-1 do for j from i+1 to nn do k:=k+1; ytemp:=evalm(wej(col(X,i),col(II,j))+wej(col(II,i),col(X,j))); y[k]:=expVec(eta,ytemp); od; od; X2:=stackmatrix(seq(y[t],t=1..d)); map(simplify,evalm(-X2)) end:"extlie(matrix) returns the representation of a Lie element on the 2nd exterior power"; wej:=proc(a,b) local ma,mb; ma:=Matrix(a); mb:=Matrix(b); map(simplify,evalm(kron(ma,mb)-kron(mb,ma))) end:"wej(vec,vec) returns the wedge product of two vectors"; ### Checkspan checkspan:=proc(BB,v) local EV,t; global ev; ev:=expVec(BB,v); EV:=evalm(add(ev[t]*BB[t],t=1..doDIM(BB))); evalm(EV),iszero(evalm(EV-v)) end:"checkspan(basis,v) checks if v is in the span of the given basis"; inv:=x->1/x:"inv defined"; ### RUNADJ runAdj:=proc(YY) local J,i,aj,Y; J:=IdentityMatrix(rowdim(YY)); aj:=adj(J-YY); Y[0]:=subs(tau=0,evalm(aj)): for i from 1 to n-1 do Y[i]:=map(simplify,(evalm(subs(tau=0,map(diff,aj,tau$i))/i!))) od: stackmatrix(seq((row(Y[b],1),b=0..d))) end:"runAdj(matrix) returns the adjoint sorted by homogeneous degree "; ### MakeIdems makeIdems:=proc(xR1,xR2,xP1,xP2) local s,t,u,i,j,S,k,l ;global r,E; r:=doDIM(xR1); s:=0: t:=0: for i in [xR1,xR2] do s:=s+1; if s=3 then s:=1 fi; for j in [xP1,xP2] do: t:=t+1; if t=3 then t:=1 fi; S:={seq(i,i=1..n)}; E[t,s]:=matrix(n,n,0): for l to r do: for k to r do if member(i[l],j[k] intersect S) then for u in j[k] do E[t,s][u,i[l]]:=1; od; S:=S minus {i[l]}; fi; od; od; od; od; print(map(matrans,[E[1,1],E[1,2],E[2,1],E[2,2]])) end:"makeidems makes idempotents for a 2-by-2 box of partitions and ranges"; ### MakeIdem makeIdem:=proc(xR1,xP1) local u,i,j,S,k,l,r,EX,nx,nn ; xx:=[]; for i to nops(xP1) do if(xP1[i] <>{}) then xx:=[op(xx),xP1[i]] fi; od; nx:=`union`(seq(xx[k],k=1..nops(xx))); nn:=nops(nx); EX:=matrix(nn,nn,0): r:=nops(xR1); for i in [xR1] do for j in [xx] do: S:={seq(i,i=1..nn)}; for l to r do: for k to r do if member(i[l],j[k] intersect S) then for u in j[k] do EX[u,i[l]]:=1; od; S:=S minus {i[l]}; fi; od; od; od; od; ## print(matrans(EX)); evalm(EX) end:"makeidem(R,P) makes an idempotent for a partition and a range"; makeFunction:=proc(xR,xP) local u,i,j,S,k,l,r,EX,nx,nn ; nx:=`union`(seq(xP[k],k=1..nops(xP))); nn:=nops(nx); EX:=matrix(nn,nn,0): for i to nops(xP) do for k to nops(xP[i]) do EX[xP[i][k],xR[i]]:=1 od; od; evalm(EX) end:"makeFunction makes a function mapping a given partition into a given range"; ### LieProduct lieproduct:=proc(v,w) local m,mm,dm; dm:=doDIM(LBAS); m:=evalm(sum(v[k]*LBAS[k],k=1..dm)); mm:=evalm(sum(w[k]*LBAS[k],k=1..dm)); expVec(LBAS,cmm(m,mm)) end:"lieproduct of two vectors returns their product with respect to the basis xi"; ### LieElement lieelement:=proc(v) local dm; dm:=doDIM(LBAS); evalm(sum(v[k]*LBAS[k],k=1..dm)) end:"lieelement returns the matrix of the corresponding element of the Lie algebra"; ### Singular Value Decomposition svd:=proc(X) local SV;global H,K,UU,VT,DD; SV:=SingularValues(Matrix(X),output=['U','S','Vt']); DD:=diag(seq(SV[2][j],j=1..rowdim(X))): UU:=SV[1]; VT:=SV[3]; H:=multiply(UU,VT); K:=multiply(transpose(VT),DD,VT); evalm(H),evalm(K) end:"svd returns the polar decomposition of the matrix"; ####### FLIPPI flip:=proc(X) global F,F2; map(simplify,multiply(F2,X,F)) end:"flip returns X in lightcone basis"; ###### ABEL LIMIT abel:=proc(X) local J; J:=IdentityMatrix(rowdim(X)); return map(limit,evalm((1-s)*inverse(J-s*X)),s=1) end:"returns the Abel limit of powers of the matrix X"; ########### SYMMETRIC POWERS sympow:=proc(X,p) local nn,d,X2,CC,i,j; nn:=rowdim(X); d:=binomial(nn,p); X2:=matrix(d,d); CC:=choose(nn,p); for i from 1 to d do for j from 1 to d do X2[i,j]:=permanent(submatrix(X,CC[i],CC[j])) od; od; map(simplify,evalm(X2)) end:"sympow(matrix,power) returns the p-th symmetric power no-diagonal of the matrix"; ############## FullySymmetric Square sym2:=proc(MM) local nn,NN,cc,a1,a2,b1,b2,y,t,x,i,j; nn:=rowdim(MM); NN:=binomial(nn,2); x:=vector(nn); y:=multiply(MM,x); t:=matrix(NN+nn,NN+nn): cc:=choose(nn,2): for i to NN do for j to NN do a1:=cc[i][1];a2:=cc[i][2]; b1:=cc[j][1];b2:=cc[j][2]; t[i,j]:=diff(y[a1]*y[a2],[x[b1],x[b2]]); od; od; for i to NN do for j from NN+1 to NN+nn do a1:=cc[i][1];a2:=cc[i][2]; t[i,j]:=1/2*diff(y[a1]*y[a2],[x[j-NN],x[j-NN]]); od; od; for i from NN+1 to NN+nn do for j to NN do b1:=cc[j][1];b2:=cc[j][2]; t[i,j]:=diff(y[i-NN]*y[i-NN],[x[b1],x[b2]]); od; od; for i from NN+1 to NN+nn do for j from NN+1 to NN+nn do t[i,j]:=1/2*diff(y[i-NN]*y[i-NN],[x[j-NN],x[j-NN]]) od; od; evalm(t) end:"sym2(matrix) returns the 2nd fully symmetric power of the matrix"; ########### MATVEC matvec:=proc(g) local nn,ii,w,kk,i,j; nn:=rowdim(convert(g,matrix)); ii:=(1+sqrt(1+8*nn))/2; w:=evalm(ZeroMatrix(ii)); kk:=1: for i from 1 to ii-1 do for j from i+1 to ii do if(i<>j) then w:=evalm(w+g[kk]*sigma[i,j]) fi ; kk:=kk+1 od ;od; evalm(w): end:"matvec converts a vector of length n choose 2 into a symmetric matrix with zero diagonal"; mexvec:=proc(g) local nn,ii,w,kk,i,j; nn:=rowdim(convert(g,matrix)); ii:=(1+sqrt(1+8*nn))/2; w:=evalm(ZeroMatrix(ii)); kk:=1: for i from 1 to ii-1 do for j from i+1 to ii do if(i<>j) then w:=evalm(w+g[kk]*rho[i,j]) fi ; kk:=kk+1 od ;od; evalm(w): end:"mexvec converts a vector of length n choose 2 into a skew-symmetric matrix"; ########### READROWS readrows:=proc(M) local S0,nn,s0,k,i,mm; S0:=[]: nn:=rowdim(M): mm:=coldim(M): for k to nn do s0:=[]; for i to mm do if(M[k,i]<>0) then s0:=[op(s0),i] fi; od; S0:=[op(S0),s0]; od: S0: end:"reads rows and makes sets according to which columns are nonzero with rows marked"; ########### READCYCLES readcycles:=proc(M) local S0,nn,s0,k,i,mm; S0:={}:s0:={}: nn:=rowdim(M): mm:=coldim(M): for k to nn do s0:={}; for i to mm do if(M[k,i]<>0) then s0:= s0 union {i} fi; od; S0:=S0 union {s0}; od: S0: end:"reads rows and makes sets according to which columns are nonzero"; ########### READ2 read2:=proc(M) local S0,nn,s0,k,i,mm,MM,sp,sn,CK,S1,sp1,sn1,s1; S0:=[]; S1:=[]; CK:=choose(n,2); MM:=tend2(M); nn:=rowdim(MM): mm:=coldim(MM): for k to nn-1 do sp:=[]; sn:=[]; sp1:=[]; sn1:=[]; for i to mm-1 do if(MM[k,i]<0) then sn:=[op(sn),i]; if(MM[k,i]=-1/2) then sn:=[op(sn),i] fi; sn1:=[op(sn1),convert(CK[i],set)]; if(MM[k,i]=-1/2) then sn1:=[op(sn1),convert(CK[i],set)] fi fi; if(MM[k,i]>0) then sp:=[op(sp),i] ; if(MM[k,i]=1/2) then sp:=[op(sp),i] fi; sp1:=[op(sp1),convert(CK[i],set)] ; if(MM[k,i]=1/2) then sp1:=[op(sp1),convert(CK[i],set)] fi fi; od; if(MM[k,mm]<0) then sn:=[op(sn),"x"] ; if(MM[k,mm]=-1/2) then sn:=[op(sn),"x"] fi; sn1:=[op(sn1),"x"] ; if(MM[k,mm]=-1/2) then sn1:=[op(sn1),"x"] fi fi; if(MM[k,mm]>0) then sp:=[op(sp),"x"] ; if(MM[k,mm]=1/2) then sp:=[op(sp),"x"] fi; sp1:=[op(sp1),"x"] ; if(MM[k,mm]=1/2) then sp1:=[op(sp1),"x"] fi fi; s0:=[[op(sp)],[op(sn)]]; s1:=[[op(sp1)],[op(sn1)]]; S0:=[op(S0),s0]; S1:=[op(S1),s1]; od: S0,S1; end:"reads rows and makes sets with plus and minus for Delta1"; tend2:=proc(X) local nn,np,bot,tob,X1,i,aa; nn:=rowdim(X); np:=nn+1; bot:=matrix(1,nn,0); tob:=matrix(np,1,0); X1:=concat(stackmatrix(X,bot),tob); for i to np do aa:=add(X1[i,k],k=1..np); X1[i,np]:=-aa; od; evalm(X1) end:"extends matrix by one row and column to rowsum zero"; ########### OMEGAK omegaK:=proc(A,B,k) local nn,n4,J4,A4,pi4,u4; nn:=rowdim(A); n4:=binomial(nn,k): A4:=evalm((1/2)*(sympow(A,k)+sympow(B,k))): J4:=evalm(IdentityMatrix(n4)): pi4:=evalm((1/x[1])*linsolve(J4-transpose(A4),vector(n4,0),'r',x)): u4:=evalm((1/x[1])*linsolve(J4-A4,vector(n4,0),'r',x)): evalm(pi4),evalm(u4),readrows(abel(sympow(A,k))),readrows(abel(sympow(B,k))) end:"prints omega level k giving pi-k and u-k separately with R and B input"; ########### piV piV:=proc(v,k) local NL,EI,i,j,va,vb,LEVEL; NL:=rowdim(convert(v,matrix)): LEVEL:=binomial(n,k-1); EI:=matrix(NL,LEVEL,0): for i to NL do for j to LEVEL do va:=convert(choose(n,k)[i],set); vb:=convert(choose(n,k-1)[j],set): if(vb subset va) then EI[i,j]:=1 fi od: od: evalm(v&*EI) end:"produces the one-level lowering operator for the pi-s"; ########### uV uV:=proc(v,k) local NL,EI,i,j,va,vb,LEVEL; NL:=rowdim(convert(v,matrix)): LEVEL:=binomial(n,k-1); EI:=matrix(LEVEL,NL,0): for i to LEVEL do for j to NL do va:=convert(choose(n,k-1)[i],set); vb:=convert(choose(n,k)[j],set): if(va subset vb) then EI[i,j]:=pi[op(1,vb minus va)] fi od: od: multiply(EI,v) end:"produces the one-level lowering operator for the u-s"; ########### REPLACE ROW/COL colreplace:=proc(X,i,v) local Y; Y:=concat(X,v); DeleteColumn(ColumnOperation(Matrix(Y),[i,-1]),-1) end:"replaces column i of X by v"; rowreplace:=proc(X,i,v) local Y; Y:=stackmatrix(X,v); DeleteRow(RowOperation(Matrix(Y),[i,-1]),-1) end:"replaces row i of X by v"; ####### GO DET EVAL MIXED DET goDET:=proc(A,B,C) local vx,i,nu; global cmat; nu:=rowdim(B): vx:=array(1..nu): for i to nu do vx[i]:=row(evalm(B),i): if(evalb(i in A)) then vx[i]:=row(evalm(C),i) fi; od; cmat:=stackmatrix(seq(vx[i],i=1..nu)); det(cmat) end:"computes the mixed determinant of arg2, such as I-A2 and arg3 such as Delta2"; ######## GO CORE CORE:=proc(g,B,C,b) local CR,idx,j,k,ng,nu; nu:=rowdim(B): CR:=array(1..nops(g)); for j to nops(g) do idx:={}; ng:=op(j,g): for k to nu do if(divide(ng,x[k])) then idx:=idx union {k} fi; od; CR[j]:=goDET(idx,evalm(B),evalm(C)); if(b) then print(idx,CR[j]) fi od: add(CR[k],k=1..nops(G)),evalm(CR) end:"computes the mixed determinants comprising the core of the system"; xtend:=proc(X) local nn,np,bot,tob,X1,i,aa; nn:=rowdim(X); np:=nn+1; bot:=matrix(1,nn,0); tob:=matrix(np,1,0); X1:=concat(stackmatrix(X,bot),tob); for i to np do aa:=add(X1[i,k],k=1..np); X1[i,np]:=1-aa; od; evalm(X1) end:"extends matrix by one row and column and makes it stochastic"; ### READVEC readVec:=proc(vv) op(1,readcycles(transpose(convert(vv,matrix)))) end:"outputs support of a vector"; ### NULLVEC nullVec:=proc(vv) local idx,rd; rd:=rowdim(convert(vv,matrix)): idx:={}; for i to rd do if(vv[i]=0) then idx:={op(idx),i} fi od; idx end:"outputs support of a vector"; ### CLEAR DENOMS clearDenoms:=proc(vv) local dd,i,dn; dd:=rowdim(convert(vv,matrix)): for i to dd do dn[i]:=denom(vv[i]) od: evalm(lcm(seq(dn[i],i=1..dd))*vv) end:"clears denominators of a vector with rational entries"; F:=proc(X) evalm(symult(A,X)+symult(Delta,X)) end:"symmetric mapping on the right"; G:=proc(X) evalm(symmult(A,X)+symmult(Delta,X)) end:"symmetric mapping on the left"; res:=proc(X,t) local jj; jj:=evalm(IdentityMatrix(rowdim(X))); inverse(t*jj-X) end:"returns the resolvent of X as a function of t"; ### GO RANGE goRange:=proc(X) local nn,nd,nul,nsn,NSN,VRN,i,j,idx,vv; nn:=coldim(X); nd:=rank(X): nul:=nn-nd; NSN:=matrix(1,nn,0): if(nullity(X)>0) then nsn:=nullspace(X): NSN:=stackmatrix(seq(nsn[k],k=1..nul)): fi: VRN:=linsolve(NSN,vector(max(1,nul),0),'rr',y); evalm(VRN) end:"returns the range of X and the corresponding partitions of the parametrized range"; goPartitions:=proc(v,nd) local i,j,idx,vv,np; np:=rowdim(convert(v,matrix)); for i to nd do idx||i:={}; vv:=map(diff,v,y[i]); for j to np do if(vv[j]=1) then idx||i:={op(idx||i),j} elif (vv[j]=-1) then idx||i:={op(idx||i),[j]} fi od: od: seq(idx||k,k=1..nd) end:"returns the partitions of a parameter vector"; ### GO PARTS goParts := proc(XY) local yvars, GR, JAC, nrows, ncols, ix, i, j, k, T, ct, q, HR, po, ne, H, GPN; yvars := [seq(y[k], k = 1 .. (linalg:-rank)(XY))]; GR := goRange(XY); JAC := evalm((linalg:-jacobian)(GR, yvars)); unassign('ps', 'ng', 'ix', 'PT'); nrows := (linalg:-rowdim)(JAC); ncols := (linalg:-coldim)(JAC); ix := {}; for i to nrows do for j to ncols do if JAC[i, j] = -1 then ix := {op(ix), i}; ps || i := []; ng || i := []; end if; end do; end do; for i in ix do for j to ncols do if JAC[i, j] = 1 then ps || i := [op(ps || i), j] end if; if JAC[i, j] = -1 then ng || i := [op(ng || i), j] end if; end do; bx || i := []; for k in (combinat:-permute)(ng || i) do for j in (combinat:-choose)(ps || i, nops(ng || i)) do bx || i := [op(bx || i), [j, k]] end do; end do; end do; TX:= (combinat:-cartprod)([seq(bx || k, k in ix)]); unassign('H'); ct := 0; while not TX[finished] do q:=TX[nextvalue](); HR := GR; for i to nops(ix) do po := q[i][1]; ne := q[i][2]; HR := subs([seq(y[ne[k]] = y[po[k]], k = 1 .. nops(ne))], evalm(HR)); end do; ct := ct + 1; H[ct] := evalm(HR); end do; unassign('GPN'); for i to ct do GPN[i] := {goPartitions(H[i], (linalg:-rank)(XY))} end do; seq(GPN[k], k = 1 .. ct); end proc:"goParts finds all compatible partitions based on N"; getParts:=proc(GP,nx,rnk,rg,RA,BA,YN) J:=evalm(IdentityMatrix(nx)); UJ:=matrix(nx,nx,1); NNY:=evalm(UJ-YN); PN:=`union`(seq({GP[k]},k=1..nops(GP))); uk:=vector(binomial(nx,rnk),1): for i to nops(PN) do tmp:=makeIdem(rg,PN[i]); EX||i:=symult(tmp,J); od: if(nops(PN)>0) then IRANK:=goRank(seq(evalm(EX||i),i=1..nops(PN))): fi; if(IRANK=1) then return [evalm(EX1)],[1/1],[],[] fi; for k in choose(nops(PN),IRANK) do UBASIS:=[seq(EX||i,i in k)]: alph:=vector(IRANK,-1); if ( goRank( seq(evalm(UBASIS[i]),i=1..IRANK) )=IRANK ) then if(checkspan(UBASIS,NNY)[2]) then alph:=expVec(UBASIS,NNY); stp:=true; else next fi; for j to IRANK do if(alph[j]<0) then stp:=false; break; fi; od; if(not stp) then next fi; if(stp) then break; fi; fi; od; PB:=[]: for i to nops(UBASIS) do if(alph[i]>0) then PB:=[op(PB),UBASIS[i]]; fi; od; DIFF:=1: PT:=`union`(seq({readcycles(PB[k])},k=1..nops(PB))); while (DIFF> 0) do q:=nops(PT): for i to q do PT:=PT union {readcycles(symult(RA,PB[i]))}; PT:=PT union {readcycles(symult(BA,PB[i]))}; od; DIFF:=nops(PT)-q; PB:=[]: for i to nops(PT) do tmp:=makeIdem(rg,PT[i]); PB:=[op(PB),symult(tmp,J)]; od: od; np:=nops(PB); APARTS:=matrix(np,np,0): RPARTS:=matrix(np,np,0): BPARTS:=matrix(np,np,0): for i to np do aa:=symult(RA,PB[i]): bb:=symult(BA,PB[i]): for j to np do if(iszero(evalm(PB[j])-aa)) then RPARTS[i,j]:=RPARTS[i,j]+1/2 fi: if( iszero(PB[j]-bb)) then BPARTS[i,j]:=BPARTS[i,j]+1/2 fi: od: od: APARTS:=evalm(RPARTS+BPARTS); jj:=evalm(IdentityMatrix(np)): nh:=nullspace(jj-transpose(APARTS)); alph:=nh[1]; alph:=evalm(alph/add(alph[k],k=1..np)); if(nops(PB)>1) then return PB, evalm(alph), readrows(evalm(RPARTS)),readrows(evalm(BPARTS)) fi; PB, evalm(alph), [], [] end:"returns the partitions, alpha and the action of R and B on the partitions";