######################
# This is code implements the algorithm to compute the canonical basis for gl(n|m) as
# described in my survey article "Representations of the general linear Lie superalgebra in the BGG category O."
#
# Load it up using GAP4. There shouldn't be any errors.
# Once successfully loaded  you set global variables n and m. It defaults to 
#
# gap> n:=2;m:=2;
#
# Then for example you call the function 
#
# gap> B([1,2,2,1]); 
#
# to compute the canonical basis b_{(1,2|2,1)}. Note this is the first example that is not multiplicity free.
# Obviously the code will only work for pretty small n and m... then it will start running out of memory or taking too long.
#
# --Jon Brundan.
######################

#######################################################################################################################################
#######################################################################################################################################

if not IsBound(IsMvpObj) then

# This code uses some library functions which I have embedded here. This is a bit gratuitous and time-wasteful
# but it is my standard technique for implementing such algorithms...

DeclareCategory("IsMvpObj", IsScalar);
DeclareCategoryCollections("IsMvpObj");
DeclareRepresentation("IsMvpRep", IsComponentObjectRep, ["coeff", "term", "normalize"]);
DeclareGlobalFunction("Mvp");


MvpStringToNumber := function(s) local n,i,ten,j;
  n := 0;
  ten := 1;
  for i in Reversed([1..Length(s)]) do
    for j in [0..9] do
      if s[i] = String(j)[1] then n := n + j*ten;fi;
    od;
    if s[i] = '-' then n := -n;fi;
    ten := ten*10;
  od;
  return(n);
end;

MvpUnpack := function(string) local i,j;
  return List(SplitString(string, ".,;|", "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz()[] "), MvpStringToNumber);
end;

MvpSet := function(e,c,s) local R,F,i,j,k,b,ee,eee,cc,new,uppers,lowers;
  eee := [];
  for i in [1..Length(e)] do
    uppers := []; lowers := [];
    for j in [1..Length(e[i])] do
      if e[i][j][1][1] < 'a' or e[i][j][1][1] >'z' then
        Add(uppers, ShallowCopy(e[i][j]));
      else
        Add(lowers, ShallowCopy(e[i][j]));
      fi;
    od;
    Sort(lowers); new := [];
    for k in [1..Length(lowers)] do
      if k=1 or lowers[k][1]<>lowers[k-1][1] then Add(new,lowers[k]);
      else new[Length(new)][2]:=new[Length(new)][2]+lowers[k][2];
      fi;
    od;
    for k in [1..Length(uppers)] do
      if k=1 or uppers[k][1]<>uppers[k-1][1] then Add(new,uppers[k]);
      else new[Length(new)][2]:=new[Length(new)][2]+uppers[k][2];
      fi;
    od;
    new:=Filtered(new,x->x[2] <> 0*x[2]);
    Add(eee, new);
  od;
  SortParallel(eee,c);
  ee := []; cc := [];
  for i in [1..Length(eee)] do
    if i=1 or eee[i]<>eee[i-1] then
      Add(ee,eee[i]);Add(cc,c[i]);
    else cc[Length(cc)]:=cc[Length(cc)]+c[i];
    fi;
  od;
  b:=List(cc,x->x <> 0*x);
  R := rec(elm:=Immutable(ListBlist(ee,b)),coeff:=Immutable(ListBlist(cc,b)),normalize:=s);
  F := NewFamily("IsMvpObj", IsMvpObj);
  return Objectify(NewType(F, IsMvpObj and IsMvpRep), R);
end;

MvpNormalize := function() end;

MvpNormalizeSub := function(thise,thiscf,normalizefunc) local i,j,lowers,uppers,flip;
  if normalizefunc <> fail then 
    lowers := []; uppers := [];
    for i in [1..Length(thise)] do 
      if thise[i][1][1] < 'a' or thise[i][1][1] > 'z' then
        Add(uppers, thise[i]);
      else
        Add(lowers,thise[i]);
      fi;
    od;
    flip := normalizefunc(uppers);
    if flip <> fail then
      return MvpNormalize(MvpSet([lowers], [thiscf], normalizefunc) * flip);
    fi;
  fi;
  return MvpSet([thise],[thiscf], normalizefunc);
