# This is a general implementation of the Groebner walk from
# the grevlex order to an elimination order as in Collart,
# Kalkbrener, and Mall -- JSC v. 24 (1997), 465-469.  Aim to
# apply this to the calculation of "bisector surfaces" and other
# computations in applied geometry, as in Elber and Kim -- "The bisector 
# surface of rational space curves" ACM Transactions on Graphics, v. 17, 
# number 1 (January 1998).)  
#
#
#  some standard definitions
#
with(linalg):
with(Groebner):
initialize:=proc(elimvars,leftvars)
   
   global A,grevlex,elim,vars,RestOfMat,grevlexMat;
   local i,j,firstrow,nextrow;

   vars:=[op(elimvars),op(leftvars)];  
   with(Ore_algebra):
   A:=poly_algebra(op(elimvars),op(leftvars)):
   grevlex:=termorder(A,tdeg(op(vars))):
   firstrow:=[];
   for i to nops(vars) do
     firstrow:=[op(firstrow),1];
   od;
   grevlexMat:=[firstrow];
   for i from nops(vars) by -1 to 2 do
     nextrow:=[];
     for j to i-1 do
       nextrow:=[op(nextrow),0]
     od;
     nextrow:=[op(nextrow),-1];
     for j from i+1 to nops(vars) do 
       nextrow:=[op(nextrow),0];
     od;
     grevlexMat:=[op(grevlexMat),nextrow];
   od;
#   elim:=termorder(A,lexdeg(elimvars,leftvars)):
   RestOfMat:=[];
   for i to nops(elimvars) do
     RestOfMat:=[op(RestOfMat),0];
   od;
   for i to nops(leftvars) do
     RestOfMat:=[op(RestOfMat),1];
   od;
   for i from nops(leftvars) by -1 to 2 do
     nextrow:=[];
     for j to nops(elimvars) do 
       nextrow:=[op(nextrow),0];
     od;
     for j to i-1 do
       nextrow:=[op(nextrow),0]
     od;
     nextrow:=[op(nextrow),-1];
     for j from i+1 to nops(leftvars) do 
       nextrow:=[op(nextrow),0];
     od;
     RestOfMat:=RestOfMat,nextrow;
   od;
   for i from nops(elimvars) by -1 to 2 do
     nextrow:=[];
     for j to i-1 do
       nextrow:=[op(nextrow),0]
     od;
     nextrow:=[op(nextrow),-1];
     for j from i+1 to nops(elimvars) do 
       nextrow:=[op(nextrow),0];
     od;
     for j to nops(leftvars) do 
       nextrow:=[op(nextrow),0];
     od;
     RestOfMat:=RestOfMat,nextrow;
   od;
   RestOfMat:=[RestOfMat];   
   end:
#
# compute the initial form of a polynomial with
# respect to a given weight vector
#
InForm:=proc(p,vars,WeightVector)
  
  local i,j,term,var,InF,weight,topweight,poly;

  poly:=p;
  InF := 0;
  topweight := 0;
  while poly <> 0 do
    term := leadmon(poly,grevlex);
    weight := 0;
    for j to nops(vars) do
       weight := weight + WeightVector[j]*degree(term[2],vars[j]);
    od;
    if weight > topweight then
       topweight := weight;
       InF := term[1]*term[2];
    elif weight = topweight then
       InF := InF + term[1]*term[2];
    fi;
    poly := poly - term[1]*term[2];
  od;
  RETURN(InF)
  end:
#
# compute the vector of exponents of a monomial
#
GetExpVector := proc(term)

  local i,vec;
  
  vec:=[];
  for i to nops(vars) do 
    vec:=[op(vec),0];
  od;
  for i to nops(vars) do
    vec[i] := degree(term,vars[i]);
  od;
  RETURN(vec)
  end:
