reduce-historical/r35/lib/wu.red

module wu;   % Simple implementation of the Wu algorithm.

% Author: Russell Bradford
%         School of Mathematical Sciences
%         University of Bath
%         Bath
%         Avon BA2 7AY
%         United Kingdom
%         E-mail: rjb@maths.bath.ac.uk

% First distributed version: 8 July 90
% Bug fixes in wupseudodivide, and misc other changes: 28 Aug 90

% This is a simple implementation of the Wu algorithm, intended to help
% myself understand the method.  As such, there is little optimization,
% and indeed, only implements the basic version from
%
% "A Zero Structure Theorem for Polynomial-Equations-Solving",
%  Wu Wen-tsun, Institute of Systems Science, Academia Sinica, Beijing

% Interface:
% much as the Groebner basis package:
%
% wu({x*y-a, x^y+y^2-b}, {x, y});
%
% uses Wu on the named polynomials with ordering on the variables x > y.
% returns a list of pairs { characteristic set, initial }
%
% { {{a^2 - b*y^2 + y^4}, y} }
%
% The zeros of the input polynomials are the the union of the zeros of
% the characteristic sets, subject to the initials being non-zero.
% Thus the zeros of {x*y-a, x^y+y^2-b} are the zeros of
% {a^2 - b*y^2 + y^4, a - x*y} subject to y neq 0.
%
% The switch
%
% on trwu;
%
% prints some tracing of the algorithm as it works, in particular the
% choice of basic sets, and the computation of characteristic sets.

% This package runs on Reduce 3.3.

% Keywords: polynomial reduction characteristic set sets initial
%           ascending
% chrstrem Wu

% All improvements and bug fixes are welcomed!!

% Possible bug fixes, improvements:
% Should use distributed polys, then class is an integer;
% rather than use union, use an insertion sort;
% return a list of {{polys},{initials}};
% fix pseudo divide for when there is a non-trivial content in the
%   remainder;
% many opportunities for reusing data from a previous iteration, e.g.,
%   when a new polynomial added into a basic set is less than all
%   current members of the basic set, and they are reduced wrt it.
% factor out monomials and numeric contents

symbolic;

fluid '(!*trwu !*trchrstrem wuvarlist!* kord!*);

switch trwu, trchrstrem;

procedure wuconstantp f;
   % A constant is a poly that does not involve any of the interesting
   % variables.
   domainp f or not memq(mvar f, wuvarlist!*);

smacro procedure wuclass f;
  if wuconstantp f then nil else mvar f;

smacro procedure wudeg f;
  if wuconstantp f then 0 else ldeg f;

smacro procedure wuinitial f;
  if wuconstantp f then f else lc f;

procedure wureducedpolysp(f, polylist);
% if f reduced wrt the polys in polylist?
  null polylist or
  (wureducedp(f, car polylist) and wureducedpolysp(f, cdr polylist));

procedure wureducedp(g, f);
% is g reduced wrt f?
  wuconstantp f or
  wuconstantp g or
  deginvar(g, wuclass f) < ldeg f;

procedure deginvar(f, x);
% the degree of x in f
  if wuconstantp f then 0
  else if mvar f = x then ldeg f
  else begin scalar kord!*;
    kord!* := list x;
    f := reorder f;
    return if mvar f = x then ldeg f else 0
  end;

% wukord* = '(x y a) means: all other symbols < x < y < a
fluid '(wukord!*);

procedure symbollessp(x, y);
% an ordering on symbols: Cambs lisp and PSL orderp differ on nils
  if null y then nil
  else if null x then t
  else if wukord!* then wuorderp(x, y)
  else not orderp(x, y);

procedure wuorderp(x, y);
% an order on the symbols has been specified
% return T if x < y
% circumlocutions abound
begin scalar kord, answ;
   if x eq y then return nil;
   kord := wukord!*;
   while kord and not answ do
     if x eq car kord
       then answ := if memq(y, cdr kord) then 'yes else 'no
     else if y eq car kord
       then answ := if memq(x, cdr kord) then 'no else 'yes
     else kord := cdr kord;
   return if answ then answ eq 'yes else not orderp(x, y)
