#
# seminormal - matrices for Young's seminormal representations of S_n
# rep_matrix - representing matrices for arbitrary permutations
#
# Calling sequence:
#  seminormal(mu);
#  rep_matrix(w,A);
#
# Parameters:
#  mu - a partition of some integer n>1
#   w - a permutation of [1,2,...,n] 
#   A - a table of matrices
#
# If mu is a partition of n>1, seminormal(mu) computes the representing
# matrices for the adjacent transpositions in Young's seminormal form for
# the irreducible representation of S_n indexed by mu.
#
# The results are returned in a table whose i-th entry is the representing
# matrix for the transposition (i,i+1). Each matrix is in 'sparse' format,
# with the rows and columns corresponding to standard Young tableaux in the
# last-letter-highest-row-order. (See the documentation for the 'tableaux'
# function elsewhere in the examples directory.)
#
# Once the matrices for the adjacent transpositions have been assigned to
# some table A, the function rep_matrix(w,A) may be used to compute the
# representing matrix for any permutation w of [1,2,...,n].
#
# Reference:
#  G. James and A. Kerber, "The Representation Theory of the Symmetric
#  Group," Section 3.3.
#
# Examples:
#  A:=seminormal([3,1]);
#  print(A[1],A[2],A[3]);
#  rep_matrix([2,4,1,3],A);

seminormal:=proc(mu) local k,m,tab,i,j,r1,r2,c1,c2,A,ad;
  tab:=tableaux(mu); m:=nops(tab); A:=table();
  for k to convert(mu,`+`)-1 do
    A[k]:=array('sparse',1..m,1..m);
    for i to m do
      if member(subs({k=k+1,k+1=k},tab[i]),tab,'j') and j<i then next fi;
      for r1 while not member(k,tab[i][r1],'c1') do od;
      for r2 while not member(k+1,tab[i][r2],'c2') do od;
      ad:=(r1+c2-r2-c1); A[k][i,i]:=1/ad;
      if ad>1 then
        A[k][j,j]:=-1/ad; A[k][i,j]:=1+1/ad; A[k][j,i]:=1-1/ad;
      fi
    od
  od; op(A)
end;

rep_matrix:=proc(w,A) local i,w0;
  for i to nops(w)-1 while w[i]<w[i+1] do od;
  if i<nops(w) then
    w0:=subsop(i=w[i+1],i+1=w[i],w);
    evalm(rep_matrix(w0,A) &* A[i])
  else
    linalg['band']([1],linalg['rowdim'](A[1]))
  fi
end;

# A copy of SF/examples/tableaux:

tableaux:=proc(mu) local nu,n,i,old,new;
  n:=convert(mu,`+`); new:=NULL;
  if n=0 then RETURN([[]]) fi;
  if nops(mu)<2 then RETURN([[[$1..n]]]) fi;
  for i to nops(mu) do
    if i<nops(mu) and mu[i]=mu[i+1] then next fi;
    if mu[i]>1 then nu:=mu[i]-1 else nu:=NULL fi;
    old:=tableaux(subsop(i=nu,mu));
    if mu[i]=1 then old:=map(x->[op(x),[]],old) fi;
    new:=new,op(map((x,j,m)->subsop(j=[op(x[j]),m],x),old,i,n))
  od;
  [new]
end;