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

module kuechl; % Walking faster, B. Amrhrein, O. Gloor, 
               % W. Kuechlin 
               % in: Calmet, Limongelli (Eds.) Design and 
               % Implementation of Symbolic Computation Systems,
               % Sept.1996

% Version 3 with a rational local solution (after letters from
%  H. M. Moeller).
% Version 4 with keeping the polynomials as DIPs converting only
%  their order mode.

fluid '(vdpsortmode!* vdpsortextension!* dipvars!* !*trgroeb
        groetime!* !*vdpinteger pcount!* !*gsugar !*groebopt
        !*groebrm !*groebdivide);

switch trgroeb;

fluid '(secondvalue!*);

put('groebner_walk,'psopfn,'groeb!-w1);

symbolic procedure groeb!-w1 u;
 begin
  if !*groebopt then
      rerror(groebner,31,
         "don't call 'groebner_walk' with 'on groebopt'");
  if null dipvars!* then
      rerror(groebner,30,"'torder' must be called before");
  groetime!* := time();        % initialization
  !*gsugar := t;               % initialization
  !*groebrm := nil;            % initialization
  u := car groeparams(u,1,1);
  groebnervars(u,nil);
  u := groeb!-list(u,'simp);
  groedomainmode();
  u := groeb!-w2 u;
  return 'list . groeb!-collect(u, 'mk!*sq); 
 end;

symbolic procedure groeb!-list (u, fcn);
% Execute the function 'fcn' for the elements of the algebriac
% list 'u'.
     <<if atom u or not(eqcar(u,'list)) then
         rerror('groebner,29, 
                   "groebner: list as argument required");
         groeb!-collect(cdr u, fcn)
      >>;

symbolic procedure groeb!-collect(l, f);
% Collect the elements of function 'f' applied to the elements of
% the symbolic list 'l'. If 'f' is a number, map 'l' only.
     for each x in l collect
       if numberp f then f else apply1(f,x);

symbolic procedure groeb!-w2 g;
  % This is (essentially) the routine Groebner_Walk.
  % G is a list of standard quotients,
  % a Groebner basis gradlex or based on a vector like [1 1 1 ...]
  % The result is the Groebner basis (standard quotients) with the
  % final term order (lex) as its main order.
  begin scalar iwv, owv, omega, gomega, gomegaplus, tt, tto, pc;
    scalar first, mx, imx, mmx, immx, nn, ll, prim;
    scalar !*vdpinteger, !*groebdivide;
    !*vdpinteger := nil;             % switch on division mode
    !*groebdivide := t;
    first := t;                     % Initialization
    pcount!* := 0;                  % Initialization
    mmx := !*i2rn 1;                % Initialization
    immx := mmx;                    % Initialization
    iwv := groeb!-collect(dipvars!*, 1); 
    omega := iwv;                   % input order vector
    owv := 1 . groeb!-collect(cdr dipvars!*, 0);
    tto := owv;                     % output order vector
    groeb!-w9('weighted,omega);     % Install omega as weighted order
    g := groeb!-collect(g, 'sq2vdp);
    pc := pcount!*;
    gbtest g;                       % Test the Groebner property
    nn := length dipvars!*;
    ll := rninv!: !*i2rn nn;        % inverse of the length
    prim := t;                      % preset
loop:
    groeb!-w9('weighted, omega);
    mx := groeb!-w6!-4 groeb!-collect(omega, 1); 
                              % Compute the maximum of \omega.
        if !*trgroeb then groebmess34 cadr mx;
    imx := rninv!: mx;
    g := if first then groeb!-collect(g, 'vdpsimpcont)
             else groeb!-w10 g;
        if !*trgroeb then groebmess29(omega);
    gomega := if first or not prim then g
         else groeb!-w3(g,omega);      % G_\omega = initials(G_\omega);
    pcount!* := pc;
        if !*trgroeb and not first then groebmess32 gomega;
    gomegaplus := if first then list gomega
        else gtraverso(gomega,nil,nil,nil);
    if cdr gomegaplus then rerror(groebner,31,
                "groebner_walk, cdr of 'groebner' must be nil")
       else gomegaplus := car gomegaplus;
        if !*trgroeb and not first then groebmess30 gomegaplus;
    if not first and prim 
      then g := groeb!-w4(gomegaplus, gomega, g)
      else if not prim then g := gomega;       
      % G = lift(G_{%omega}{plus},<{plus},G_{%omega), G, <)
    if not first
         then g := for each x in g collect gsetsugar (x, nil);
        if !*trgroeb and not first then groebmess31 g;
    if groeb!-w5(omega,imx,tto,immx) then goto ret;
       % stop if tt has been 1 once.
    if not first and rnonep!: tt then goto ret; % Secodary abort crit.
    tt := groeb!-w6!-6(g,tto,immx,omega,imx,ll); % determine_border
        if !*trgroeb then groebmess36 tt;
    if null tt then go to ret;
        % criterion: take primary only if tt neq 1
    prim := not rnonep!: tt;
        if !*trgroeb then groebmess37 prim;
          %\omega = (1 - t)*\omega + t*tau 
    omega := groeb!-w7(tt, omega, imx, tto, immx);
        if !*trgroeb then groebmess35 omega;
    first := nil;
    go to loop;
ret:
        if !*trgroeb then groebmess33 g;
    g := groeb!-collect(g, 'vdpsimpcont);
    g := groeb!-collect(g, 'vdp2sq);
    return g 
end;

symbolic procedure groeb!-w3(g, omega);
% Extract head terms of g corresponding to omega
 begin scalar x, y, gg, ff;
   gg := for each f in g collect
   <<ff := vdpfmon(vdplbc f,vdpevlmon f);
     gsetsugar(ff, nil);
     x := evweightedcomp2(0, vdpevlmon ff, omega);
     y := x;
     f := vdpred f;
     while not vdpzero!? f and y=x do 
     <<y := evweightedcomp2(0, vdpevlmon f,omega);
       if y=x then 
           ff := vdpsum(ff,vdpfmon(vdplbc f,vdpevlmon f)); 
       f := vdpred f>>;
    ff>>;
   return gg;
end; 
        
symbolic procedure groeb!-w4(gb, gomega, g);          
   %gb     Groebner basis of gomega,
   %gomega head term system g_\omega of g,
   %g      full (original) system of polynomials.
begin scalar x;   
   for each y in gb do gsetsugar(y,nil);
   x := for each y in gomega collect groeb!-w8(y, gb);
   x := for each z in x collect groeb!-w4!-1(z, g);
   return x
end;

symbolic procedure groeb!-w4!-1(pl, fs);
% pl is a list of polynomials corresponding to the full system fs.
% Compute the sum of pl * fs.
% Result is the sum.
  begin scalar z;
   z := vdpzero();
   gsetsugar(z,0);
   for each p in pair(pl,fs) do
     if car p then
      z := vdpsum(z,vdpprod(car p , cdr p));
   z := vdpsimpcont z;
   return z
  end;

symbolic procedure groeb!-w5(ev1, x1, ev2, x2);
% ev1 = ev2 equality test.
      groeb!-w5!-1(x1, ev1, x2, ev2);

symbolic procedure groeb!-w5!-1(x1, ev1, x2, ev2);
    (null ev1 and null ev2) or
    (rntimes!:(!*i2rn car ev1 , x1) = rntimes!:(!*i2rn car ev2 , x2)
         and groeb!-w5!-1(x1 ,cdr ev1 ,x2 , cdr ev2));

symbolic procedure groeb!-w6!-4 omega;
% Compute the weighted length of \omega.
    groeb!-w6!-5 (omega, vdpsortextension!*, 0);

symbolic procedure groeb!-w6!-5(omega, v, m);
     if null omega then !*i2rn m
     else if 0 = car omega
           then groeb!-w6!-5(cdr omega, cdr v, m)
     else if 1 = car omega
           then groeb!-w6!-5(cdr omega, cdr v, m #+ car v) 
     else groeb!-w6!-5(cdr omega, cdr v, m #+ car omega #* car v);

symbolic procedure groeb!-w6!-6(gb, tt, ifactt, tp, ifactp, ll);
% Compute the weight border (minimum over all polynomials of gb).
    begin 
      scalar mn, x, zero, one;
        zero := !*i2rn 0;
        one := !*i2rn 1;
        while not null gb do
        << x := groeb!-w6!-7(car gb, tt, ifactt, tp, ifactp,
                            zero, one, ll);
           if null mn or (x and rnminusp!: rndifference!:(x, mn)) 
                 then mn := x;
           gb := cdr gb;>>;
        return mn
    end;

symbolic procedure groeb!-w6!-7(pol, tt, ifactt, tp, ifactp,
                                  zero, one, ll);
% Compute the minimal weight for one polynomial; the iea is,
% that the polynomial has a degree greater than 0.
    begin
       scalar a,b,ev1,ev2,x,y,z,mn;
       ev1 := vdpevlmon pol;
       a := evweightedcomp2(0, ev1, vdpsortextension!*);
       y := groeb!-w6!-8(ev1,tt, ifactt, tp, ifactp,
                        zero, zero, one, ll);
       y := (rnminus!: car y) . (rnminus!: cdr y);
       pol := vdpred pol;
       while not (vdpzero!? pol) do
       <<ev2 := vdpevlmon pol;
         pol := vdpred pol;
         b := evweightedcomp2(0, ev2, vdpsortextension!*);
         if not (a = b) then
         <<x := groeb!-w6!-9(ev2, tt, ifactt, tp, ifactp,
                             car y, cdr y, one, ll, nil);
           if x then 
           <<z := rndifference!:(x, one);
             if  rnminusp!: rndifference!:(zero, x) and
                (rnminusp!: z or rnzerop!: z)
               and
                (null mn or rnminusp!: rndifference!:(x, mn))
            then mn := x>>
         >>
       >>;
       return mn
     end;

symbolic procedure groeb!-w6!-8(ev, tt, ifactt, tp, ifactp,
                                   sum1, sum2, m, dm);
 begin
  scalar x, y, z;
  if ev then <<x := rntimes!:(!*i2rn car ev, m); 
               y := rntimes!:(!*i2rn car tp, ifactp);
               z := rntimes!:(!*i2rn car tt, ifactt)>>;
  return
     if null ev then sum1 . sum2 else
       groeb!-w6!-8(cdr ev, cdr tt, ifactt, cdr tp, ifactp,
           rnplus!:(sum1,rntimes!:(y, x)),
          rnplus!:(sum2, rntimes!:( rndifference!:(z, y),x )),
         rndifference!:(m, dm),
         dm)
  end;

symbolic procedure groeb!-w6!-9(ev, tt, ifactt, tp, ifactp,
                                y1 , y2, m, dm, done);
% Compute the rational solution s: 
% (tp + s * (tt - tp)) * ev1 = (tp + s * (tt - tp)) * evn.
% The sum with ev1 is collected already in y1 and y2
% (with negative sign).
% This routine collects the sum with evn and computes the solution.
 begin
  scalar x, y, z;
  if ev then <<x := rntimes!:(!*i2rn car ev, m);
               y := rntimes!:(!*i2rn car tp, ifactp),
               z := rntimes!:(!*i2rn car tt, ifactt)>>;
  return
     if null ev then
          if null done then nil
              else
                rnquotient!:(rnminus!: y1, y2)
        else
         groeb!-w6!-9( cdr ev, cdr tt, ifactt, cdr tp, ifactp,
              rnplus!:(y1, rntimes!:(y, x)),
              rnplus!:(y2, rntimes!:(rndifference!:(z, y),x )),
              rndifference!:(m, dm),
              dm,
              done or not(car ev = 0))
   end;

symbolic procedure groeb!-w7 (tt, omega, x, tto, y);  
% Compute omega*x * (1-tt) + tto*y *tt.
% tt is a rational number.
% x and y are rational numbers (inverses of the legths of omega/tt).
    begin scalar n,z;
      n := !*i2rn 1;
      omega := for each g in omega collect
        <<  z := rnplus!:(
                  rntimes!:(
                   rntimes!:(!*i2rn g, x),
                   rndifference!:(!*i2rn 1, tt)),
                  rntimes!:(
                   rntimes!:(!*i2rn car tto, y),
                   tt));
            tto := cdr tto;
            n := groeb!-w7!-1(n, rninv!: z);
            z>>;
      omega := for each a in omega collect 
            rnequiv rntimes!:(a , !*i2rn n);
      return omega;
    end;

symbolic procedure groeb!-w7!-1(n, m);
% Compute lcm of n and m. N and m are rational numbers.
% Return the lcm.
% Ignore the denominators of n and m.
     begin
       scalar x, y, z;
       if atom n then x := n else
       <<x := rnprep!: n;
         if not atom x then x := cadr x>>;
       if atom m then y := m else
       <<y := rnprep!: m;
         if not atom y then y := cadr y>>;
       z := lcm( x, y);
       return z
    end;

symbolic procedure groeb!-w8(p, gb);
% Computes the cofactor of p wrt gb.
% Result is a list of cofactors corresponding to g.
% The cofactor 0 is represented as nil.
    begin
      scalar x,y;
        x := groeb!-w8!-1(p,gb);
        p := secondvalue!*;
        while not vdpzero!? p do
        <<y := groeb!-w8!-1(p,gb);
          p := secondvalue!*;
          x := for each pp in pair(x,y) do
            if null car pp then cdr pp
            else if null cdr pp then car pp
            else vdpsum(car pp, cdr pp) >>;
         return x
     end;

symbolic procedure groeb!-w8!-1(p, gb);
% Search in groebner basis gb the polynomial which divides the
% head monomial of the polynomial p. The walk version of
% groebSearchInList 
% Result: the sequence corresponding to g with the monomial
% factor inserted.
  begin
    scalar e, cc, r, done, pp;
    pp := vdpevlmon p;
    cc := vdplbc p;
    r := for each poly in gb collect
        if done then nil else
        if vevdivides!?(vdpevlmon poly,pp) then
            <<done := t;
              e := poly; 
              cc := vbcquot(cc,vdplbc poly);   
              pp := vevdif(pp, vdpevlmon poly); 
              secondvalue!* := 
                vdpsum(
                  vdpprod(gsetsugar(vdpfmon(vbcneg cc, pp),nil),poly),
                  p);
              vdpfmon(cc, pp)>>
        else nil;
    if null e then 
       <<print p;
         print "-----------------";
         print gb;
         rerror(groebner,28,"groeb-w8-1 illegal structure")>>;
   return r
   end;

symbolic procedure groeb!-w9(mode, ext);
% Switch on vdp order mode "mode" with extension "ext".
% Result is the previous extension.
    begin scalar x;
       x := vdpsortextension!*;
       vdpsortextension!* := ext;
       torder2 mode;
       return x
   end;

symbolic procedure groeb!-w10 s;
% Convert the dips in s corresponding to the actual order.
    groeb!-collect(s, 'groeb!-w10!-1);

symbolic procedure groeb!-w10!-1 p;
% Convert the dip p corresponding to the actal order.
  begin scalar x;
     x := vdpfmon(vdplbc p, vdpevlmon p);
     x := gsetsugar(vdpenumerate x, nil);
     p := vdpred p;
     while not vdpzero!? p do
     <<x := vdpsum(vdpfmon(vdplbc p, vdpevlmon p), x);
       p := vdpred p>>;
     return x
  end;

symbolic procedure rninv!: x;
% Return inverse of a (rational) number x: 1/x.
    <<if atom x then x := !*i2rn x;
      x := cdr x;
      if car x < 0 then mkrn(- cdr x, - car x)
            else mkrn(cdr x, car x) >>;

symbolic procedure sq2vdp s;
% Standard quotient to VDP.
    begin
       scalar x,y;
       x := f2vdp numr s;
       gsetsugar(x, nil);
       y := f2vdp denr s;
       gsetsugar(y, 0);
       s := vdpdivmon(x,vdplbc y,vdpevlmon y);
       return s;
    end;

symbolic procedure vdp2sq v;
% Conversion VDP to standard quotient
    begin
     scalar x, y, z, one;
     one := 1 ./ 1;
     x := nil ./ 1;
     while not vdpzero!? v do
     <<y := vdp2f vdpfmon(one, vdpevlmon v) ./ 1;
       z := vdplbc v;
       x := addsq(x, multsq(z,y));
       v := vdpred v>>;
      return x;
    end;

endmodule;

end;