#
#   
# Determine a weight vector which represents a given
# order on a given ideal.  G should be a Groebner basis
# with respect to that order.  This version finds
# the most "economical" epsilon.  (The motivation for
# is the limitation that weight vectors must be 
# integer vectors for the Groebner package(!), so big
# denominators such as those produced by GetWeight2 
# caused big problems.) 
#
GetWeight3:=proc(G,TO,TOmat)

   global vars;
   local n,i,j,epsilon,maxtdeg,rowmax,total,vec,p,q,r,Done,SmallEnough;
   
   print(`in GetWeight3`);
   print(TOmat);
   maxtdeg := degree(G[1]);
   for i from 2 to nops(G) do
     if degree(G[i]) > maxtdeg then
        maxtdeg := degree(G[i]);
     fi;
   od;
   total := 1;
   for i from 1 to nops(vars) do
     rowmax := TOmat[i,1];
     for j from 2 to nops(vars) do 
       if abs(TOmat[i,j]) > rowmax then
          rowmax := abs(TOmat[i,j]);
       fi;
       total := total + rowmax;
       od;
     od;
#   print(total*maxtdeg);
   n:=2;
   Done := false;
   while evalb((n <= total*maxtdeg) and (Done = false)) do
     SmallEnough := true;
     epsilon:=1/n;
     vec:=TOmat[1];
     for i from 2 to nops(vars) do
        vec:=vec + epsilon^(i-1)*TOmat[i];
     od;
#     print(vec);
     j:=1;
     while j <= nops(G) do  
#        print(j);
        q:= leadmon(G[j],TO); 
#        print(q);
        r:= InForm(G[j],vars,vec);
#        print(r);
        if evalb(simplify(r - q[1]*q[2]) <> 0) then 
           SmallEnough := false;
           break;
        fi;
        j:=j+1; 
     od;
     if SmallEnough then
       RETURN(eval(vec));
       Done := true;
     else 
       print(n,vec);
       n := n+1;
     fi;
   od;
   end:
    
#
# test output from previous procedure(s)
# by comparing initial forms and leading terms
#             
TestVec := proc(G,vec,TO)
  
  local q,r,j;

  for j to nops(G) do 
     q:= leadmon(G[j],TO); 
     r:= InForm(G[j],vars,vec);
     if evalb(simplify(r - q[1]*q[2]) <> 0) then 
       print(j,q,r) fi; 
     od;
  end: 
#
#  Find first time we leave current Groebner cone
#  along a particular path -- G is a Groebner basis
#  with respect to TO
#
LeaveCone2 := proc(G,weight,TO)

  local i,p,pleft,lead,leadvec,term,termvec,firsttime,ttime,el,w,dif;
 
  #print(`in LeaveCone2`);
  firsttime := 1;
  for p in G do
    lead := leadmon(p,TO);
    leadvec:=GetExpVector(lead[2]);
    pleft := simplify(p - lead[1]*lead[2]);
    while pleft <> 0 do
       term := leadmon(pleft,TO);
       termvec:=GetExpVector(term[2]);
       dif:=leadvec-termvec;
       w:=0;
       for i to nops(vars) do 
          w:=w+dif[i]*weight[i];
       od;
       el:=0;
       for i to nops(elimvars) do
         el:=el+dif[i];
       od;
       if w <> el then
         ttime:=w/(w-el);
   #      print(ttime);
         if (ttime > 0) and (ttime < firsttime) then
           firsttime:=ttime;
         fi;
       fi;
       pleft := simplify(pleft - term[1]*term[2]);
       od;
    od;
  RETURN( firsttime )
  end:
#
#  Division algorithm with quotients:
# 
#  Alternate implementation of:
#
#  Division algorithm with quotients:
#
DivQuots:=proc(f,G,TO)

  local quots,rem,i,j,p,s,divocc,lmG,lcG,lmp,lmpe,q;

  s:=nops(G);
  lmG := [];
  lcG := [];
  quots := [];
  for i to s do
    lmG := [op(lmG),leadterm(G[i],TO)];
    lcG := [op(lcG),leadcoeff(G[i],TO)];
    quots := [op(quots),0];
  od;
  rem := 0;
  p := f;
  while p <> 0 do
    i := 1;
    divocc:=false;
    lmp := leadmon(p,TO);
    while i <= s do
        if denom(lmp[2]/lmG[i]) = 1 then
           q := simplify(lmp[1]*lmp[2]/(lcG[i]*lmG[i]));
           quots[i] := quots[i] + q;
           p := simplify(p - q*G[i]);
           divocc:=true;
           break;
        else
           i := i+1;
        fi;
    od;
    if divocc = false then
       rem := rem + lmp[1]*lmp[2];
       p := simplify(p - lmp[1]*lmp[2]);
    fi;
  od;
  RETURN(quots,rem)
  end:     
