# wreduce.txt
#
# Copyright (c) 1998 Andy Stubbs <agstubbs@liverpool.ac.uk>
#
# routines to help reduce the "hyperlliptic weight" of a rational
# polynomial... hardwired to genus 3 at present, but easily
# generalised.

gpolyutil[g3weight]:=proc(poly)
  local pnum,pdem, w;
  pnum:=collect(expand(subs(x[1]=w^2*x[1],y[1]=w^7*y[1],numer(poly))),w);
  pdem:=collect(expand(subs(x[1]=w^2*x[1],y[1]=w^7*y[1],denom(poly))),w);
### WARNING: degree(0,x) now returns -infinity
### WARNING: degree(0,x) now returns -infinity
  degree(pnum,w)-degree(pdem,w);
end:

gpolyutil[get_heaviest_terms]:=proc(p,g)
  local m, n, terms, poly;
  poly:=expand(x[1]*y[1]*p):
  terms:=0:
### WARNING: degree(0,x) now returns -infinity
  for m from 1 to degree(poly,y[1]) do
    expand(coeff(collect(poly,y[1]),y[1]^m));
### WARNING: degree(0,x) now returns -infinity
    n:=degree(%,x[1]);
    if n>0 then
      coeff(collect(expand(%%),x[1]),x[1]^n);
      if gpolyutil[g3weight](terms)=((2*g+1)*(m-1)+2*(n-1))
        then terms:=terms+x[1]^(n-1)*y[1]^(m-1)*% fi;
      if gpolyutil[g3weight](terms)<((2*g+1)*(m-1)+2*(n-1))
        then terms:=x[1]^(n-1)*y[1]^(m-1)*% fi;
    fi;
  od:
  expand(terms);
end:

gpolyutil[permute_vars]:=proc(poly)
  local s,t,basicpoly,i;
  basicpoly := subs(seq(x[i] = t[i],i = 1 .. nargs-1),
                    seq(y[i] = s[i],i = 1 .. nargs-1),poly);
  subs( seq(t[i] = x[args[1+i]],i = 1 .. nargs-1),
        seq(s[i] = y[args[1+i]],i = 1 .. nargs-1),
        basicpoly )
end:

gpolyutil[get_orbit]:=proc(poly,genus)
  local perms, terms, newterms, i, j;
  perms:=combinat[permute](genus):
  if type(expand(poly),`+`)
    then terms:=[op(expand(poly))]
    else terms:=[poly] fi;
  newterms:={};
  for i from 1 to nops(terms) do
    newterms:={ op(newterms),
                seq(gpolyutil[permute_vars](terms[i],op(perms[j])),
                    j=1..nops(perms))};
  od;
  sum('newterms[i]','i'=1..nops(newterms));
end:

gpolyutil[try_simple_reduce_g3]:=proc(poly)
  local heavy_numer, heavy_denom, weight_quot;
  heavy_numer := gpolyutil[get_heaviest_terms](numer(poly),3);
  heavy_denom := gpolyutil[get_heaviest_terms](denom(poly),3);
  simplify(%%/%);
  if not type(%,`polynom`) then 0 else % fi;
### WARNING: degree(0,x) now returns -infinity
### WARNING: degree(0,x) now returns -infinity
  if degree(collect(expand(%),x[1]),x[1])<degree(collect(expand(%),x[2]),x[2])
    then 0
    else % fi;
  gpolyutil[get_orbit](%,3);
end:

gpolyutil[naive_reduce_g3]:=proc(poly)
  local result, delta_poly;
  result := poly;
  if type(simplify(result),`polynom`) then RETURN(result) fi;
  delta_poly := 1;
  while ( delta_poly <> 0 ) do
    delta_poly := gpolyutil[try_simple_reduce_g3]( result );
    result := simplify(expand(result-delta_poly));
  od:
  result;
end:

# extra routines 28/11/98 to automate the process

# return all the terms of weight greater than wght
#
gpolyutil[get_weight_greaterthan]:=proc(wght::integer,p::polynom,genus::integer)
  local symp, tempval, retval;
  if (genus <> 3) then ERROR(`only genus 3 implemented - weight function`) fi;
  symp:=get_orbit(p,genus);
  tempval := 0;
  retval := 0;
  while (g3weight(symp) > wght) do
    tempval := get_orbit(get_heaviest_terms(symp,genus),genus);
    retval := retval + tempval;
    symp := simplify(symp-tempval);
  od;
  retval;
end: