#
# Compute scalars for transforming an orthogonal hereditary matrix model
#  for the i-th irrep of W(R) into an optimal rational seminormal form.
#  Requires the orthogonal model to be built first.
# Reference: Sections 3E and 3F of "Explicit matrices for irreducible..."
#
# Set infolevel[seminormal]:=2 to see info during the optimization step.

seminormal:=proc(i,n) local inds,eq,N,j,k,v,ch,sc,T,r; global Semi;
  if not assigned(Semi[n]) then Semi[n]:={} fi;
  if n=0 or member(i,Semi[n]) then RETURN(`Done`)
    elif not assigned(Rep[i,n]) then ERROR(`orthogonal data required`)
    elif i<0 then seminormal(-i,n);
    for j in Branch[i,n] do Semi[j,i,n]:=Semi[j,-i,n] od;
    Semi[n]:={op(Semi[n]),i}; RETURN(`Done`)
  fi;
  for j in Branch[i,n] do seminormal(j,n-1) od;
  inds:=sort(map(op,[indices(Rep[i,n])]));
  eq:=NULL; v:=[seq(cat('e',j),j=Branch[i,n])];
  if nops(v)>1 then for j in inds do
    ch:=Rep[i,n][j]; N:=nops(ch[1]);
    sc:=[seq(cat('e',T[N-1]),T=ch)];
    for k to N-2 do r:=n+k-N+1;
      sc:=zip((x,y)->x*y,[seq(Semi[T[k],T[k+1],r],T=ch)],sc)
    od;
    sc:=simplify(sc,'sqrt');
    eq:=eq,op(subs({seq(e[op(ch[k])]=e[op(ch[k])]/sc[k],k=1..nops(ch))},
      [seq(sc[k]*Rep[ch[k],n],k=1..nops(ch))]));
  od fi;
  eq:=map(coeffs,expand([eq]),indets([eq],'indexed'));
  eq:=map(proc(x) if not type(x,'constant') then x fi end,eq);
  sc:=`seminormal/rat`(eq,v);
  for j to nops(Branch[i,n]) do Semi[Branch[i,n][j],i,n]:=sc[j] od;
  Semi[n]:={op(Semi[n]),i}; `Done`;
end;

# Given a list a of terms of the form  r*x/y , where r^2 is rational and
# x,y are variables in the list v, find a rescaling of the variables in v
# such that all of the scalars r are rational (given that one exists).

`seminormal/rat`:=proc(a,v) local b,res,so;
  if nops(v)=1 then RETURN([1]) fi;
  b:=a; res:=v;
  while nops(b)>0 do
    if type(b[1],'constant') then ERROR(`this cannot happen`) fi;
    so:=solve(b[1]=sign(evalf(lcoeff(b[1]))));
    b:=map(proc(x) if not type(x,'rational') then x fi end,
      simplify(subs(so,b),'sqrt'));
    res:=subs(so,res);
  od;
  res:=map(lcoeff,res);
  b:=subs(zip((x,y)->x=y*x,v,res),a);
  b:={op(simplify(b,'sqrt'))};
  res:=zip((x,y)->x*y,res,`seminormal/min`(b,v));
  res:=map((x,y)->simplify(x/y,'sqrt'),res,res[1]);
  if min(op(evalf(res)))>0 then res
    else ERROR(`negative rescaling cannot occur`) fi;
end:

# Given a set a of terms of the form  r*x/y , where r is rational and
# x,y are variables in the list v, find a rescaling of the variables
# in v that minimizes the least common denominator. Unfortunately, there
# are many optimal solutions; choosing a "bad" optimum means that the
# solutions for W[n+1] will not be as good.

`seminormal/min`:=proc(a,v) local f,b,res,u,sat,U,p;
  b:=a; res:=[1$nops(v)]; sat:=0;
  while sat<nops(v) do
    sat:=sat+1; u:=v[sat];
    f:=`seminormal/line`(b,u,1);
    f:=f/`seminormal/line`(b,u,-1);
    f:=`seminormal/sqpt`(f);  
    if f=1 then next fi;
    userinfo(2,seminormal,cat(`rescaling `,u,` by`),f);
    res:=subsop(sat=res[sat]/f,res);
    b:=subs(u=u/f,b); sat:=0;
  od;
  f:=ilcm(op(map(x->denom(lcoeff(x)),b)));
  for p in numtheory['factorset'](f) do
    U:=`seminormal/opt`(b,v,p);
    if {op(U)}<>{1} then userinfo(2,seminormal,`improved denom`,U) fi;
    b:=subs(zip((x,y)->x=y*x,v,U),b);
    res:=zip((x,y)->y*x,res,U);
    U:=`seminormal/n_opt`(b,v,p);
    if {op(U)}<>{1} then userinfo(2,seminormal,`improved numer`,U) fi;
    b:=subs(zip((x,y)->x=y*x,v,U),b);
    res:=zip((x,y)->y*x,res,U);
  od; res
end:

# Compute the "rational gcd" of all leading terms involving u^c in b.

`seminormal/line`:=proc(b,u,c) local f;
  f:=map((x,y)->lcoeff(coeff(x,y[1],y[2])),b,[u,c]);
  igcd(op(map(numer,f)))/ilcm(op(map(denom,f)));
end:

# Extract the square part of nonzero rational r.

`seminormal/sqpt`:=proc(r) local res,q,p,a;
  res:=1; q:=numer(r);
  for p in numtheory['factorset'](q) do
    while irem(q,p^2,'a')=0 do q:=a; res:=res*p od;
  od;
  q:=denom(r);
  for p in numtheory['factorset'](q) do
    while irem(q,p^2,'a')=0 do q:=a; res:=res/p od;
  od; res
end:

# Get the power of prime p involved in nonzero rational r.

`seminormal/getpow`:=proc(r,p) local q,res,c;
  q:=numer(r); res:=0;
  while irem(q,p,'c')=0 do q:=c; res:=res+1 od;
  q:=denom(r);
  while irem(q,p,'c')=0 do q:=c; res:=res-1 od;
  res
end:

# Use linear programming to minimize the highest power of p occurring
# in any denominator.

`seminormal/opt`:=proc(a,v,p) local c,neq,c0,f,g,j,k,so;
  neq:=NULL;
  for j to nops(v) do
    g:=map(coeff,a,v[j]);
    for k to nops(v) do
      f:=subs(0=NULL,map(coeff,g,v[k],-1));
      if nops(f)>0 then
        f:=min(op(map(`seminormal/getpow`,f,p)));
        neq:=neq,c[k]-c[j]<=f-c[0];
      fi
    od
  od;
  neq:=subs(c[1]=0,{neq});
  so:=simplex['maximize'](c[0],neq); c0:=subs(so,c[0]);
  if not type(c0,'integer') then c0:=floor(c0);
    so:=simplex['maximize'](c[0],{op(neq),c[0]=c0});
  fi;
  so:=subs(so,[0,seq(c[j],j=2..nops(v))]);
  if not type(so,'list(integer)') then ERROR(
    `simplex[maximize] incorrectly returned a non-integer solution`) fi;
  map((x,y)->y^x,so,p);
end:

# Use linear programming to minimize the highest power of p occurring
# in any numerator, subject to not increasing the highest power of p
# occurring among the denominators.

`seminormal/n_opt`:=proc(a,v,p) local c,neq,m1,c0,c1,f,g,j,k,so;
  neq:=NULL; m1:=NULL;
  for j to nops(v) do
    g:=map(coeff,a,v[j]);
    for k to nops(v) do
      f:=subs(0=NULL,map(coeff,g,v[k],-1));
      if nops(f)>0 then
        f:=op(map(`seminormal/getpow`,f,p)); m1:=min(m1,f);
        neq:=neq,c[k]-c[j]<=min(f)-c1,c[j]-c[k]<=c[0]-max(f);
      fi
    od
  od;
  neq:=subs({c1=m1,c[1]=0},{neq});
  so:=simplex['minimize'](c[0],neq); c0:=subs(so,c[0]);
  if not type(c0,'integer') then c0:=ceil(c0);
    so:=simplex['minimize'](c[0],{op(neq),c[0]=c0});
  fi;
  so:=subs(so,[0,seq(c[j],j=2..nops(v))]);
  if not type(so,'list(integer)') then ERROR(
    `simplex[minimize] incorrectly returned a non-integer solution`) fi;
  map((x,y)->y^x,so,p);
end: