reduce-historical/r38/packages/groebner/groebtra.red

module groebtra;
 
% Calculation of a Groebner base with the Buchberger algorithm
% including the backtracking information which denotes the
% dependency between base and input polynomials .
 
% Authors: H. Melenk, H.M. Moeller, W. Neun;date : August 2000 
 
switch groebopt,groebfac,trgroeb,trgroebs,trgroeb1,
       trgroebr,groebstat,groebprot;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interface
 
symbolic procedure groebnertraeval u;
% Backtracking Groebner calculation .
begin integer n;scalar !*groebfac,!*groebrm,!*groebprot,!*gsugar;
 n:=length u;if n=1 then return groebnertra1(reval car u,nil,nil)
  else if n neq 2 then rerror(groebnr2,10,
		    "groebnert called with wrong number of arguments")
  else return groebnertra1(reval car u,reval cadr u,nil)end;

put('groebnert,'psopfn,'groebnertraeval);

symbolic procedure groebnertra1(u,v,mod1);
% 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,y,z,x,np,oldorder,groetags!*,tagvars;
 integer pcount!*,nmod;
 u:=for each j in getrlist u collect 
 <<if not eqcar(j,'equal)
       or not(idp(w:=cadr j)or(pairp w and w=reval w and
                          get(car w,'simpfn)='simpiden))
        then rerror(groebnr2,11,
	   "groebnert parameter not {...,name=polynomial,...}")
        else cadr j . caddr j>>;
    if null u then rerror(groebnr2,12,"empty list in groebner")
     else if null cdr u then return'list . {'equal,cdar u,caar u};
    w:=for each x in u collect cdr x;
    if mod1 then
    <<z:=nil;
       for each vect in w do 
       <<if not eqcar(vect,'list)then
           rerror(groebnr2,13,"non list given to groebner");
          if nmod=0 then nmod:=length cdr vect else
           if nmod <(x:=length cdr vect)then nmod=x;
          z:=append(cdr vect,z)>>;
       if not eqcar(mod1,'list) then
        rerror(groebnr2,14,"illegal column weights specified");
       vars:=groebnervars(z,v);
       tagvars:=for i:=1 : nmod collect mkid('! col,i);
       w:=for each vect in w collect
       <<z:=tagvars;y:=cdr mod1;
       'plus . for each p in cdr vect collect
        'times . {'expt,car z,car y} . 
           <<z:=cdr z;y:=cdr y;
              if null y then y:='(1);{p}>>>>;
       groetags!*:=length vars + 1;
       vars:=append(vars,tagvars)>>
      else vars:=groebnervars(w,v);
      groedomainmode();
      if null vars then vdperr'groebner;
      oldorder:=vdpinit vars;
                  % Cancel common denominators .
      u:=pair(for each x in u collect car x,w);
      u:=for each x in u collect
      <<z:=simp cdr x;
         multsq(simp car x,denr z ./ 1). reorder numr z>>;
      w:=for each j in u collect cdr j;
                  % Optimize varable sequence if desired .
      if !*groebopt then <<w:=vdpvordopt(w,vars);vars:=cdr w;
                           w:=car w;vdpinit vars>>;
      w:=pair(for each x in u collect car x,w);
      w:=for each j in w collect
      <<z:= f2vdp cdr j;vdpputprop(z,'cofact,car 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>>>>;
      w:=groebtra2 w;
      w:=if mod1 then groebnermodres(w,nmod,tagvars)else
                       groebnertrares w;
      setkorder oldorder;
      gvarslast:='list . vars;return w end;
 
symbolic procedure groebnertrares w;
begin scalar c,u;
 return'list . for each j in w collect
 <<c:=prepsq vdpgetprop(j,'cofact);
    u:=vdp2a j;if c=0 then u else {'equal,u,c}>> end;

symbolic procedure groebnermodres(w,n,tagvars);
begin scalar x,c,oldorder;
 c:=for each u in w collect prepsq vdpgetprop(u,'cofact);
 oldorder:=setkorder tagvars;
 w:=for each u in w collect 
'list . <<u:=numr simp vdp2a u;
          for i:=1 : n collect
           prepf if not domainp u and mvar u=nth(tagvars,i)
             then <<x:=lc u;u:=red u;x>> else nil>>;
 setkorder oldorder;
         % Reestablish term order for output .
 w:=for each u in w collect vdp2a a2vdp u;
 w:=pair(w,c);
 return'list . for each p in w collect 
  if cdr p=0 then car p else
   {'equal,car p,cdr p} end;

symbolic procedure preduceteval pars;
% Trace version of PREDUCE;
% parameters:
%     1      expression to be reduced,
%        (formula or equation)
%     2      polynomials or equations;base for reduction;
%            must be equations with atomic lhs;
%     3      optional: list of variables .
begin scalar vars,x,y,u,v,w,z,oldorder,!*factor,!*exp,
 !*gsugar;integer pcount!*;!*exp:=t;
 pars:=groeparams(pars,2,3);
 y:=car pars;u:=cadr pars;v:= caddr pars;
 u:=for each j in getrlist u collect
 <<if not eqcar(j,'equal)then
            j . j else cadr j . caddr j>>;
 if null u then rerror(groebnr2,15,"empty list in preducet");
 w:=for each p in u collect cdr p;% The polynomials .
 groedomainmode();
 vars:=if null v then
  for each j in gvarlis w collect !*a2k j else  getrlist v;
 if not vars then vdperr'preducet;
 oldorder:=vdpinit vars;
 u:=for each x in u collect
 <<z:=simp cdr x;
    multsq(simp car x,denr z ./ 1). reorder numr z>>;
 w:=for each j in u collect
 <<z:= f2vdp cdr j;vdpputprop(z,'cofact,car j)>>;
 if not eqcar(y,'equal)then y:={'equal,y,y};
 x:=a2vdp caddr y;          % The expression .
 vdpputprop(x,'cofact,simp cadr y);% The lhs(name etc.) .
 w:=tranormalform(x,w,'sort,'f);
 u:={'equal,vdp2a w,prepsq vdpgetprop(w,'cofact)};
 setkorder oldorder;return u end;

put('preducet,'psopfn,'preduceteval);

symbolic procedure groebnermodeval u;
% Groebner for moduli calculation .
( if n=0 or n > 3 then 
   rerror(groebnr2,16,
    "groebnerm called with wrong number of arguments")
  else groebnertra1(reval car u,
        if n >= 2 then reval cadr u else nil,
        if n >= 3 then reval caddr u else'(list 1))
  )where n=length u;

put('groebnerm,'psopfn,'groebnermodeval);
 
symbolic procedure groebtra2 p;
% Setup all global variables for the Buchberger algorithm;
% printing of statistics .
begin scalar groetime!*,tim1,spac,spac1,p1,
 pairsdone!*,factorlvevel!*;integer factortime!*;
 groetime!*:=time();
 vdponepol();% We construct dynamically .
 hcount!*:=pcount!*:=mcount!*:=fcount!*:=
 bcount!*:=b4count!*:=hzerocount!*:=basecount!*:=0;
 if !*trgroeb then
 <<prin2 "Groebner calculation with backtracking starts ";terprit 2>>;
 spac:=gctime();p1:= groebtra3 p;
 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  ";
    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;
 
symbolic procedure groebtra3 g0;
begin scalar x,g,d,d1,d2,p,p1,s,h,g99,one;
 x:=for each fj in g0 collect vdpenumerate trasimpcont fj;
 for each fj in x do g:=vdplsortin(fj,g0);
 g0:=g;g:=nil;
% 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:={nil,h,h}>>
 else
 <<       % Take next poly from pairs .
    p:=car d;d:=cdr d;
    s:=traspolynom(cadr p, caddr p);tramess3(p,s);
    h:=groebnormalform(s,g99,'tree);% Piloting wo cofact .
    if vdpzero!? h then groebmess4(p,d)
     else h:=trasimpcont tranormalform(s,g99,'tree,'h)>>;
 if vdpzero!? h then goto bott;
 if vevzero!? vdpevlmon h then % Base 1 found .
 <<       tramess5(p,h);
   g0:=d:=nil;g:={h};goto bott>>;
 s:= nil;
                  % h polynomial is accepted now .
 h:=vdpenumerate h;
                        tramess5(p,h);
                              % Construct new critical pairs .
 d1:=nil;
 for each f in g do
  if groebmoducrit(f,h)then
  <<d1:=groebcplistsortin( {tt(f,h),f,h},d1);
     if tt(f,h)=vdpevlmon f then
     <<g:=delete(f,g);groebmess2 f>>>>;
                  groebmess51 d1;
 d2:=nil;
 while d1 do
 <<d1:=groebinvokecritf d1;
    p1:=car d1;
    d1:=cdr d1;
    if groebbuchcrit4(cadr p1,caddr p1,car p1)
     then d2:=append(d2,{p1});
    d1:=groebinvokecritm(p1,d1)>>;
 d:=groebinvokecritb(h,d);
 d:=groebcplistmerge(d,d2);
 g:=h . g;
 g99:=groebstreeadd(h,g99);
                            groebmess8(g,d);
bott: end;
 return groebtra3post g end;
 
symbolic procedure groebtra3post g;
% Final reduction .
begin scalar r,p;
 g:=vdplsort g;
 while g do
 <<p:=tranormalform(car g,cdr g,'sort,'f);
    if not vdpzero!? p then r:=trasimpcont p . r;g:=cdr g>>;
 return reversip r end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Reduction of polynomials .
%
 
symbolic procedure tranormalform(f,g,type,mode);
% 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,break;
 while not vdpzero!? f and not break do
 begin vev:=vdpevlmon f;c:=vdplbc f;
  divisor:=groebsearchinlist(vev,g);
  if divisor and !*trgroebs then
    <<prin2 "//-";prin2 vdpnumber divisor>>;
  if divisor then if !*vdpinteger  then
   f:=trareduceonestepint(f,nil,c,vev,divisor)
  else f:=trareduceonesteprat(f,nil,c,vev,divisor)
  else break:=t end;
 if mode='f then f:=tranormalform1(f,g,type,mode);
 return f end;

symbolic procedure tranormalform1(f,g,type,mode);
% Reduction of subsequent terms .
begin scalar c,vev,divisor,break,f1;
 mode:=nil;f1:=f;type:=nil;
 while not vdpzero!? f and not vdpzero!? f1 do
 <<f1:=f;break:=nil;
    while not vdpzero!? f1 and not break do
    <<vev:=vdpevlmon f1;c:=vdplbc f1;f1:=vdpred f1;
       divisor:=groebsearchinlist(vev,g);
       if divisor and !*trgroebs then
       <<prin2 "//-";prin2 vdpnumber divisor>>;
       if divisor then
       <<if !*vdpInteger  then
            f:=trareduceonestepint(f,nil,c,vev,divisor)
           else
            f:=trareduceonesteprat(f,nil,c,vev,divisor);
           break:=t>>>>>>;
 return f end;
 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
%  special reduction procedures
 
symbolic procedure trareduceonestepint(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;
 dummy:=nil;fcofa:=vdpgetprop(f,'cofact);
 if null fcofa then fcofa:=nil ./ 1;
 gcofa:=vdpgetprop(g1,'cofact);
 if null gcofa then gcofa:=nil ./ 1;
 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(f,a,vevzero(),g1,vbcneg b,vevlcm);
 x:=vdpilcomb1tra(fcofa,a,vevzero(),gcofa,vbcneg b,vevlcm);
 vdpputprop(f,'cofact,x);return f end;
 
symbolic procedure trareduceonesteprat(f,dummy,c,vev,g1);
% Reduction step for rational case :
% calculate f=f - g / vdpLbc(f)
begin scalar x,fcofa,gcofa,vev;
 dummy:=nil;fcofa:=vdpgetprop(f,'cofact);
 gcofa:=vdpgetprop(g1,'cofact);
 vev:=vevdif(vev,vdpevlmon g1);
 x:=vbcneg vbcquot(c,vdplbc g1);
 f:=vdpilcomb1(f,a2vbc 1,vevzero(),g1,x,vev);
 x:=vdpilcomb1tra(fcofa,a2vbc 1,vevzero(),gcofa,x,vev);
 vdpputprop(f,'cofact,x);return f end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Calculation of an S-polynomial .
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
symbolic procedure traspolynom(p1,p2);
begin scalar s,ep1,ep2,ep,rp1,rp2,db1,db2,x,
 cofac1,cofac2;
 if vdpzero!? p1  then return p1;if vdpzero!? p1  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:=vdpilcomb1tra(cofac2,db1,ep2,cofac1,db2,ep1);
 vdpputprop(s,'cofact,x);return s end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Normalisation with cofactors taken into accounta .
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
symbolic procedure trasimpcont p;
 if !*vdpinteger then trasimpconti p else trasimpcontr p;
 
% Routines for integer coefficient case :
% calculation of contents and dividing all coefficients by it .
 
symbolic procedure trasimpconti 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 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);
 cofac:=vdpreducecontitra(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 trasimpcontr p;
% Calculate the contents of p and divide all coefficients by it .
begin scalar res,cofac;
 cofac:=vdpgetprop(p,'cofact);
 if vdpzero!? p or vdplbc p then return p;
 res:=vdpreduceconti(p,vdplbc p,nil);
 cofac:=vdpreducecontitra(cofac,vdplbc p,nil);
 res:=vdpputprop(res,'cofact,cofac);return res end;
 
symbolic procedure vdpilcomb1tra(cofac1,db1,ep1,cofac2,db2,ep2);
% The linear combination, here done for the cofactors(standard quotients);
 addsq(multsq(cofac1,vdp2f vdpfmon(db1,ep1)./ 1),
         multsq(cofac2,vdp2f vdpfmon(db2,ep2)./ 1));
 
symbolic procedure vdpreducecontitra(cofac,num,dummy);
% Divide the cofactor by a number .
<<dummy:=nil;quotsq(cofac,simp vbc2a num)>>;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Special handling of moduli .
%

symbolic procedure groebmoducrit(p1,p2);
 null groetags!*
  or pnth(vdpevlmon p1,groetags!*)=pnth(vdpevlmon p2,groetags!*);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Trace messages .
%
 
symbolic procedure tramess0 x;
 if !*trgroeb then
  <<prin2t "adding member to intermediate quotient base:";vdpprint x;terpri()>>;
 
symbolic procedure tramess1(x,p1,p2);
 if !*trgroeb then
 <<prin2 "pair(";prin2 vdpnumber p1;prin2 ",";prin2 vdpnumber p2;
    prin2t ") adding member to intermediate quotient base:";vdpprint x;
    terpri()>>;
 
symbolic procedure tramess5(p,pp);
 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 pp;prin2t "with cofactor";
        writepri(mkquote prepsq vdpgetprop(pp,'cofact),'only);
        groetimeprint()>>>>
  else if !*trgroeb then             % print for input polys
  <<prin2t "candidate from input:";vdpprint pp;groetimeprint()>>;
 
symbolic procedure tramess3(p,s);
 if !*trgroebs then <<
   prin2 "S-polynomial from ";groebpairprint p;vdpprint s;prin2t "with cofactor";
   writepri(mkquote prepsq  vdpgetprop(s,'cofact),'only);
   groetimeprint();terprit 3>>;
 
endmodule;;end;