reduce-historical/r37/packages/groebner/buchbg.red

module buchbg;  % Central Groebner base code: Buchberger algorithm.

%   Authors:   H. Melenk, H.M. Moeller,  W. Neun
%   ZIB Berlin 
%   July 1988 

fluid '(!*gcd);

fluid '(                           % switches from the user interface
        !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
        !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
        !*groebfullreduction !*groebstat !*groebdivide !*groebprot
        !*groebheufact !*groebweak !*groelterms
 
        vdpvars!*                  % external names of the variables
        !*vdpinteger !*vdpmodular  % indicating type of algorithm
        vdpsortmode!*              % term ordering mode
        secondvalue!* thirdvalue!* % auxiliary: multiple return values
        fourthvalue!*
        groebdomain!*              % domain mode if selected explicitly
        factortime!*               % computing time spent in factoring
        factorlevel!*              % bookkeeping of factor tree
        groefeedback!*             % sideeffect factorization
        groesfactors!*             % data structure for act. fact.
        pairsdone!*                % list of pairs already calculated
        probcount!*                % counting subproblems
        vbccurrentmode!*           % current domain for base coeffs.
        vbcmodule!*                % for modular calculation: 
                                   % current prime
        groetime!*
        !*gsugar                   % enable Traverso's sugar technique
        groecontcount!*            % control of content reduction     
        gmodule!*                  % internal module basis
        groebabort!*               % input abort conditions
       );
 
global '(groebrestriction          % interface:
                                   %  name of restriction function
         groebresmax            % maximum number of internal results
         gvarslast                 % output: variable list 
         groebmonfac               % minimum exponent for reduction of
                                   % monomial factors
         groebprotfile             % reduction protocol
         glterms                   % list for lterms collection
        );
 
flag ('(groebrestriction groebresmax gvarslast groebmonfac
         groebprotfile glterms), 'share);
 
groebrestriction := nil;
groebresmax := 300;
groebmonfac := 1;
groecontcount!* := 10;
!*groebfullreduction := t;
!*gsugar := t;
!*groelterms := t;
 
switch groebopt,groebres,trgroeb,trgroebs,trgroeb1,
       trgroebr,groebfullreduction,groebstat,groebprot;
 
% variables for counting and numbering
fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!* 
        basecount!* hzerocount!*);
 
%      option ' groebopt'    "optimizes" the given input
%                            polynomial set ( variable
%                            ordering )
%      option ' trgroeb'     prints intermediate
%                            results on the output file
%      option ' trgroeb1'    prints internal representation
%                            of critical pair list d
%      option ' trgroebs '   prints S - polynomials
%                            on the output file
%      option   trgroebr     prints (intermediate) results and
%                            computation statistics
% groebstat  the statistics are printed
% groebres   the H- polynomials are optimised using resultant
%            and factorisation method
%
% groebrm    multiplicities of factors in h-polynomials are reduced
%            to simple factors.
%
% groebdivide
%            the algorithm avoids all divisions (only for modular
%            calculation) , if this switch is set off;
%
% groebprot  Write a protocol to the variable "groebprotfile";
 
!*groebfullreduction := t;
%!*groebPreReduce := T;
!*groebdivide := t;
 

% the code for checkpointing is factored out
% This version: NO CHECKPOINT FACILITY

smacro procedure groebcheckpoint1(); list nil;
smacro procedure groebcheckpoint2(); list nil;
smacro procedure groebcheckpoint2a(); list nil;
smacro procedure groebcheckpoint3(); list nil;
smacro procedure groebcheckpoint4(); list nil;
smacro procedure groebcheckpoint5(); list nil;

 
symbolic procedure buch!-vevdivides!? (vev1,vev2);
    % test: vev1 divides vev2 ? for exponent vectors vev1,vev2
    vevmtest!? (vev2,vev1) and 
     (null gmodule!* or gevcompatible1(vev1,vev2,gmodule!*));
 
