File r36/src/groebner.red artifact 20644b6448 part of check-in 255e9d69e6


module groebner;

% Author: Herbert Melenk
% in cooperation with Winfried Neun, Michael Moeller.

% module structure:
%
%     GROEBNER  package + GROEBNR2 package
%        polynomial arithmetic:
%         VDP2DIP package included in DIPOLY package
%        factorization:
%         GROEBFAC package is also included here, but could be loaded
%         on demand.

% Modifications:
%
% 17-Sep-1994 (HM):  The ideal variables are now declared in the TORDER
% statement.  The calling conventions can be still used, but are removed
% from the documents.

fluid '(currentvdpmodule!*);

currentvdpmodule!* := 'vdp2dip;

imports vdp2dip;

create!-package('(groebner grinterf grinter2 buchbg groebcri
         groesolv groebopt groebsea groebsor groebspa groebfac
         groebidq kredelw traverso),
         '(contrib groebner));

put('groebner,'version,2.9);

% Other packages needed.

load!-package 'dipoly;

if(null v or v<3.7) where v=get('dipoly,'version)
  then rederr "wrong dipoly module";

smacro procedure tt(s1,s2);
  % lcm of leading terms of s1 and s2
       vevlcm(vdpevlmon s1,vdpevlmon s2);

endmodule;


module grinterf;  % Interface of Groebner package to REDUCE.
                  % Entry points to the main module and general
                  % interface support.


fluid '(!*factor !*complex !*exp asymplis!* current!-modulus);

global '(!*match powlis!* largest!-small!-modulus);

fluid '(                           % switches from the user interface
        !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
        !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
        !*fullreduction !*groebstat !*groebprot !*gltbasis
        !*groebsubs

        !*grmod!*                  % indicating modular coefficients
        !*vdpinteger !*vdpmodular  % indicating type of algorithm
        vdpsortmode!*              % term ordering mode
        secondvalue!* thirdvalue!* % auxiliary: multiple return values
        fourthvalue!*
        factortime!*               % computing time spent in factoring
        factorlvevel!*             % bookkeeping of factor tree
        pairsdone!*                % list of pairs already calculated
        probcount!*                % counting subproblems
        vbccurrentmode!*           % current domain for base coeffs.
        vbcmodule!*                % for modular calculation:
                                   %     current prime
        bczerodivl!*               % coefficient zero divisors (list of
                                   % standard forms)
        gmodule                    % external module basis
        gmodule!*                  % internal module basis
        global!-dipvars!*          % predefined variable list
        !*gsugar                   % enable sugar strategy
       );

global '(groebrestriction          % interface: name of function
         groebresmax               % maximum number of internal results
         gvarslast                 % output: variable list
         groebprotfile
         gltb
         glterms                   % list for lterms collection
        );

flag ('(groebrestriction groebresmax gvarslast groebprotfile
        gltb glterms gmodule),'share);

switch groebopt,groebres,trgroeb,trgroebs,trgroeb1,
       trgroebr,groebstat,gltbasis,gsugar;

% variables for counting and numbering
fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
        basecount!* hzerocount!*);

% control of the polynomial arithmetic actually loaded
fluid '(currentvdpmodule!*);

vdpsortmode!* := 'lex;     % initial mode
gltb := '(list);           % initially empty

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    interface functions

symbolic procedure groebnereval u;
  % non factorizing Groebner calculation
   begin integer n; scalar !*groebfac,!*groebrm,!*factor,
                !*exp; !*exp := t;
      n := length u;
      if n=1 then return cadr groebner1(reval car u,nil,nil)
       else if n neq 2
        then rerror(groebner,1,
                    "GROEBNER called with wrong number of arguments");
       u:= groebner1(reval car u,reval cadr u,nil);
       if !*gltbasis then gltb := cadr gltb;
       return cadr u;
      end;

put('groebner,'psopfn,'groebnereval);

symbolic procedure groebnerfeval u;
  % non factorizing Groebner calculation
   begin integer n; scalar !*groebfac,!*groebrm,!*factor,
         !*exp,!*ezgcd,s,r,q; !*exp := t;
      if null getd 'groebfactorize then load!-package 'groebfac;
      !*groebrm := !*groebfac := t;
      groebrestriction!* := reval groebrestriction;
      if null dmode!* then !*ezgcd:=t;
      n := length u;
      r:=  if n=1 then groebner1(reval car u,nil,nil) else
           if n=2 then groebner1(reval car u,reval cadr u,nil) else
           if n neq 3
        then rerror(groebner,2,
                    "GROEBNER called with wrong number of arguments")
             else  groebner1(reval car u,reval cadr u,reval caddr u);
      q := r;
       % remove duplicates.
      while q do <<s := car q; q := cdr q;
                    if member(s,q) then r := delete(s,r)>>;
      return r;
    end;

put('groebnerf,'psopfn,'groebnerfeval);

symbolic procedure idquotienteval u;
   begin integer n; scalar !*factor,!*exp; !*exp := t;
      n := length u;
      if n=2 then return groebidq(reval car u,reval cadr u,nil)
       else if n neq 3
        then rerror(groebner,3,
                    "IDQUOTIENT called with wrong number of arguments")
       else return groebidq(reval car u,reval cadr u,reval caddr u)
      end;

put('idealquotient,'psopfn,'idquotienteval);

smacro procedure vdpnumber f;
     vdpgetprop(f,'number) ;


symbolic procedure groebner1(u,v,r);
%    Buchberger algorithm system driver. u is a list of expressions
%    and v a list of variables or NIL in which case the variables in u
%    are used.
   begin scalar vars,w,np,oldorder,!*grmod!*;
         integer pcount!*;
      w := for each j in groerevlist u
              collect if eqexpr j then !*eqn2a j else j;
      if null w then rerror(groebner,4,"Empty list in Groebner");
      vars := groebnervars(w,v);
      if r then r := groerevlist r;
      groedomainmode();
      if vars then goto notempty;
      u:=0; for each p in w do if p neq 0 then u:=1;
      return {'list,{'list,u}};
   notempty:
      if dmode!* eq '!:mod!: and null setdiff(gvarlis w,vars)
         and current!-modulus < largest!-small!-modulus
        then !*grmod!* := t;
      oldorder := vdpinit vars;
                  % cancel common denominators
      w := for each j in w collect reorder numr simp j;
                  % optimize variable sequence if desired.
      if !*groebopt and vdpsortmode!* memq '(lex gradlex revgradlex)
        then  << w:=vdpvordopt (w,vars); vars := cdr w;
                 w := car w; vdpinit vars>>;
      w := for each j in w collect f2vdp j;
      if not !*vdpinteger then
      <<np := t;
        for each p in w do
          np := if np then vdpcoeffcientsfromdomain!? p else nil;
        if not np then <<!*vdpmodular:= nil; !*vdpinteger := t>>;
      >>;
      if !*groebprot then
          <<groebprotfile := list 'list>>;
      if r then r := for each p in r collect
           vdpsimpcont f2vdp numr simp p;
      w := groebner2(w,r);

      if cdr w then  % Remove redundant partial bases.
      begin scalar !*gsugar;
        for each b in w do
         for each c in w do
          if b and b neq c then
           <<v:=t; for each p in c do
               v:=v and vdpzero!? groebnormalform(p,b,'list);
             if v then <<w:=delete(b,w); b:=nil>>;
           >>;
      end;
      if !*gltbasis then
        gltb :=
           'list . for each base in w collect
                   'list . for each j in base collect
                      vdp2a vdpfmon(a2vbc 1,vdpevlmon j);
      w := 'list . for each base in w collect
                   'list . for each j in base collect vdp2a j;
      vdpcleanup();
      gvarslast := 'list . vars;
      return w;
   end;

symbolic procedure groebnervars(w,v);
    begin scalar z,dv,gdv,vars;
      if v='(list) then v:=nil;
      v:=v or (gdv:=cdr global!-dipvars!*) and global!-dipvars!*;
      vars :=
       if null v then
                for each j in gvarlis w collect !*a2k j
          else                      % test, if vars are really used
             << z := gvarlis w;
                groebnerzerobc setdiff(z,v:= groerevlist v);
                for each j in v do
                    if member(j,z) then dv := !*a2k j . dv;
                dv := reversip dv;
                if not (length v = length dv) and !*trgroeb then
                      << prin2 " Groebner: ";
                         prin2 (length v - length dv);
                         prin2t " of the variables not used";
                         terpri () >>;
                dv>>;
       return gdv or vars;
    end;

symbolic procedure groebnerzerobc u;
% u is the list of parameters in a Groebner job. Extract the
% corresponding rules from !*match and powlis!*.
 if u then
  begin scalar w,m,p;
    bczerodivl!* := nil;
    m:=!*match; !*match:=nil;
    p:=powlis!*; powlis!*:=nil;
    for each r in m do if cadr r='(nil . t) then
    <<w:=(numr simp {'difference,'times.for each q in car r collect
{'expt,car q,cdr q}
                    ,caddr r});
      for each x in kernels w do if not member(x,u) then w:=nil;
      if w then bczerodivl!* := w . bczerodivl!*;
    >>;
    for each r in p do if member(car r,u) and caddr r='(nil . t) then
    <<w:=(numr simp {'difference, {'expt,car r,cadr r} ,cadddr r});
      bczerodivl!* := w . bczerodivl!*;
    >>;
    for each r in asymplis!* do if member(car r,u) then
      bczerodivl!* := (r .*1 .+nil) . bczerodivl!*;
    !*match:=m; powlis!*:=p;
  end;

% symbolic procedure maklist pl;
% make list of polynomials. pl is a list of polynomials.
% maklist pl returns a list of distributive polynomials.
%  for each p in pl collect f2vdp car p;

symbolic procedure gvarlis u;
   % Finds variables (kernels) in the list of expressions u.
   sort(gvarlis1(u,nil),function ordop);

symbolic procedure gvarlis1(u,v);
   if null u then v
    else union(gvar1(car u,v),gvarlis1(cdr u,v));

symbolic procedure gvar1(u,v);
   if null u or numberp u or (u eq 'i and !*complex) then v
    else if atom u then if u member v then v else u . v
    else if get(car u,'dname) then v
    else if car u memq '(plus times expt difference minus)
     then gvarlis1(cdr u,v)
    else if car u eq 'quotient then gvar1(cadr u,v)
    else if u member v then v
    else u . v;

symbolic procedure groebidq(u,f,v);
%  Ideal quotient. u is a list of expressions (Gbasis), f a polynomial
%  and v a list of variables or NIL
   begin scalar vars,w,np,oldorder,!*factor,!*exp;
         integer pcount!*; !*exp := t;
      w := for each j in groerevlist u
              collect if eqexpr j then !*eqn2a j else j;
      if null w then rerror(groebner,5,"Empty list in IDEALQUOTIENT");
      if eqexpr f then f := !*eqn2a f;
      vars := groebnervars(w,v);
      groedomainmode();
      if null vars then vdperr 'idealquotient;
      oldorder := vdpinit vars;
                  % cancel common denominators
      w := for each j in w collect numr simp j;
      f :=  numr simp f;
      w := for each j in w collect f2vdp j;
      f := f2vdp f; % now do the conversions
      if not !*vdpinteger then
      <<np := t;
        for each p in f.w do
          np := if np then vdpcoeffcientsfromdomain!? p
                else nil;
        if not np then <<!*vdpmodular:= nil; !*vdpinteger := t>>;
      >>;
      w := groebidq2 (w,f);
      w := 'list . for each j in w collect vdp2a j;
      setkorder oldorder;
      return w;
   end;

fluid '(!*backtrace);

symbolic procedure vdperr name;
  % case that no variables were found
    <<prin2 "**** Groebner illegal parmeter in "; prin2 name;
      if !*backtrace then backtrace();
      rerror(groebner,6," ,e.g. no relevant variables found")>>;

symbolic procedure groeparams(u,nmin,nmax);
   % u is a list of psopfn-parameters; they are given to REVAL and
   % the number of parameters is controlled to be between nmin, nmax
   % result is the list of evaluated parameters padded with NILs
     begin integer n; scalar w; n:= length u;
       if n<nmin or n>nmax then rerror(groebner,7,
         "Illegal number of parameters in call to Groebner package");
       u:= for each v in u collect
         <<w := reval v;
           if eqcar(w,'list) then 'list . groerevlist w else w>>;
       while length u < nmax do u := append(u,'(nil));
       return u;
     end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialization of the distributive plynomial arithmetic
%

symbolic procedure vdpinit (vars);
     begin scalar r,gm;
          % Eventually set up module basis.
        if eqcar(gmodule,'list) and cdr gmodule then
         gm := for each y in cdr gmodule collect
          <<y := reval y;
            if not member(y,vars) then vars:=append(vars,{y});
            y >>;
        r:=vdpinit2(vars);
          % convert an eventual module basis.
        gmodule!* := if gm then vdpevlmon a2vdp ('times . gm);
        return r;
    end;

symbolic procedure groedomainmode();
     <<!*vdpinteger := !*vdpmodular := nil;
      if not flagp(dmode!*,'field) then !*vdpinteger := t
           else
      if !*modular then !*vdpmodular := t>>;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% some lisp functions which are not member of standard lisp
%

symbolic procedure groedeletip(a,b);
   begin scalar q;
       while b and  a= car b do b:= cdr b;
       if null b then return nil;
       q := b;
       while cdr b do if a=cadr b then cdr b := cddr b else b:= cdr b;
       return q;
   end;

symbolic procedure groerevlist u;
      <<if idp u then u := reval u;
       for each p in getrlist u collect reval p>>;

endmodule;


module grinter2;  % Interface of Groebner package to REDUCE:
                  % autoloadint entry points to operators of groebnr2;


symbolic procedure groebnr2entry(fn,u);
       <<groebnr2load();
         apply(get(fn,'psopfn),list u)>>;

symbolic procedure groebnr2load();
       load!-package 'groebnr2;

put('gsort,'psopfn,
      '(lambda(u)(groebnr2entry 'gsort u)));

put('gsplit,'psopfn,
     '(lambda(u)(groebnr2entry 'gsplit u)));


put('gspoly,'psopfn,
     '(lambda(u)(groebnr2entry 'gspoly u)));

put('gvars,'psopfn,
     '(lambda(u)(groebnr2entry 'gvars u)));

put('greduce,'psopfn,
     '(lambda(u)(groebnr2entry 'greduce u)));

put('preduce,'psopfn,
     '(lambda(u)(groebnr2entry 'preduce u)));

put('groebnert,'psopfn,
     '(lambda(u)(groebnr2entry 'groebnert u)));

