#
# orthogonal(i,R) builds hereditary orthogonal matrices for the i-th irrep
#  of W(R) and stores the results in a table.
# By default, assume that the simple roots are ordered gracefully.
# Add the optional argument 'slow' to indicate an ungraceful order.
# Note that this code is not stand alone: it requires quadsolve, sparseops,
#  clone_test, unpack, orthopt, and the coxeter package to be loaded first.

orthogonal:=proc(i,R) local S,S0,j,m,cl,n,ch,tr,b; global Rep,Branch;
  S:=coxeter['base'](R); n:=nops(S);
  if n=0 or assigned(Rep[i,n]) then RETURN(`Done`)
    elif i<0 then RETURN(build_neg(-i,n,R)) fi;
  m:=coxeter['cox_matrix'](S);
  m:=[seq(m[j,n],j=1..n)];
  if not assigned(Rep[n]) then Rep[n]:=coxeter['irr_chars'](S) fi;
  if not assigned(Rep[0,n]) then
    for j while m[j]=2 do od;
    if m[j]=1 then j:=j+1 fi; Rep[0,n]:=j-1;
  fi;
  S0:=subsop(n=NULL,S);
  if not assigned(Rep[n-1]) then Rep[n-1]:=coxeter['irr_chars'](S0) fi;
  cl:=clone_test(i,S0,S); b:=Branch[i,n];
  if nops(b)<>nops({op(b)}) then
    b:=fast_prod(Branch[n],Rep[n][i]);
    if max(op(b))>2 then ERROR(`multiplicities are too high`) fi;
    b:=map(proc(x,y) if y[x]=1 then x elif y[x]=2 then x,-x fi end,
      [$1..nops(b)],b);
    Branch[i,n]:=b;
  fi;
  for j in b do orthogonal(j,S0) od;
  ch:=listchains(i,Rep[0,n],n);
  b:=map(x->op(1,x),{op(ch)});
  if member('slow',[args]) then tr:=slow_trace else tr:=cheap_trace fi;
  for j in b do Rep[i,n][j]:=map((x,y)->if x[1]=y then x fi,ch,j) od;
  if traperror(build_op(i,n,m,cl,tr))=lasterror then
    unassign(Rep[i,n]); ERROR(lasterror) else `Done` fi;
end;

# The main engine: generates the defining equations; solves them;
#  verifies; installs the results in tables.
 
build_op:=proc(i,n,m,cl,tr) local a,b,c,eq,j,k,l,T,sol,ch,ch0,ct,st;
  global Rep;
  c:=table('symmetric'); ct:=0; st:=time();
  for j in [indices(Rep[i,n])] do
    ch:=Rep[i,n][op(j)];
    for k to nops(ch) do
      Rep[ch[k],n]:=convert(
        [seq(c[ct+k,ct+l]*e[op(ch[l])],l=1..nops(ch))],`+`);
    od;
    ct:=ct+nops(ch);
  od;
  ch:=map(op,map(op,[entries(Rep[i,n])]));
  a:=restore(n,n,ch);
  eq:=[seq(cat('e',j),j=1..nops(a))];
  eq:=braid_eqns(1,fast_prod(a,a),eq);
  for k from Rep[0,n]+1 to n-1 do
    ch0:=listchains(i,min(Rep[0,n],Rep[0,k]),n);
    b:=restore(k,n,ch0);
    eq:={op(eq),op(braid_eqns(m[k],restore(n,n,ch0),b))};
  od;
  eq:={op(eq),seq(tr(k[1],i,n)-k[2],k=cl),pick_zero(i,n,ch)};
  lprint(cat(nops(eq),` equations, `,nops(indets(eq)),` vars`),time()-st);
  st:=time(); sol:=quadsolve(eq);
  if nops([sol])=0 then ERROR(`no solution found`) fi;
  lprint(cat(`found solution `),time()-st); st:=time();
  sol:=tri_solve(sol);
  if simplify(subs(sol,eq),'sqrt')<>{0}
    then ERROR(`sanity check failed`) fi;
  for T in ch do Rep[T,n]:=subs(sol,Rep[T,n]) od;
  lprint(`verification took `,time()-st);
  if min(op(Branch[i,n]))<0 then
    st:=time(); orthopt(i,n);
    lprint(`orthogonal optimization took`,time()-st)
  fi; NULL
end;

# Get branching data for restricting the i-th irrep of W(S) to W(S0).
# Branch[n] holds the restriction matrix between W(S) and W(S0):
#  given a character ch for W(S), fast_prod(Branch[n],ch) is the list of
#  irreducible multiplicities in the restriction of ch to W(S0).

get_branch:=proc(i,n,S0,S) local j,f; global Branch;
  if assigned(Branch[i,n]) then RETURN(Branch[i,n]) fi;
  if not assigned(Branch[n]) then
    f:=[seq(cat('e',j),j=1..nops(Rep[n]))];
    f:=coxeter['restrict'](f,S0,S);
    Branch[n]:=map(coxeter['cprod'],Rep[n-1],f,S0);
  fi;
  f:=fast_prod(Branch[n],Rep[n][i]);
  Branch[i,n]:=[seq(j$f[j],j=1..nops(f))]
end;

# For non-free reps, install into the database a second copy of each
#  W[n-1]-irrep that occurs with multiplicity two, using a negative index
#  to distinguish it from the first copy.

build_neg:=proc(i,n,R) local j,ch,chn,k; global Rep,Branch;
  orthogonal(i,R);
  Branch[-i,n]:=Branch[i,n];
  for j in [indices(Rep[i,n])] do
    ch:=Rep[i,n][op(j)];
    chn:=map((x,y)->subsop(nops(x)=y,x),ch,-i);
    Rep[-i,n][op(j)]:=chn;
    for k to nops(ch) do
      Rep[chn[k],n]:=subs(
        zip((x,y)->e[op(x)]=e[op(y)],ch,chn),Rep[ch[k],n]);
    od;
  od; `Done`
end;
 
# For non-free reps, identify a set of matrix entries for s[n] whose 
#  vanishing may be used to specify a unique orthogonal model of rational 
#  type. This could fail if we are unlucky and pick matrix entries that
#  vanish generically by accident. (A more robust implementation would set 
#  error traps for failures and make alternative choices.)
# Reference: Section 4A of "Explicit matrices for irreducible..."
 
pick_zero:=proc(i,n,ch) local k,j,u,v,res,neg;
  k:=n-Rep[0,n]; res:=NULL;
  neg:=map(proc(x) if x<0 then -x fi end,Branch[i,n]);
  for j in neg do
    for u in ch do 
      if u[k]=-j then v:=Rep[i,n][u[1]];
        v:=map(proc(x,y,z) if not member(abs(x[z]),y) then x fi end,v,neg,k);
        if nops(v)>0 then res:=res,coeff(Rep[u,n],e[op(v[1])]); break fi
      fi
    od
  od; res
end;