#
# Special purpose implementation of division algorithm;
# list of divisors is segregated into separate lists of 
# monomials and polynomials
#
DivMP:=proc(f,Mons,Polys,TO)

  local rem,i,j,p,sM,sP,divocc,lmP,lcP,lmp,q;

  sM:=nops(Mons); sP:=nops(Polys);
  lmP := [];
  lcP := [];
  for i to sP do
    lmP := [op(lmP),leadterm(Polys[i],TO)];
    lcP := [op(lcP),leadcoeff(Polys[i],TO)];
  od;
  rem := 0;
  p := f;
  while p <> 0 do
    divocc:=false;
    lmp := leadmon(p,TO);
    i:=1;
    while i <= sM do
        if denom(op(1,Mons[i])*lmp[2]/Mons[i]) = 1 then
           p:=simplify(p-lmp[1]*lmp[2]);
           divocc:=true;
           break;
        else
           i := i+1;
        fi;
    od;
    if divocc = false then
      i := 1;
      while i <= sP do
        if denom(lmp[2]/lmP[i]) = 1 then
           q := simplify(lmp[1]*lmp[2]/(lcP[i]*lmP[i]));
           p := simplify(p - q*Polys[i]);
           divocc:=true;
           break;
        else
           i := i+1;
        fi;
      od;
    fi; 
    if divocc = false then
       rem := rem + lmp[1]*lmp[2];
       p := simplify(p - lmp[1]*lmp[2]);
    fi;
  od;
  RETURN(rem)
  end:     
#
# Special purpose implementation of Buchberger's Algorithm
# tailored for the special case encountered in the Walk:
# some (most?) of the polynomials actually monomials; all with 
# few terms.

  SPBuchberger:=proc(F,TO)

    local G,B,m,p,i,j,LT1,LT2,S,Mons,Polys;
    
    
    G:=F;
    B:=[];
    Mons:=[]; Polys:=[];
    for i to nops(G)-1 do                                
      if type(G[i],`+`)=true then
        Polys:=[op(Polys),G[i]];
        for j from i+1 to nops(G) do 
          B:=[op(B),[i,j]]; 
        od;
      else
        Mons:=[op(Mons),G[i]];
        for j from i+1 to nops(G) do
          if type(G[j],`+`) = true then
            B:=[op(B),[i,j]];
          fi;                                                      
        od;
      fi;                                                        
    od;
    if type(G[nops(G)],`+`) = true then
      Mons:=[op(Mons),G[nops(G)]];
    else
      Polys:=[op(Polys),G[nops(G)]];
    fi;
#   print(nops(B));
    while B <> [] do
#     print(nops(B));
      LT1:=leadterm(G[B[1][1]],TO); LT2:=leadterm(G[B[1][2]],TO);
      if lcm(LT1,LT2) <> LT1*LT2 then
        S := DivMP(spoly(G[B[1][1]],G[B[1][2]],TO),Mons,Polys,TO);
        if S <> 0 then
          G:=[op(G),S];
          if type(S,`+`) = true then
             Polys:=[op(Polys),S];        # S is a poly; add all pairs
             for i to nops(G) - 1 do
               B:=[op(B),[i,nops(G)]];
             od
          else
             Mons:=[S,op(Mons)];
             for i to nops(G) - 1 do     # S is a mon; add only pairs with mons
                if type(G[i],`+`) = true then
                   B:=[op(B),[i,nops(G)]];
                fi;
             od;
          fi;
        fi;
      fi;
      B:=[op(2..nops(B),B)];
    od;
    RETURN(G)
    end:
      

