RCOUT.mw

 > #####################################################################################

 > #################### !!!!!!!!!!!!!!!!!! START RUNS HERE !!!!!!!!!!!!! ######################################

 > #####################################################################################

 > print(matrans(R),matrans(B),submatrix(Q,1..n-d,1..n),pi);

 > "rank Delta"=rank(Delta);"rank A"=rank(A);"min-poly",factor(minpoly(Delta,x));unassign('phi','ee');print("pi",pi);

 > ########################## RANK RUNS ###############################

 > d:=rank(Delta);n:=rowdim(A):zeta:=vector(n,0):ND:=nullspace(Delta);P:=concat(seq(ND[k],k=1..nops(ND)),seq(zeta,j=1..d)):nullspace(J-P);

 > ee:=[-1,1,1,-1,1,1,-1,1,1,1]:phi:=diag(seq(ee[q],q=1..n)):Detta:=phi&*Delta:R:=evalm(A+Detta):B:=evalm(A-Detta):matrans(R),matrans(B);factor(minpoly(R,s)),factor(minpoly(B,s));nullspace(transpose(Detta));minpoly(Detta,x),factor(minpoly(Detta,x));

 > rterm:={}:v:=evalm(multiply(e,R)):for i to n do if (v[i]=0) then rterm:=rterm union {i} fi od:bterm:={}:v:=evalm(multiply(e,B)):for i to n do if (v[i]=0) then bterm:=bterm union {i} fi od:print("terminal points of R",rterm,"terminal points of B",bterm);

 > RCYC:=map(limit,evalm((1-s)*inverse(J-s*R)),s=1):BCYC:=map(limit,evalm((1-s)*inverse(J-s*B)),s=1):

 > unassign('S0','T0','s0','t0');RC:={}:S0:={}:s0:={}:k:=0:a:=1:for IX to n do for i to n do if(RCYC[IX,i]<>0) then k:=k+1;s0:={i}union s0 fi; if( RCYC[IX,i]<0 and a=1 ) then a:=-1 fi; od;S0:=S0 union {s0};RC:=RC union s0;s0:={}: od:BC:={}:j:=0:b:=1:T0:={}:t0:={}:for IX to n do for i to n do if(BCYC[IX,i]<>0) then j:=j+1;t0:={i} union t0 fi; if(BCYC[IX,i]<0 and b=1) then b:=-1 fi od;T0:=T0 union {t0};BC:=BC union t0;t0:={} od:rcyc:=a*k*RCYC:bcyc:=evalm(b*j*BCYC):print("cycles in R",S0,"cycles in B",T0);

 > for i to n do if iszero(evalm((R^2)^i-R^i)) then print ("order of R is",i);break fi od;N:=i:

 > vv:=multiply(Delta,R^N,e):P0:={}:for i to n do if vv[i]=0 then P0:=P0 union {i} fi od: print("periodic points",P0);print("PN points",P0 minus RC);

 > for i to n do if iszero(evalm((B^2)^i-B^i)) then print ("order of B is",i);break fi od;N:=i:

 > vv:=multiply(Delta,B^N,e):Q0:={}:for i to n do if vv[i]=0 then Q0:=Q0 union {i} fi od: print("periodic points",Q0);print("PN points",Q0 minus BC);

 > antR:=(rterm minus RC) minus P0;antB:=(bterm minus BC) minus Q0;

 > Y0[0]:=Omega:for i from 1 to d do Y0[i]:=map(simplify,multiply(Y0[i-1],Detta)) od:K0:=stackmatrix(seq(row(Y0[b],1),b=1..d)): print(rank(K0),K0);

 > NK:=nullspace(transpose(K0));dd:=d+1:

 > for k from 1 to n  do phi[k,k]:=-ee[k]; Y0[0]:=Omega:for i from 1 to d do Y0[i]:=map(simplify,multiply(Y0[i-1],Detta)) od:KA[k]:=stackmatrix(seq(row(Y0[b],1),b=1..d)); print(k,rank(KA[k]));phi[k,k]:=ee[k]; od:

 > if ( rank(K0)=d) then mm:=minpoly(Detta,s):NKN:=CoefficientList(mm,s);dd:=degree(mm)+1 ;