end;

smacro procedure classlessp(c1, c2);
% an order on classes, which are symbols in this implementation
  symbollessp(c1, c2);

procedure wulessp(f, g);
% standard forms f and g
% a partial order
  classlessp(wuclass f, wuclass g) or
  (wuclass f = wuclass g and wudeg f < wudeg g);

procedure wulessp!*(f, g);
% as above, but use some arbitrary means to complete to a total order
  if wulessp(f, g) then t
  else if wulessp(g, f) then nil
  else totallessp(f, g);

smacro procedure nil2zero f;
  f or 0;

procedure totallessp(f, g);
% a total order on polynomials
  totalcompare(f, g) = 'less;

procedure totalcompare(f, g);
% order f and g
% horrid bit of code
  if f = g then 'equal
  else if wulessp(f, g) then 'less
  else if wulessp(g, f) then 'greater
  else if wuconstantp f then                % and so wuconstantp g
    totalcompareconstants(f, g)
  else begin scalar answ;
    answ := totalcompare(lc f, lc g);
    if answ neq 'equal then return answ;
    return totalcompare(red f, red g)
  end;

procedure totalcompareconstants(f, g);
% order the constants f and g
  if f = g then 'equal
  else if domainp f then
    if domainp g then      % Assumption of ints
      if nil2zero f < nil2zero g then 'less else 'greater
    else 'less
  else if domainp g then 'greater
  else begin scalar wukord!*, wuvarlist!*, answ;
    if symbollessp(mvar f, mvar g) then return 'less
    else if symbollessp(mvar g, mvar f) then return 'greater
    else answ := totalcompareconstants(lc f, lc g);
    if answ neq 'equal then return answ;
    return totalcompareconstants(red f, red g)
  end;