end;

MvpNormalize := function(x) local normalized,i;
  if IsCyc(x) or IsRationalFunction(x) then return(x);fi;
  normalized := Zero(x);
  for i in [1..Length(x!.elm)] do
    if x!.coeff[i] <> 0*x!.coeff[i] then
      normalized := normalized + MvpNormalizeSub(x!.elm[i],x!.coeff[i],x!.normalize);
    fi;
  od;
  return(normalized);
end;

# Useful routines
# The first truncates all powers of var in pol > deg

MvpTruncateAbove := function(poly,deg,var) local e,c,i,j,valid,t;
  if IsCyc(poly) or IsRationalFunction(poly) then
    if deg < 0 then return Zero(poly);fi;
    return poly;
  fi;

  c := []; e := [];
  for i in [1..Length(poly!.coeff)] do
    valid := true;
    if poly!.elm[i] = [] then
       if deg < 0 then valid := false;fi;
    else
      for t in poly!.elm[i] do
        if t[1] = var and t[2] > deg then valid := false; fi;
      od;
    fi;
    if valid then
      Add(c, poly!.coeff[i]);
      Add(e, ShallowCopy(poly!.elm[i]));
    fi;
  od;
  return MvpSet(e,c,poly!.normalize);
end;

MvpTruncateBelow := function(poly,deg,var) local e,c,i,j,valid,t;
  if IsCyc(poly) or IsRationalFunction(poly) then
    if deg > 0 then return Zero(poly);fi;
    return poly;
  fi;
  c := []; e := [];
  for i in [1..Length(poly!.coeff)] do
    valid := true;
    if poly!.elm[i] = [] then
       if deg > 0 then valid := false;fi;
    else
      for t in poly!.elm[i] do
        if t[1] = var and t[2] < deg then valid := false; fi;
      od;
    fi;
    if valid then
      Add(c, poly!.coeff[i]);
      Add(e, ShallowCopy(poly!.elm[i]));
    fi;
  od;
  return MvpSet(e,c,poly!.normalize);
end;

# Return the deg coefficient of a polynomial
# which is of course usually another mvp.
# However if it is truly an integer it returns it as such!

MvpCoefficient := function(poly,var,deg) local i,j,thisdeg,c,e,t,newe;
  if IsCyc(poly) or IsRationalFunction(poly) then
    if deg = 0 then return poly;fi;
    return 0;
  fi;
  c := []; e := [];
  for i in [1..Length(poly!.coeff)] do
    thisdeg := 0;
    for t in poly!.elm[i] do
      if t[1] = var then thisdeg := t[2]; fi;
    od;
    if thisdeg = deg then
      newe := [];
      for t in poly!.elm[i] do
        if t[1] <> var then Add(newe, ShallowCopy(t)); fi;
      od;
      Add(e, newe);
      Add(c, poly!.coeff[i]);
    fi;
  od;
  if e = [[]] then return c[1];fi;
  return MvpSet(e,c,poly!.normalize);
end;

# Change all degrees to their negatives

MvpBar := function(poly, var) local c,e,i,j,this;
  if IsCyc(poly) or IsRationalFunction(poly) then return(poly);fi;

  c := ShallowCopy(poly!.coeff); 
  e := [];
  for i in [1..Length(poly!.elm)] do
    this := [];
    for j in [1..Length(poly!.elm[i])] do
      if poly!.elm[i][j][1] = var then
        Add(this, [poly!.elm[i][j][1], -poly!.elm[i][j][2]]);
      else
        Add(this, [poly!.elm[i][j][1], poly!.elm[i][j][2]]);
      fi;
    od;
    Add(e,this);
  od;
  return MvpSet(e,c,poly!.normalize);
