# # 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 satdenom(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: