RCOUT.mw

> ############ START HERE ######################

> Delta:=evalm((R-B)/2):A:=evalm((R+B)/2): `R`=matrans(R);`B`=matrans(B);

> unassign('e','a','w'):pi:=evalm(1/w[1]*linsolve(J-transpose(A),vector(n,0),'r',w));d:=rank(Delta):`rank of Delta`=d,`for n equals`=n;`rank of A`=rank(A);NA:=NullSpace(Matrix(A)):NTA:=NullSpace(Matrix(transpose(A))):DA:=evalm(Delta&*A):(Delta,DA):`RankCheck1`=is(0=1-n+rank(stackmatrix(Delta,DA)));nu2:=nullity(evalm(Delta^2));nu:=nullity(Delta);`RankCheck2`=is((nu2-nu)<2);zeta:=vector(n,0):u=vector(n,1):print(evalm(R-B));

R = vector([4, 3, 1, 2])

B = vector([3, 4, 4, 3])

pi := vector([1, 1, 2, 2])

`rank of Delta` = 3, `for n equals` = 4

`rank of A` = 3

RankCheck1 = true

nu2 := 1

nu := 1

RankCheck2 = true

matrix([[0, 0, -1, 1], [0, 0, 1, -1], [1, 0, 0, -1], [0, 1, -1, 0]])

> Eigenvectors(Matrix(Delta)):TT:=map(fnormal,map(evalf,%[2]),4);factor(minpoly(Delta,x));map(abs,map(fnormal,multiply(pi,TT),3));

TT := Matrix([[0., 1., .5000-1.323*I, .5000+1.323*I], [0., 1., -.5000+1.323*I, -.5000-1.323*I], [1., 1., -1., -1.], [1., 1., 1., 1.]])

1/4*x*(2*x+1)*(2*x^2-x+1)

vector([4., 6., 0., 0.])

> FF:=FrobeniusForm(Matrix(Delta),output=['F','Q']);QD:=FF[2]:DQD:=FF[1]:

FF := Matrix([[0, 0, 0, 0], [1, 0, 0, (-1)/4], [0, 1, 0, (-1)/4], [0, 0, 1, 0]]), Matrix([[1, 0, (-1)/4, (-1)/8], [0, 0, 1/4, 1/8], [0, 1/2, 0, 0], [0, 0, (-1)/4, 1/8]])

> ED:=Eigenvalues(Matrix(Delta)):EDD:=map(fnormal,map(evalf,ED),3):print(EDD,map(fnormal,map(abs,EDD),3));

Vector[column]([[0.], [-.500], [.250+.661*I], [.250-.661*I]]), Vector[column]([[0.], [.500], [.707], [.707]])

> evalm(A^9);trace(A);

matrix([[85/512, 85/512, 171/512, 171/512], [85/512, 85/512, 171/512, 171/512], [43/256, 85/512, 85/256, 171/512], [85/512, 43/256, 171/512, 85/256]])

0

> phi:=diag(seq(1,i=1..n)):

> phi[n,n]:=-1:detta:=multiply(phi,Delta):

> Eigenvectors(Matrix(detta)):TT:=map(fnormal,map(evalf,%[2]),4);

TT := Matrix([[1., 1., -1.000-0.*I, -1.000+0.*I], [-1., 1., 1.000+0.*I, 1.000-0.*I], [0., 1., .5000+1.323*I, .5000-1.323*I], [1., 1., 1., 1.]])

> for i to d do map(fnormal,multiply(A^i,TT),3);kol[i]:=col(%,2) od:stackmatrix(seq(kol[i],i=1..d));rank(%);

matrix([[1.00, 1.00, 1.00, 1.00], [1.00, 1.00, 1.00, 1.00], [1.00, 1.00, 1.00, 1.00]])

1

> unassign('tau');e:=vector(n):phi:=matrix(n,n,(i,j)->if i=j then e[i] else 0 fi):

> deta:=factor(det(A+tau*phi&*Delta));