end;  

# Returns a list [bot,top] giving the bottom and top degrees of
# variable var in poly. If it involves no var's at all, returns false.

MvpRange := function(poly,var) local bot,top,i,j,thisdeg;
  if IsCyc(poly) or IsRationalFunction(poly) then
    return false;
  fi;

  bot := false;
  top := false;
  for i in [1..Length(poly!.elm)] do
    thisdeg := false;
    for j in [1..Length(poly!.elm[i])] do
      if poly!.elm[i][j][1] = var then 
        if thisdeg = false then thisdeg := 0;fi;
        thisdeg := thisdeg + poly!.elm[i][j][2];
      fi;
    od;
    if thisdeg <> false then
      if bot = false or thisdeg < bot then bot := thisdeg;fi;
      if top = false or thisdeg > top then top := thisdeg;fi;
    fi;
  od;
  if bot = false then return false;fi;
  return [bot,top];
end;

# Here are the methods

InstallGlobalFunction(Mvp, 
  "mvp",
  function(arg) local newelm;
    if not IsString(arg[1]) or arg[1] = [] then
      newelm := [arg[1]];
    else
      newelm := [[[arg[1],1]]];
    fi;
    if Length(arg) = 1 then 
      return MvpSet(newelm,[1],fail);
    fi;
    if Length(arg) = 2 then
      if IsFunction(arg[2]) then
        return MvpSet(newelm,[1],arg[2]);
      fi;
        return MvpSet(newelm, [arg[2]],fail);
    fi;
    if Length(arg) = 3 then
      return MvpSet(newelm, [arg[2]], arg[3]);
    fi;
  end
);

InstallMethod(PrintObj,
  "mvp",
  [IsMvpObj and IsMvpRep], 
  function(h) local i,j,coeff,elm,s,res;
    res:="";
    if h!.elm=[] then Append(res,"0");fi;
    for i in [1..Length(h!.elm)] do
      elm:="";
      for s in h!.elm[i] do Append(elm,s[1]);
        if s[2]<>1 then Append(elm,"^");Append(elm,String(s[2]));fi;
      od;
      coeff:=String(h!.coeff[i]);
      if elm<>"" then
        if coeff="1" then coeff:="";
	elif coeff="-1" then coeff:="-";
	fi;
      fi;
      if '-' in coeff{[2..Length(coeff)]} 
          or '+' in coeff{[2..Length(coeff)]} then
        coeff:=Concatenation("(",coeff,")");
      fi;
      if Position(coeff,'-')<>1 and i<>1 then Append(res,"+");fi;
      Append(res,coeff);
      Append(res,elm);
    od;
    Print(String(res));
  end 
);

InstallMethod(\=,
  "mvp=mvp",
  [IsMvpObj and IsMvpRep, IsMvpObj and IsMvpRep],
  function( x, y ) local diff;
    diff := MvpNormalize(x-y);
    return(diff!.elm = []);
  end
);

InstallMethod(\=,
  "mvp=cyc",
  [IsMvpObj and IsMvpRep, IsCyc],
  function( x, y ) local diff;
    diff := MvpNormalize(x-y);
    return(diff!.elm = []);
  end
);

InstallMethod(\=,
  "mvp=cyc",
  [IsMvpObj and IsMvpRep, IsRationalFunction],
  function( x, y ) local diff;
    diff := MvpNormalize(x-y);
    return(diff!.elm = []);
  end
);

InstallMethod(\=,
  "cyc=mvp",
  [IsCyc,IsMvpObj and IsMvpRep],
  function( x, y ) 
    return y=x;
  end
);

InstallMethod(\=,
  "cyc=mvp",
  [IsRationalFunction,IsMvpObj and IsMvpRep],
  function( x, y ) 
    return y=x;
  end
);

