load("power-6.lisp");
lt(poly):=block([rlt:lterm(poly),glt],
   glt:ratdisrep(rlt),glt/numfactor(glt));
lc(poly):=block([rp:rat(poly),rlt:lterm(poly),glt],
   glt:ratdisrep(rlt),numfactor(glt)/ratdenom(rp));
p1:x1*x4+x3-x1*x2;
p2:2*x4^2-2*x3*x4+5*x1*x2*x4-5*x1*x2*x3;
g:x1*x4^2+x4^2-x1*x2*x4-x2*x4+x1*x2+3*x2;
h1:x4^2-x3*x4-x2*x4+x1*x2+3*x2;
/* This is a quick hack to check if a division produces a negitive power */ 
NegitivePowerP(expr):=block([lov:listofvars(expr),flag:false],
   for v in lov do (
      if lopow(expr,v) < 0 then return(flag:true)
   ),
   flag
)$
/* This function subtracts a possibly non constant multiple of p from f,
   and returns the difference, and a flag.
   If more subtraction is possible, then return true, otherwise return false */
OneReduction(f,p):=block(
   [leadingTermOfP:lt(p),leadingCoeffOfP:lc(p),
    leadingTermOfF:lt(f),leadingCoeffOfF:lc(f)],
   /* no subtraction was possible */
   if orderlessp(leadingTermOfF,leadingTermOfP) then ([f,false])
   else ( /* here we may have more to subtract */
      if ordergreatp(leadingTermOfF,leadingTermOfP) then (
         if NegitivePowerP(leadingTermOfF/leadingTermOfP) then
             [f,false]
	 else
         [expand(f - (
            (leadingCoeffOfF/leadingCoeffOfP)*
               (leadingTermOfF/leadingTermOfP)*p)),true]
      )
      else ( /* monomials are the same, so dont neet to call this again */
         [expand(f - ((leadingCoeffOfF/leadingCoeffOfP)*p)),false]
      )
   )
)$

NormalForm(f,_ps):=block(
   [ps:sort(_ps,lambda([x,y],if ordergreatp(lt(x),lt(y)) then true else false)),
    last,flag],
   for p in ps do (
      flag:true,
      while flag do (
         last:OneReduction(f,p),
	 f:last[1],
	 flag:last[2]
      )
   ),
   f
)$
trace(NormalForm,OneReduction);
NormalForm(g,[p1,p2]);