deta := -1/8*(1+tau*e[3])*(1+tau*e[4])*tau*(e[2]+e[1])

> #ee:=[1,1,1,1,1,1,1,1]:#r:=ecliffsubsmat(evalm(A+phi&*Delta),ee);b:=ecliffsubsmat(evalm(A-phi&*Delta),ee);iszero(r+b-2*A);

> #factor(ecliffsubs(deta,ee));

> ISX:=NullSpace(Matrix(transpose(evalm(Delta)))):ND:=NullSpace(Matrix(Delta)):

> CA:=r->if r>0 then concat(seq(NA[q],q=1..r))  else "" fi:CTA:=r->if r>0 then concat(seq(NTA[q],q=1..r))  else "" fi:

> "ker A",CA(nullity(A)),"ker Tr A",CTA(nullity(A)),"ker Delta",concat(seq(ND[q],q=1..nullity(Delta))),"ker  tr Delta",concat(seq(ISX[q],q=1..nullity(Delta)));

> P:=concat(seq(ND[q],q=1..nu),seq(zeta,q=1..d)):Q:=transpose(concat(seq(ISX[q],q=1..nu),seq(zeta,q=1..d))):Q:=submatrix(Q,1..n-d,1..n):print(Q);iszero(evalm(Delta&*P)),iszero(evalm(Q&*Delta));evalm(2*u&*Delta);

matrix([[1, 1, 0, 0]])

true, true

vector([1, 1, -1, -1])

> `charpoly of A`=factor(collect(charpoly(A,s),s)):`rank of Delta`=rank(Delta);`charpoly of Delta`=factor(charpoly(Delta,s));`minpoly`=factor(minpoly(Delta,s));`invcharpoly`=collect(det(J-tau*Delta),tau):`d-th deriv of adjoint of id minus tau Delta`=row(map(diff,adjoint(J-tau*Delta),tau$d),1);ee:=vector(n,1):"rank of Omega union Delta"=rank(stackmatrix(pi,Delta));

`rank of Delta` = 3

`charpoly of Delta` = 1/4*s*(1+2*s)*(2*s^2-s+1)

`linalg:-minpoly` = 1/4*s*(1+2*s)*(2*s^2-s+1)

`d-th deriv of adjoint of id minus tau Delta` = vector([3/4, 3/4, 0, 0])

> sigma:=multiply(pi,u):spi:=evalm(pi/sigma);Omega:=stackmatrix(seq(spi,q=1..n)):#Hpi:=multiply(spi,F);Fpi:=multiply(spi,QD);

spi := vector([1/6, 1/6, 1/3, 1/3])

> factor(minpoly(R,s)),factor(minpoly(B,s));JordanForm(Matrix(R),output=['J']),JordanForm(Matrix(B),output=['J']);

(s-1)*(s+1)*(s^2+1), s*(s-1)*(s+1)

Matrix([[-1, 0, 0, 0], [0, 1, 0, 0], [0, 0, I, 0], [0, 0, 0, -I]]), Matrix([[0, 0, 0, 0], [0, -1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]])

> ###################### V CALCULATED HERE #################

> #Eta:=evalm(J-(A-Omega)):EE:=evalm(Eta^(-1)):V:=evalm(Delta&*Eta^(-1)):print(V);N:=nullity(Delta):

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

> UT(Matrix(R),Matrix(B));

Error, (in UT) invalid subscript selector

> "n"=n;"rank Delta"=d;matrans(R),matrans(B);print(pi);OMEGA0:=row(adj(J-R),1);OMEGA1:=row(adj(J-B),1);VV:=multiply(Delta,inverse(J-B+adj(J-B))):VVV:=multiply(Delta,inverse(J-B+Omega)):zaz:=iszero(evalm((VV-VVV)));

vector([4, 3, 1, 2]), vector([3, 4, 4, 3])

vector([1, 1, 2, 2])

OMEGA0 := vector([1, 1, 1, 1])

OMEGA1 := vector([0, 0, 1, 1])

zaz := true

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

> ################################# MAIN RUN ############################

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

> e:=vector(n);phi:=diag(seq(e[q],q=1..n)):deta:=multiply(phi,Delta);

e := array(1 .. 4, [])

deta := matrix([[0, 0, -1/2*e[1], 1/2*e[1]], [0, 0, 1/2*e[2], -1/2*e[2]], [1/2*e[3], 0, 0, -1/2*e[3]], [0, 1/2*e[4], -1/2*e[4], 0]])

> #unassign('Y','YY');Y[0]:=evalm(spi):for i from 1 to d do Y[i]:=map(simplify,multiply(Y[i-1],deta)) od:KU:=((stackmatrix(seq(spi+Y[b],b=0..d))));

> #KUU:=(map(simplify,multiply(KU,transpose(KU))));

> #KUV:=ecliffe(KUU);

> ee:=vector(n,1):print(ee),ecliffsubsmat(KUU,ee);det(%);

vector([1, 1, 1, 1])

matrix([[10/9, 4/9, 11/18, 19/36], [4/9, 5/18, 7/36, 17/72], [11/18, 7/36, 13/36, 5/18], [19/36, 17/72, 5/18, 37/144]])

0

> #z1:=((det(KUU)));

> z:=factor(ecliffit(z1));

z := 1/10368*e[4]*e[1]-11/93312*e[3]*e[1]*e[4]*e[2]-1/10368*e[4]*e[2]+7/93312*e[4]*e[3]-25/93312*e[1]*e[2]-1/10368*e[3]*e[1]+1/10368*e[3]*e[2]+29/93312

> ZVal:=ecliffsubs(z1,vector(n,0));ifactor(%);

ZVal := 0

0

> ZVal:=ecliffsubs(z,vector(n,0));ifactor(%);

ZVal := 29/93312

``(29)/(``(2)^7*``(3)^6)

> ee:=vector(n,1):ZVal:=ecliffsubs(z,ee);ee[2]:=-1:ee[6]:=1:ZVal:=ecliffsubs(z,ee);

ZVal := 0

Error, 1st index, 6, larger than upper array bound 4

ZVal := 1/1296

> gZ:=grad(z,e):ecliffsubsmat(gZ,ee);

vector([1/2592, (-1)/2592, 0, 1/2592])

> c:=8;ee[c]:=-ee[c]:print("ee=",ee);Zval:=ecliffsubs(z,ee);

c := 8

Error, 1st index, 8, larger than upper array bound 4

Zval := 1/1296

> KK:=ecliffsubsmat(KUV,vector(n,1));

KK := matrix([[10/9, 4/9, 11/18, 19/36], [4/9, 5/18, 7/36, 17/72], [11/18, 7/36, 13/36, 5/18], [19/36, 17/72, 5/18, 37/144]])

> map(Re,map(evalf,Eigenvalues(Matrix(KK))));

Vector[column]([[0.], [0.], [1.879470605], [.1274738387]])

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

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

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

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

vector([4, 3, 1, 2]), vector([3, 4, 4, 3]), matrix([[1, 1, 0, 0]]), vector([1, 1, 2, 2])

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

> OR:=multiply(R,R,B);OB:=multiply(R,B,R);checkspan(KBAS,OR),checkspan(KBAS,OB);expVec(KBAS,OR),expVec(KBAS,OB);

OR := matrix([[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 0, 1]])

OB := matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0]])

matrix([[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), true, matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0]]), true

vector([1, 0, 1, 0, 0, 0]), vector([0, 1, 0, 0, 0, 1])

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

d := 3

ND := {vector([1, 1, 1, 1])}

{vector([1, 1, 1, 1])}

> ee:=[1,1,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));

vector([4, 3, 1, 2]), vector([3, 4, 4, 3])

(s-1)*(s+1)*(s^2+1), s*(s-1)*(s+1)

{vector([1, 1, 0, 0])}

1/4*x+1/4*x^2+x^4, 1/4*x*(2*x+1)*(2*x^2-x+1)