InstallMethod(\+,
  "mvp+mvp",
  [IsMvpObj and IsMvpRep, IsMvpObj and IsMvpRep],
  function( x, y )
    return MvpSet(Concatenation(x!.elm,y!.elm),Concatenation(x!.coeff,y!.coeff),x!.normalize);
  end
);

InstallMethod(\+,
  "mvp+cyc",
  [IsMvpObj and IsMvpRep, IsCyc],
  function( x, y )
    return MvpSet(Concatenation(x!.elm,[[]]),Concatenation(x!.coeff,[y]),x!.normalize);
  end
);

InstallMethod(\+,
  "mvp+cyc",
  [IsMvpObj and IsMvpRep, IsRationalFunction],
  function( x, y )
    return MvpSet(Concatenation(x!.elm,[[]]),Concatenation(x!.coeff,[y]),x!.normalize);
  end
);

InstallMethod(\+,
  "cyc+mvp",
  [IsCyc, IsMvpObj and IsMvpRep],
  function( x, y )
    return y+x;
  end
);

InstallMethod(\+,
  "cyc+mvp",
  [IsRationalFunction, IsMvpObj and IsMvpRep],
  function( x, y )
    return y+x;
  end
);

InstallMethod(ZeroOp,
  "mvp",
  [IsMvpObj],
  x -> MvpSet([[]],[0],x!.normalize));

InstallMethod(AdditiveInverseOp,
  "mvp",
  [IsMvpObj and IsMvpRep],
  x -> MvpSet(ShallowCopy(x!.elm),-ShallowCopy(x!.coeff),x!.normalize));

InstallMethod(\*,
  "mvp*mvp",
  [IsMvpObj and IsMvpRep, IsMvpObj and IsMvpRep],
  function( x, y ) local relm,rcoeff,i,j;
    relm:=[];rcoeff:=[];
    for i in [1..Length(x!.elm)] do
      for j in [1..Length(y!.elm)] do
	Add(relm,Concatenation(x!.elm[i],y!.elm[j]));
        Add(rcoeff,x!.coeff[i]*y!.coeff[j]);
      od;
    od;
    return MvpSet(relm,rcoeff,x!.normalize);
  end
);

InstallMethod(\*,
  "mvp*cyc",
  [IsMvpObj and IsMvpRep, IsCyc],
  function( x, y ) local relm,rcoeff,i;
    relm:=[];rcoeff:=[];
    for i in [1..Length(x!.elm)] do
      Add(relm,x!.elm[i]);
      Add(rcoeff,x!.coeff[i]*y);
    od;
    return MvpSet(relm,rcoeff,x!.normalize);
  end
);

InstallMethod(\*,
  "mvp*cyc",
  [IsMvpObj and IsMvpRep, IsRationalFunction],
  function( x, y ) local relm,rcoeff,i;
    relm:=[];rcoeff:=[];
    for i in [1..Length(x!.elm)] do
      Add(relm,x!.elm[i]);
      Add(rcoeff,x!.coeff[i]*y);
    od;
    return MvpSet(relm,rcoeff,x!.normalize);
  end
);

InstallMethod(\*,
  "cyc*mvp",
  [IsCyc,IsMvpObj and IsMvpRep],
  function( x, y ) 
    return y*x;
  end
);

InstallMethod(\*,
  "cyc*mvp",
  [IsRationalFunction,IsMvpObj and IsMvpRep],
  function( x, y ) 
    return y*x;
  end
);



InstallMethod(OneOp,
  "mvp",
  [IsMvpObj],
  x -> MvpSet([[]],[1],x!.normalize));

