RCOUT.html

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

> #################### !!!!!!!!!!!!!!!!!! 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 ;

> p:=add(NKN[j]*s^(j-1),j=1..dd):q:=add(abs(NKN[j])*t^(j-1),j=1..dd):print(p,q);

> PP:=evalm(subs(s=Detta,p)):QQ:=evalm(subs(t=A,q)):SA:=evalm(PP+QQ):rowsum:=multiply(SA,vector(n,1)):ST:=evalm(SA/rowsum[1]):SB:=evalm(-PP+QQ):rowsum:=multiply(SB,vector(n,1)):SU:=map(simplify,evalm(SB/rowsum[1])):print(PP,QQ,ST,SU,"cmm check",cmm(ST,Omega),cmm(SU,Omega));print( map(limit,evalm((1-s)*inverse(evalm(J-s*ST))),s=1),map(limit,evalm((1-s)*inverse(evalm(J-s*SU))),s=1)); fi;

> if(rank(K0)<d) then for i to nullity(K0) do NKN:=[0,seq(NK[i][k],k=1..d)];

> p:=add(NKN[j]*s^(j-1),j=1..dd):q:=add(abs(NKN[j])*t^(j-1),j=1..dd):print(p,q);

> PP:=evalm(subs(s=Detta,p)):QQ:=evalm(subs(t=A,q)):SA:=evalm(PP+QQ):rowsum:=multiply(SA,vector(n,1)):ST:=map(simplify,evalm(SA/rowsum[1])):SB:=evalm(-PP+QQ):rowsum:=multiply(SB,vector(n,1)):SU:=map(simplify,evalm(SB/rowsum[1])):print(PP,QQ,ST,SU,"cmm check",cmm(ST,Omega),cmm(SU,Omega));print( map(limit,evalm((1-s)*inverse(evalm(J-s*ST))),s=1),map(limit,evalm((1-s)*inverse(evalm(J-s*SU))),s=1)); od;fi;