# 
#       
# The Groebner Walk with "target order" the elimination order
# on the first set of nops(elimvars) variables.  (This current version 
# is a "smarter" implementation, using "lifting" rather than the full 
# Buchberger algorithm to do individual conversion steps.) 
#
# OCTOBER 21, 1998 -- This version also incorporates user-defined term
# orders (using the ``fixed'' version of the Groebner package routines
# in MyTermOrder.map and MatOrder.map
#  
read `/home/fac/little/MathWork/MyTermOrder.map`;
read `/home/fac/little/MathWork/MatOrder.map`;
#
IGenWalkX:=proc(G,Weight,TOStart,SuddenDeathOption,VerboseOutputOption)
  
  global M;
  local TOCurrent,TONew,MatNew,GCurrent,GNew,WeightCurrent,WeightNew,
        WNPrime,g,h,LMList,ltime,Done,i,InGCurrent,InGNew,dq,ltd,
        path,evec,gs,v;
   
  GCurrent:=G;
  TOCurrent:=TOStart;
  WeightCurrent:=Weight;
  evec:=[];
  for i to nops(elimvars) do
    evec:=[op(evec),1];
  od;
  for i to nops(leftvars) do
    evec:=[op(evec),0];
  od;
  Done := false;
  while evalb(Done = false) do
     LMList := [];
     for g in GCurrent do 
        LMList:=[op(LMList),leadmon(g,TOCurrent)[2]]
     od:
     if VerboseOutputOption = true then
       print(LMList);
     fi; 
     path:=evalm((1-t)*WeightCurrent+t*evec);
     ltime:=LeaveCone2(GCurrent,WeightCurrent,TOCurrent);
    # print(ltime);
     ltd:=denom(ltime);
     if ltime = 1 then
        Done := true;
     fi;
     WeightNew := [];
     for i to nops(vars) do 
        WeightNew:=[op(WeightNew),ltd*subs(t=ltime,path[i])] 
     od:
     print(WeightNew);
  #   WNPrime:=WeightNew+evec;          #!!!!!!!!!!!CAUTION!!!!!!!!!!!
                                       #This is a short-cut that
                                       #can cause problems if the 
                                       #the Groebner fan for the ideal
                                       #contains any extremely narrow
                                       #cones.  If so, "scaling up" all
                                       #weights sufficiently should solve 
                                       #the problem.  There is a limit
                                       #to the size of the integers that
                                       #can appear in weight vectors, however,
                                       #so this scaling can also fail!
     MatNew:=[WeightNew,op(RestOfMat)]:
     M:=MatNew;
     TONew:=MyTermOrder(A,MatrixCompare,vars);
#    TONew:=termorder(A,'matrix'(MatNew,vars));
     InGCurrent:=[];
     for g in GCurrent do
        InGCurrent:=[op(InGCurrent),InForm(g,vars,WeightNew)];
     od;
     # print(InGCurrent);
     InGNew:=SPBuchberger(InGCurrent,TONew);
     # print(InGNew);
     GNew:=[];
     for g in InGNew do
        dq:=DivQuots(g,InGCurrent,TOCurrent);
        h:=0;
        for i to nops(InGCurrent) do
           h := h + dq[1][i]*GCurrent[i]; 
        od;
        GNew:=[op(GNew),simplify(h)];
     od;        
     GNew:=inter_reduce(GNew,TONew);
     if VerboseOutputOption = true then
       print(`next Groebner basis finished`);
     fi;
     if SuddenDeathOption = true then
       for g in GNew do
          gs := g;
          for v in elimvars do
             gs := subs(v=0,gs);
          od;
          if gs = g then
            print(`sudden death`);
            print(g);
            Done := true;
          fi;
       od;
     fi;
     if Done = false then
        WeightCurrent := WeightNew;
        GCurrent:=GNew;
        TOCurrent:=TONew;
     fi;
     od;
  RETURN(GNew)
  end:
#
#