InstallMethod(InverseOp,
  "mvp",
  [IsMvpObj and IsMvpRep],
  function(x) local h,newelm;
    if Length(x!.elm) <> 1 or x!.coeff[1] = 0 then return fail; fi;
    newelm := [];
    for h in [1..Length(x!.elm[1])] do
      Add(newelm, ShallowCopy(x!.elm[1][h]));
      newelm[Length(newelm)][2] := -newelm[Length(newelm)][2];
    od;
    return MvpSet([Reversed(newelm)], [x!.coeff[1]^-1],x!.normalize);
  end
);

fi;

#######################################################################################################################################
#######################################################################################################################################
# That is the end of the library stuff, here is the code proper...

if not IsBound(n) then  
  n := 2; m := 2;
  Indeterminate(Rationals,"q");
  q := LaurentPolynomialByCoefficients(FamilyObj(1), [1], 1);
fi;

# The basic data type is a term, which is an n|m-tuple, and then a vec, which
# is a linear combination of V's stored in an MVP

V := function(term) local string,i;
  string := "V(";
  for i in [1..n] do
    Append(string, String(term[i]));
    if i < n then Append(string,",");fi;
  od;
  if m <> 0 and n <> 0 then Append(string,"|");fi;
  for i in [n+1..n+m] do
    Append(string, String(term[i]));
    if i < n+m then Append(string, ",");fi;
  od;
  Append(string, ")");
  return q^0*Mvp(string);
end;

###########################
# This is the Bruhat ordering
# Returns true if term1 >= term2

Dominates := function(term1,term2) local i,psi1,psi2,a;
  for a in [Minimum(Minimum(term1), Minimum(term2))..Maximum(Maximum(term1), Maximum(term2))] do
    psi1 := 0; psi2 := 0;
    for i in [1..n] do
      if term1[i] >= a then
        psi1 := psi1 + 1;
      fi;
      if term2[i] >= a then
        psi2 := psi2 + 1;
      fi;
      if psi1 < psi2 then return(false);fi;
    od;
    for i in [1..m] do
      if term1[n+i] >= a then
        psi1 := psi1 - 1;
      fi;
      if term2[n+i] >= a then
        psi2 := psi2 - 1;
      fi;
      if psi1 < psi2 then return(false);fi;
    od;
    if psi1 <> psi2 then return(false);fi;
  od;
  return(true);
end;

###########################
# Here is the correction routine
# You pass it the leading guy which must indeed be leading with coefficient 1

B := function() end;	      # Forward reference

Correct := function(vec,leading) local nextleading,p,leadingok,nextcf,t,poly,i,this;

# Scan through all the terms.
# Check that the leading coefficient is 1, and find the 
# next lowest term in dominance that involves 1+(stuff) as its coeff.

  nextleading := false;
  leadingok := false;
  for t in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[t][1][1]);
    if not Dominates(this, leading) then Assert(0,"Not lower!");fi;
    if this = leading then
      if vec!.coeff[t] <> q^0 then Assert(0,"Bad leading coefficient!");fi;
      leadingok := true;
    elif nextleading = false or (this <> nextleading and Dominates(nextleading,this)) then
      poly := CoefficientsOfLaurentPolynomial(vec!.coeff[t]);	# See if it has coefficients in constant or negative powers of q
      if poly[2] <= 0 then
        nextleading := this;
        nextcf := poly;
      fi;
    fi;
  od;

  if leadingok = false then Assert(0,"Wrong leading term!");fi;
  if nextleading = false then return(vec);fi;

  p := 0*q;
  for i in [nextcf[2]..0] do
    if i < 0 and IsBound(nextcf[1][i+1-nextcf[2]]) then
      p := p + (q^i + q^(-i))*(nextcf[1][i+1-nextcf[2]]);
    fi;
    if i = 0 and IsBound(nextcf[1][i+1-nextcf[2]]) then
      p := p + q^0*(nextcf[1][i+1-nextcf[2]]);
    fi;
  od;
  return Correct(vec - p*B(nextleading),leading);
end;

###########################
# You can act on vecs...

