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

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;
      dummy := nil;
      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;
      dummy := nil;
      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;


end;