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

module glexconv;  % Newbase-Algorithm:
   % Faugere, Gianni, Lazard, Mora


fluid '(!*factor !*complex !*exp !*gcd !*ezgcd); % REDUCE switch
fluid '(                           % switches from the user interface
        !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
        !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
        !*fullreduction !*groebstat !*groebprot !*gltbasis
        !*groebsubs

        !*vdpinteger !*vdpmodular  % indicating type of algorithm
        !*gsugar                   % sugar mode
        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
        glexdomain!*               % information wrt current domain
        !*gsugar
       );

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

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

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

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

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

fluid '(glexmat!*); % matrix for the indirect lex ordering

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

% parameters;
%      glexconvert(basis,[vars],[maxdeg=n],[newvars={x,y,..}])

symbolic procedure glexconverteval u;
   begin integer n; scalar !*groebfac,!*groebrm,!*factor,!*gsugar,
          v, bas,vars,maxdeg,newvars,!*exp; !*exp := t;
      u := for each p in u collect reval p;
      bas := car u ; u := cdr u;
      while u do
       << v := car u; u := cdr u;
          if eqcar(v,'list) and null vars then vars := v
          else if eqcar(v,'equal) then
           if(v := cdr v)and eqcar(v,'maxdeg) then maxdeg := cadr v
             else if eqcar(v,'newvars) then newvars := cadr v
             else <<prin2 (car v);
                    rerror(groebnr2,4,"Glexconvert, keyword unknown")>>
            else rerror(groebnr2,5,
                       "Glexconvert, too many positional parameters")>>;
       return glexbase1(bas,vars,maxdeg,newvars);
      end;

put('glexconvert,'psopfn,'glexconverteval);

symbolic procedure glexbase1(u,v,maxdeg,nv);
   begin scalar vars,w,nd,oldorder,!*gcd,!*ezgcd,!*gsugar;
         integer pcount!*;
      !*gcd := t;
      w := for each j in groerevlist u
              collect if eqexpr j then !*eqn2a j else j;
      if null w then rerror(groebnr2,6,"Empty list in Groebner");
      vars := groebnervars(w,v);
      !*vdpinteger := !*vdpmodular := nil;
      if not flagp(dmode!*,'field) then !*vdpinteger := t
           else
      if !*modular then !*vdpmodular := t;

      if null vars then vdperr 'groebner;
      oldorder := vdpinit vars;
                  % cancel common denominators
      w := for each j in w collect reorder numr simp j;
                  % optimize varable sequence if desired
      w := for each j in w collect f2vdp j;
      for each p in w do
          nd := nd or not vdpcoeffcientsfromdomain!? p;
      if nd then <<!*vdpmodular:= nil;
                    !*vdpinteger := t;
                    glexdomain!* := 2>>
          else  glexdomain!* := 1;
      if glexdomain!*=1 and not !*vdpmodular then !*ezgcd:=t;
      if null maxdeg then maxdeg := 200;
      if nv then nv := groerevlist nv;
      if null nv then nv := vars else
          for each x in nv do if not member(x,vars) then
           <<rerror(groebnr2,7,list("new variable ",x,
                           " is not a basis variable"))>>;
      u := for each v in nv collect a2vdp v;
      gbtest w;
      w := glexbase2(w,u,maxdeg);

      w := 'list . for each j in w collect prepf j;
      setkorder oldorder;
      gvarslast := 'list . vars;
      return w;
   end;

fluid '(glexeqsys!* glexvars!* glexcount!* glexsub!*);

symbolic procedure glexbase2(oldbase,vars,maxdeg);
   % in contrast to documented algorithm monbase ist a list of
   %  triplets  (mon.cof.vect)
   %   such that cof * mon == vect modulo oldbase
   %       (cof is needed because of the integer algoritm)
   begin scalar lexbase, staircase, monbase;
         scalar monom, listofnexts, vect,q,glexeqsys!*,glexvars!*,
                glexsub!*;
         integer n;
      if not groezerodim!?(oldbase,length vars) then
        prin2t "####### warning: ideal is not zerodimensional ######";
             % prepare matrix for the indirect lex ordering
      glexmat!* := for each u in vars collect vdpevlmon u;
      monbase := staircase := lexbase := nil;
      monom := a2vdp 1;
      listofnexts := nil;
      while not(monom = nil) do
      <<
        if not glexmultipletest(monom,staircase) then
        << vect := glexnormalform(monom,oldbase);
           q := glexlinrel(monom,vect,monbase);
           if q then
           <<lexbase := q . lexbase; maxdeg := nil;
             staircase := monom . staircase;
           >>
           else
           <<monbase := glexaddtomonbase(monom,vect,monbase);
             n := n #+1;
             if maxdeg and n#> maxdeg then
                  rerror(groebnr2,8,
                         "No univar. polynomial within degree bound");
             listofnexts :=
                 glexinsernexts(monom,listofnexts,vars)>>;
        >>;
        if null listofnexts then monom := nil
           else << monom := car listofnexts ;
                   listofnexts := cdr listofnexts >>;
      >>;
      return lexbase;
    end;

symbolic procedure glexinsernexts(monom,l,vars);
   begin scalar x;
     for each v in vars do
     << x := vdpprod(monom,v);
        if not vdpmember(x,l) then
        <<vdpputprop(x,'factor,monom); vdpputprop(x,'monfac,v);
          l := glexinsernexts1(x,l);
        >>;
     >>;
     return l;
   end;

symbolic procedure glexmultipletest(monom,staircase);
    if null staircase then nil
    else if vevmtest!?(vdpevlmon monom, vdpevlmon car staircase)
          then t
    else glexmultipletest(monom, cdr staircase);

symbolic procedure glexinsernexts1(m,l);
    if null l then list m
    else if glexcomp(vdpevlmon m,vdpevlmon car l) then m.l
    else car l . glexinsernexts1(m,cdr l);

symbolic procedure glexcomp(ev1,ev2);
  % true if ev1 is greater than ev2
  % we use an indirect ordering here (mapping via newbase variables)
     glexcomp0(glexcompmap(ev1,glexmat!*),
               glexcompmap(ev2,glexmat!*));

symbolic procedure glexcomp0(ev1,ev2);
     if null ev1 then nil
     else if null ev2 then glexcomp0(ev1,'(0))
     else if (car ev1 #- car ev2) = 0
                  then glexcomp0(cdr ev1,cdr ev2)
     else if car ev1 #< car ev2 then t
     else nil;

symbolic procedure glexcompmap (ev,ma);
      if null ma then nil
      else glexcompmap1(ev,car ma) . glexcompmap(ev,cdr ma);

symbolic procedure glexcompmap1(ev1,ev2);
   % the dot product of two vectors
      if null ev1 or null ev2 then 0
      else (car ev1 #* car ev2) #+ glexcompmap1(cdr ev1,cdr ev2);

symbolic procedure glexaddtomonbase(monom,vect,monbase);
  % primary effect: (monom . vect) . monbase
  % secondary effect: builds the equation system
   begin scalar x;
       if null glexeqsys!* then
        <<glexeqsys!* := a2vdp 0; glexcount!*:=-1>>;
       x := mkid('gunivar,glexcount!*:=glexcount!*+1);
       glexeqsys!* := vdpsum(glexeqsys!*,vdpprod(a2vdp x,cdr vect));
       glexsub!* := (x .(monom . vect)) . glexsub!*;
       glexvars!* := x . glexvars!*;
       return (monom . vect) . monbase;
   end;

symbolic procedure glexlinrelold(monom,vect,monbase);
    if monbase then
    begin scalar sys,sub,auxvars,r,v,x;
       integer n;
         v := cdr vect;
         for each b in reverse monbase do
         <<x := mkid('gunivar,n); n := n+1;
           v := vdpsum(v,vdpprod(a2vdp x,cddr b));
           sub := (x . b) . sub;
           auxvars := x . auxvars>>;
       while not vdpzero!? v do
       <<sys := vdp2f vdpfmon(vdplbc v,nil) . sys; v := vdpred v>>;
       x := sys;
       sys := groelinsolve(sys,auxvars);
       if null sys then return nil;
          % construct the lex polynomial
       if !*trgroeb
         then prin2t " ======= constructing new basis polynomial";
       r := vdp2f vdpprod(monom,car vect) ./ 1;
       for each s in sub do
         r:= addsq(r,multsq(vdp2f vdpprod(cadr s,caddr s) ./ 1,
                            cdr assoc(car s,sys)));
       r := vdp2f vdpsimpcont f2vdp numr r;
       return r;
    end;

symbolic procedure glexlinrel(monom,vect,monbase);
    if monbase then
    begin scalar sys,r,v,x;
       v := vdpsum(cdr vect,glexeqsys!*);
       while not vdpzero!? v do
       <<sys := vdp2f vdpfmon(vdplbc v,nil) . sys; v := vdpred v>>;
       x := sys;
       sys := groelinsolve(sys,glexvars!*);
       if null sys then return nil;
          % construct the lex polynomial
       r := vdp2f vdpprod(monom,car vect) ./ 1;
       for each s in glexsub!* do
         r:= addsq(r,multsq(vdp2f vdpprod(cadr s,caddr s) ./ 1,
                            cdr assoc(car s,sys)));
       r := vdp2f vdpsimpcont f2vdp numr r;
       return r;
    end;



symbolic procedure glexnormalform(m,g);
  % reduce m wrt basis G;
  % the reduction product is preserved in m for later usage
  begin scalar cof,vect,r,f,fac1;
      if !*trgroeb then prin2t " ======= reducing ";
      fac1 := vdpgetprop(m,'factor);
      if fac1 then vect := vdpgetprop(fac1,'vector);
      if vect then
      << f := vdpprod(cdr vect,vdpgetprop(m,'monfac));
	 cof := car vect>>
      else
      << f := m; cof :=  a2vdp 1>>;
      r := glexnormalform1(f,g,cof);
      vdpputprop(m,'vector,r);
      if !*trgroeb then
      <<vdpprint vdpprod(car r,m); prin2t " =====> ";
        vdpprint cdr r>>;
      return r;
   end;


symbolic procedure glexnormalform1(f,g,cof);
  begin scalar f1,c,vev,divisor,done,fold,a,b;
      fold := f; f1 := vdpzero(); a:= a2vdp 1;
      while not vdpzero!? f do
          begin
            vev:=vdpevlmon f; c:=vdplbc f;
            divisor := groebsearchinlist (vev,g);
            if divisor then done := t;
            if divisor then
                 if !*vdpinteger then
                  << f:=groebreduceonestepint(f,a,c,vev,divisor);
                     b := secondvalue!*;
                     cof := vdpprod(b,cof);
                     if not vdpzero!? f1 then
                          f1 := vdpprod(b,f1);
                  >>
                  else
                   f := groebreduceonesteprat(f,nil,c,vev,divisor)
             else
                  <<f1 := vdpappendmon (f1,vdplbc f,vdpevlmon f);
                    f := vdpred f;
                    >>;

          end;
      if not done then return cof . fold;
      f := groebsimpcont2(f1,cof); cof := secondvalue!*;
      return cof . f;
   end;


symbolic procedure groelinsolve(equations,xvars);
    (begin scalar r,q,test,oldmod,oldmodulus;
       if !*trgroeb then prin2t " ======= testing linear dependency ";
       r := t;
       if not !*modular and glexdomain!*=1 then
       <<oldmod := dmode!*;
         if oldmod then setdmode(get(oldmod,'dname),nil);
         oldmodulus := current!-modulus;
         setmod list 16381;  % = 2**14-3
         setdmode('modular,t);
         r := groelinsolve1(for each u in equations collect
                                numr simp prepf u, xvars);
         setdmode('modular,nil);
         setmod list oldmodulus;
         if oldmod then setdmode(get(oldmod,'dname),t);
       >> where !*ezgcd=nil;
       if null r then return nil;
       r := groelinsolve1(equations,xvars);
       if null r then return nil;
         % divide out the common content
       for each s in r do
         if not(denr cdr s = 1) then test := t;
       if test then return r;
       q := numr cdr car r;
   %   for each s in cdr r do
   %     if q neq 1 then
   %        q := gcdf!*(q,numr cdr s);
   %   if q=1 then return r;
   %   r := for each s in r collect
   %      car s . (quotf(numr cdr s, q) ./ 1);
       return r;
     end) where !*ezgcd=!*ezgcd;    % stack old value

symbolic procedure groelinsolve1(equations,xvars);
  % Gaussian elimination in integer mode
  %     free of unexact divisions (see Davenport et al,CA, pp86-87
  % special cases: trivial equations are ruled out early
  % INPUT:
  % equations:     list of standard forms
  % xvars:         variables for the solution
  % OUTPUT:
  % list of pairs (var . solu) where solu is a standard quotient
  %
  % internal data structure: standard forms as polynomials in xvars
   begin scalar oldorder,x,p,solutions,val,later,break,gc,field;
    oldorder := setkorder xvars;
    field := dmode!* and flagp(dmode!*,'field);
    equations := for each eqa in equations collect reorder eqa;
    for each eqa in equations do
         if eqa and domainp eqa then break:= t;
    if break then goto empty;
    equations := sort(equations,function grloelinord);

  again:
    break := nil;
    for each eqa in equations do if not break then
     % car step: eliminate equations of type 23 = 0
     %                                    and 17 * u = 0
     %                                    and 17 * u + 22 = 0;
     << if null eqa then equations := delete(eqa,equations)
       else if domainp eqa then break := t  % inconsistent system
       else if not member(mvar eqa,xvars) then break := t
       else if domainp red eqa or not member(mvar red eqa,xvars) then
       <<equations := delete (eqa,equations);
         x := mvar eqa;
         val := if lc eqa = 1 then negf red eqa ./ 1
                      else multsq(negf red eqa ./ 1, 1 ./lc eqa);
         solutions := (x . val) . solutions;
         equations := for each q in equations collect
                               groelinsub(q,list(x.val));
         later :=
              for each q in later collect groelinsub(q,list(x.val));
         break := 0;
        >>;
     >>;
     if break = 0 then goto again else if break then goto empty;

     % perform an elimination loop
     if null equations then goto ready;
     equations := sort(equations,function grloelinord);
     p := car equations; x := mvar p;
     equations := for each eqa in cdr equations collect
       if mvar eqa = x then
          << if field then
            eqa := addf(eqa, negf multf(quotf(lc eqa,lc p),p))
            else
            <<gc := gcdf(lc p,lc eqa);
              eqa := addf(multf(quotf(lc p,gc),eqa),
                        negf multf(quotf(lc eqa,gc),p))>>;
            if not domainp eqa then
                    eqa := numr multsq(eqa ./ 1, 1 ./ lc eqa);
            %%%%%%eqa := groelinscont(eqa,xvars);
            eqa>>
                  else eqa;
     later := p . later;
     goto again;


ready:   % do backsubstitutions
     while later do
     <<p := car later; later := cdr later;
       p := groelinsub(p,solutions);
       if domainp p or not member(mvar p,xvars) or
            (not domainp red p and member(mvar red p,xvars)) then
                               <<break := t; later := nil>>;
       x := mvar p;
       val := if lc p = 1 then negf red p ./ 1
             else quotsq(negf red p ./ 1 , lc p ./ 1);
       solutions := (x . val) . solutions;
     >>;
     if break then goto empty else goto finis;
 empty: solutions := nil;
 finis: setkorder oldorder;
     solutions := for each s in solutions collect
            car s . (reorder numr cdr s ./ reorder denr cdr s);
      return solutions;
  end;


symbolic procedure grloelinord(u,v);
    % apply ordop to the mainvars of u and v
         ordop(mvar u, mvar v);

symbolic procedure groelinscont(f,vars);
   % reduce content from standard form f
  if domainp f then f else
  begin scalar c;
    c := groelinscont1(lc f,red f,vars);
    if c=1 then return f;
    prin2 "*************content: ";print c;
    return quotf(f,c);
  end;

symbolic procedure groelinscont1(q,f,vars);
  % calculate the contents of standard form f
   if null f or q=1 then q
   else if domainp f or not member(mvar f,vars) then gcdf!*(q,f)
   else groelinscont1(gcdf!*(q,lc f),red f,vars);

symbolic procedure groelinsub(s,a);
  % s is a standard form linear in the top level variables
  % a is an assiciation list (variable . sq) . ...
  % The value is the standard form, where all substitutions
  % from a are done in s (common denominator ignored).
         numr groelinsub1(s,a);

symbolic procedure groelinsub1(s,a);
    if domainp s then s  ./ 1
    else (if x then addsq(multsq(cdr x,lc s ./ 1),y)
               else addsq(lt s .+ nil ./ 1,y))
         where x=assoc(mvar s,a),y=groelinsub1(red s,a);

endmodule;

end;