> 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:for IX to n do for i to n do if(RCYC[IX,i]<>0) then k:=k+1;s0:={i}union s0 fi; od;S0:=S0 union {s0};RC:=RC union s0;s0:={}: od:BC:={}:j:=0: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; od;T0:=T0 union {t0};BC:=BC union t0;t0:={} od:rcyc:=k*RCYC:bcyc:=evalm(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;

antR := {}

antB := {1, 2}

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

matrix([[1/6, 1/6, 1/3, 1/3], [1/3, 1/3, 1/6, 1/6], [1/12, 1/12, 5/12, 5/12], [5/24, 5/24, 7/24, 7/24]]), 2,

matrix([[1/6, 1/6, 1/3, 1/3], [0, 0, 1/2, 1/2], [1/4, 1/4, 1/4, 1/4], [1/8, 1/8, 3/8, 3/8]]), 2,

matrix([[1/6, 1/6, (-1)/6, (-1)/6], [(-1)/12, (-1)/12, 1/12, 1/12], [1/24, 1/24, (-1)/24, (-1)/24]])

> K2:=stackmatrix(Y0[0],row(K0,2),row(K1,3),row(K0,4));

K2 := matrix([[1/6, 1/6, 1/3, 1/3], [1/3, 1/3, 1/6, 1/6], [1/4, 1/4, 1/4, 1/4], [5/24, 5/24, 7/24, 7/24]])

> charpoly(Detta,t);factor(%);

t^4+1/4*t+1/4*t^2

1/4*t*(2*t+1)*(2*t^2-t+1)

> #map(simplify,multiply(spi,inverse(Matrix(evalm(J-t*Detta)))));

> factor(1/simplify(det(J-t*Detta)));

4/((t+2)*(t^2-t+2))

> Series:=map(simplify,evalm(multiply(spi,add(evalm(psi[n]*Detta^n),n=0..d)))):Resolvent:=map(simplify,multiply(spi,inverse(Matrix(evalm(J-t*Detta))))):MM:=evalm(Series-Resolvent):

> sys:=solve({seq(MM[i],i=1..n)},{seq(psi[n],n=0..d)});

sys := {psi[3] = psi[3], psi[2] = psi[2], psi[0] = 1, psi[1] = -1/4*(2*psi[3]+psi[3]*t-2*psi[2]*t-4*psi[2]-8*t)/(t+2)}

> #RowSpace(Matrix(K0));RowSpace(Matrix(K1));RowSpace(Matrix(K2));

> factor(det(multiply(K2,transpose(K2))));

0

> map(fnormal,map(evalf,K2),3);

matrix([[.167, .167, .333, .333], [.333, .333, .167, .167], [.250, .250, .250, .250], [.208, .208, .292, .292]])

> KK:=stackmatrix(seq(Y0[b],b=1..d));

KK := matrix([[1/6, 1/6, (-1)/6, (-1)/6], [(-1)/12, (-1)/12, 1/12, 1/12], [1/24, 1/24, (-1)/24, (-1)/24]])

> NK:=nullspace(transpose(KK));

NK := {vector([0, 1, 2]), vector([1, 0, -4])}

> 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=0..d)); print(k,rank(KA[k]));phi[k,k]:=ee[k]; od:

1, 4

2, 4

3, 4

4, 4

> #for k from 1 to n  do phi[k,k]:=-ee[k];

> #Resolvent:=map(simplify,multiply(spi,inverse(Matrix(evalm(J-t*Detta))))):

> #print(k,Resolvent);

> #phi[k,k]:=ee[k]; od:

> "minpoly",minpoly(Detta,s),factor(minpoly(Detta,s));

> if ( rank(KK)=d) then mm:=minpoly(Detta,s):NKN:=CoefficientList(mm,s):ddd:=degree(mm)+1 :

> p:=add(NKN[j]*s^(j-1),j=1..ddd):q:=add(abs(NKN[j])*t^(j-1),j=1..ddd):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(KK)<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)); print(checkspan(KBAS,SU),checkspan(KBAS,ST),expVec(KBAS,SU),expVec(KBAS,ST)) od;fi;