symbolic procedure gevcompatible1(v1,v2,g);
  % test whether v1 and v2 belong to the same vector column.
   if null g then t 
   else if null v1 then (null v2 or gevcompatible1('(0),v2,g))
   else if null v2 then gevcompatible1(v1,'(0),g) else
    (car g = 0 or car v1 = car v2) and
         gevcompatible1(cdr v1,cdr v2,cdr g);
      
symbolic procedure gcompatible(f,h);
   null gmodule!* or 
     gevcompatible1(vdpevlmon f,vdpevlmon h,gmodule!*);

symbolic procedure groebmakepair(f,h);
     % construct a pair from polynomials f and h
   begin scalar ttt,sf,sh;
      ttt:=tt(f,h);
      return if !*gsugar then
      <<sf:=gsugar(f) #+ vevtdeg vevdif(ttt,vdpevlmon f);
        sh:=gsugar(h) #+ vevtdeg vevdif(ttt,vdpevlmon h);
        list(ttt,f,h,max(sf,sh))>>
       else list(ttt,f,h);
   end;

% the 1-polynomial will be constructed at run time
% because the length of the vev is not known in advance
fluid '(vdpone!*);
 
symbolic procedure vdponepol;
    % construct the polynomial=1
    vdpone!* := vdpfmon(a2vbc 1,vevzero());
 
symbolic procedure groebner2(p,r);
  % setup all global variables for the Buchberger algorithm
  % printing of statistics
 begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*,
              pairsdone!*,factorlevel!*,groesfactors!*,!*gcd;
       factortime!* := 0;
       groetime!*   := time();
       vdponepol(); % we construct dynamically
       hcount!* := 0; mcount!* := 0; fcount!* := 0;
       bcount!* := 0; b4count!* := 0; hzerocount!* := 0;
       basecount!* := 0; !*gcd := t; glterms := list('list);
       groecontcount!* := 10;
       if !*trgroeb then
       << prin2 "Groebner Calculation starting ";
          terprit 2;
          prin2 " groebopt: "; print !*groebopt;
       >>;
       spac := gctime();
 
       p1:= if !*groebfac or null !*gsugar
          then groebbasein (p,!*groebfac,!*groebres,r) 
                 where !*gsugar=nil
            else gtraverso(p,nil,nil,nil);
 
       if !*trgroeb or !*trgroebr or !*groebstat then
       <<
          spac1 := gctime() - spac;
          terpri();
          prin2t "statistics for GROEBNER calculation";
          prin2t "===================================";
          prin2 " total computing time (including gc): ";
          prin2((tim1 := time()) - groetime!*);
          prin2t "          milliseconds  ";
       if factortime!* neq 0 then
        <<prin2 " (time spent in FACTOR (excl. gc):    ";
          prin2 factortime!*;
          prin2t "          milliseconds)";
        >>;
          prin2 " (time spent for garbage collection:  ";
          prin2 spac1;
          prin2t "          milliseconds)";
          terprit 1;
          prin2  "H-polynomials total: "; prin2t hcount!*;
          prin2  "H-polynomials zero : "; prin2t hzerocount!*;
          prin2  "Crit M hits: "; prin2t mcount!*;
          prin2  "Crit F hits: "; prin2t fcount!*;
          prin2  "Crit B hits: "; prin2t bcount!*;
          prin2  "Crit B4 hits: "; prin2t b4count!*;
       >>;
       return p1;
 end;
 
 
smacro procedure testabort h;
     vdpmember (h,abort1) or
      'cancel = (abort2 := groebtestabort(h,abort2));
 
symbolic procedure groebenumerate f;
  % f is a temporary result. Prepare it for medium range storage
  % and ssign a number
  if vdpzero!? f then f else
   << vdpcondense f;
      if not vdpgetprop(f,'number) then
      <<f := vdpputprop(f,'number,(pcount!* := pcount!* #+ 1));
        if !*groebprot then 
          << groebprotsetq(mkid('poly,pcount!*),'candidate);
             vdpputprop(f,'groebprot,mkid('poly,pcount!*));
          >>;
      >>;
     f>>;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Buchberger's Algorithm
%
% INPUT :   G0 = { f1, ... , fr}  set of nonzero polynomials
% OUTPUT:   groebner basis (list of nonzero polynomials)
%
% internal variables:
%
%    problems  list of problems to be processed. problems is non nil,
%              if the inital problem was split by a successful factoring
%    results   collection of results from problems
%    G         basis under construction
%    G1        local pointer to G
%    D         list of critical pairs during algorithm
%    D1,D2     local lists of pairs during update of D
%    f,fj      polynomials
%    p,p1,p2   pairs
%    s,h       polynomials in the central part of the algorithm
%              (the "s-poly" and the "h-poly" selon Buchberger
%    G99       set of polynomials used for reduction. Held as
%              a search tree
%    abort1    list of polynomials in factorization context.
%              a calculation branch can be cancelled if one of
%              these polys is detected
%    abort2    list of sets of polynomials. If a new h polynomial
%              is calculated, it should be removed from the sets.
%              if one set becomes null, the set restriction is
%              fulfilled and the branch can be cancelled
 
symbolic procedure groebbasein (g0,fact,groebres,abort1);
   begin scalar problems, results, abort2,x,
                g,g1,d,d1,d2,p,p1,s,h,g99,hlist,lv,lasth;
         integer gvbc, probcount!*;
         groebabort!* := abort1; 
         lv := length vdpvars!*;
         groebcheckpoint1();
         for each p in g0 do 
            if vdpzero!? p then g0 := delete(p,g0);
          if !*groebprereduce then g0 := groebprereduce g0;
          x := for each fj in g0 collect
          <<groebsavelterm fj;
            gsetsugar(vdpenumerate vdpsimpcont fj,nil)>>;
          if !*groebprot then
           for each f in x do
           << groebprotsetq(mkid('poly,h:= vdpnumber f), vdp2a f);
              vdpputprop(f,'groebprot,mkid('poly,h))>>;
          g0 := x;
                                       % establish the initial problem
        problems := list (list (nil,nil,nil,g0,abort1,nil,nil,
                          vbccurrentmode!*, nil,nil) );
        !*trgroeb and groebmess1(g,d);
        goto macroloop;
        groebcheckpoint2();
macroloop:
  while problems and gvbc < groebresmax do
  begin
     groebcheckpoint2a();
                                       % pick up next problem
     x := car problems;
     groebcheckpoint3();
     d := car x;
     g := cadr x;
     g99 := groebstreereconstruct caddr x;
     g0 := cadddr x; 
     abort1 := nth(x,5);
     abort2 := nth(x,6);
     pairsdone!* := nth(x,7);
     h := nth(x,8);                     % vbcCurrentMode!*
     factorlevel!* := nth(x,9);
     groesfactors!* := nth(x,10);
     groebcheckpoint4();
     problems := cdr problems;
     g0 := % sort G0, but keep factor in first position
       if factorlevel!* and  g0 and cdr g0 then car g0 . vdplsort cdr g0
        else vdplsort g0;
     x := nil;
     lasth := nil;
                !*trgroeb and groebmess23 (g0,abort1,abort2);
      while d or g0 do
       begin
          if groebfasttest(g0,g,d,g99) then goto stop;
                 !*trgroeb and groebmess50(g);
          if g0 then
          <<    h := car g0;
               g0 := cdr g0;
               gsetsugar(h,nil);
               groebsavelterm h;
               p := list(nil,h,h); >>
             else
          <<   p := car d;
               d := delete (p,d);
               s := groebspolynom (cadr p, caddr p);
               if fact then
                    pairsdone!* := (vdpnumber cadr p
                                      . vdpnumber caddr p)
                              . pairsdone!*;
                           !*trgroeb and groebmess3 (p,s);
               h:=groebnormalform(s,g99,'tree);
               groebsavelterm h;
               h:=groebsimpcontnormalform h;
                 if vdpzero!? h then !*trgroeb and groebmess4(p,d);
                                      % test for possible chains
               if not vdpzero!? h then  % only for real h's
                   << s := groebchain (h,cadr p,g99);
                      if s = h then
                        h := groebchain (h,caddr p,g99);
                        if secondvalue!* then
                            g := delete(secondvalue!*,g)>>; >>;
          if vdpzero!? h then goto bott;
          if vevzero!? vdpevlmon h then % base 1 found
                  <<   !*trgroeb and groebmess5(p,h); goto stop>>;
          if testabort(h) then
                  <<   !*trgroeb and groebmess19(h,abort1,abort2);
                   goto stop>>;
          s:= nil;
                    % look for implicit or explicit factorization
          hlist := nil;
          if groebrestriction!* then
                    hlist := groebtestrestriction(h,abort1);
          if not hlist and fact then 
                    hlist := groebfactorize(h,abort1,g,g99);
          if hlist = 'zero then goto bott;
          if groefeedback!* then g0 := append(groefeedback!*,g0);
          groefeedback!* := nil;
                   % factorisation found but only one factor survived
          if hlist and length hlist = 2 then
                   <<h := car cadr hlist; hlist := nil>>;
          if not hlist and groebres then
             <<hlist := groebtestresultant(lasth,h,lv);
               if hlist then groebres := nil>>;
          if hlist then
                <<if hlist neq 'cancel then
                   problems:=
                    groebencapsulate(hlist,d,g0,g,g99,
                           abort1,abort2,problems,fact);
                  go stop>>;
                  % h polynomial is accepted now
               h := groebenumerate h; !*trgroeb and groebmess5(p,h);
                              % construct new critical pairs
               d1 := nil;
                   !*trgroeb and groebmess50(g);
              for each f in g do
                  if(car p or % that means "not an input polynomial"
                     not member (vdpnumber h . vdpnumber f, pairsdone!*)
                    ) and gcompatible(f,h) then
                    <<d1 := groebcplistsortin(groebmakepair(f,h),d1);
                      if tt(f,h) = vdpevlmon(f) then 
                              <<g := delete (f,g);
                                !*trgroeb and groebmess2 f>> >>;
                  !*trgroeb and groebmess51(d1);
               d2 := nil;
               while d1 do
                  <<d1 := groebinvokecritf d1;         
                    p1 := car d1; d1 := cdr d1;
                    d2 := groebinvokecritbuch4 (p1,d2);
                    d1 := groebinvokecritm (p1,d1) >>;
                 d := groebinvokecritb (h,d);  
                 d := groebcplistmerge(d,d2);   
                         % monomials and binomials
                 if vdplength h < 3 and car p then 
                   <<g := groebsecondaryreduction (h,g,g99,d,nil,t);
                     if g = 'abort then goto stop;
                     g99 := secondvalue!*;
                     d := thirdvalue!*>>;
                 g := h . g;
                 lasth := h;
                 g99 := groebstreeadd(h, g99);
                            !*trgroeb and groebmess8(g,d);
                 goto bott;
    stop:        d := g := g0 := nil;
    bott: groebcheckpoint5();
   end;
   g := vdplsort g;       % such that T descending
   x := groebbasein2(g,g99,problems,abort1,abort2,fact);
   g1 := car x; problems := cdr x;
   if g1 then <<results := g1 . results; gvbc := gvbc+1>>;
   !*trgroeb and groebmess13(g1,problems);
  end;
  if gvbc  >= groebresmax then
      lpriw("########","warning: GROEBRESMAX limit reached");
  return groebbasein3 results;
end;
 
symbolic procedure groebfasttest(g0,g,d,g99); 
      <<g := g0 := d := g99 := nil; nil>>;
           % a hook for special techniques

symbolic procedure groebbasein2(g,g99,problems,abort1,abort2,fact);
  % final reduction for a base G: reduce each polynomial with the
  % other members; here again support of factorization
  begin scalar !*groebfullreduction,!*groebheufact;     % saving value
   scalar g1,f,h,hlist,x,!*gsugar; integer cnt;
   !*groebfullreduction := t;
   g1 := nil;
   while g do
      <<h := car g;
        g := cdr g;
        if !*groebprot then
          groebprotsetq('candidate,mkid('poly,vdpnumber h));
        h := groebnormalform (h,g,'sort);
        f := groebsimpcontnormalform h;
                           !*trgroeb and groebmess26(h,f);
        if !*groebprot then 
          groebprotsetq({'gb,cnt:=cnt+1},'candidate);
        if vdpone!? f then <<g1 := g := nil>>;  % base {1} found
                                                % very late now
        if fact and not vdpzero!? f then
            << hlist := groebfactorize (f,abort1,nil,nil);
               if not null hlist then
                 <<          % lift structure
                  hlist := for each x in cdr hlist collect car x;
                             % discard superfluous factors
                  for each h in hlist do
                       if vdpmember(h,abort1) then
                           <<hlist := delete(h,hlist);
                               !*trgroeb and 
                               groebmess19(h,abort1,abort2)>>;
                                       % generate new subproblems
                  x := 0;
                  for each h in hlist do
                    <<hlist := delete(h,hlist);
                      h := groebenumerate h;
                      problems:=
                         list(nil,                 % null D
                              append(g1,g),        % base
                              g99,                 % G99
                              list h,              % G0
                              append(hlist,abort1),
                              abort2,
                              pairsdone!*,
                              vbccurrentmode!*,
                              (x := x+1) . factorlevel!*,
                              groesfactors!*
                             ) . problems;
                    >>;
                  g := g1 := nil;     % cancel actual final reduction
                  f := nil;
                >>  
              >>;
        if f and vdpevlmon h neq vdpevlmon f then
        <<g:= vdplsort append(g,f . g1); g1 := nil>> else
        if f and not vdpzero!? f then g1 := append (g1 ,list f);
        >>;
     return g1.problems;
   end;

symbolic procedure groebbasein3 results;
 % final postprocessing : remove multiple bases from the result
 begin scalar x,g,f,p1,p2;
  x := nil; g := results; p1 := p2 := 0;
  while results do
      <<if vdpone!? car car results      % exclude multiple {1}
              then p1 := p1 + 1              % count ones
          else
      <<f := for each p in car results     % delete props for member 
                     collect vdpremallprops p;
        if  member (f,x)                     % each base only once
              then p2 := p2 + 1              % count multiples
        else if not groeb!-abort!-id(f,groebabort!*) 
           then x := f . x;
         results := cdr results>> >>;
  results := if null x then list list vdpone!* else x;
  return results;
end;

fluid '( !*vbccompress);
 
 
symbolic procedure groebchain(h,f,g99);
    % test if a chain from h-plynomials can be computed from the h
   begin scalar h1,h2,h3,count,found;
         secondvalue!* := nil;
       return h; % erst einmal
         if not buch!-vevdivides!? (vdpevlmon h, vdpevlmon f) 
              then return h;
         h2 := h;
         h1 := f;
         found := t;
         count := 0;
         while found do
                  <<h3 := groebspolynom(h1,h2);
                    h3 := groebnormalform(h3,g99,'tree);
                    h3 := vdpsimpcont h3;
                    if vdpzero!? h3 then
                           <<found := nil;
                            prin2t "chain---------------------------";
                             vdpprint h1;
                             vdpprint h2;
                             vdpprint h3;
                             secondvalue!* := f;
                             count := 9999>>
                      else
                   if vdpone!? h3 then
                           <<found := nil;
                            prin2t "chain---------------------------";
                             vdpprint h1;
                             vdpprint h2;
                             vdpprint h3;
                             h2 := h3;
                             count := 9999>>
 
                      else
                    if buch!-vevdivides!?(vdpevlmon h3, vdpevlmon h2) 
		      then <<found := t;
                             prin2t "chain---------------------------";
                             vdpprint h1;
                             vdpprint h2;
                             vdpprint h3;
                             h1 := h2;
                             h2 := h3;
                             count := count+1>>
                      else
                           found := nil;
                  >>;
       return if count > 0 then
             << prin2 "CHAIN :"; prin2t count; h2>>
            else h;
    end;
 
 
symbolic procedure groebencapsulate(hlist,d,g0,g,g99,
                   abort1,abort2,problems,fact);
% hlist is a factorized h poly. This procedure has the job to
% form new problems from hlist and to add them to problems.
% Result is problems.
% Standard procedure: only creation of subproblems
    begin scalar factl,      % list of factorizations under way
                 u,y,z;
          integer fc;
      if length vdpvars!*>10 or car hlist neq 'factor then
          return groebencapsulatehardcase(hlist,d,g0,g,g99,
                   abort1,abort2,problems,fact);
                            % encapsulate for each factor
      factl := groebrecfactl list hlist;;
                   !*trgroeb and groebmess22 (factl,abort1,abort2);
      for each x in reverse factl do
           <<y := append (car x, g0);
             z := vdpunion (cadr x,abort1);
             u := append(caddr x,abort2);
                  problems :=  list(
                            d,
                            g,
                            g,                       % future G99
                            y,                       % as new problem
                            z,                       % abort1
                            u,                       % abort2
                            pairsdone!*,             % pairsdone!*
                            vbccurrentmode!*,
                            (fc:= fc+1) . factorlevel!*,
                            groesfactors!*
                           ) . problems;
            >>;
        return problems;
end;
 
 
 
 
symbolic procedure groebencapsulatehardcase(hlist,d,g0,g,g99,
                   abort1,abort2,problems,fact);
% hlist is a factorized h poly. This procedure has the job to
% form new problems from hlist and to add them to problems.
% Result is problems.
% First the procedure tries to compute new h-polynomials from the
% remaining pairs which are not affected by the factors in hlist.
% purpose is to find further factorizations and to do calculations
% in common for all factors in order to shorten the separate later
% branches;
    begin scalar factl,      % list of factorizations under way
                 factr,      % variables under factorization
                 u,h,d1,d2,p1,y,z,p,s,f,gc,pd,break,fl1;
          integer fc;
      factl := list hlist;
      factr := vdpspace car cadr hlist;
      for each x in cdr hlist do
         for each p in x do
             factr := vevunion(factr,vdpspace p);
 
% ITER:
                                 % now process additional pairs
      while d or g0 do
        begin
          break := nil;
          if g0 then
          <<               % next poly from input
               s := car g0; g0 := cdr g0; p := list(nil,s,s);
          >>
             else
          <<               % next poly fropm pairs
               p := car d; d := delete (p,d);
               if not vdporthspacep(car p,factr) then
                 s:= nil else
               <<s := groebspolynom (cadr p, caddr p);
                        !*trgroeb and groebmess3 (p,s);>>;
          >>;
 
          if null s or not vdporthspacep(vdpevlmon s,factr) then
              <<                     % throw away s polynomial
                 f := cadr p;
                 if not vdpmember3(f,g0,g,gc)
                      then gc := f . gc;
                 f := caddr p;
                 if car p and not vdpmember3 (f,g0,g,gc)
                      then gc := f . gc;
                  goto bott>>;
 
          h := groebnormalform(s,g99,'tree);
                  if vdpzero!? h and car p then
                               !*trgroeb and groebmess4(p,d);
 
          if not vdporthspacep(vdpevlmon h,factr) then
              <<                     % throw away h polynomial
                 f := cadr p;
                 if not vdpmember3(f,g0,g,gc)
                      then gc := f . gc;
                 f := caddr p;
                 if car p and not vdpmember3 (f,g0,g,gc)
                      then gc := f . gc;
                 goto bott>>;
 
 %%%         if car p then
 %%%               pairsdone!* := (vdpnumber cadr p . vdpnumber caddr p)
 %%%                             . pairsdone!*;
 
          if vdpzero!? h then goto bott;
 
          if vevzero!? vdpevlmon h then     % base 1 found
                   goto stop;
 
          h := groebsimpcontnormalform h;   %  coefficients normalized
 
          if testabort h then
                 <<!*trgroeb and groebmess19(h,abort1,abort2);
                   goto stop>>;
 
          s:= nil;
 
          hlist := nil;
          if groebrestriction!* then
                    hlist := groebtestrestriction(h,abort1);
          if hlist = 'cancel then goto stop;
 
          if not hlist and fact then
                  hlist := groebfactorize(h,abort1,g,g99);
          if groefeedback!* then g0 := append(groefeedback!*,g0);
          groefeedback!* := nil;
          if hlist and length hlist = 2 then
                   <<h := car cadr hlist; hlist := nil>>;
          if hlist then
                << for each x in cdr hlist do
                      for each h in x do
                           factr := vevunion(factr,vdpspace h);
                  factl := hlist . factl; % add to factors
                  goto bott>>;
 
          h := groebenumerate h;             % ready now
 
                        !*trgroeb and groebmess5(p,h);
 
                              % construct new critical pairs
          d1 := nil;
 
 
          for each f in g do
             if tt(f,h) = vdpevlmon(f) and gcompatible(f,h) then
              <<g := delete (f,g);
                d1 := groebcplistsortin( groebmakepair(f,h) , d1);
                             !*trgroeb and   groebmess2 f;
              >>;
                     !*trgroeb and groebmess51(d1);
          d2 := nil;
          while d1 do
                  <<d1 := groebinvokecritf d1;         
                    p1 := car d1; d1 := cdr d1;
                    d2 := groebinvokecritbuch4 (p1,d2); 
                    d1 := groebinvokecritm (p1,d1);  
                  >>;
 
          d := groebinvokecritb (h,d);            
          d := groebcplistmerge(d,d2);
 
          if vdplength h < 3 then
                    <<g := groebsecondaryreduction (h,g,g99,d,gc,t);
                      if g = 'abort then goto stop;
                      g99 := secondvalue!*;
                      d := thirdvalue!*;
                      gc := fourthvalue!*>>;
          g := h . g;
          g99 := groebstreeadd(h, g99);
 
                            !*trgroeb and groebmess8(g,d);
          goto bott;
 
  stop:
         d := g := g0 := gc := factl := nil;
  bott:
         end;     %ITER
 
                            % now collect all relvevant polys
        g0 := vdpunion(g0,vdpunion(g,gc));
 
                            % encapsulate for each factor
        if factl then
           << factl := groebrecfactl factl;
              !*trgroeb and groebmess22 (factl,abort1,abort2);
           >>;
        for each x in reverse factl do
           <<fl1 := (fc := fc+1) . factorlevel!*;
             break:= nil;
             y := append (car x, g0);
             z := vdpunion (cadr x,abort1);
             u := append(caddr x,abort2);
             if vdpmember(vdpone!*,y) then break:=vdpone!*;
 
                        % inspect the unreduced list first 
             if not break then for each p in z do 
                  if vdpmember(p,y) then break := p; 
 
                        % now prepare the reduced list
            if not break then
            <<y := append (car x,groebreducefromfactors(g0,car x));
              pd := secondvalue!*;
              if vdpmember(vdpone!*,y) then break := vdpone!* else
                 for each p in z do if vdpmember(p,y) then break := p;
              pd := subla(pd,pairsdone!*);
            >>;
 
            if not break then
                  problems :=  list(
                            nil,                     % new D
                            nil,                     % new G
                            nil,                     % future G99
                            y,                       % as new problem
                            z,                       % abort1
                            u,                       % abort2
                            nil,                     % pairsdone!*
                            vbccurrentmode!*,
                            fl1,                     % factorlevel!*,
                            groesfactors!*           % factor db
                           ) . problems
                else !*trgroeb and groebmess19a(break,fl1);
            >>;
        return problems;
end;
 
symbolic procedure groebrecfactl (factl);
% factl is a list of factorizations:a list of lists of vdps
% generate calculation pathes from them
   begin scalar rf,res,type;
     if null factl then return list list(nil,nil,nil);
     rf := groebrecfactl (cdr factl);
     factl := car factl;
     type := car factl;   % FACTOR or RESTRICT
     factl := cdr factl;
     while factl do
        <<for each y in rf do
           if vdpdisjoint!?(car factl,cadr y) then
             res := list( vdpunion(car factl,car y),
                     (if type = 'factor then
                         append (for each x in cdr factl collect car x,
                                                                cadr y)
                        else
                      if type = 'resultant then
                         append (for each x in cdr factl collect cadr x,
                                                                cadr y)
                        else cadr y),
                     (if type neq 'factor and type neq 'resultant then
                         append(cdr factl,caddr y)
                        else caddr y)
                    ) . res;
          factl := cdr factl>>;
     return res;
   end;
 
 
symbolic procedure groebtestabort (h,abort2);
   % tests if h is member of one of the sets in abort2.
   % if yes, it is deleted. If one wet becomes null, the message
   % "CANCEL is returned, otherwise the updated abort2.
    begin scalar x,break,res;
              % car test the occurence
       x := abort2;
       while x and not break do
        << if vdpmember(h,car x) then break := t;
           x := cdr x>>;
       if not break then return abort2;  % not relvevant
       break := nil;
       while abort2 and not break do
        <<x := vdpdeletemember(h,car abort2);
          if null x then break := t;
          res := x . res;
          abort2 := cdr abort2;
        >>;
       !*trgroeb and groebmess25(h,res);
       if break then return 'cancel;
       return res;
    end;
 
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Reduction of polynomials
%
 
symbolic procedure groebnormalform(f,g,type);
     % general procedure for reduction of one polynomial from a set
     %  f is a polynomial, G is a Set of polynomials either in
     %     a search tree or in a sorted list
     % type describes the ordering of the set G:
     %    'TREE     G is a search tree
     %    'SORT     G is a sorted list
     %    'LIST     G is a list, but not sorted
     %  f has to be reduced modulo G
  begin scalar f0,f1,c,vev,divisor,break,done,fold,a;
        integer n,s1,s2;
      if !*groebweak and !*vdpinteger 
            and groebweakzerotest(f,g,type) then return f2vdp nil;
      fold := f; f1 := vdpzero(); a:= vbcfi 1;
      gsetsugar(f1,gsugar f);
      while not vdpzero!? f do
          begin
            vev:=vdpevlmon f; c:=vdplbc f;
 
            if not !*groebfullreduction and not vdpzero!? f1 then
                    g := nil;                 % cut off

            if type = 'sort then
                 while g and vevcompless!?(vev,vdpevlmon car g)
                            do g := cdr g;
            if null g then
                  <<f1:=vdpsum (f1,f); f:=vdpzero();
                    break := t; divisor := nil; goto ready>>;
 
            divisor :=
               if type = 'tree then groebsearchinstree(vev,g)
                               else groebsearchinlist (vev,g);
 
            if divisor then done := t;    % true action indicator
 
            if divisor and !*trgroebs then
                       << prin2 "//-";
                          prin2 vdpnumber divisor;
                          >>;
            if divisor then
                      if vdplength divisor = 1 then
                            f := vdpcancelmvev(f,vdpevlmon divisor)
                          else
                      if !*vdpinteger or not !*groebdivide then
                       << f:=groebreduceonestepint(f,f1,c,vev,divisor);
                          f1 := secondvalue!*; n := n+1;
                          if  not vdpzero!? f and 
                              n #> groecontcount!* then
                          <<f0 := f;
                            f:=groebsimpcont2(f,f1);
                            groecontentcontrol(f neq f0);
                            f1 := secondvalue!*; n := 0 >>;
                       >>
                       else
                        f := groebreduceonesteprat(f,nil,c,vev,divisor)
               else
                  <<!*gsugar and <<s1:=gsugar(f);s2:=gsugar(f1)>>;
                    f1 := vdpappendmon (f1,vdplbc f,vdpevlmon f);
                    f := vdpred f;
                    !*gsugar and <<gsetsugar(f,s1); 
                                   gsetsugar(f1,max(s2,s1))>>; 
                    >>;
            ready:
          end;
      if !*groebprot then groebreductionprotocolborder();
      if not done then f1 := fold;
      return f1;
   end;
 
symbolic procedure groecontentcontrol u;
   % u indicates, that a substantial content reduction was done;
   % update content reduction limit from u.
  groecontcount!*:= if not numberp groecontcount!* then 10 else
                    if u then max(0,groecontcount!*-1)
                         else min(10,groecontcount!*+1);

symbolic procedure groebvbcbig!? a;
    % test if a is a "big" coefficient
     (if numberp x then (x > 1000000000000 or x < -1000000000000)
                   else t) 
             where x=vbcnumber a;

symbolic procedure groebsimpcontnormalform h;
  % simpCont version preserving the property SUGAR 
     if vdpzero!? h then h else
     begin scalar sugar,c;
         sugar :=gsugar h; c:=vdplbc h; 
         h := vdpsimpcont h;
         gsetsugar(h,sugar);
         if !*groebprot and not(c=vdplbc h)then
          groebreductionprotocol2 
             reval list('quotient,vbc2a vdplbc h,vbc2a c);
         return h;
     end;

symbolic procedure groebsimpcont2(f,f1);
  % simplify two polynomials with the gcd of their contents;
    begin scalar c,s1,s2;
         s1:=gsugar f; s2:=gsugar f1;
         c := vdpcontent f;
         if vbcone!? vbcabs c then goto ready;
         if not vdpzero!? f1 then 
         << c := vdpcontent1(f1,c);
            if vbcone!? vbcabs c then goto ready;
            f1 :=  vdpdivmon(f1,c,nil)>>;
         f := vdpdivmon(f,c,nil);
                 !*trgroeb and groebmess28 c;
         groebsaveltermbc c;
         gsetsugar(f,s1); gsetsugar(f1,s2);
      ready:
         secondvalue!* := f1; return f;
    end;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  special case reductions
%
 
symbolic procedure groebprereduce g;
% reduce the polynomials in g with themselves.
% the reduction is continued until headterms are stable
% is possible;
  begin scalar res,work,oldvev,f,oldf,!*groebweak,!*groebfullreduction;
        integer count;
     if !*trgroebs then
     <<  g := for each p in g collect vdpenumerate p;
         for each p in g do vdpprint p>>;
     res := nil; % delete zero polynomials from G
     for each f in g do if not vdpzero!? f then res := f . res;
     work := g := res := reversip res;
     while work do
      << g := vdplsort res;   % sort prvevious result
         if !*trgroebs then prin2t "Starting cycle in prereduction.";
         res := nil;
         count := count + 1;
         work := nil;
         while g do
         << oldf := f:= car g; g := cdr g;
            oldvev := vdpevlmon f;
            f := vdpsimpcont groebnormalform (f,g,'sort);
            if (!*trgroebs or !*groebprot) and not vdpequal(f,oldf) then
            <<f := vdpenumerate f;
              if !*groebprot then
                if not vdpzero!? f then
                 groebprotsetq(mkid('poly,vdpnumber f), vdp2a f)
                 else groebprotval 0;
              if !*trgroebs then
                 <<prin2t "reducing"; vdpprint oldf; prin2t "to";
                                                    vdpprint f>>;
            >>;
            if not vdpzero!? f then
            <<if oldvev neq vdpevlmon f then work := t;
              res := f . res>>;
         >>;
      >>;
     return for each f in res collect vdpsimpcont f;
     end;
 
symbolic procedure groebreducefromfactors (g,facts);
% reduce the polynomials in G from those in facts.
   begin scalar new,gnew,f,nold,nnew,numbers;
     if !*trgroebs then <<
                  prin2t "replacing from factors:";
                  for each x in facts do vdpprin2t x
                        >>;
     while g do
          <<f := car g;
            g := cdr g;
            nold := vdpnumber(f);
            new :=  groebnormalform(f,facts,'list);
            if vdpzero!? new then
               << if !*trgroebs then <<prin2 "vanishes ";
                                       prin2 vdpnumber f
                                     >>;
               >>
            else
            if vevzero!? vdpevlmon new then
               << if !*trgroebs then <<prin2 "ONEPOL ";
                                       prin2 vdpnumber f
                                     >>;
                 g := nil;
                 gnew := list vdpone!*;
               >>
            else
            <<
              if new neq f then
               <<new := vdpenumerate vdpsimpcont new;
                 nnew := vdpnumber new;
                 numbers := (nold . nnew) . numbers;
                 if !*trgroebs then <<prin2 "replacing ";
                                      prin2 vdpnumber f;
                                      prin2 " by ";
                                      prin2t vdpnumber new
                                      >>;
               >>;
              gnew := new . gnew;
            >>;
          >>;
    secondvalue!* := numbers;
    return gnew;
    end;
 
 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  support for Reduction by "simple" polynomials
 
symbolic procedure groebnormalform1(f,p);
    % short version; reduce f by p
    % special case: p is a nomomial
     if vdplength p = 1 then vdpcancelmvev(f,vdpevlmon p)
            else groebnormalform(f,list p,nil);
 
symbolic procedure groebprofitsfromvev(p,vev);
    % tests, if at least one monomial from p would be reduced by vev
    if vdpzero!? p then nil
        else
    if buch!-vevdivides!?(vev, vdpevlmon p) then t
        else
            groebprofitsfromvev(vdpred p,vev);
 
 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  special reduction procedures
 
 
 
symbolic procedure groebreduceonestepint(f,f1,c,vev,g1);
      % reduction step for integer case:
      % calculate f= a*f - b*g a,b such that leading term vanishes
      %                        (vev of lvbc g divides vev of lvbc f)
      %
      % and  calculate f1 = a*f1
      % return value=f, secondvalue=f1
 begin scalar vevlcm,a,b,cg,x,rg1;
            % trivial case: g1 single monomial
      if vdpzero!? (rg1:= vdpred g1)
          then return <<f := vdpred f; secondvalue!* := f1; f>>;
      vevlcm := vevdif(vev,vdpevlmon g1);
      cg := vdplbc g1;
            % calculate coefficient factors
      x := if not !*groebdivide then vbcfi 1 else vbcgcd(c,cg);
      a := vbcquot(cg,x);
      b := vbcquot(c,x);
            % multiply relvevant parts from f and f1 by a (vbc)
      if f1 and not vdpzero!? f1 then f1 := vdpvbcprod(f1,a);
      if !*groebprot then groebreductionprotocol(a,vbcneg b,vevlcm,g1);
      f:= vdpilcomb1 (vdpred f, a, vevzero(),
                     rg1,vbcneg b, vevlcm);
            % return with f and f1
      secondvalue!*:= f1;
      thirdvalue!* := a;
      return f;
 end;
 
 
symbolic procedure groebreduceonesteprat(f,dummy,c,vev,g1);
      % reduction step for rational case:
      % calculate f= f - g/vdpLbc(f)
      %
 begin scalar x,rg1,vevlcm;
            % trivial case: g1 single monomial
      dummy := nil;
      if vdpzero!? (rg1 := vdpred g1) then return vdpred f ;
            % calculate coefficient factors
      x := vbcneg vbcquot(c,vdplbc g1);
      vevlcm := vevdif(vev,vdpevlmon g1);
      if !*groebprot then 
         groebreductionprotocol(a2vbc 1,x,vevlcm,g1);
      return vdpilcomb1(vdpred f,a2vbc 1,vevzero(),
                            rg1,x,vevlcm);
 end;
 
 
symbolic procedure groebreductionprotocol(a,b,vevlcm,g1);
    if !*groebprot then
        groebprotfile :=
        if not vbcone!? a then
         append(groebprotfile,
                list(
                   list('equal,
                        'candidate,
                         list('times,'candidate,vbc2a a)),
                   list('equal,
                        'candidate,
                        list('plus,
                             'candidate,
                             list('times,
                                  vdp2a vdpfmon(b,vevlcm),
                                  mkid('poly,vdpnumber g1) )))
               )    )
        else
        append(groebprotfile,
                list(
                   list('equal,
                        'candidate,
                        list('plus,
                             'candidate,
                             list('times,
                                  vdp2a vdpfmon(b,vevlcm),
                                  mkid('poly,vdpnumber g1) )))
 
               )    ) ;
 
symbolic procedure groebreductionprotocol2 a;
      if !*groebprot then
        groebprotfile :=
	if not(a=1) then
         append(groebprotfile,
                list(
                   list('equal,
                        'candidate,
                         list('times,'candidate,a))));

symbolic procedure groebreductionprotocolborder();
        append(groebprotfile,'!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+ . nil);
 
 
 
symbolic procedure groebprotsetq(a,b);
    groebprotfile :=
      append (groebprotfile,
              list (list ('equal,a,b)));
 
symbolic procedure groebprotval a;
    groebprotfile :=
      append (groebprotfile,
              list (list ('equal,'intermediateresult,a)));
 
 

symbolic procedure subset!?(s1,s2);
    % tests if s1 is a subset of s2
       if null s1 then t
          else
       if member(car s1,s2) then subset!?(cdr s1,s2)
          else
            nil;

symbolic procedure vevsplit (vev);
      % split vev such that each exponent vector has only one 1
      begin scalar n,vp,e;
        n := 0;
        for each x in vev do
            <<n := n+1;
              if x neq 0 then
                <<e := append (vdpevlmon vdpone!*,nil);
                  rplaca(pnth(e,n),1);
                  vp := e . vp;
                >>;
            >>;
        return vp;
     end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   calculation of an S-polynomial
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%general strategy:
%
%  groebSpolynom4 calculates the traditional s-Polynomial from p1,p2
%    (linear combination such that the highest term vanishes).
%  groebSpolynom2 subtracts multiples of p2 from the s-polynomial such
%    that head terms are eliminated early.
 
symbolic procedure groebspolynom (p1,p2);
     groebspolynom2(p1,p2);
 
 
symbolic procedure groebspolynom2 (p1,p2);
     if vdpzero!? p1 then p2 else if vdpzero!? p2 then p1 else
   begin scalar s,tp1,tp2,ts,cand;
     s := groebspolynom3(p1,p2);
     if vdpzero!? s or vdpone!? s or !*groebprot then return s;
     tp1 := vdpevlmon p1; tp2 := vdpevlmon p2;
     while not vdpzero!? s
       and ((buch!-vevdivides!?(tp2,(ts := vdpevlmon s))
		and (cand := p2))
	   or
	   (buch!-vevdivides!?(tp1,(ts := vdpevlmon s))
		and (cand := p1)))
	do << if !*vdpinteger then
                    s := % vdpsimpcont 
                          groebreduceonestepint(s,nil,vdplbc s,ts,cand)
                       else
                          % rational, float and modular case
                        s := groebreduceonesteprat
                                  (s,nil,vdplbc s,ts,cand);
           >>;
     return s;
   end;
 
symbolic procedure groebspolynom3(p,q);
   begin scalar r;
     r:=groebspolynom4(p,q);
     groebsavelterm r;
     return r;
   end;

symbolic procedure groebspolynom4 (p1,p2);
  begin scalar ep1,ep2,ep,rp1,rp2,db1,db2,x,r;
    ep1 := vdpevlmon p1; ep2 := vdpevlmon p2;
    ep := vevlcm(ep1, ep2);
    rp1 := vdpred p1; rp2 := vdpred p2;
    gsetsugar(rp1,gsugar p1); gsetsugar(rp2,gsugar p2);
    r:=     ( if vdpzero!? rp1 and vdpzero!? rp2 then rp1
              else ( if vdpzero!? rp1 then
                      <<db2:=a2vbc 0;
                        vdpprod(rp2,
                                vdpfmon(db1:=a2vbc 1,
                                       vevdif(ep, ep2) ) )
                      >>
                     else if vdpzero!? rp2 then
                      <<db1:=a2vbc 0;
                        vdpprod(rp1,
                                vdpfmon(db2:=a2vbc 1,
                                        vevdif(ep, ep1) ) )
                      >>
                          else
                             <<db1 := vdplbc p1;
                               db2 := vdplbc p2;
                               if !*vdpinteger then
                                  << x:= vbcgcd (db1,db2);
                                     if not vbcone!? x then
                                     << db1 := vbcquot (db1,x);
                                        db2 := vbcquot (db2,x);
                                  >> >>;
                                vdpilcomb1 (rp2,db1,vevdif(ep,ep2),
                                          rp1,vbcneg db2,vevdif(ep,ep1))
                             >>
                   )
            );
     if !*groebprot then
        groebprotsetq('candidate, 
         {'difference,
           {'times,vdp2a vdpfmon(db2,vevdif(ep,ep2)), 
                   mkid('poly,vdpnumber p1)},
           {'times,vdp2a vdpfmon(db1,vevdif(ep,ep1)),
                      mkid('poly,vdpnumber p2)}} );
        
     return r;
end;
 
symbolic procedure groebsavelterm r;
     if !*groelterms and not vdpzero!? r then
        groebsaveltermbc vdplbc r;

symbolic procedure groebsaveltermbc r;
     <<r:=vbc2a r;
       if not numberp r and not constant_exprp r then
         for each p in cdr fctrf numr simp r do
         <<p:=prepf car p;
           if not member(p,glterms) then nconc(glterms,list p);
      >> >>;

symbolic procedure sfcont f;
   % Calculate the integer content of standard form f.
      if domainp f then f else
      gcdf(sfcont lc f, sfcont red f);

symbolic procedure vdplmon u; vdpfmon (vdplbc u,vdplbc u);
 
symbolic procedure vdpmember3(p,g1,g2,g3);
     % test membership of p in one of then lists g1,g2,g3
       vdpmember(p,g1) or vdpmember(p,g2) or vdpmember(p,g3);

symbolic procedure groeb!-abort!-id(base,abort1);
   % Test whether one of the elements in abort1 is
   % member of the ideal described by base. Definite
   % test here.
  if null abort1 then nil else 
   vdpzero!?(groebnormalform(car abort1,base,'list))
     or groeb!-abort!-id(base,cdr abort1);

endmodule;

end;