ActC := function(vec,i) local ans,this,j,new;
  ans := 0;
  for j in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[j][1][1]);
    if i >= 1 and i < n then
      if this[i] = this[i+1] then
        ans := ans + vec!.coeff[j]*V(this)*(q+q^(-1));
      elif this[i] < this[i+1] then
        new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i];
	ans := ans + vec!.coeff[j]*V(new) + vec!.coeff[j]*V(this)*q^(-1);
      else
        new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i];
	ans := ans + vec!.coeff[j]*V(new) + vec!.coeff[j]*V(this)*q;
      fi;
    elif i > n and i <= n+m-1 then
      if this[i] = this[i+1] then
        ans := ans + vec!.coeff[j]*V(this)*(q+q^(-1));
      elif this[i] < this[i+1] then
        new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i];
	ans := ans + vec!.coeff[j]*V(new) + vec!.coeff[j]*V(this)*q;
      else
        new := ShallowCopy(this); new[i] := this[i+1]; new[i+1] := this[i];
	ans := ans + vec!.coeff[j]*V(new) + vec!.coeff[j]*V(this)*q^(-1);
      fi;
    else Assert(0,"Bad index!");fi;
  od;
  return ans;
end;

# This acts by K_i in the rth position

ActKp := function(vec,i,r,power) local ans,j,cf,this;
  if vec = 0*q then return 0*q;fi;
  ans := 0*q;
  for j in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[j][1][1]);
    cf := vec!.coeff[j];
    if r <= n then
      if this[r]=i then cf := cf * (q^power);fi;
      if this[r]=i+1 then cf := cf / (q^power);fi;
    else
      if this[r]=i then cf := cf / (q^power);fi;
      if this[r]=i+1 then cf := cf * (q^power);fi;
    fi;
    ans := ans + cf * V(this);
  od;
  return ans;
end;

# This acts by K_i^{-1} on the rth position

ActKm := function(vec,i,r,power) local ans,j,cf,this;
  if vec = 0*q then return 0*q;fi;
  ans := 0*q;
  for j in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[j][1][1]);
    cf := vec!.coeff[j];
    if r <= n then
      if this[r]=i then cf := cf / (q^power);fi;
      if this[r]=i+1 then cf := cf * (q^power);fi;
    else
      if this[r]=i then cf := cf * (q^power);fi;
      if this[r]=i+1 then cf := cf / (q^power);fi;
    fi;
    ans := ans + cf * V(this);
  od;
  return ans;
end;

 # This acts by E_i on the rth position

ActE := function(vec,i,r) local ans,j,this;
  if vec = 0*q then return 0*q;fi;
  ans := 0*q;
  for j in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[j][1][1]);
    if r <= n and this[r] = i+1 then
      this[r] := i;
      ans := ans + vec!.coeff[j]*V(this);
    elif r >= n+1 and this[r] = i then
      this[r] := i+1;
      ans := ans + vec!.coeff[j]*V(this);
    fi;
  od;
  return ans;
end;

# This acts by F_i on the rth position

ActF := function(vec,i,r) local ans,j,this;
  if vec = 0*q then return 0*q;fi;
  ans := 0*q;
  for j in [1..Length(vec!.coeff)] do
    this := MvpUnpack(vec!.elm[j][1][1]);
    if r <= n and this[r] = i then
      this[r] := i+1;
      ans := ans + vec!.coeff[j]*V(this);
    elif r >= n+1 and this[r] = i+1 then
      this[r] := i;
      ans := ans + vec!.coeff[j]*V(this);
    fi;
  od;
  return ans;
end;

# Given a vector in T^{n|m} and a list of r positions, this acts
# K_i^{-r}\otimes...\otimes K_i^{-r}\otimes q^r K_i^{-r} E_i\otimes K_i^{1-r}\otimes ...\otimes K_i^{1-r}\otimes q^{r-1} K_i^{1-r} E_i
# \otimes...\otimes E_i \otimes 1 \otimes... \otimes 1 