s^2+2*s^3, t^2+2*t^3

matrix([[(-1)/2, 1/2, 0, 0], [1/2, (-1)/2, 0, 0], [0, 0, (-1)/2, 1/2], [0, 0, 1/2, (-1)/2]]), matrix([[1/2, 1/2, 1, 1], [1/2, 1/2, 1, 1], [1/2, 1/2, 1, 1], [1/2, 1/2, 1, 1]]), matrix([[0, 1/3, 1/3, 1/...

matrix([[1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3]]), matrix([[1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3]])

matrix([[1/3, 0, 1/3, 1/3], [0, 1/3, 1/3, 1/3], [1/6, 1/6, 1/2, 1/6], [1/6, 1/6, 1/6, 1/2]]), true, matrix([[0, 1/3, 1/3, 1/3], [1/3, 0, 1/3, 1/3], [1/6, 1/6, 1/6, 1/2], [1/6, 1/6, 1/2, 1/6]]), true, ...

s-4*s^3, t+4*t^3

matrix([[1/2, (-1)/2, -1, 1], [(-1)/2, 1/2, 1, -1], [1/2, (-1)/2, 1, -1], [(-1)/2, 1/2, -1, 1]]), matrix([[1/2, 1/2, 2, 2], [1/2, 1/2, 2, 2], [3/2, 1/2, 1, 2], [1/2, 3/2, 2, 1]]), matrix([[1/5, 0, 1/5...

matrix([[1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3]]), matrix([[1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3]])

matrix([[1/20, 1/4, 11/20, 3/20], [1/4, 1/20, 3/20, 11/20], [3/20, 3/20, 1/20, 13/20], [3/20, 3/20, 13/20, 1/20]]), false, matrix([[1/4, 1/20, 3/20, 11/20], [1/20, 1/4, 11/20, 3/20], [7/20, (-1)/20, 9...

> print(PP);for i from 1 to doDIM(IDEMS) do evalm(IDEMS[i]&*PP) od;

matrix([[1/2, (-1)/2, -1, 1], [(-1)/2, 1/2, 1, -1], [1/2, (-1)/2, 1, -1], [(-1)/2, 1/2, -1, 1]])

matrix([[(-1)/2, 1/2, -1, 1], [1/2, (-1)/2, 1, -1], [1/2, (-1)/2, 1, -1], [(-1)/2, 1/2, -1, 1]])

matrix([[1/2, (-1)/2, 1, -1], [(-1)/2, 1/2, -1, 1], [1/2, (-1)/2, 1, -1], [(-1)/2, 1/2, -1, 1]])

matrix([[1/2, (-1)/2, -1, 1], [(-1)/2, 1/2, 1, -1], [(-1)/2, 1/2, 1, -1], [1/2, (-1)/2, -1, 1]])

matrix([[1/2, (-1)/2, -1, 1], [(-1)/2, 1/2, 1, -1], [1/2, (-1)/2, -1, 1], [(-1)/2, 1/2, 1, -1]])

> cmm(ST,A),cmm(SU,A);

matrix([[(-1)/6, 0, 0, (-1)/6, 1/6, 1/6], [(-1)/6, 0, 0, (-1)/6, 1/6, 1/6], [0, 1/6, 1/6, 0, (-1)/6, (-1)/6], [1/6, (-1)/6, (-1)/6, 1/6, 0, 0], [1/6, (-1)/6, (-1)/6, 1/6, 0, 0], [0, 1/6, 1/6, 0, (-1)/...

> goLie(SA,A);

3,

> JordanForm(Matrix(evalm(ST))),FrobeniusForm(Matrix(ST)),cmm(SA,Omega),det(SA);

Matrix([[0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0], [0, 0, (-1)/3, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, (-1)/3, 0], [0, 0, 0, 0, 0, 0]]), Matrix([[0, 0, 0, 0, 0, 0], [1, 0, 1/3, 0, 0, 0], [0, 1, 2/3...

> unassign('IDEMS'):read "/home/ph/GAP/idems";

IDEMS[1] := [[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 0, 1]]

IDEMS[2] := [[0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]]

IDEMS[3] := [[1, 0, 0, 0], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0]]

IDEMS[4] := [[1, 0, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0]]

> #for i to doDIM(IDEMS) do evalm(IDEMS[i]) od;

> for i to doDIM(IDEMS) do matrans(IDEMS[i]) od;

vector([4, 3, 3, 4])

vector([3, 4, 3, 4])

vector([1, 2, 2, 1])

vector([1, 2, 1, 2])

> #for i to doDIM(IDEMS) do multiply(R,IDEMS[i]),multiply(B,IDEMS[i]),multiply(2*Delta,IDEMS[i]) od;

> #evalm(2*Delta);

> unassign('KROSPS'):

> for i to doDIM(IDEMS) do KRISPS[i]:=multiply(Omega,IDEMS[i]):KROSPS[i]:=multiply(spi,IDEMS[i]):print(matrans(IDEMS[i]),KROSPS[i]) od:

vector([4, 3, 3, 4]), vector([0, 0, 1/2, 1/2])

vector([3, 4, 3, 4]), vector([0, 0, 1/2, 1/2])

vector([1, 2, 2, 1]), vector([1/2, 1/2, 0, 0])

vector([1, 2, 1, 2]), vector([1/2, 1/2, 0, 0])

> uu:=matrix(n,1,1):

> KRSPS:=stackmatrix(seq(KROSPS[i],i=1..doDIM(IDEMS))):RS:=rowspace(KRSPS);for i to doDIM(RS) do VS[i]:=evalm(uu&*transpose(convert(RS[i],matrix))) od;

RS := {vector([1, 1, 0, 0]), vector([0, 0, 1, 1])}

VS[1] := matrix([[1, 1, 0, 0], [1, 1, 0, 0], [1, 1, 0, 0], [1, 1, 0, 0]])

VS[2] := matrix([[0, 0, 1, 1], [0, 0, 1, 1], [0, 0, 1, 1], [0, 0, 1, 1]])

> checkspan(VS,Omega),expVec(VS,Omega);

matrix([[1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3], [1/6, 1/6, 1/3, 1/3]]), true, vector([1/6, 1/3])

> doDIM(IDEMS);

4

> MK:=stackmatrix(seq(KROSPS[i],i=1..doDIM(IDEMS)));rank(MK);

MK := matrix([[0, 0, 1/2, 1/2], [0, 0, 1/2, 1/2], [1/2, 1/2, 0, 0], [1/2, 1/2, 0, 0]])

2

> #unassign('Q'):read "/home/ph/GAP/Q":for i to doDIM(Q) do matrans(evalm(Q[i])),multiply(spi,Q[i]) od;

> unassign('KBAS'):read "/home/ph/GAP/kbas":for i to doDIM(KBAS) do evalm(KBAS[i]) od;

matrix([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, -1], [0, 0, -1, 1]])

matrix([[0, 0, 0, 0], [0, 0, 0, 0], [1, -1, 0, 0], [-1, 1, 0, 0]])

matrix([[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 1, 0]])

matrix([[0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 0, 1], [0, 0, 1, 0]])

matrix([[0, 1, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0]])

matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0]])

> #for i to doDIM(Q) do checkspan(KBAS,Q[i]);expVec(KBAS,Q[i]) od;

> for i to doDIM(IDEMS) do checkspan(KBAS,IDEMS[i]);expVec(KBAS,IDEMS[i]) od;

matrix([[0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 1, 0], [0, 0, 0, 1]]), true

vector([1, 0, 1, 0, 0, 0])

matrix([[0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]]), true

vector([1, 0, 0, 1, 0, 0])

matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0]]), true

vector([0, 0, 0, 0, 0, 1])

matrix([[1, 0, 0, 0], [0, 1, 0, 0], [1, 0, 0, 0], [0, 1, 0, 0]]), true

vector([0, 1, 0, 0, 0, 1])