put('preducet,'psopfn,
     '(lambda(u)(groebnr2entry 'preducet u)));

put('groebnerm,'psopfn,
     '(lambda(u)(groebnr2entry 'groebnerm u)));

put('glexconvert,'psopfn,
     '(lambda(u)(groebnr2entry 'glexconvert u)));

put('hilbertpolynomial,'psopfn,
     '(lambda(u)(groebnr2entry 'hilbertpolynomial u)));

put('gzerodim!?,'psopfn,
     '(lambda(u)(groebnr2entry 'gzerodim!? u)));

put('dd_groebner,'psopfn,
     '(lambda(u)(groebnr2entry 'dd_groebner u)));

put('trgroeb,'simpfg,'((t (groebnr2load))));
put('trgroebs,'simpfg,'((t (groebnr2load)(setq !*trgroeb t))
                        (nil (setq !*trgroeb nil))));
put('trgroebr,'simpfg,'((t (groebnr2load))));
put('groebstat,'simpfg,'((t ( groebnr2load))));
put('groebweak,'simpfg,'((t (groebnr2load))));
put('groebres,'simpfg,'((t (groebnr2load))));

symbolic procedure groebtestrestriction (a1,a2);
   <<groebnr2load(); apply('groebtestrestriction,list(a1,a2))>>;
endmodule;


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;
         integer lastcheckpoint,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;


% RESTART:
        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);

% ITERATION :
      while d or g0 do
       begin
          if groebfasttest(g0,g,d,g99) then goto stop;
                 !*trgroeb and groebmess50(g);
          if g0 then
          <<          % take next poly from input
               h := car g0;
               g0 := cdr g0;
               gsetsugar(h,nil);
               groebsavelterm h;
               p := list(nil,h,h);
          >>
             else
          <<          % take next poly from pairs
               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;      %  ITERATION

   g := vdplsort g;       % such that T descending

% FINALREDUCTION:
   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;     % MACROLOOP

  if gvbc  >= groebresmax then
      lpriw("########","warning: GROEBRESMAX limit reached");

  return groebbasein3 results;
end;

symbolic procedure groebfasttest(g0,g,d,g99); 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
      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;


module groebcri;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  criteria for the Buchberger algorithm
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fluid '(bcount!* b4count!* mcount!* fcount!*);


smacro procedure tt(s1,s2);
  % lcm of leading terms of s1 and s2
       vevlcm(vdpevlmon s1,vdpevlmon s2);

smacro procedure atleast2elementsin (u);
    % test if u has at least a cadr element
    u and cdr u;

symbolic procedure groebbuchcrit4(p1,p2,e);
% Buchberger criterion 4. p1 and p2 are distributive
% polynomials. e is the least common multiple of
% the leading exponent vectors of the distributive
% polynomials p1 and p2. groebBuchcrit4(p1,p2,e) returns a
% boolean expression. True if the reduction of the
% distributive polynomials p1 and p2 is necessary
% else false.
% orig:
%     e neq vevsum( vdpevlmon p1, vdpevlmon p2);
     groebbuchcrit4t(vdpevlmon p1,vdpevlmon p2);

symbolic procedure groebbuchcrit4t(e1,e2);
% nonconstructive test of lcm(e1,e2) = e1 + e2
% equivalent: no matches of nonzero elements.
  if null e1 or null e2 then nil
     else
  if (car e1 neq 0) and (car e2 neq 0) then t
     else groebbuchcrit4t(cdr e1,cdr e2);

symbolic procedure groebinvokecritbuch4 (p,d2);
% Buchberger's criterion 4 is tested on the pair p and the list
% D2 of critical pairs is updated with respect to that crit.
% Result is the updated D2;
    begin scalar p1,p2,vev1,vev2,f1,f2,fd,b4;
          p1 := cadr p;
          p2 := caddr p;
          vev1 := vdpevlmon p1;
          vev2 := vdpevlmon p2;
          f1 := vdpgetprop(p1,'monfac);
          f2 := vdpgetprop(p2,'monfac);
                   % discard known common factors first
          if f1 and f2 then
             <<fd := vevmin (f1,f2);
               b4 := groebbuchcrit4t(vevdif(vev1,fd),vevdif(vev2,fd));
               if b4 and    % is the body itself a common factor?
                     vevdif(vev1,f1) = vevdif(vev2,f2)
                         % test if the polys reduced by their monom.
                         % factor are equal
                     and groebbuchcrit4compatible(p1,f1,p2,f2)
                       then b4 := nil;
             >>
           else
               b4 := groebbuchcrit4t(vev1,vev2);

          if b4 then
                d2 := append (d2, list p)
          else
                b4count!* := b4count!* + 1;
          return d2;
     end;

symbolic procedure groebbuchcrit4compatible (p1,f1,p2,f2);
% p1,p2 polys, f1,f2 exponent vectors (monomials), which are known to
% be factors of their f;
% tests, if p1/f1 = p2/f2
   if vdpzero!? p1 then vdpzero!? p2
       else
   if vdplbc p1 = vdplbc p2
              and
      groebbuchcrit4compatiblevev(vdpevlmon p1,f1,vdpevlmon p2,f2)
       then groebbuchcrit4compatible (vdpred p1,f1,vdpred p2,f2)
      else
          nil;


symbolic procedure groebbuchcrit4compatiblevev (vev1,f1,vev2,f2);
     if null vev1 then null vev2
          else
     if (if f1 then car vev1 - car f1 else car vev1)
                         =
        (if f2 then car vev2 - car f2 else car vev2)
          then
            groebbuchcrit4compatiblevev (cdr vev1,
                                   if f1 then cdr f1 else nil,
                                   cdr vev2,
                                   if f2 then cdr f2 else nil)
         else
            nil;
symbolic procedure groebinvokecritf d1;
% groebInvokeCritF tests a list D1 of critical pairs. It cancels all
% critical pairs but one in D1 having the same lcm (i.e. car
% component) as car(D1). This only one is chosen, if possible,
% such that it doesn't satisfy groebBuchcrit4.
% Version: moeller upgraded 5.7.87

begin
  scalar tp1,p2,active;
  tp1 := car(car(d1));
  active := atleast2elementsin d1;
  while active do
    << p2 := cadr d1;
       if car(p2) = tp1 then
          << fcount!* := fcount!* +1;
             if not groebbuchcrit4t(cadr p2, caddr p2)
                 then d1 := cdr(d1)
                 else d1 := groedeletip(p2,d1);
             active := atleast2elementsin d1;
          >>
         else active := nil
    >>;
  return d1;
end;

symbolic procedure groebinvokecritm (p1,d1);
% D1 is a list of critical pairs, p1 is a critical pair.
% crit M tests, if the lcm of p1 divides one of the lcm's in D1.
% If so, this object is eliminated.
% Result is the updated D1;
   << for each p3 in d1 do
         if buch!-vevdivides!?(car(p1), car(p3)) then
             <<mcount!* := mcount!*+1;
               d1 := groedeletip (p3,d1)>>;       %  Criterion M
      d1>>;



symbolic procedure groebinvokecritb (fj,d);
% D is a list of critical pairs, fj is a polynomial.
% Crit B allows to eliminate a pair from D, if the leading monomial
% of fj divides the lcm of the pair, but the lcm of fj with each of
% the members of the pair is not the lcm of the pair itself
% Result is the updated D;
  << for each p in d do
         if buch!-vevdivides!?(vdpevlmon(fj), car(p)) and
            tt(fj,cadr(p)) neq car(p) and % Criterion B
            tt(fj,caddr(p)) neq car(p) then
                     <<bcount!* := bcount!* +1;
                       d:= delete  (p,d)>>;
    d>>;

endmodule;


module groesolv; % Tools for solving systems of polynomials (and poly-
%               % nomial equations) based on Groebner basis techniques.

%   Authors:     H.Melenk (ZIB Berlin)
%                H.M Moeller (Fernunversitaet Hagen)
%                W.Neun (ZIB Berlin)
%
%    Aug 1992:   accepting external solutions for univariate pols.
%
%    March 1994: access to roots/multroot.
%
%    Feb 1995: assumptions, requirements added.

% operators:
%
%      GROESOLVE      does the whole job of solving a nonlinear
%                     set of expressions and / or equations
%
%      GROEPOSTPROC   expects that its car parameter is a
%                     lexical groebner base already


fluid '(groesolvelvevel!* !*groebprereduce);
fluid '(groesoldb!* groesoldmode!* !*groebnumval !*groebcomplex
        !*groebopt !*groebprot !*groesolrecurs
        !*groesolgarbage denominators!* variables!*
        groebroots!*    % a-list for predefined root expressions
        depl!*          % the reduce dependency list
        !*vdpinteger    % coefficient mode
        !*compxroots    % switch for multroot
        gmodule
        !*arbvars
       );

fluid '(!*convert !*ezgcd !*msg !*precise dmode!* !*varopt);
global '(groebprotfile gvarslast !*trgroesolv glterms assumptions
         requirements);

groesolvelvevel!* := 0;

symbolic procedure groesolveeval u;
    begin scalar !*ezgcd,gblist,oldtorder,res,!*groebopt,!*precise,
            y,fail,problems,denominators!*,variables!*,gmodule,at;
      if null dmode!* then !*ezgcd := t;
      !*groebopt := !*varopt; gvarslast :='(list);
      oldtorder := apply1('torder,'(lex));
      groesoldmode!* := get(dmode!*,'dname);
      !*groebcomplex := !*complex;
      groesetdmode(groesoldmode!*,nil);
      problems:={u};
      while problems and not fail do
      <<u:=car problems; problems:=cdr problems;
        gblist := cdr groebnerfeval u;
        at := union(at,cdr glterms);
        !*groebopt := nil;    % 29.8.88
        groesetdmode(groesoldmode!*,t);
        if null variables!* then variables!*:=gvarslast;
        if not(gblist = '((list 1))) then
         for each gb in gblist do
         <<if !*trgroesolv then
             <<writepri("starting from basis",'first);
               writepri(mkquote gb,'last)>>;
           for each r in res do if gb
                % do not compare with the mother problem;
             and not subsetp(car r,car u)
              then
           if groesolvidsubset!?(gb,car r,variables!*) then
              res:=delete(r,res) else
           if groesolvidsubset!?(car r,gb,variables!*) then
              <<gb:=nil;
                if !*trgroesolv then writepri("redundant",'only)
              >>;
           if gb then
           <<y:=groesolvearb(groesolve0(gb,variables!*),variables!*);
             if y neq 'failed then res:=(gb.y).res else fail:=t;
             if !*trgroesolv then
             <<writepri("partial result:",'first);
               writepri(mkquote('list.cdar res),'last)
             >>;
             for each d in denominators!* do
              problems:={append(gb,{d}),variables!*}.problems;
             denominators!*:=nil;
           >>;
          >>;
        >>;
      apply1('torder,{oldtorder});
      problems:=nil; if fail then res:=nil;
      if null res then requirements:=append(requirements,at)
        else assumptions := append(assumptions,at);
      for each r in res do problems:=union(cdr r,problems);
      return 'list . groesolve!-redun2 problems
     end;

symbolic procedure groesolve!-redun2 sol;
  % Sol is a list of solutions; remove redundancies, now not
  % by ideal theory but by simple substitution.
 begin scalar b;
   for each s in sol do
     if s memq sol then
     <<b:=nil;
       for each r in sol do
         if not(r eq s) then
           if not b and groesolve!-redun2a(r,s) then b:=r;
       if b then sol:=delete(s,sol);
     >>;
   return sol;
 end;

symbolic procedure groesolve!-redun2a(r,s);
  % redundancy test: if sub(s,r)=> trivial then t because s
  % is a special case of r.
    if smemq('root_of,s) then nil else
  begin scalar q,!*evallhseqp,!*protfg;
    !*evallhseqp:=t; !*protfg:=t;
    q:=errorset({'subeval,mkquote {s,r}},nil,nil);
    if errorp q then <<erfg!*:=nil;return nil>>;
    q:=cdar q;
    while q and 0=reval{'difference,cadar q,caddar q} do q:=cdr q;
    return null q;
  end;

symbolic procedure groesolvidsubset!?(b1,b2,vars);
   % test if ideal(b1) is a subset of ideal(b2)
   % (b2 is a specialization of b1 wrt zeros).
  null b1 or (car b1='list or 0=preduceeval{car b1,b2,vars}) and
      groesolvidsubset!?(cdr b1,b2,vars);

symbolic procedure groesolvearb(r,vars);
 % Cover unmentioned variables.
  if atom r or not !*arbvars then r else
  for each s in r collect
   <<for each x in cdr vars do
       if not smember(x,s) then
            s:=append(s,{{'equal,x,prepf makearbcomplex()}});
     s>>;

%------------------- driver for the postprocessor ----------------

symbolic procedure groesolve0(a,vars);
  begin scalar r,ids,newvars,newa;
    if(r:=groepostnumsolve(a,vars))then return r;
    if(r:=groepostfastsolve(a,vars))then return r;
    r:=groepostsolveeval{a,vars};
    if r neq 'failed then return cdr r;
    ids:=cdr gindependent_seteval{a,vars};
    if null ids then goto nullerr;
    ids:=car ids;
    newvars:='list.for each x in cdr vars join
                 if not(x memq ids) then {x};
    newa:=groebnereval{a,newvars};
    denominators!*:=cdr glterms;
    if newa='(list 1) then rerror(groebner,24,
        "recomputation for dim=0 failed");
    r:=groepostfastsolve(newa,newvars);
    if r then return r;
    r:=groepostsolveeval{a,vars};
    if r neq 'failed then return cdr r;
nullerr:
    rerror(groebner,23,
        "Moeller ideal decomposition failed with 0 dim ideal");
    end;

symbolic procedure groepostnumsolve(gb,vars);
  if not errorp errorset('(load!-package 'roots2),nil,nil)
  and getd 'multroot0
  and get(dmode!*,'dname) member '(rounded complex!-rounded)
  and length gb = length vars and groepostnumsolve1(gb,vars)
  then (cdr reval multroot0(precision 0,gb)) where !*compxroots=t;

symbolic procedure groepostnumsolve1(gb,vars);
  if null gb then t else
  groepostnumsolve1(cdr gb,cdr vars) and
  <<for each x in kernels numr simp car gb do q:=q and x member vars;
    q>> where q=t;

symbolic procedure groepostfastsolve(gb,vars);
   % try to find a fast solution.
  begin scalar u,p1,p2,fail,kl,res;
    if !*trgroesolv then prin2t "fast solve attempt";
    groesoldmode!* := get(dmode!*,'dname);
    !*groebnumval := member(groesoldmode!* ,
                  '(rounded complex!-rounded));
    groesetdmode(groesoldmode!*,'nil);
    u:=kl:=for each p in cdr gb collect
         <<p:=numr simp p;
           intersection(vars,kernels p).p>>;
    if u='((nil)) then goto trivial;
    while u and cdr u do
    <<p1:=car u; p2:=cadr u; u:= cdr u;
      car p1:=setdiff(car p1,car p2);
      fail:=fail or null car p1>>;
    if fail then goto exit;
    res:= for each r in groepostfastsolve1(reverse kl,nil,0)
        collect 'list.reverse r;
    goto exit;
  trivial:
    res:= {'list.for each x in cdr vars collect
                {'equal,x,mvar makearbcomplex()}};
  exit:
    groesetdmode(groesoldmode!*,t);
    return res;
   end;

fluid '(f);

symbolic procedure groepostfastsolve1(fl,sub,n);
   if null fl then '(nil) else
   begin scalar u,f,v,sub1;
     n:=n+1;
     f:=car fl; v:=car f;  f:=numr subf(cdr f,sub);
     if null f then return groepostfastsolve1(cdr fl,sub,n);
     % v:=car sort(v,function(lambda(x,y);degr(f,x)>degr(f,y)));
     v := car v;
     (f:=reorder f) where kord!*={v};
     if not domainp lc f then groepostcollectden reorder lc f;
     u:=groesolvepolyv(prepf f,v);
     return for each s in u join
     <<sub1:=if smemq('root_of,s) then sub else
         (v . caddar s) . sub;
       for each q in groepostfastsolve1(cdr fl,sub1,n) collect
         car s.q
     >>;
   end;

unfluid '(f);

symbolic procedure groepostcollectden d;
  % d is a non trivial denominator (standard form);
  % collect its factors.
 for each p in cdr fctrf d do
  if not member(p:=prepf car p,denominators!*) then
      denominators!*:=p.denominators!*;

put('groesolve,'psopfn,'groesolveeval);

symbolic procedure groepostsolveeval u;
    begin scalar a,b,vars,oldorder;
          scalar groesoldb!*;
          scalar !*groebprereduce,!*groebopt,!*groesolgarbage;
        groesoldmode!* := get(dmode!*,'dname);
        groesetdmode(groesoldmode!*,'nil);
        !*groebnumval := member(groesoldmode!* ,
                          '(rounded complex!-rounded));
        if vdpsortmode!* = 'lex then t else
          rerror(groebner,8,
                 "Groepostproc, illegal torder; (only LEX allowed)");
        a := groerevlist reval car u;
        vars := cdr u and groerevlist cadr u or groebnervars(a,nil);
        oldorder := setkorder vars;
        b :=  groesolve1(a,a,vars);
        a := nil;
        if b eq 'failed then a:=b else
        <<for each sol in b do % delete multiplicities
            if not member(sol,a) then a := sol . a;
          a := 'list . for each sol in a collect 'list . sol;
        >>;
        setkorder oldorder;
        groesetdmode(groesoldmode!*,t);
        return a;
    end;

put('groepostproc ,'psopfn,'groepostsolveeval);


% DATA structure:
%
%  all polynomials are held in prefix form (expressions)
%  transformation to standard quotients/ standard forms is done locally
%  only; distributive form is not used here.
%
%  a zero is a set of equations, if possible with a variable on the
%  lhs each
%       e.g. {y=17,z=b+8};
%        internally: ((equal y 17)(equal z (plus b 8)))
%  a zeroset is a list of zeros
%       elgl {{y=17,z=b+8},{y=17,z=b-8}}
%  internally the sets (lists) are kept untagged as lists; the
%  tag 'list is only added to the results and to those lists which
%  are parameters to algebraic operators not in this package.

symbolic procedure groesolve1 (a,fulla,vars);
 %  A        lex Groebner basis or tail of lex Groebner basis
 %  fullA    the complete lex Groebner basis to A
 %  vars     the list of variables
  if null a or a='(1) then nil else
  <<begin scalar f1,a1,res,q,gi,ng1,ng2,ngq,qg,mv,test;
      res := assoc (a,groesoldb!*); if res then return cdr res;
      groesolvelvevel!* := groesolvelvevel!* + 1;
      %%tr prin2t  "=================================================";
      %%tr prin2l list( " groesolvelvevel!*: ",groesolvelvevel!*);
      %%tr           prin2t " Problem:";
      %%tr                   writepri (mkquote ('list . a),'only);;
      if member(a,!*groesolrecurs) then return 'failed;
%        <<vars := append(cdr vars,{car vars});
%          if member(vars,!*groesolrecurs) then
%          <<!*groesolgarbage := T;
%            return list for each p in A collect list('EQUAL,p,0) >>;
%          !*groesolrecurs := vars . !*groesolrecurs;
%          A := cdr Groebnereval{'list.A,'list . vars};
%        >>;
      !*groesolrecurs := a . !*groesolrecurs;
      if length a = 1 then
       << %%tr print  "single polynomial; solve it";
          res := groesolvepoly car a; goto ready>>;
         % step 1
      f1 := car a;
      a1 := cdr a;
      test := nil;
      mv := intersection(vars,ltermvariables f1);    % test Buchcrit 4
      for each p in a1 do
         if intersection (mv,ltermvariables p) then test := t;
      if not test then
      << %%tr print "f1 has unique main variable";
         ngq := groesolve1(a1,a1,vars);
         if ngq eq 'failed then <<res:='failed;goto ready>>;
         res := zerosetintersection(ngq,f1,vars);
         goto ready>>;
%     Q := cdr GroebIdq('list . A1,f1,'list . vars);    % A1:f1
      q := groesolvidq(a1,f1,vars);
        %%tr print "A1:f1"; %%tr writepri (mkquote ('list . Q),'only);;
      if q='(1) then     % f1 already member of A1; skip it
        <<%%tr print "f1 already in A1; ignore";
          res:= groesolve1(a1,fulla,vars);
          goto ready>>;
      ngq := groesolve1(q,q,vars);
      if ngq eq 'failed then <<res:='failed;goto ready>>;
      ng1 := zerosetintersection (ngq,f1,vars);
         % step 4
      if groesolvidqtest(a1,f1,vars) then
      << while q do
         << gi := car q; q := cdr q;
            gi := preduceeval list (gi , 'list . a, 'list . vars);
            if gi = 0 then q:= nil else
               a := cdr groebidq('list . a ,gi,'list . vars);
         >>;
         ng2 := groesolve1(a,a,vars);
         if ng2 eq 'failed then <<res:='failed;goto ready>>;
      >> else
      <<ng2 := ();
         if length q = 1 then
         << gi := preduceeval list (car q, 'list . fulla, 'list . vars);
           if gi neq 0 then
           << qg := cdr groebidq('list . fulla,gi,'list . vars); % A1:gi
              ng2 := groesolve1(qg,qg,vars);
              if ng2 eq 'failed then <<res:='failed;goto ready>>;
         >> >>
         else
            <<ng2 := groesolve2(a1,q,vars);
              if ng2 eq 'failed then <<res:='failed;goto ready>>
            >>
      >>;
      res:= zerosetunion (ng1,ng2);
    ready:
      %%tr print list( "Result, level!*: ",groesolvelvevel!*);
      %%tr writeres res;
      %%tr print "...................................................";
      groesolvelvevel!* := groesolvelvevel!* -1;
      groesoldb!* := (a . res) . groesoldb!*;
      return res;
   end
   >> where !*groesolrecurs = !*groesolrecurs; % recursive fluid!

symbolic procedure groesolvidqtest(a1,f1,vars);
      not(deg(f1,car vars) eq deg(car a1,car vars));

symbolic procedure groesolvidq(a1,f1,vars);
      begin scalar x,temp;
         x := car  vars;
         if not groesolvidqtest(a1,f1,vars)
             then return cdr groebidq('list . a1,f1,'list . vars);
         temp :=
           for each p in a1 collect
             reval car reverse coeffeval list(p,x);
         return cdr groebnereval list ('list . temp,'list . vars);
       end;

symbolic procedure groesolve2(a,q,vars);
  % Calculation of the zeroset A1:(g1,g2,,,,gs),
  % the gi given as members of Q.
     groesolvetree (a,q,q,vars);

symbolic procedure groesolvetree(a,t1,phi,vars);
   begin scalar q,ngtemp,ngall,t2,h,g,a2,phi2;
      %%tr prin2t ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
      %%tr prin2t "groesolvetree; A:";
      %%tr writepri(mkquote ('list . a),'only);
      %%tr writepri( "T1 :",'car);
      %%tr writepri(mkquote ('list . T1) ,'last);
      %%tr writepri( "phi:",'car);
      %%tr writepri(mkquote ('list .  phi),'last);
       if null phi then return nil else
       if null t1 then
          <<q := cdr groebidq('list.a, 'times.phi,'list.vars);
            %%tr prin2t "<< << << << << << << << << << << << << <<<";
            return if car q = 1 then nil
                                else groesolve1(q,q,vars) >>;
       for each g in t1 do
       <<h:=preduceeval list(g,'list.a,'list.vars);
         phi := delete(g,phi);
         if not(h=0) then <<t2:=h.t2; phi:=h.phi>>;
       >>;
       if null phi then return nil;  % 29.8.88
       t1 := t2;
       q := cdr groebidq('list.a, 'times.phi,'list.vars);
        %%tr writepri( "PHI :",'car);
        %%tr writepri(mkquote ('TIMES.phi) ,'last);
        %%tr writepri( "A:PHI :",'car);
        %%tr  writepri(mkquote ('LIST.Q) ,'last);
       if not(car q = 1) then
       << ngall := groesolve1(q,q,vars);
          if ngall eq 'failed then return 'failed;
       >>;
       if !*groesolgarbage then
               return groesolverestruct(q,phi,vars,ngall);
        while t1 do
          <<g:=car t1; t1:=cdr t1;
            phi2 := delete(g,phi);
            if phi2 then
            <<a2 := cdr groebnereval list('list . g . a,'list . vars);
              if not(car a2 =1) then
               <<ngtemp := groesolvetree(a2,t1,phi2,vars);
                 ngall := zerosetunion(ngtemp,ngall)>>;
           >>>>;
       %%tr prin2t "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<";
       return ngall;
   end;

symbolic procedure groesolverestruct(a,phi,vars,ngall);
     % there was a problem with an embedded solution in phi such that
     % A:phi = A
     % we try a heuristic by making one variable a formal parameter
    begin scalar newa,newvars,mv,oldorder,solu;
        mv := ltermvariables ('times.phi);
        mv := car mv;
         %%tr prin2 "groesolverestruct:taking variable ";prin2t mv;
        newvars := delete(mv,vars);
        oldorder := setkorder newvars;
        newa := cdr groebnereval list('list . a,'list . newvars);
         %%tr    prin2t "new Groebner basis:";
         %%tr    writepri(mkquote ('LIST . newA),'ONLY);
        !*groesolgarbage := nil;
        solu := groesolve1(newa,newa,newvars);
         %%tr if !*groesolgarbage then prin2t "Heuristics failed"
         %%tr                     else prin2t "better result";
        setkorder oldorder;
        return if !*groesolgarbage then ngall else solu;
    end;



%%tr symbolic procedure writeres r;
%%tr writepri (mkquote ('list.for each x in r collect 'list.x),'ONLY);

symbolic procedure ltermvariables u;
  % extract variables of leading term in u
    begin scalar v;
      u := numr simp u;
      while not domainp u do
      <<v := mvar u . v;
        u := lc u>>;
      return reversip v;
    end;

symbolic procedure zerosetintersection (ng,poly,vars);
     % NG is a zeroset
     % poly is a polynomial
     % the routine maps the zeros in NG by the polynomial:
     %   each zero is substituted into the polynomial;
     %   that gives a univariate
     %   solved by SOLVE or numerical techniques.
     % the result is the solution NG, including the solutions of the
     % polynomial.
   begin scalar res,ns,testpoly,ppoly,sol,s,var,dom;
      %%tr print "Intersection ";
      %%tr writepri (mkquote ('list . NG),'only);
      %%tr         print " with ";
       %%tr writepri(mkquote poly,'only);
      res := ();
      poly := simp poly;
      var:=if not domainp numr poly then groesolmvar(numr poly,vars)
             else 'constant;
    loop:
      if ng=() then goto finish;
      ns := car ng; ng := cdr ng;
      testpoly := poly;
      dom := groesoldmode!* or 'rational;
      groesetdmode(dom,t);
      testpoly := simp prepsq testpoly;
      for each u in ns do
       if idp lhs u and not smemq('root_of,rhs u) then
       <<s := rhs u;
         testpoly := subsq(testpoly,list (lhs u . s));
       >>;
      groesetdmode(dom,nil);
      ppoly := prepf numr testpoly;
      sol := groesolvepolyv(ppoly,var);
      res := append(res , for each r in sol collect append(r,ns) ) ;
      goto loop;
    finish:
      %%tr print  "Schnittresultat: ";
      %%tr writeres res;
      return res;
   end;

symbolic procedure groesolmvar(poly,vars);
  % select main variable wrt vars sequence.
   <<while vars and not smember(car vars,poly) do
      vars:=cdr vars;
     if null vars then rerror('groebner,27,"illegal poly");
     car vars>>;


% solving a single polynomial with respect to its main variable

symbolic procedure groesolvepoly p; groesolvepolyv(p,mainvar p);

symbolic procedure groesolvepolyv(p,var);
   % find the zeros for one polynomial p in the
   % variable "var".
   % current dmode is NIL.
  (begin scalar res,u,!*convert,y,z;
      if (u:=assoc(var,depl!*)) then depl!*:=delete(u,depl!*);
      if !*trgroesolv then
      <<writepri("   solving univariate with respect to ",'first);
        writepri(mkquote var,'last);
        writepri(mkquote p,'only);
      >>;
      for each s in groebroots!* do
       if 0=reval{'difference,p,car s} then res:=cdr s;
      if res then return res;
      groesetdmode(groesoldmode!*,t);
      u := numr simp p;
      res := if !*groebnumval and univariatepolynomial!? u
               then groeroots(p,var)
              else (solveeval list (p,var))
                where kord!*=nil,
                alglist!* = nil . nil;
      res := cdr res;
       % Collect nontrivial denominator factors.
       % Reorder for different local order during solveeval.
      for each x in res do
      <<y:=prepf (z:=reorder denr simp caddr x);
        if dependsl(y,variables!*) then groepostcollectden z;
      >>;
      res := for each x in res collect list x;
      groesetdmode(groesoldmode!*,nil);
      return res;
  end) where depl!*=depl!*;

symbolic procedure univariatepolynomial!? fm;
      domainp fm or univariatepolynomial!?1 (fm,mvar fm);

symbolic procedure univariatepolynomial!?1 (fm,v);
       domainp fm or
         domainp lc fm and v = mvar fm and
                           univariatepolynomial!?1(red fm,v);

symbolic procedure predecessor (r,l);
  % looks for the  predecessor of r in l
   if not pairp l or not pairp cdr l or r = car l
          then rerror(groebner,9,"No predecessor available") else
   if r = cadr l then car l else predecessor (r,cdr l);

symbolic procedure zerosetunion(ng1,ng2);
     <<%print list( "union von ",ng1,ng2;
       ng1 := zerosetunion1(ng1,ng2);
     % print list( "-->",ng1);
       ng1>>;

symbolic procedure zerosetunion1(ng1,ng2);
    % Vereinigung von Nullstellengebilden
      if ng1 = () then ng2
            else
      if zerosetmember(car ng1,ng2) then zerosetunion1(cdr ng1,ng2)
            else
      car ng1 . zerosetunion1(cdr ng1,ng2);

symbolic procedure zerosetmember (ns,ng);
      if ng = () then nil else
      if zeroequal(ns,car ng) then ng else
         zerosetmember (ns,cdr ng);

symbolic procedure zeroequal(ns1,ns2);
      if zerosubset(ns1,ns2) then zerosubset(ns2,ns1) else nil;

symbolic procedure zerosubset(ns1,ns2);
      if null ns1 then t else
      if member(car ns1,ns2) then zerosubset(cdr ns1,ns2)
         else nil;

symbolic procedure groesetdmode(dmode,dir);
  % Interface for switching an arbitrary domain on/off.
  % Preserve complex mode. Turn on EZGCD whenever possible.
   if null dmode then nil else
   begin scalar !*msg,x,y;
   if null dir then
   << if !*complex then y:=setdmode('complex,nil);
      !*complex := nil;
      if dmode!* = '!:rd!: then !*rounded :=nil;
      if dmode!* then y:=setdmode(get(dmode!*,'dname),nil);
      if !*groebcomplex then
      <<setdmode('complex,t); !*complex:=t>>;
   >>
   else
   << if memq(dmode,'(complex complex!-rounded complex!-rational))
        then <<!*complex := t; y:=setdmode('complex,t);
                 if (x:=atsoc(dmode,'((complex!-rounded . rounded)
                           (complex!-rational . rational)) ))
                 then y:=setdmode(cdr x,t)>>
           else y:=setdmode(dmode,t);
      if memq(dmode,'(rounded complex!-rounded)) then !*rounded :=t;
    >>;
    !*ezgcd := null dmode!*;
    return y;
  end;

symbolic procedure preduceeval pars;
%  Polynomial reduction driver. u is an expression and v a list of
% expressions. Preduce calculates the polynomial u reduced wrt the list
% of expressions v.
% parameters:
%     1      expression to be reduced
%     2      polynomials or equations; base for reduction
%     3      optional: list of variables
   begin scalar vars,x,u,v,w,oldorder;
         scalar !*factor,!*exp,!*gsugar,!*vdpinteger;
         integer n,pcount!*; !*exp := t;
      if !*groebprot then groebprotfile := list 'list;
      n := length pars;
      x := reval car pars;
      u := reval cadr pars;
      v := if n>2 then reval caddr pars else nil;
      w := for each j in groerevlist u
              collect if eqexpr j then !*eqn2a j else j;
      if null w then rerror(groebnr2,3,"Empty list in Preduce");
      vars := groebnervars(w,v);
      if not vars then vdperr 'preduce;
      oldorder := vdpinit vars;
      w := for each j in w collect a2vdp j;
      x := a2vdp x;
      if !*groebprot then
          <<w := for each j in w collect vdpenumerate j;
            groebprotsetq('candidate,vdp2a x);
            for each j in w do groebprotsetq(mkid('poly,vdpnumber j),
                                             vdp2a j)>>;
      w := groebnormalform(x,w, 'sort);
      w := vdp2a w;

      setkorder oldorder;
      return if w then w else 0;
   end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% the following code is the interface to Stan's rootfinder
%

symbolic procedure groeroots(p,x);
   begin scalar r;
      r := reval{'roots,p};
         % re-evaluate rhs in order to get prefix form
      r := for each e in cdr r collect
            list('equal,cadr e,reval caddr e);
      return 'list . r;
   end;

endmodule;


module groebopt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   optimization of the sequence of variables
%

% Optimization of variable sequence;
% the theoretical background can be found
% in  Boege/Gebauer/Kredel , J.Symb.Comp(1986)I,83-98
% Techniques modfied to the following algorithm
%
% x > y if
%  x appears in a higher power than y
%       or
%  the highest powers are equal, but x appears more often with
% that power.
%
% An explicit dependency DEPENDS X,Y will supersede the optimality.

fluid '(vdpvars!*     % list of intersting variables (and korder)
        intvdpvars!*
        !*trgroeb     % control of print messages
       );

symbolic procedure vdpvordopt (w,vars);
  % w : list of polynomials (standard forms)
  % vars: list of variables
  % returns (w . vars), both reorderdered
  begin; scalar c;
        vars := sort(vars,'ordop);
        c := for each x in vars collect x . 0 . 0;
        for each poly in w do vdpvordopt1 (poly,vars,c);
        c := sort(c,function vdpvordopt2);
        intvdpvars!* := for each v in c collect car v;
        vars := sort(intvdpvars!*,function vdpvordopt3);
        if !*trgroeb then
        <<prin2 "optimized sequence of kernels: "; prin2t vars>>;
        return (for each poly in w collect reorder poly) . vars;
  end;


symbolic procedure vdpvordopt1 (p,vl,c);
   if null p then 0
      else
   if domainp p or null vl then 1
      else
   if mvar p neq car vl then vdpvordopt1(p,cdr vl,c)
      else
   begin scalar var,pow,slot; integer n;
      n := vdpvordopt1 (lc p,cdr vl,c);
      var := mvar p;
      pow := ldeg p;
      slot := assoc(var,c);
      if pow #> cadr slot then
          <<rplaca(cdr slot,pow); rplacd(cdr slot,n)>>
        else
          rplacd(cdr slot,n #+ cddr slot);
      return n #+ vdpvordopt1 (red p,vl,c);
   end;

symbolic procedure vdpvordopt2(sl1,sl2);
  % compare two slots from the power table
    <<sl1 := cdr sl1;
      sl2 := cdr sl2;
      car sl1 #< car sl2 or car sl1 = car sl2 and cdr sl1 #< cdr sl2
    >>;

symbolic procedure vdpvordopt3(x,y);
  % two level compare: 1. dependency, 2. optimal order.
     depends(x,y) or member(y,member(x,intvdpvars!*));

endmodule;


module groebsea;

% support of search for reduction polynomials

fluid '(thirdvalue!* fourthvalue!* hcount!* !*groebweak);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  search for reduction candidates in a list

symbolic procedure groebsearchinlist (vev,g);
  % search for a polynomial in the list G, such that the lcm divides
  % vev; G is expected to be sorted in descending sequence
    if null g then nil
    else if buch!-vevdivides!?(vdpevlmon car g, vev) then car g
    else groebsearchinlist (vev,cdr g);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  search tree for polynomials
%  simple variant: mapped to search list
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure groebstreeadd (poly,stru);
% add one polynomial to the tree
% if this is a simple polynomial (mono or bino), reform
% the tree
   if hcount!* #< 5000 then vdplsortin(poly,stru)
     else vdplsortinreplacing(poly,stru);

symbolic procedure groebsearchinstree (vev,stru);
% search a polynomial corresponding to the exponent vector vev
     groebsearchinlist (vev,stru);

symbolic procedure groebstreeextract stru;
% gives a linear list of all polynomials in the tree
                 stru;

symbolic procedure groebstreereconstruct u;
% reconstructs a tree from a linear list of polynomials
               vdplsort u;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% reforming G, D and G99 when a very simple polynomial was
% found (e.g. a monomial, a binomial)

symbolic procedure groebsecondaryreduction(poly,g,g99,d,gc,mode);
    % if poly is a simple polynomial, the polynomials in G and G99
    % are reduced in a second pass. Result is G, secondvalue is G99.
    % mode says, that G99 has to be modified in place.
  begin scalar vev,p,pl,x,rep,first,rpoly,break;   % y;
    secondvalue!* := g99; thirdvalue!* := d; fourthvalue!* := gc;
    vev := vdpevlmon poly; rpoly := vdpred poly;
        % cancel redundant elements in G99
    for each p in g99 do if buch!-vevdivides!?(vev,vdpevlmon p)
                  then g99:=delete(p,g99);
    if vdplength poly > 2 or  vevzero!? vev then return g;
    if !*groebweak and not vdpzero!? rpoly
       and (groebweaktestbranch!=1(poly,g,d)) then return 'abort;
        !*trgroeb and groebmess50 g;
    pl := union(g,g99);
    first := t;

    while pl and not break do
      << p:= car pl; pl := cdr pl;
         if groebprofitsfromvev(p,vev) then
                          % replace by simplified version
            <<x := groebnormalform1(p,poly);
              x := groebsimpcontnormalform x;
              x := vdpenumerate x;
                 if first then !*trgroeb and groebmess20(poly);
                      first := nil;
                      !*trgroeb and groebmess21(p,x);
                  rep := (p . x) . rep;
              if not vdpzero!? x and
                 vevzero!? vdpevlmon x then break := t; % 1 found
             >> >>;

    if break then return 'abort;

                           % reform G99
    g99 := for each p in g99 collect groebsecondaryreplace(p,rep);
    secondvalue!* :=  groebsecondaryremovemultiples g99;
                           % reform D
 %  thirdvalue!* := nil;
 %  for each p in D do
 %  <<x := groebSecondaryReplace(cadr p,rep);
 %    y := groebSecondaryReplace(caddr p,rep);
 %    thirdvalue!*:=groebcplistSortIn(groebmakepair(x,y),thirdvalue!*);
 %  >>;
    thirdvalue!* := d;
                          % reform Gc
    fourthvalue!* :=
      groebsecondaryremovemultiples
        for each y in gc collect groebsecondaryreplace(y,rep);

    g:=for each y in g collect groebsecondaryreplace(y,rep);
         !*trgroeb and groebmess50 g;
      % print('groebSecondaryReduction .
      %          for each p in g collect vdpnumber p);
    return groebsecondaryremovemultiples g;
  end;

symbolic procedure groebsecondaryremovemultiples g;
    if null g then nil else
    if vdpzero!? car g or member(car g,cdr g) then
                groebsecondaryremovemultiples cdr g else
    car g  . groebsecondaryremovemultiples cdr g;

symbolic procedure groebsecondaryreplace(x,rep);
   (if y then cdr y else x) where y = atsoc(x,rep);

endmodule;


module groebsor;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  maintenance of lists of critical pairs (sorting etc.)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure groebcplistsortin (p,pl);
% Distributive polynomial critical pair list sort. pl is a
% special list for Groebner calculation, p is a pair.
% returns the updated list pl (p sorted into);
% if null pl then list p
%   else
% if vevcompless!? (car car pl, car p)
%     then car pl  . groebcplistsortin (p, cdr pl)
%   else
%      p . pl;
%  non copying version:
  if null pl then list p
    else
  <<groebcplistsortin1 (p,pl); pl>>;

symbolic procedure groebcplistsortin1(p,pl);
% destructive insert of p into nonnull pl
   if not groebcpcompless!?(car pl, p)
        then <<rplacd(pl,car pl . cdr pl); rplaca(pl,p)>>
      else
   if null cdr pl then rplacd(pl,list p)
      else
        groebcplistsortin1(p,cdr pl);


symbolic procedure groebcplistsort g;
  <<for each p in g do gg:=groebcplistsortin(p,gg); gg>>
      where gg=nil;

symbolic procedure groebcplistmerge(pl1,pl2);
% Distributive polynomial critical pair list merge. pl1 and pl2
% are critical pair lists used in the Groebner calculation.
% groebcplistmerge(pl1,pl2) returns the merged list.
begin scalar cpl1,cpl2,sl;
      if null pl1 then return pl2;
      if null pl2 then return pl1;
      cpl1 := car pl1; cpl2 := car pl2;
      sl := groebcpcompless!?(cpl1, cpl2);
      return
           (if sl then cpl1 . groebcplistmerge(cdr pl1,pl2)
                  else cpl2 . groebcplistmerge(pl1,cdr pl2) )
end;

symbolic procedure groebcpcompless!?(p1,p2);
% compare 2 pairs srt their sugar(=cadddr) or their lcm (=car).
    if !*gsugar then
         (if not(d=0) then d<0 else
          if not(q=0) then q<0 else
            vdpnumber(caddr p1)<vdpnumber(caddr p2)
         )  where d= cadddr p1 - cadddr p2, q=vevcomp(car p1,car p2)
       else vevcompless!?(car p1,car p2);

endmodule;


module groebspa;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% manipulation of subspaces.
% A subspace among the variables is described by an exponent vector
%  with only zeroes and ones. It terminates with the last
% one. It may be null (nil);

expr procedure vevunion(e1,e2);
    begin scalar x,y;
       y := vevunion1(e1,e2);
       x := reverse y;
       if car x = 1 then return y;
       while x and car x = 0 do x := cdr x;
       return reversip x;
    end;

expr procedure vevunion1(e1,e2);
    if vdpsubspacep(e1,e2) then e2
       else
    if vdpsubspacep(e2,e1) then e1
       else
    if car e1 neq 0 or car e2 neq 0 then 1 . vevunion1(cdr e1,cdr e2)
       else
          0 . vevunion1(cdr e1,cdr e2);

expr procedure vdpsubspacep(e1,e2);
% test if e1 describes a subspace from e2
    if null e1 then t
      else
    if null e2 then vdpspacenullp(e1)
      else
    if car e1 > car e2 then nil
      else
    if e1 = e2 then t
      else
         vdpsubspacep(cdr e1,cdr e2);

expr procedure vdporthspacep(e1,e2);
% test if e1 and e2 describe orthogonal spaces(no intersection);
    if null e1 or null e2 then t
       else
    if car e2 = 0 or car e1 = 0
          then vdporthspacep(cdr e1,cdr e2)
       else nil;


expr procedure vdpspacenullp(e1);
% test if e1 describes an null space
     if null e1 then t
        else
     if car e1 = 0 then vdpspacenullp(cdr e1)
        else nil;


expr procedure vdpspace(p);
  % determine the variables of the polynomial.
    begin scalar x,y;
        if vdpzero!? p then return nil;
        x := vdpgetprop(p,'subroom);
        if x then return x;
        x := vevunion(nil,vdpevlmon p);
        y := vdpred p;
        while not vdpzero!? y do
         <<x := vevunion(x,vdpevlmon y);
           y := vdpred y>>;
        vdpputprop (p,'subroom,x);
        return x;
    end;

symbolic procedure vdpunivariate!?(p);
    if vdpgetprop(p,'univariate) then t
     else begin scalar ev; integer n;
         ev := vdpevlmon p;
         for each x in ev do
           if not(x=0) then n := n#+1;
         if not(n=1) then return nil;
         ev := vdpspace(p);
         for each x in ev do
           if not(x=0) then n := n#+1;
         if not(n=1) then return nil;
         vdpputprop(p,'univariate,t);
         return t;
     end;

endmodule;


module groebfac; % Factorization of polynomials during Groebner calc'n.

% create!-package('(groebfac),'(contrib groebner));

% Other packages needed.

% load!-package 'vdp2dip;

imports groebner,vdp2dip,factor;

fluid '(groebactualg99!* groebactualg!* factortime!* !*vdpmodular
          vdpone!* groebfabort!* !*factor !*groebrm current!-modulus
          !*trgroeb !*gsugar);

symbolic procedure groebfactorize (h,abort1,g,g99);
  begin scalar r,tim,gctim,h1,groebactualg99!*,groebfabort!*,test;
    scalar s;
    s:=!*gsugar and gsugar h;
    groebactualg99!* := g99; groebactualg!* := g;
    groebfabort!* := abort1;
    if vdpgetprop(h,'irreducible) then return groebfactorize3 h;
    tim := time();
    gctim := gctime();
         !*trgroeb and groebmess7 h;
    r := if r := vdpgetprop(h,'factors) then r
     else if !*groebrm then groebfactorize1 h
     else if not !*vdpmodular then groebfactorize2 h
     else nil;
    factortime!* := factortime!* + time() - tim -(gctime()-gctim);
    if null r then <<vdpputprop(h,'irreducible,t);
                     return groebfactorize3 h>>;
               if cdr r then !*trgroeb and groebmess14 (h,r);
    vdpputprop(h,'factors,r);
    for each p in r do
       if vdpmember(car p,g) then test:= car p;
    if test then
       <<!*trgroeb and groebmess27a(h,test); return 'zero>>;
    h1 := car r;
    for each p in r do
       if vdpmember(car p,abort1) then
            <<r := delete(p,r); !*trgroeb and groebmess27 car p >>
       else vdpputprop(car p,'irreducible,t);
    if null r then r := list h1;     % at least one
    if null cdr r then groebfactorize3 caar r;
       % inherit sugar if no substantial factor.
    if !*gsugar then
     if null cdr r then gsetsugar(caar r,s) else
       for each p in r do gsetsugar(car p,vdptdeg car p);
    return 'factor . r;
  end;


symbolic procedure groebfactorize1 h;
  % factorize: separate monomial factors which were detected already;
   begin scalar monf,vp,n,e,h1,h2,vp2;
        monf := vdpgetprop(h,'monfac);
        if null monf then
             return if not !*vdpmodular then groebfactorize2 h
                                       else nil;     % no factor
        h2 := vdpdivmon (h,vbcfi 1,monf);

       if groebmonfac neq 0 then
       <<                                 % now build a polynomial from
        n := 0;                           % each variable in MONFAC
        for each x in monf do
            <<n := n#+1;
              if x neq 0 then
                <<e := list x;
                  for i:=2:n do e := 0 . e; % prefix with n-1 zeros.
                  vp := vdpfmon(a2vbc 1,e) . vp;
                >>;
            >>;
       >>
        else
        !*trgroeb and groebmess15 monf;
                  % append body of orig. poly, factorized
        if not vdpzero!? h2 and
           not vevzero!? vdpevlmon h2 then
               <<if not !*vdpmodular then vp2 := groebfactorize2 h2;
                 vp2 := if not vp2 then list h2
                     else for each v in vp2 collect car v;
                 vp := nconc(vp,vp2)>>;
                   % ascending sorting
    %   if length vp = 1 then return nil;
        h1 := vp;
        return
               reverse  for each x in h1 collect list vdpenumerate x;
   end;

symbolic procedure groebfactorize2 h;
  % tries to factorize a h-polynomial via REDUCE factorizer
   begin scalar h1,h2,!*factor; !*factor := t;
%       h1 := vdp2a vdprectoint (h,vdplcm h);
%       h1 := fctrf !*q2f simp h1;   % factorf
        h1 := groefctrf vdp2f h;
        if null cdr h1 then return nil;
        if null cddr h1      % only one element in factorization list
           and cdr cadr h1 = 1     % and multiplicity = 1
                then return nil;
        h2 := for each l in cdr h1 join
            for i:=1:cdr l collect car l;
        h2 := vdplsort for each p in h2 collect vdpsimpcont f2vdp p;
        return for each x in h2 collect list vdpenumerate x;
   end;

symbolic procedure groefctrf p;
   (fctrf p) where !*factor=t,current!-modulus = current!-modulus;

symbolic procedure groebfactorize3 h;
  % additional efforts to factor something.
       nil;

endmodule;


module groebidq;

% calculation of ideal quotient using a modified Buchberger algorithm.

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

fluid '(                           % switches from the user interface
        !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
        !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
        !*groebstat !*groebdivide !*groebprot
        !*groebidqbasis

        vdpvars!*                  % external names of the variables
        !*vdpinteger !*vdpmodular  % indicating type of algorithm
        vdpsortmode!*              % term ordering mode
        secondvalue!* thirdvalue!* % auxiliary: multiple return values
        fourthvalue!*
        factortime!*               % computing time spent in factoring
        factorlvevel!*              % bookkeeping of factor tree
        pairsdone!*                % list of pairs already calculated
        probcount!*                % counting subproblems
        vbccurrentmode!*            % current domain for base coeffs.
        groetime!*
        !*gsugar

       );


switch groebopt,groebfac,groebrm,groebres,trgroeb,trgroebs,trgroeb1,
       trgroebr,groebstat,groebprot;
!*groebidqbasis := t;              %default: basis from idq

% variables for counting and numbering
fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
        basecount!* hzerocount!*);


symbolic procedure groebidq2 (p,f);
  % setup all global variables for the Buchberger algorithm
  % printing of statistics
 begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*,
       pairsdone!*, factorlvevel!*,!*gsugar;
       factortime!* := 0;
       groetime!* := time();
       vdponepol(); % we construct dynamically
       hcount!* := 0;
       pcount!* := 0;
       mcount!* := 0;
       fcount!* := 0;
       bcount!* := 0;
       b4count!* := 0;
       hzerocount!* := 0;
       basecount!* := 0;
       if !*trgroeb then
       <<
          prin2 "IDQ Calculation starting ";
          terprit 2;
       >>;
       spac := gctime();

       p1:=  groebidq3 (p, f);

       if !*trgroeb or !*trgroebr or !*groebstat then
       <<
          spac1 := gctime() - spac;
          terpri();
          prin2t "statistics for IDQ calculation";
          prin2t "==============================";
          prin2 " total computing time (including gc): ";
          prin2((tim1 := time()) - groetime!*);
          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 if !*groebidqbasis then car groebner2(p1,nil) else p1;
 end;


symbolic procedure groebidq3 (g0,fff);
   begin scalar result,x,g,d,d1,d2,p,p1,s,h,g99,one,gi;
         gi := g0;
         fff := vdpsimpcont fff;
         vdpputprop(fff,'number,0);         % assign number 0
         vdpputprop(fff,'cofact,a2vdp 1);   % assign cofactor 1
         x := for each fj in g0 collect
           << fj:=vdpenumerate vdpsimpcont fj;
              vdpputprop(fj,'cofact,a2vdp 0);   % assign cofactor 0
              fj>>;
         g0 := list fff;
         for each fj in x do g0 := vdplsortin(fj,g0);

% ITERATION :
      while (d or g0) and not one do
       begin
          if g0 then
          <<          % take next poly from input
               h := car g0;
               g0 := cdr g0;
               p := list(nil,h,h);
          >>
             else
          <<          % take next poly from pairs
               p := car d;
               d := delete (p,d);
               s := idqspolynom (cadr p, caddr p);
                          idqmess3(p,s);
               h := idqsimpcont idqnormalform(s,g99,'tree);
               if vdpzero!? h then
               <<      !*trgroeb and groebmess4(p,d);
                 x := vdpgetprop(h,'cofact);
                 if not vdpzero!? x then
                    if vevzero!? vdpevlmon x then one:= t else
                     << result := idqtoresult(x , result);
                        idqmess0 x>>;
               >>
          >>;
          if vdpzero!? h then goto bott;
          if vevzero!? vdpevlmon h then % base 1 found
                  <<          idqmess5(p,h);
                 result:=gi; d:=g0:=nil;
                 goto bott;>>;
          s:= nil;
                  % h polynomial is accepted now
               h := vdpenumerate h;
                        idqmess5(p,h);
                              % construct new critical pairs
               d1 := nil;
               for each f in g do
                    <<d1 := groebcplistsortin( list(tt(f,h),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;
                    if groebbuchcrit4t(cadr p1,caddr p1)
                       then d2 := append (d2, list p1)
                           else
                    <<x := idqdirectelement(cadr p1,caddr p1);
                      if not vdpzero!? x then
                         if vevzero!? vdpevlmon x then one:= t else
                              <<idqmess1 (x,cadr p1,caddr p1);
                                result := idqtoresult(x,result);
                                >> >>;

                    d1 := groebinvokecritm (p1,d1);
                  >>;
             %   D := groebInvokeCritB (h,D);
                 d := groebcplistmerge(d,d2);
                 g := h . g;
                 g99 := groebstreeadd(h, g99);
                            !*trgroeb and groebmess8(g,d);
    bott:
   end;      %  ITERATION
                      % now calculate groebner base from quotient base
   if one then result := list vdpfmon(a2vbc 1,vevzero());
        idqmess2 result;
   return result;
end;     % MACROLOOP

symbolic procedure idqtoresult(x,r);
    % x is a new element for the quotient r;
    % is is reduced by r and then added.
      <<x := groebsimpcontnormalform groebnormalform(x,r,'sort);
        if vdpzero!? x then r else vdplsortin(x,r)>>;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Reduction of polynomials
%

symbolic procedure idqnormalform(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
     % version for idealQuotient: doing side effect calculations for
     % the cofactors; only headterm reduction
  begin scalar c,vev,divisor,done,fold;
      fold := f;
      while not vdpzero!? f and g do
          begin
            vev:=vdpevlmon f; c:=vdplbc f;
            if type = 'sort then
                 while g
                       and vevcompless!? (vev,vdpevlmon (car g))
                            do g := cdr g;
            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 !*vdpinteger  then
                          f := idqreduceonestepint(f,nil,c,vev,divisor)
                       else
                          f := idqreduceonesteprat(f,nil,c,vev,divisor)
               else
                  g := nil;
          end;
      return if done then f else fold; %in order to preserve history
   end;



% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  special reduction procedures



symbolic procedure idqreduceonestepint(f,dummy,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,fcofa,gcofa;
      fcofa := vdpgetprop(f,'cofact);
      gcofa := vdpgetprop(g1,'cofact);
      vevlcm := vevdif(vev,vdpevlmon g1);
      cg := vdplbc g1;
            % calculate coefficient factors
      x := vbcgcd(c,cg);
      a := vbcquot(cg,x);
      b := vbcquot(c,x);
      f:= vdpilcomb1 (vdpred f, a, vevzero(),
                      vdpred g1,vbcneg b, vevlcm);
      x := vdpilcomb1 (fcofa, a, vevzero(),
                       gcofa ,vbcneg b, vevlcm);
      vdpputprop(f,'cofact,x);
      return f;
 end;


symbolic procedure idqreduceonesteprat(f,dummy,c,vev,g1);
      % reduction step for rational case:
      % calculate f= f - g/vdpLbc(f)
      %
 begin scalar x,fcofa,gcofa,vev;
      fcofa := vdpgetprop(f,'cofact);
      gcofa := vdpgetprop(g1,'cofact);
      vev := vevdif(vev,vdpevlmon g1);
      x := vbcneg vbcquot(c,vdplbc g1);
      f := vdpilcomb1(vdpred f,a2vbc 1,vevzero(),
                      vdpred g1,x,vev);
      x := vdpilcomb1 (fcofa,a2vbc 1 , vevzero(),
                       gcofa,x,vev);
      vdpputprop(f,'cofact,x);
      return f;
 end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   calculation of an S-polynomial and related things
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure idqspolynom (p1,p2);
  begin scalar s,ep1,ep2,ep,rp1,rp2,db1,db2,x,cofac1,cofac2;
    if vdpzero!? p1  then return p1; if vdpzero!? p2  then return p2;
    cofac1 := vdpgetprop(p1,'cofact); cofac2 := vdpgetprop(p2,'cofact);
    ep1 := vdpevlmon p1;             ep2 := vdpevlmon p2;
    ep := vevlcm(ep1, ep2);
    rp1 := vdpred p1;               rp2 := vdpred p2;
    db1 := vdplbc p1;                db2 := vdplbc p2;
    if !*vdpinteger then
     <<x:=vbcgcd(db1,db2); db1:=vbcquot(db1,x); db2:=vbcquot(db2,x)>>;
    ep1 := vevdif(ep,ep1); ep2 := vevdif(ep,ep2); db2 := vbcneg db2;
    s := vdpilcomb1 (rp2,db1,ep2,rp1,db2,ep1);
    x := vdpilcomb1 (cofac2,db1,ep2,cofac1,db2,ep1);
    vdpputprop(s,'cofact,x);
    return s;
   end;


symbolic procedure idqdirectelement(p1,p2);
  % the s-Polynomial is reducable to zero because of
  % buchcrit 4. So we can calculate the corresponing cofactor
  % directly
    (if vdpzero!? c1 and vdpzero!? c2 then c1
     else vdpdif(vdpprod(p1,c2),vdpprod(p2,c1))
                        )  where c1 = vdpgetprop(p1,'cofact),
                                 c2 = vdpgetprop(p2,'cofact);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Normailsation with cofactors taken into account
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


symbolic procedure idqsimpcont(p);
    if !*vdpinteger then idqsimpconti p else idqsimpcontr p;

% routines for integer coefficient case:
% calculation of contents and dividing all coefficients by it

symbolic procedure idqsimpconti (p);
%   calculate the contents of p and divide all coefficients by it
  begin scalar res,num,cofac;
    if vdpzero!?  p then return p;
    cofac := vdpgetprop(p,'cofact);
    num := car vdpcontenti p;
    if not vdpzero!? cofac then num:=vbcgcd(num,car vdpcontenti cofac);
    if not vbcplus!? num then num := vbcneg num;
    if not vbcplus!? vdplbc p then num:=vbcneg num;
    if vbcone!? num then return p;
    res := vdpreduceconti (p,num,nil);
    if not vdpzero!? cofac then cofac:=vdpreduceconti(cofac,num,nil);
    res := vdpputprop(res,'cofact,cofac);
    return res;
    end;


% routines for rational coefficient case:
% calculation of contents and dividing all coefficients by it

symbolic procedure idqsimpcontr (p);
%   calculate the contents of p and divide all coefficients by it
  begin scalar res,cofac;
    cofac := vdpgetprop(p,'cofact);
    if vdpzero!?  p then return p;
    if vbcone!? vdplbc p then return p;
    res := vdpreduceconti (p,vdplbc p,nil);
    if not vdpzero!? cofac then
                 cofac:=vdpreduceconti(cofac,vdplbc p,nil);
    res := vdpputprop(res,'cofact,cofac);
    return res;
    end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  trace messages
%

symbolic procedure idqmess0 x;
    if !*trgroeb then
    <<prin2t "adding member to intermediate quotient basis:";
      vdpprint x;
      terpri()>>;

symbolic procedure idqmess1 (x,p1,p2);
    if !*trgroeb then
    <<prin2 "pair ("; prin2 vdpnumber p1; prin2 ",";
      prin2 vdpnumber p2;
      prin2t ") adding member to intermediate quotient basis:";
      vdpprint x;
      terpri()>>;

symbolic procedure idqmess2 b;
    if !*trgroeb then
    <<prin2t "---------------------------------------------------";
      prin2 "the full intermediate base of the ideal quotient is:";
      for each x in b do vdpprin3t x;
      prin2t "---------------------------------------------------";
      terpri()>>;


symbolic procedure idqmess5(p,h);
   if car p then                  % print for true h-Polys
  << hcount!* := hcount!* + 1;
     if !*trgroeb then <<
                       terpri();
                       prin2  "H-polynomial ";
                       prin2 pcount!*;
                       groebmessff(" from pair (",cadr p,nil);
                       groebmessff(",", caddr p,")");
                       vdpprint h;
                       prin2t "with cofactor";
                       vdpprint vdpgetprop(h,'cofact);
                       groetimeprint() >> >>
     else
     if !*trgroeb then <<          % print for input polys
                       prin2t "candidate from input:";
                       vdpprint h;
                       groetimeprint() >>;



symbolic procedure idqmess3(p,s);
     if !*trgroebs then <<
                       prin2 "S-polynomial  ";
                       prin2 " from ";
                       groebpairprint p;
                       vdpprint s;
                       prin2t "with cofactor";
                       vdpprint vdpgetprop(s,'cofact);
                            groetimeprint();
                            terprit 3 >>;

endmodule;


module kredelw; % Kredel Weispfenning algorithm.

%   Author:     H.Melenk (ZIB Berlin)


fluid '(vdpsortmode!* !*groebopt !*gsugar);

symbolic procedure gdimension_eval u;
  begin integer n,m;
   for each s in cdr gindependent_seteval u
    do if (m:=length cdr s)>n then n:=m;
   return n;
  end;

put('gdimension,'psopfn,'gdimension_eval);

symbolic procedure gindependent_seteval pars;
% Independent set algorithm (Kredel/Weispfenning).
% parameters:
%     1      Groebner basis
%     2      optional: list of variables
   begin scalar vars,u,v,w,oldorder,!*factor,!*exp,!*gsugar,
            !*groebopt;
      !*exp := t;
      u := reval car pars;
      v := if cdr pars then reval cadr pars else nil;
      w := for each j in groerevlist u
              collect if eqexpr j then !*eqn2a j else j;
      if null w then rerror(groebnr2,3,"Empty list");
      vars :=  if null v then
                for each j in gvarlis w collect !*a2k j
            else  groerevlist v;
      if not vars then return '(list);
      oldorder := vdpinit vars;
      w:=for each j in w collect vdpevlmon a2vdp j;
      vars:=for each y in vars collect y.vdpevlmon a2vdp y;
      w:=groebkwprec(vars,nil,w,nil);
      return 'list. for each s in w collect
         'list. reversip for each x in s collect car x;
   end;

put('gindependent_sets,'psopfn,'gindependent_seteval);

symbolic procedure groebkwprec(vars,s,lt,m);
  % recursive Kredel Weispfennig algorithm.
  % vars: unprocessed variables,
  % s:    current subset of s,
  % lt:   leading term basis,
  % m:    collection of independent sets so far.
  % returns: updated m.
  begin scalar x,s1,bool;
    s1:=for each y in s collect cdr y;
    while vars do
    <<x:=car vars; vars:= cdr vars;
      if groebkwprec1(cdr x.s1,lt) then
         m:=groebkwprec(vars,x.s,lt,m);
    >>;
    bool:=t;
    for each y in m do % bool and not subsetp(s,y);
        bool:=bool and not(length s=length intersection(s,y));
    return if bool then s.m else m;
  end;

symbolic procedure groebkwprec1(s,lt);
  % t if intersection of T(s) and lt is empty.
   if null lt then t else
   groebkwprec2(s,car lt) and groebkwprec1(s,cdr lt);

symbolic procedure groebkwprec2(s,mon);
  % t if monomial not in T(s).
 <<for each m in s do mon:=vevcan0(m,mon);
   not vevzero!? mon>>;

symbolic procedure vevcan0(m,mon);
  % divide multiples of m1 out of mon.
  if vevzero!? m then mon else
  if vevzero!? mon then nil else
  (if car m neq 0 then 0 else car mon).vevcan0(cdr m,cdr mon);

endmodule;


module traverso;

% Buchberger algorithm base on "sugar" strategy
%  see Giovini-Mora-Niesi-Robbiano-Traverso:
%   One sugar gube, please. ISSAC 91 proceddings, pp 49-54

fluid '(!*gtraverso!-sloppy !*gsugar);

!*gtraverso!-sloppy := t;

symbolic procedure gtraverso (g0,fact,groebres,abort1);
  begin scalar g,d,s,h,p,!*gsugar;

    !*gsugar := t;

    g0:=for each fj in g0 join
          if not vdpzero!? fj then
          <<groebsavelterm fj;
               {gsetsugar(vdpenumerate vdpsimpcont fj,nil)}>>;

  main_loop:
    if null g0 and null d then return gtraverso!-final g;

    if g0 then
         <<h:=car g0;g0:=cdr g0;
           p := list(nil,h,h)
         >>
       else
         <<p := car d;
           d := cdr d;
           s := groebspolynom (cadr p, caddr p);
                  !*trgroeb and groebmess3 (p,s);
           h:=groebsimpcontnormalform groebnormalform(s,g,'list);
           if vdpzero!? h then
           <<!*trgroeb and groebmess4(p,d); goto main_loop>>;

            if vevzero!? vdpevlmon h then % base 1 found
                  <<   !*trgroeb and groebmess5(p,h);
                       d:=g:=g0:=nil;
                  >>;
         >>;

         h := groebenumerate h; !*trgroeb and groebmess5(p,h);
         groebsavelterm h;

          % new pair list
         d := gtraverso!-pairlist(h,g,d);

          % new basis
         g := nconc(g,{h});
         goto main_loop;
   end;

symbolic procedure gtraverso!-pairlist(gk,g,d);
  % gk: new polynomial,
  % g:  current basis,
  % d:  old pair list.
  begin scalar ev,r,n,nn,q;
     % delete triange relations from old pair list.
    d := gtraverso!-pairs!-discard1(gk,d);

     % build new pair list.
    ev := vdpevlmon gk;
    for each p in g do
     if not groebbuchcrit4t(ev,vdpevlmon p)
       then r := vevlcm(ev,vdpevlmon p).r
       else n := groebmakepair(p,gk) . n;

     % delete from new pairs equivalents to coprime lcm.
    for each q in r do
     for each p in n do
      if car p=q then n:=delete(p,n);

     % discard multiples: collect survivers in n
    if !*gtraverso!-sloppy then !*gsugar:=nil;
    n := groebcplistsort(n);
    !*gsugar := t;

    nn := n; n:=nil;
    for each p in nn do
    <<q:=nil;
      for each r in n do
        q:=q or vevdivides!?(car r,car p);
      if not q then n:=groebcplistsortin(p,n);
    >>;

    return groebcplistmerge(d,reversip n);
  end;

symbolic procedure gtraverso!-pairs!-discard1(gk,d);
  % crit B
  begin scalar gi,gj,tij,evk;
   evk:=vdpevlmon gk;
   for each pij in d do
   <<tij := car pij; gi:=cadr pij; gj:=caddr pij;
    if vevstrictlydivides!?(tt(gi,gk),tij)
       and vevstrictlydivides!?(tt(gj,gk),tij)
      then d:=delete(pij,d);
   >>;
   return d;
  end;

symbolic procedure vevstrictlydivides!?(ev1,ev2);
   not(ev1=ev2) and vevdivides!?(ev1,ev2);

symbolic procedure gtraverso!-final g;
  % final reduction and sorting;
  begin scalar r,p,!*gsugar;
   g:=vdplsort g; % descending
   while g do
   <<p:=car g; g:=cdr g;
     if not groebsearchinlist(vdpevlmon p,g) then
       r := groebsimpcontnormalform groebnormalform(p,g,'list) . r;
   >>;
   return list reversip r;
  end;

% copyd('groebbasein,'traverso);

endmodule;


end;


REDUCE Historical
REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]