# Tools for manipulating sparse matrices.

# A matrix is a list of linear polynomials in the variables e1,e2,...
# The i-th polynomial represents the i-th row of the matrix.
# Example: e3-2*e5 represents the row vector [0,0,1,0,-2,...]

# Compute the matrix product a * b.

fast_prod:=proc(a,b) local i;
  map(expand,subs({seq(cat('e',i)=b[i],i=1..nops(b))},a))
end;

# Note: If a and b are matrices in this format, then versions of Maple
# starting with V.4 automatically "do the right thing" with linear
# combinations such as 3*a-2*b.

########################

# transpose a *square* matrix in sparse format

fast_transpose:=proc(a) local i,j,N,b,res,co,tm,v,q;
  N:=nops(a); res:=table([seq(i=0,i=1..N)]);
  b:=subs({seq(cat('e',i)=q^i,i=1..N)},a);
  for i to N do
    co:=[coeffs(b[i],[q],'tm')];
    tm:=map(degree,[tm],q); v:=cat('e',i);
    for j to nops(tm) do res[tm[j]]:=res[tm[j]]+co[j]*v od;
  od;
  res:=[seq(res[i],i=1..N)];
end;

########################

# test whether sparse matrix a is symmetric

is_sym:=proc(a) local b,j,q,k,co,tm;
  b:=subs({seq(cat('e',j)=q^j,j=1..nops(a))},a);
  for j to nops(a) do
    co:=[coeffs(b[j],[q],'tm')]; tm:=[tm];
    co:={seq(coeff(b[degree(tm[k],q)],q,j)-co[k],k=1..nops(co))};
    if simplify(co,'sqrt')<>{0} then RETURN(false) fi
  od; true
end;

########################

# Given symmetric matrices a and b and an integer m>=1, compute the
# set of matrix entries of the braid difference (a*b*a*...)-(b*a*b*...)
# (each product having m terms). Since the difference is symmetric (m odd)
# or skew-symmetric (m even), we only extract coefficients on or above the
# main diagonal.

# If there is a fourth argument (i.e., a flag), don't assume that a and b
# are symmetric, and return *all* coefficients of the braid difference.

# Note: to build hereditary irreps relative to the standard chains,
# we never need to apply the braid relations with m=4 for matrices larger
# than 3x3, and m=6 cases are never larger than 2x2.

braid_eqns:=proc(m,a,b) local u,co,i,vars;
  u:=[a,b]; co:=table();
  for i to m-1 do u:=[fast_prod(u[2],a),fast_prod(u[1],b)] od;
  u:=zip((x,y)->x-y,u[1],u[2]);
  vars:={seq(cat('e',i),i=1..nops(a))};
  if nargs>3 then RETURN(map(coeffs,{op(u)},vars)) fi;
  for i to nops(a) do
    co[i]:={coeffs(u[1],vars)};
    vars:=subs(cat('e',i)=NULL,vars);
    u:=subs(cat('e',i)=0,subsop(1=NULL,u));
  od;
  {seq(op(co[i]),i=1..nops(a))}
end;

#####################

# Given a group element w=[i_1,...,i_l] of para-Coxeter type (i.e., a
# product of *distinct* simple reflections), evaluate the trace of w in
# the i-th irrep of the Weyl group of rank n in the database.
#
# Requires the diagonal of the matrix representing w to be the product
# of the diagonals of its factors. (True for graceful parabolic chains.)

cheap_trace:=proc(w,i,n) local tr,rg,k,T;
  tr:=listchains(i,0,n);
  rg:=[seq(Rep[0,k]+1..k+1,k=1..n)];
  tr:=[seq(convert([seq(coeff(
    Rep[[op(rg[k],T)],k],e[op(rg[k],T)]),k=w)],`*`),T=tr)];
  convert(tr,`+`);
end;

# Compute the trace of w in the i-th irrep in an honest way.
# We unpack a full matrix for each term appearing in w, so if w has
# repetitions and the dimension is large, this is grossly inefficient.

slow_trace:=proc(w,i,n) local a,j,ch;
  if nops(w)=0 then RETURN(Rep[n][i][1]) fi;
  ch:=listchains(i,0,n); a:=restore(w[1],n,ch);
  for j from 2 to nops(w) do a:=fast_prod(a,restore(w[j],n,ch)) od;
  convert([seq(coeff(a[j],cat('e',j)),j=1..nops(a))],`+`)
end;