ActMultipleE := function(vec,i,positions) local mult,r,j;
  if vec = 0*q then return 0*q;fi;
  mult := 0;
  for r in Reversed([1..n+m]) do
    if r in positions then
      vec := ActE(vec,i,r);
      vec := ActKm((q^mult)*vec,i,r,mult);
      mult := mult + 1;
    else
      vec := ActKm(vec,i,r,mult);
    fi;
  od;
  return(vec);
end;

ActMultipleF := function(vec,i,positions) local mult,r,j;
  if vec = 0*q then return 0*q;fi;
  mult := 0;
  for r in [1..n+m] do
    if r in positions then
      vec := ActF(vec,i,r);
      vec := ActKp((q^mult)*vec,i,r,mult);
      mult := mult + 1;
    else
      vec := ActKp(vec,i,r,mult);
    fi;
  od;
  return(vec);
end;

ActDividedE := function(vec, i, s) local v,positions;
  if vec = 0*q then return 0*q;fi;
  v := 0*q;
  for positions in Combinations([1..n+m],s) do
    v := v+ActMultipleE(vec,i,positions);
  od;
  return(v);
end;

ActDividedF := function(vec, i, s) local v,positions;
  if vec = 0*q then return 0*q;fi;
  v := 0*q;
  for positions in Combinations([1..n+m],s) do
    v := v+ActMultipleF(vec, i,positions);
  od;
  return(v);
end;

############################
# This is a step in the bumping algorithm

# Given a word thisside and i, it finds i...j all consecutive and adds one to them all,
# tacking on [i,side*mult of i],...,[j,side*mult of j] to monomial and returning j+1

Bump := function(thisside,side,i,monomial) local oldside,j,mult,this;
  j := i;
  oldside := ShallowCopy(thisside);
  while j in thisside do j := j+1;od;
  while i < j do
    mult := 0;
    for this in [1..Length(oldside)] do 
      if oldside[this] = i then thisside[this]:=thisside[this]+1;mult:=mult+1;fi;
    od;
    Add(monomial,[i,side*mult]);
    i := i+1;
  od;
  return j;
end;

############################
# This is the main routine

startside := 1;

B := function(term) local i,new,intersect,top,side,thisside,left,right,vec,monomial;

# Find biggest match

  left := []; right := [];
  for i in [1..n] do
    Add(left,term[i]);
  od;
  for i in [n+1..n+m] do
    Add(right,term[i]);
  od;
  intersect := Intersection(left,right);
  if intersect <> [] then
    top := Maximum(intersect);
    monomial := [];	# This is stored as a list of [i,m] where m > 0 is multiplicitiy of e and m < 0 is minus multiplicicity of f
    side := startside;	# 1 for left hand side, -1 for right hand side
    if side = 1 then thisside := left; else thisside := right;fi;
    while top in thisside do
      top := Bump(thisside,side,top,monomial);
      side := -side;
      if side = 1 then thisside := left; else thisside := right;fi;
    od; 
    Append(left,right);
    vec := B(left);
    for i in monomial do
      if i[2] > 0 then
        vec := ActDividedE(vec,i[1],i[2]);
      else
        vec := ActDividedF(vec,i[1],-i[2]);
      fi;
    od;
    return Correct(vec,term);
  fi;

# Here if typical

  for i in [1..n-1] do
    if term[i] < term[i+1] then
      new := ShallowCopy(term);
      new[i] := term[i+1]; new[i+1] := term[i];
      return Correct(ActC(B(new),i), term);
    fi;
  od;
  for i in [n+1..n+m-1] do
    if term[i] > term[i+1] then
      new := ShallowCopy(term);
      new[i] := term[i+1]; new[i+1] := term[i];
      return Correct(ActC(B(new),i), term);
    fi;
  od;

# Here if typical and weakly dominant

  return V(term);
end;