procedure wusort polylist;
% sort a list of polys into Wu order
  sort(polylist, 'wulessp!*); 

procedure collectvars polylist;
% make a list of the variables appearing in the list of polys
begin scalar varlist;
  varlist := for each poly in polylist conc collectpolyvars poly;
  return sort(union(varlist, nil), 'symbollessp)
end;

procedure collectpolyvars poly;
  collectpolyvarsaux(poly, nil);

procedure collectpolyvarsaux(poly, sofar);
  if domainp poly then sofar
  else union(
         union(sofar, list mvar poly),
         union(collectpolyvarsaux(lc poly, nil),
               collectpolyvarsaux(red poly, nil)));

procedure pickbasicset polylist;
% find a basic set from the ordered list of polys
begin scalar basicset;
  foreach var in wuvarlist!* do <<
    while polylist and symbollessp(mvar car polylist, var) do
      polylist := cdr polylist;
    while polylist and var = mvar car polylist and
          not wureducedpolysp(car polylist, basicset) do
      polylist := cdr polylist;
    if polylist and var = mvar car polylist then <<
      basicset := car polylist . basicset;
      polylist := cdr polylist
    >>
  >>;
  return reversip basicset
end;

procedure wupseudodivide(f, g, x);
% not a true pseudo divide---multiply f by the smallest power
% of lc g necessary to make a fraction-free division
begin scalar origf, oldkord, lcoeff, degf, degg, answ, fudge;
  origf := f;
  oldkord := setkorder list x;
  f := reorder f;
  if wuconstantp f or mvar f neq x then <<
    setkorder oldkord;
    return nil . origf
  >>;
  g := reorder g;
  if wuconstantp g or mvar g neq x then <<
    f := multf(f, quotf(g, gcdf!*(lc f, g)));
    setkorder oldkord;
    return reorder f . nil
  >>;
  degf := ldeg f;
  degg := ldeg g;
  if degf - degg + 1 < 0 then <<
    setkorder oldkord;
    return nil . origf
  >>;
  lcoeff := lc g;
  lcoeff := exptf(lcoeff, degf - degg + 1);
  answ := qremf(multf(lcoeff, f), g);
  fudge := gcdf!*(gcdf!*(lcoeff, cdr answ), car answ);
  answ := quotf(car answ, fudge) . quotf(cdr answ, fudge);
  setkorder oldkord;
  return reorder car answ . reorder cdr answ;
end;

procedure simpwupseudodivide u;
begin scalar f, g, x, answ;
  f := !*a2f car u;
  g := !*a2f cadr u;
  x := if cddr u then !*a2k caddr u else mvar f;
  answ := wupseudodivide(f, g, x);
  return list('list, mk!*sq !*f2q car answ,
                     mk!*sq !*f2q cdr answ)
end;

put('wudiv, 'psopfn, 'simpwupseudodivide);

procedure findremainder(f, polylist);
% form the Wu-remainder of f wrt those polys in polylist
<< foreach poly in polylist do
     f := cdr wupseudodivide(f, poly, mvar poly);
   f
>>;

procedure prin2t!* u;
% a useful procedure
<< prin2!* u;
   terpri!* t
>>;

procedure chrstrem polylist;
% polylist a list of polynomials, to be Wu'd
% horrible circumlocutions here
begin scalar revbasicset, pols, rem, remainders;

  if !*trwu or !*trchrstrem then <<
    terpri!* t;
    prin2t!* "--------------------------------------------------------";
  >>;

  repeat <<
    polylist := wusort polylist;

    if !*trwu or !*trchrstrem then <<
      prin2t!* "The new pol-set in ascending order is";
      foreach poly in polylist do printsf poly;
      terpri!* t;
    >>;

    if wuconstantp car polylist then <<
      if !*trwu then prin2t!* "which is trivially trivial";
      remainders := 'inconsistent;
      revbasicset := list 1;
    >>
    else <<
      remainders := nil;
      % Keep in reverse order.
      revbasicset := reversip pickbasicset polylist;
    >>;

    if !*trwu and null remainders then <<
      prin2t!* "A basic set is";
      foreach poly in reverse revbasicset do printsf poly;
      terpri!* t;
    >>;

    pols := setdiff(polylist, revbasicset);
    foreach poly in pols do
      if remainders neq 'inconsistent then <<

        if !*trwu then <<
          prin2!* "The remainder of ";
          printsf poly;
          prin2!* "wrt the basic set is "
        >>;

        rem := findremainder(poly, revbasicset);

        if !*trwu then <<
          printsf rem;
        >>;

        if rem then
          if wuconstantp rem then <<
            remainders := 'inconsistent;

            if !*trwu then <<
              prin2t "which is a non-zero constant, and so";
              prin2t "the equations are inconsistent."
            >>

          >>
          else remainders := union(list absf rem, remainders);
      >>;
    if remainders and remainders neq 'inconsistent then
      polylist := append(polylist, remainders)
  >> until null remainders or remainders = 'inconsistent;

  if remainders = 'inconsistent then revbasicset := list 1;

  if !*trwu or !*trchrstrem then <<
    terpri!* t;terpri!* t;
    prin2t!* "The final characteristic set is:";
    foreach poly in reverse revbasicset do printsf poly
  >>;

  return reversip foreach poly in revbasicset collect absf poly
end;

procedure simpchrstrem u;
begin scalar answ, polylist, wuvarlist!*;
  polylist := foreach f in u collect !*a2f f;
  wuvarlist!* := colectvars polylist;
  answ := chrstrem polylist;
  return 'list . foreach f in answ collect mk!*sq !*f2q f;
end;

put('chrstrem, 'psopfn, 'simpchrstrem);

procedure wu(polylist, varlist);
% Do the Wu algorithm.
% Vars in varlist arranged in increasing order.
% Return (((poly, poly, ... ) . initial) ... ), a list of characteristic
% sets dotted onto the product of their initials.
% Very parallelizable.
  begin scalar stufftodo, answ, polset, chrset, initialset, initial,
        wuvarlist!*;
  stufftodo := list delete(nil,
                       union(foreach poly in polylist collect absf poly,
                             nil));
  if null car stufftodo then <<
    if !*trwu then prin2t!* "trivial CHS";
    return list(list nil . 1);
  >>;
  if null varlist then <<
    if !*trwu then prin2t!* "trivial CHS";
    return list(list 1 . 1);
  >>;
  wuvarlist!* := varlist;
  while stufftodo do <<
     polset := wusort car stufftodo;
     stufftodo := cdr stufftodo;
     chrset := chrstrem polset;
     if chrset neq '(1) then <<
       initialset := foreach pol in chrset collect wuinitial pol;
       initial := 1;
       foreach pol in initialset do initial := multf(initial, pol);

       if !*trwu then <<
         prin2!* "with initial ";
         printsf initial;
       >>;

       if member(initial, chrset) then <<
         if !*trwu then prin2t!*
           "which we discard, as the initial is a member of the CHS";
       >>
       else answ := union(list(chrset . initial), answ);

       foreach initial in initialset do
         if not wuconstantp initial then <<
           if member(initial, polset) then <<
             prin2t!*
            "*** Something awry: the initial is a member of the polset";
             answ := union(list(polset . 1), answ) % unsure of this one.
           >>
           else stufftodo := union(list wusort(initial . polset),
                                   stufftodo)
         >>
     >>
  >>;
  if null answ then answ := list(list 1 . 1);

  if !*trwu then <<
    terpri!* t;terpri!* t;
    prin2t!* "--------------------------------------------------------";
    prin2t!* "Final result:";
    foreach zset in answ do <<
      prin2t!* "Ascending set";
      foreach f in car zset do printsf f;
      prin2!* "with initial ";
      printsf cdr zset;
      terpri!* t
    >>
  >>;

  return answ;
end;

procedure simpwu u;
% rebind kord* to reflect the wu order of kernels
begin scalar pols, vars, oldkord, answ, nargs;
  nargs := length u;
  if nargs = 0 or nargs > 2 then
    rederr "Wu called with wrong number of arguments";

  pols := aeval car u;
  if nargs = 2 then vars := aeval cadr u;
  if (nargs = 1 and not eqcar(pols, 'list)) or
     (nargs = 2 and not eqcar(vars, 'list)) then
    rederr "Wu: syntax wu({poly, ...}) or wu({poly, ...}, {var, ...})";

  oldkord := kord!*;

  if nargs = 1 then
    begin scalar kord!*, polset, vars;
      kord!* := if wukord!* then reverse wukord!* else oldkord;
      polset := foreach f in cdr pols collect reorder !*a2f f;
      vars := collectvars polset;

      if !*trwu then <<
        terpri!* t;
        prin2!* "Wu variables in decreasing order: ";
        foreach id in reverse vars do <<
          prin2!* id;
          prin2!* " "
        >>;
        terpri!* t
      >>;

      answ := wu(polset, vars)
    end
  else     % nargs = 2
    begin scalar kord!*, polset, wukord!*;
      kord!* := foreach k in cdr vars collect !*a2k k;
      wukord!* := reverse kord!*;
      polset := foreach f in cdr pols collect reorder !*a2f f;
      answ := wu(polset, wukord!*)
    end;

    return 'list . foreach zset in answ collect
             'list . list('list . foreach f in car zset collect
                                  mk!*sq !*f2q absf reorder f,
                          mk!*sq !*f2q absf reorder cdr zset)
end;

put('wu, 'psopfn, 'simpwu);
remprop('wu, 'number!-of!-args);

%procedure wukord u;
%% hack to specify order of kernels in Wu
%% wukord a,y,x => other kernels < a < y < x
%  wukord!* := if u = '(nil) then nil
%              else foreach x in u collect !*a2k x;
%
%rlistat '(wukord);

algebraic;

endmodule;

end;