File r37/packages/cgb/gb.red artifact bc0910d1d1 part of check-in f2fda60abd


% ----------------------------------------------------------------------
% $Id: gb.red,v 1.31 1999/04/15 07:08:01 dolzmann Exp $
% ----------------------------------------------------------------------
% Copyright (c) 1999 Andreas Dolzmann and Thomas Sturm
% ----------------------------------------------------------------------
% $Log: gb.red,v $
% Revision 1.31  1999/04/15 07:08:01  dolzmann
% Do not load rltools during compilation.
%
% Revision 1.30  1999/04/13 20:57:04  dolzmann
% Renamed switches to cgb...
% Removed !*gsugar.
% Sort the input system.
%
% Revision 1.29  1999/04/13 13:51:07  dolzmann
% gb_gb was called with three arguments instead of one argument.
%
% Revision 1.28  1999/04/11 11:30:55  dolzmann
% Added standard form interface gb_gbf. Added a procedure gb_gsys!$ for
% computing non-parametric Groebner systems.
%
% Revision 1.27  1999/04/11 09:50:40  dolzmann
% Completely rewritten the interface code for the AM.
% Moved the module vdp from dp.red into this file.
% Adappted the code to the dip_init procedure.
%
% Revision 1.26  1999/04/06 11:53:57  dolzmann
% Removed switches trgroeb, trgroebr, trgreobs, and related code.
%
% Revision 1.25  1999/04/04 16:46:19  sturm
% Added copyright and CVS fluids.
%
% Revision 1.24  1999/04/04 14:50:38  sturm
% Implemented switch tdusetorder.
%
% Revision 1.23  1999/04/04 14:09:33  sturm
% Moved dip_ilcomb and dip_ilcombr from cgb.red to dp.red.
% Created vdp_ilcomb and vdp_ilcombr for gb.red.
%
% Revision 1.22  1999/04/04 12:23:16  dolzmann
% Procedure gb_spolynomial expexts a critical pair instead of two polynomials.
%
% Revision 1.21  1999/04/03 13:38:37  sturm
% Adapted to new dip_init/dip_cleanup.
% Fixed bug in gb_strange!-reduction.
%
% Revision 1.20  1999/03/31 14:09:56  sturm
% Fixed numerous details encountered during CGB reimplementation.
%
% Revision 1.19  1999/03/30 11:29:00  dolzmann
% gb_a2s returns now a list of SF's.
% Reimplemented procedure gb_vars.
%
% Revision 1.18  1999/03/30 09:36:43  dolzmann
% Removed unused procedure gb_vdpvordopt3.
%
% Revision 1.17  1999/03/30 09:34:44  dolzmann
% Procedure gb_groebner1 binds all locally used fluids.
%
% Revision 1.16  1999/03/17 12:34:06  dolzmann
% Removed fluids !*vdpinteger, !*vdpmodular, !*grmod  and their bindings.
% Use vdp_monp instead of vdp_length.
% Added the lost procedure min!# under the new name gb_min!#.
%
% Revision 1.15  1999/03/12 13:29:41  sturm
% Added a procedure gb_updbase. Currently used in the style of groebner.red
% incontrast to the Asir variant. Switch gbupdb.
%
% Revision 1.14  1999/03/12 10:58:47  dolzmann
% Introduced verbose output for the length of g.
% Use procedure ev_sdivp of package dp instead of gb_vevsdivp.
% Removed test for zero polynomials in gb_spolynomial.
% Moved procedure vdp_ilcomb1 and vdp_ilcomb1r into package dp and
% renamed the procedures accordingly.
%
% Revision 1.13  1999/03/12 08:57:59  dolzmann
% Introduced switch gbcheckg and related code.
% All statistics code is guarded by the respective switches.
% Replaced all calls of vdp_putprop and vdp_getprop by access functions
% for the properties of the dp package.
%
% Revision 1.12  1999/03/05 10:35:41  dolzmann
% Adapted to newly created dp package.
%
% Revision 1.11  1999/03/02 15:52:11  dolzmann
% Added switch gbcounthf for counting reducible H-polynomials with groebstat.
%
% Revision 1.10  1999/03/02 15:20:02  dolzmann
% Use private gb_vdpilcomb1 and gb_vdpilcomb1r. Only slightly faster.
%
% Revision 1.9  1999/02/25 16:04:56  sturm
% Added switch gbverbose: output pairs left, and frequency of heuristic
% content reduction when altered.
% Added switch gbcontred and fluid gb_mincontred!*.
% Switch groebstat works independently from trgroeb.
% Use ioto for printing statistics.
%
% Revision 1.8  1999/02/25 09:24:36  sturm
% Fixed memq to member in gb_tr2crit.
%
% Revision 1.7  1999/02/25 08:27:45  sturm
% Initialize gb_strangecount!* and spac in groebner2.
%
% Revision 1.6  1999/02/25 06:11:21  dolzmann
% Fixed a bug in gb_tr3crit. gb is as fast as grobener on p6_sc wrt.
% revgradlex and only 10 times slower as groebner wrt. lex.
%
% Revision 1.5  1999/02/24 10:48:43  sturm
% Added switch and fluid declarations and binding. In particular
% !*groebfullreduction, !*gtraverso!-sloppy, !*ezgcd, !*groebdivide.
% gb is faster than groebner on p6_sc now.
%
% Revision 1.4  1999/02/24 08:44:28  sturm
% Added gb_searchinlist. Successfully computes gb({x*y+1,x*z+1}).
%
% Revision 1.3  1999/02/24 08:36:13  sturm
% Added gb_vdpvordopt including patch code. Fixed calls to tt.
%
% Revision 1.2  1999/02/23 16:43:27  sturm
% Checked and corrected by AD. Compiles now.
%
% Revision 1.1  1999/02/23 16:41:38  sturm
% Initial check-in. Obtained from REDUCE 3.6 Groebner by selecting and
% reformatting.
%
% ----------------------------------------------------------------------
lisp <<
   fluid '(gb_rcsid!* gb_copyright!*);
   gb_rcsid!* := "$Id: gb.red,v 1.31 1999/04/15 07:08:01 dolzmann Exp $";
   gb_copyright!* := "Copyright (c) 1999 by A. Dolzmann and T. Sturm"
>>;

module gb;

%load!-package 'dp;
load!-package 'ezgcd;
remflag('(load!-package),'eval);  % for bootstrapping
load!-package 'rltools;
flag('(load!-package),'eval);

switch gbltbasis,groebopt,cgbstat,cgbfullred,cgbverbose,cgbcontred,
   cgbcounthf,cgbcheckg;

fluid '(!*gltbasis !*groebopt !*cgbstat !*cgbfullred !*cgbverbose
   !*cgbcontred !*cgbcounthf !*cgbcheckg !*cgbsloppy);

off1 'cgbstat;
on1 'cgbfullred;
off1 'cgbverbose;
off1 'cgbcontred;
off1 'cgbcounthf;
off1 'cgbcheckg;
!*cgbsloppy := T;

switch cgbupdb;
fluid '(cgb_updbcount!* cgb_updbcountp!* cgb_updbcalls!* !*cgbupdb);
off1 'cgbupdb;

fluid '(intvdpvars!*);

fluid '(!*gcd !*ezgcd !*factor !*exp dmode!* depl!* !*backtrace);

fluid '(secondvalue!* thirdvalue!*);

fluid '(dip_vars!* vdp_pcount!*);

fluid '(cgb_hcount!* cgb_hzerocount!* cgb_tr1count!* cgb_tr2count!*
   cgb_tr3count!* cgb_b4count!* cgb_strangecount!* cgb_paircount!*
   cgb_hfaccount!* cgb_gcount!* cgb_gstat!* cgb_gbcount!*);

global '(gvarslast gltb);

gltb := {'list};

fluid '(cgb_contcount!* cgb_mincontred!*);
cgb_mincontred!* := 20;  % originally 10

put('gb,'psopfn,'gb_gb!$);

smacro procedure gb_tt(s1,s2);
   % lcm of leading terms of s1 and s2
   ev_lcm(vdp_evlmon s1,vdp_evlmon s2);

procedure gb_gb!$(u);
   % Groebner bases GB computation for the AM. [u] is a list of length
   % 1 and [car u] evaluates to an AMPSYS. Returns an AMPSYS. Computes
   % a GB of the ideal generated by all polynomials in [u] wrt. to the
   % current polynomial ring.
   begin scalar w,v,odipenv;
      if length u neq 1 then
      	 rederr {"gb called with",length u,"arguments instead of 1 argument"};
      gb_domainchk();
      u := gb_a2s!-psys reval car u;
      v := gb_vars(u,td_vars());
      gvarslast := 'list . car v;
      odipenv := vdp_init(car v,td_sortmode(),td_sortextension());
      w := errorset({'gb_gb1!$,mkquote u,mkquote car v,mkquote cdr v},
      	 T,!*backtrace);
      vdp_cleanup odipenv;
      if errorp w then
      	 rederr "Error during gb computation";
      return car w
   end;

procedure gb_gb1!$(u,m,p);
   % Groebner bases GB computation for the AM subroutine. [u] is a
   % list of SF's; [m] and [p] are list of variables. Returns an
   % AMPSYS. Computes a GB of the ideal generated by all polynomials
   % in [u] wrt. to the current polynomial ring. [m] is the list of
   % the main variables and [p] is the list of parameters.
   begin scalar w;
      u := for each x in u collect vdp_f2vdp x;
      w := gb_gb u;
      if !*gltbasis then
         gltb := gb_gb2gltb w;
      return gb_s2a!-gb w
   end;

procedure gb_gsys!$(u);
   % Groebner bases Groebner system computation for the AM. [u] is a
   % list of length 1 and [car u] evaluates to an AMPSYS. Returns an
   % AMPSYS. Computes a Groebner system of the ideal generated by all
   % polynomials in [u] wrt. to the current polynomial ring.
   begin scalar w,v,odipenv;
      if length u neq 1 then
      	 rederr {"gsys called with",length u,"arguments instead of 1 argument"};
      gb_domainchk();
      u := gb_a2s!-psys reval car u;
      v := gb_vars(u,td_vars());
      gvarslast := 'list . car v;
      odipenv := vdp_init(car v,td_sortmode(),td_sortextension());
      w := errorset({'gb_gsys1!$,mkquote u,mkquote car v,mkquote cdr v},
      	 T,!*backtrace);
      vdp_cleanup odipenv;
      if errorp w then
      	 rederr "Error during gb computation";
      return car w
   end;

procedure gb_gsys1!$(u,m,p);
   <<
      u := for each x in u collect vdp_f2vdp x;
      {'list,{'list,rl_mk!*fof 'true,gb_s2a!-gb gb_gb u}}
   >>;

procedure gb_gbf(u,m,sm,sx);
   % Groebner bases GB computation for forms. [u] is a list of SF's;
   % [m] and [p] are list of variables, [sm] is a term order and [sx]
   % is a list of additional arguments to [sm]. Returns a list of
   % SF's.
   begin scalar w,v,odipenv;
      v := gb_vars(u,m);
      odipenv := vdp_init(car v,sm,sx);      
      w := errorset({'gb_gbf1,mkquote u,mkquote car v,mkquote cdr v},
      	 T,!*backtrace);
      vdp_cleanup odipenv;
      if errorp w then
      	 rederr "Error during gb computation";
      return car w
   end;

procedure gb_gbf1(u,m,p);
   % Groebner bases GB computation form subroutine. [u] is a list of
   % SF's; [m] and [p] are list of variables. Returns a list of SF's.
   <<
      u := for each x in u collect vdp_f2vdp x;
      gb_gb!-sfl gb_gb u
   >>;

procedure gb_gb!-sfl(u);
   % Groebner bases GB to SF list. [u] is a list of CGP's. Returns a
   % list of SF's.
   for each p in u collect vdp_2f p;

procedure gb_a2s!-psys(l);
   % Groebner bases algebraic mode to symbolic mode polynomial system.
   % [l] is an AMPSYS. Returns an FPSYS.
   begin scalar w,resl;
      for each j in getrlist l do <<
      	 w := numr simp j;
	 if w and not(w member resl) then
	    resl := w . resl
      >>;
      return sort(resl,'ordp)
   end;

procedure gb_s2a!-gb(u);
   % Groebner bases symbolic mode to algebraic mode GB. [u] is a list
   % of CGP's. Returns an AMPSYS.
   'list . for each x in u collect vdp_2a x;

procedure gb_domainchk();
   % Groebner bases domain check. No argument. Return value not
   % defined. Raises an error if the current domain is not valid for
   % GB computations.
   if not memq(dmode!*,'(nil)) then
      rederr bldmsg("gb does not support domain: %w",get(dmode!*,'dname));
      
procedure gb_vars(l,vl);     %DROPPED: depend,rules,zero divisors.
   % Groebner bases variables. [l] is a list of SF's; [vl] is the list
   % of main variables. Returns a pair $(m . p)$ where $m$ and $p$ are
   % list of variables. $m$ is the list of used main variables and $p$
   % is the list of used parameters. An emty [vl] means use all
   % variables as main variables.
   begin scalar w,m,p;
      for each f in l do
	 w := union(w,kernels f);
      if vl then <<
      	 m := gb_intersection(vl,w);
      	 p := setdiff(w,vl)
      >> else
	 m := w;
      return gb_varsopt(l,m) . p
   end;

procedure gb_intersection(a,b);
   % Groebner bases intersection. [a] and [b] are lists. Returns a
   % list. The returned list contains all elements occuring in [a] and
   % in [b]. The order of the elements is the same as in [a].
   for each x in a join
      if x member b then
	 {x};

procedure gb_varsopt(l,vl);
   % Groebner bases variables optimize. [l] is a list of SF's; [vl] is
   % the list of main variables. Returns a possibly reorderd list of
   % main variables.
   if !*groebopt and td_sortmode() memq '(lex gradlex revgradlex) then
      gb_vdpvordopt(l,vl)
   else
      vl;

procedure gb_gb2gltb(base);
   'list . for each j in base collect
      vdp_2a vdp_fmon(bc_a2bc 1,vdp_evlmon j);

procedure gb_gb(p);
   begin scalar spac,p1,savetime,!*factor,!*exp,intvdpvars!*,!*gcd,!*ezgcd,
	 dip_vars!*,secondvalue!*,thirdvalue!*,cgb_gstat!*;
      integer vdp_pcount!*,cgb_contcount!*,cgb_hcount!*,cgb_hzerocount!*,
	 cgb_tr1count!*,cgb_tr2count!*,cgb_tr3count!*,cgb_b4count!*,
	 cgb_strangecount!*,cgb_paircount!*,cgb_hfaccount!*,cgb_gcount!*,
	 cgb_gbcount!*,cgb_updbcount!*,cgb_updbcountp!*,cgb_updbcalls!*;
      !*exp := !*gcd := !*ezgcd := T;
      if !*cgbstat then
      	 savetime := time();
      if !*cgbcheckg then	 
      	 cgb_gstat!* := nil;
      cgb_contcount!* := cgb_mincontred!*;
      if !*cgbstat then
      	 spac := gctime();
      p1 := if !*cgbupdb then
 	 gb_traverso!-sturm!-experimental p
      else
 	 gb_traverso p;
      if !*cgbstat then <<
	 ioto_tprin2t "Statistics for GB computation:";
	 ioto_prin2t {"Time: ",time() - savetime," ms plus GC time: ",
	    gctime() - spac," ms"};
	 ioto_prin2t {"H-polynomials total: ",cgb_hcount!*};
	 ioto_prin2t {"H-polynomials zero: ",cgb_hzerocount!*};
	 if !*cgbcounthf then
	    ioto_prin2t {"H-polynomials reducible: ",cgb_hfaccount!*};
	 if !*cgbcheckg then
	    ioto_prin2t {"H-polynomials gaussible: ",cgb_gcount!*,
	       " ",cgb_gstat!*};
	 ioto_prin2t {"Crit Tr1 hits: ",cgb_tr1count!*};
	 ioto_prin2t {"Crit B4 hits: ",cgb_b4count!*," (Buchberger 1)"};
	 ioto_prin2t {"Crit Tr2 hits: ",cgb_tr2count!*};
	 ioto_prin2t {"Crit Tr3 hits: ",cgb_tr3count!*};
	 if !*cgbupdb then
	    ioto_prin2t {"updbase: calls ",cgb_updbcalls!*,", del ",
	       cgb_updbcountp!*,"/",cgb_updbcount!*};
	 ioto_prin2t {"Strange reductions: ",cgb_strangecount!*}
      >>;
      return p1
   end;

procedure gb_traverso!-sturm!-experimental(g0);
   begin scalar gall,g,d,s,h,p;
      g0 := for each fj in g0 join
	 if not vdp_zero!? fj then
	    {vdp_setsugar(vdp_enumerate vdp_simpcont fj,vdp_tdeg fj)};
      for each h in g0 do <<  % create initial critical pairs
	 p := {nil,h,h};
      	 h := gb_enumerate h;
      	 d := gb_traverso!-pairlist(h,g,d);
	 g := nconc(g,{h});
	 gall := nconc(gall,{h});
      >>;
      while d do <<  % critical pairs left
	 if !*cgbverbose then <<
	    ioto_prin2 {"[",cgb_paircount!*,"] "};
	    cgb_paircount!* := cgb_paircount!* #- 1
	 >>;
	 p := car d;
	 d := cdr d;
	 s := gb_spolynomial(p);
	 h := gb_simpcontnormalform gb_normalform(s,gall);
	 cgb_hcount!* := cgb_hcount!* #+ 1;
	 if vdp_zero!? h then <<
	    cgb_hzerocount!* := cgb_hzerocount!* #+ 1;
	 >> else if ev_zero!? vdp_evlmon h then <<  % base 1 found
	    h := gb_enumerate h;
	    d := nil;
	    g := {h}
	 >> else <<
	    if !*cgbcounthf then
	       if cddr fctrf vdp_2f h then
		  cgb_hfaccount!* := cgb_hfaccount!* #+ 1;
      	    h := gb_enumerate h;
      	    d := gb_traverso!-pairlist(h,g,d);
      	    gall := gb_updbase(gall,h);
      	    g := nconc(g,{h})
	 >>
      >>;
      return gb_traverso!-final g
   end;

procedure gb_traverso(g0);
   begin scalar g,d,s,h,p,gstat;
      g0 := for each fj in g0 join
	 if not vdp_zero!? fj then
	    {vdp_setsugar(vdp_enumerate vdp_simpcont fj,vdp_tdeg fj)};
      for each h in g0 do <<  % create initial critical pairs
	 p := {nil,h,h};
      	 h := gb_enumerate h;
      	 d := gb_traverso!-pairlist(h,g,d);
	 g := nconc(g,{h})
      >>;
      if !*cgbverbose then
	 cgb_gbcount!* := length g;
      while d do <<  % critical pairs left
	 if !*cgbverbose then <<
	    ioto_prin2 {"[",cgb_paircount!*,"] "};
	    cgb_paircount!* := cgb_paircount!* #- 1
	 >>;
	 p := car d;
	 d := cdr d;
	 s := gb_spolynomial(p);
	 h := gb_simpcontnormalform gb_normalform(s,g);
	 if !*cgbstat then
	    cgb_hcount!* := cgb_hcount!* #+ 1;
	 if vdp_zero!? h then
	    (if !*cgbstat then
	       cgb_hzerocount!* := cgb_hzerocount!* #+ 1)	    
	 else if vdp_unit!? h then <<
	    h := gb_enumerate h;
	    d := nil;
	    g := {h}
	 >> else <<
	    if !*cgbcounthf then
	       if cddr fctrf vdp_2f h then
		  cgb_hfaccount!* := cgb_hfaccount!* #+ 1;
	    if !*cgbcheckg then <<
	       gstat := gb_chkgauss vdp_poly h;
	       if 1 member gstat then <<
		  cgb_gcount!* := cgb_gcount!* #+ 1;
	       	  cgb_gstat!* := gb_chkgauss!-stat2vl gstat . cgb_gstat!*
	       >>
	    >>;
      	    h := gb_enumerate h;
      	    d := gb_traverso!-pairlist(h,g,d);
      	    g := nconc(g,{h});
	    if !*cgbverbose then
	       cgb_gbcount!* := cgb_gbcount!* #+ 1
	 >>
      >>;
      return gb_traverso!-final g
   end;

procedure gb_updbase(g,h);
   begin scalar hev,oc;
      if !*cgbstat then
      	 oc := cgb_updbcountp!*;
      hev := vdp_evlmon h;
      g := for each p in g join
      	 if not ev_divides!?(hev,vdp_evlmon p) then
 	    {p}
	 else <<
	    if !*cgbverbose then
	       ioto_prin2 "#";
	    if !*cgbstat then
	       cgb_updbcountp!* := cgb_updbcountp!* #+ 1;
	    nil
	 >>;
      if !*cgbstat then <<
      	 if not (oc #= cgb_updbcountp!*) then
	    cgb_updbcount!* := cgb_updbcount!* #+ 1;
	 cgb_updbcalls!* := cgb_updbcalls!* + 1
      >>;
      return nconc(g,{h})
   end;

procedure gb_chkgauss(p);
   % [p] is a dipoly.
   begin scalar stat;
      stat := for each x in dip_vars!* collect 0;
      while p do <<
	 stat := gb_chkgauss1(dip_evlmon p,stat);
   	 p := dip_mred p
      >>;
      return stat
   end;

procedure gb_chkgauss1(ev,stat);
   begin scalar nstat,e,s,td;
      td := ev_tdeg ev;
      if td = 0 then
	 return stat;
      td := if td #= 1 then 1 else -1;
      while ev do <<
	 e := car ev;
	 ev := cdr ev;
	 s := car stat;
	 stat := cdr stat;
	 if e #=0 then
	    nstat := s . nstat
	 else if e #> 1 or s #= -1 then
	    nstat := (-1) . nstat
	 else if e #= 1 and s #= 0 then
	    nstat := td . nstat
	 else if s #=0 and e #=0 then
	    nstat := 0 . nstat
	 else
	    rederr "ich sehe es anders"
      >>;
      return reversip nstat
   end;

procedure gb_chkgauss!-stat2vl(gstat);
   begin scalar scdv,r;
      scdv := dip_vars!*;
      while gstat do <<
	 if eqcar(gstat,1) then
 	    r := car scdv . r;
	 gstat := cdr gstat;
	 scdv := cdr scdv
      >>;
      return reversip r
   end;

procedure gb_enumerate(f);
   % f is a temporary result. Prepare it for medium range storage, and
   % assign a number.
   if vdp_zero!? f then
      f
   else <<
      vdp_condense f;
      if vdp_number f #= 0 then
	 f := vdp_setnumber(f,vdp_pcount!* := vdp_pcount!* #+ 1);
      f
   >>;

procedure gb_traverso!-pairlist(gk,g,d);
   % gk: new polynomial, g: current basis, d: old pair list.
   begin scalar ev,r,n;
      d := gb_traverso!-pairs!-discard1(gk,d);
      % build new pair list:
      ev := vdp_evlmon gk;
      for each p in g do
     	 if not gb_buchcrit4t(ev,vdp_evlmon p) then <<
	    if !*cgbstat then
	       cgb_b4count!* := cgb_b4count!* #+ 1;
 	    r := ev_lcm(ev,vdp_evlmon p) . r
       	 >> else
 	    n := gb_makepair(p,gk) . n;
      n := gb_tr2crit(n,r);
      n := gb_cplistsort(n,!*cgbsloppy);
      n := gb_tr3crit n;
      if !*cgbverbose and n then <<
	 cgb_paircount!* := cgb_paircount!* #+ length n;
	 ioto_cterpri();
	 ioto_prin2 {"(",cgb_gbcount!*,") "}
      >>;
      return gb_cplistmerge(d,reversip n)
   end;

procedure gb_tr2crit(n,r);
   % delete equivalents to coprime lcm
   for each p in n join 
      if ev_member(car p,r) then <<
	 if !*cgbstat then
	    cgb_tr2count!* := cgb_tr2count!* #+ 1;
	 nil
      >> else
	 {p};

procedure gb_tr3crit(n);
   begin scalar newn,scannewn,q;
      for each p in n do <<
	 scannewn := newn;
	 q := nil;
 	 while scannewn do
	    if ev_divides!?(caar scannewn,car p) then <<
	       q := t;
	       scannewn := nil;
	       if !*cgbstat then
	       	  cgb_tr3count!* := cgb_tr3count!* #+ 1
	    >> else
	       scannewn := cdr scannewn;
	 if not q then
	    newn := gb_cplistsortin(p,newn,nil)
      >>;
      return newn
   end;

procedure gb_traverso!-pairs!-discard1(gk,d);
   % crit B. Delete triange relations.
   for each pij in d join
      if gb_traverso!-trianglep(cadr pij,caddr pij,gk,car pij) then <<
	 if !*cgbstat then
	    cgb_tr1count!* := cgb_tr1count!* #+ 1;
	 if !*cgbverbose then
	    cgb_paircount!* := cgb_paircount!* #- 1;
	 nil
      >> else
	 {pij};

procedure gb_traverso!-trianglep(gi,gj,gk,tij);
   ev_sdivp(gb_tt(gi,gk),tij) and ev_sdivp(gb_tt(gj,gk),tij);

procedure gb_traverso!-final(g);
   % Final reduction and sorting.
   for each rg on vdp_lsort g join
      if not gb_searchinlist(vdp_evlmon car rg,cdr rg) then
	 {gb_simpcontnormalform gb_normalform(car rg,cdr rg)};

procedure gb_buchcrit4t(e1,e2);
   % nonconstructive test of lcm(e1,e2) = e1 + e2 equivalent: no
   % matches of nonzero elements.
   not ev_disjointp(e1,e2);
   
procedure gb_cplistsort(g,sloppy);
   begin scalar gg;
      for each p in g do
 	 gg := gb_cplistsortin(p,gg,sloppy);
      return gg
   end;

procedure gb_cplistsortin(p,pl,sloppy);
   % Distributive polynomial critical pair list sort. pl is a special
   % list for Groebner calculation, p is a pair. Returns the updated
   % list pl (p sorted into).
   if null pl then
      {p}
   else <<
      gb_cplistsortin1(p,pl,sloppy);
      pl
   >>;

procedure gb_cplistsortin1(p,pl,sloppy);
   % Destructive insert of p into nonnull pl.
   if not gb_cpcompless!?(car pl,p,sloppy) then <<
      rplacd(pl,car pl . cdr pl);
      rplaca(pl,p)
   >> else if null cdr pl then
      rplacd(pl,{p})
   else
      gb_cplistsortin1(p,cdr pl,sloppy);

procedure gb_cpcompless!?(p1,p2,sloppy);
   % Compare 2 pairs wrt. their sugar (=cadddr) or their lcm (=car).
   if sloppy then
      ev_compless!?(car p1,car p2)
   else
      gb_cpcompless!?s(p1,p2,cadddr p1 #- cadddr p2,ev_comp(car p1,car p2));
	 
procedure gb_cpcompless!?s(p1,p2,d,q);
   if d neq 0 then
      d #< 0
   else if q neq 0 then
      q < 0
   else
      vdp_number(caddr p1) #< vdp_number(caddr p2);

procedure gb_cplistmerge(pl1,pl2);
   % Distributive polynomial critical pair list merge. pl1 and pl2 are
   % critical pair lists used in the Groebner calculation.
   % groebcplistmerge(pl1,pl2) returns the merged list.
   begin scalar cpl1,cpl2;
      if null pl1 then
 	 return pl2;
      if null pl2 then
 	 return pl1;
      cpl1 := car pl1;
      cpl2 := car pl2;
      return if gb_cpcompless!?(cpl1,cpl2,nil) then
 	 cpl1 . gb_cplistmerge(cdr pl1,pl2)
      else
 	 cpl2 . gb_cplistmerge(pl1,cdr pl2)
   end;

procedure gb_makepair(f,h);
   % Construct a pair from polynomials f and h.
   begin scalar ttt,sf,sh;
      ttt := gb_tt(f,h);
      sf := vdp_sugar(f) #+ ev_tdeg ev_dif(ttt,vdp_evlmon f);
      sh := vdp_sugar(h) #+ ev_tdeg ev_dif(ttt,vdp_evlmon h);
      return {ttt,f,h,ev_max!#(sf,sh)}
   end;

procedure gb_spolynomial(pr);
   begin scalar p1,p2,s;
      p1 := cadr pr;
      p2 := caddr pr;
      s := gb_spolynomial1(p1,p2);   % TODO: Switch for strange reduction
      if vdp_zero!? s or vdp_unit!? s then
 	 return s;
%      return vdp_setsugar(gb_strange!-reduction(s,p1,p2),cadddr pr)
      return gb_strange!-reduction(s,p1,p2)  % TODO: normal suger for
					     % special cases. 
   end;

procedure gb_spolynomial1(p1,p2);
   begin scalar ep1,ep2,ep,rp1,rp2,db1,db2,x;
      ep1 := vdp_evlmon p1;
      ep2 := vdp_evlmon p2;
      ep := ev_lcm(ep1,ep2);
      rp1 := vdp_mred p1;
      rp2 := vdp_mred p2;
      if vdp_zero!? rp1 and vdp_zero!? rp2 then
 	 return rp1;
      if vdp_zero!? rp1 then
	 return vdp_prod(rp2,vdp_fmon(bc_a2bc 1,ev_dif(ep,ep2)));
      if vdp_zero!? rp2 then
	 return vdp_prod(rp1,vdp_fmon(bc_a2bc 1,ev_dif(ep,ep1)));
      db1 := vdp_lbc p1;
      db2 := vdp_lbc p2;
      x := bc_gcd(db1,db2);
      if not bc_one!? x then <<
	 db1 := bc_quot(db1,x);
	 db2 := bc_quot(db2,x)
      >>;
      return vdp_ilcomb(rp2,db1,ev_dif(ep,ep2),rp1,bc_neg db2,ev_dif(ep,ep1))
   end;

procedure gb_strange!-reduction(s,p1,p2);
   % Subtracts multiples of p2 from the s-polynomial such that head
   % terms are eliminated early.
   begin scalar tp1,tp2,ts,c,saves;
      saves := s;
      tp1 := vdp_evlmon p1;
      tp2 := vdp_evlmon p2;
      c := T; while c and not vdp_zero!? s do <<
	 ts := vdp_evlmon s;
	 if gb_buch!-ev_divides!?(tp2,ts) then
	    s := gb_reduceonestepint(s,vdp_zero(),vdp_lbc s,ts,p2)
	 else if gb_buch!-ev_divides!?(tp1,ts) then
	    s := gb_reduceonestepint(s,vdp_zero(),vdp_lbc s,ts,p1)
	 else
	    c := nil
      >>;
      if !*cgbstat and not (s eq saves) then
	 cgb_strangecount!* := cgb_strangecount!* #+ 1;
      return s
   end;

procedure gb_buch!-ev_divides!?(vev1,vev2);
   % Test if vev1 divides vev2 for exponent vectors vev1 and vev2.
   ev_mtest!?(vev2,vev1);

procedure gb_normalform(f,g);
   % General procedure for the reduction of one polynomial modulo a
   % set. f is a polynomial, g is a set of polynomials. Procedure
   % behaves like type='list. f is reduced modulo g.
   begin scalar f1,c,vev,divisor,tai,fold; integer n;
      fold := f;
      f1 := vdp_setsugar(vdp_zero(),vdp_sugar f);
      while not vdp_zero!? f do <<
	 vev := vdp_evlmon f;
 	 c := vdp_lbc f;
	 divisor := gb_searchinlist(vev,g);
	 if divisor then <<
	    tai := T;
	    if vdp_monp divisor then
	       f := vdp_cancelmev(f,vdp_evlmon divisor)
	    else <<
	       f := gb_reduceonestepint(f,f1,c,vev,divisor);
	       f1 := secondvalue!*;
	       if !*cgbcontred then <<
	       	  f := gb_adtssimpcont(f,f1,n);
	       	  f1 := secondvalue!*;
	       	  n := thirdvalue!*
	       >>
	    >>
	 >> else if !*cgbfullred then <<
	    f := gb_shift(f,f1);
	    f1 := secondvalue!*
	 >> else <<
	    f1 := vdp_sum(f1,f);
 	    f := vdp_zero()
	 >>
      >>;
      return if tai then f1 else fold
   end;

procedure gb_searchinlist(vev,g);
   % search for a polynomial in the list G, such that the lcm divides
   % vev; G is expected to be sorted in descending sequence.
   if null g then
      nil
   else if gb_buch!-ev_divides!?(vdp_evlmon car g, vev) then
      car g
   else
      gb_searchinlist(vev,cdr g);

procedure gb_adtssimpcont(f,f1,n);
   begin scalar f0;
      if vdp_zero!? f then <<
      	 secondvalue!* := f1;
      	 thirdvalue!* := 0;
      	 return f
      >>;
      n := n + 1;
      if n #> cgb_contcount!* then <<
      	 f0 := f;
      	 f := gb_simpcont2(f,f1);
      	 f1 := secondvalue!*;
      	 gb_contentcontrol(f neq f0);
	 n := 0
      >>;
      secondvalue!* := f1;
      thirdvalue!* := n;
      return f
   end;

procedure gb_simpcont2(f,f1);
   % Simplify two polynomials with the gcd of their contents. [f] is
   % not zero.
   begin scalar c,s1,s2;
      c := vdp_content f;
      if bc_one!? bc_abs c then <<
	 secondvalue!* := f1;
	 return f
      >>;
      s1 := vdp_sugar f;
      s2 := vdp_sugar f1;
      if not vdp_zero!? f1 then <<
 	 c := vdp_content1(f1,c);
	 if bc_one!? bc_abs c then <<
	    secondvalue!* := f1;
	    return f
      	 >>;
	 f1 := vdp_bcquot(f1,c)
      >>;
      f := vdp_bcquot(f,c);
      vdp_setsugar(f,s1);
      vdp_setsugar(f1,s2);
      secondvalue!* := f1;
      return f
    end;

procedure gb_contentcontrol(u);
   % u indicates, that a substantial content reduction was done;
   % update content reduction limit from u.
   <<
      cgb_contcount!* := if u then
      	 ev_max!#(0,cgb_contcount!* #- 1)
      else
      	 gb_min!#(cgb_mincontred!*,cgb_contcount!* #+ 1);
      if !*cgbverbose then
	 ioto_prin2 {"<",cgb_contcount!*,"> "}
   >>;

procedure gb_min!#(a,b);
   if a #< b then a else b;

procedure gb_shift(f,f1);
   begin scalar s,s1;
      s := vdp_sugar f;
      s1 := vdp_sugar f1;
      f1 := vdp_nconcmon(f1,vdp_lbc f,vdp_evlmon f);
      f := vdp_mred f;
      vdp_setsugar(f,s);
      vdp_setsugar(f1,ev_max!#(s1,s));
      secondvalue!* := f1;
      return f
   end;

procedure gb_reduceonestepint(f,f1,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 vevcof,a,b,cg,x,rg1;
      rg1 := vdp_mred g1;
      if vdp_zero!? rg1 then <<  % g1 is monomial
	 f := vdp_mred f;
	 secondvalue!* := f1;
	 return f
      >>;
      vevcof := ev_dif(vev,vdp_evlmon g1);  % nix lcm
      cg := vdp_lbc g1;
      x := bc_gcd(c,cg);
      a := bc_quot(cg,x);
      b := bc_quot(c,x);
      % multiply relevant parts of f and f1 by a (vbc)
      if not vdp_zero!? f1 then
 	 f1 := vdp_bcprod(f1,a);
      f := vdp_ilcombr(vdp_mred f,a,rg1,bc_neg b,vevcof);
      secondvalue!*:= f1;
      return f
   end;

procedure gb_simpcontnormalform(h);
   % simpCont version preserving the property SUGAR.
   if vdp_zero!? h then
      h
   else
      vdp_setsugar(vdp_simpcont h,vdp_sugar h);

procedure gb_vdpvordopt(w,vars);
   % w : list of polynomials (standard forms)
   % vars: list of variables
   % return vars
   begin scalar c;
      vars := sort(vars,'ordop);
      c := for each x in vars collect x . 0 . 0;
      for each poly in w do
 	 gb_vdpvordopt1(poly,vars,c);
      c := sort(c,function gb_vdpvordopt2);
      intvdpvars!* := for each v in c collect car v;
      vars := gb_vdpvordopt31 intvdpvars!*;
      return vars
   end;
 
procedure gb_vdpvordopt31(u);
   begin scalar v,y;
      if null u then
 	 return nil;
      v := for each x in u join <<
 	 y := assoc(x,depl!*);
	 if null y or null xnp(cdr y,u) then
 	    {x}
      >>;
      return nconc(gb_vdpvordopt31 setdiff(u,v),v)
   end;

procedure gb_vdpvordopt1(p,vl,c);
   if null p then
      0
   else if domainp p or null vl then
      1
   else if mvar p neq car vl then
      gb_vdpvordopt1(p,cdr vl,c)
   else begin scalar var,pow,slot; integer n;
      n := gb_vdpvordopt1 (lc p,cdr vl,c);
      var := mvar p;
      pow := ldeg p;
      slot := assoc(var,c);
      if pow #> cadr slot then <<
	 rplaca(cdr slot,pow);
 	 rplacd(cdr slot,n)
      >> else
	 rplacd(cdr slot,n #+ cddr slot);
      return n #+ gb_vdpvordopt1 (red p,vl,c)
   end;

procedure gb_vdpvordopt2(sl1,sl2);
   % compare two slots from the power table
   <<
      sl1 := cdr sl1;
      sl2 := cdr sl2;
      car sl1 #< car sl2 or car sl1 = car sl2 and cdr sl1 #< cdr sl2
   >>;

endmodule;  % gb

module vdp;

procedure vdp_lbc(u);
   caddr u;

procedure vdp_evlmon(u);
   cadr u;

procedure vdp_poly(u);
   car cdddr u;

procedure vdp_zero!?(u);
   null vdp_poly u;

procedure vdp_plist(u);
   cadr cdddr u;

procedure vdp_number(f);
   vdp_getprop(f,'number) or 0;

procedure vdp_sugar(f);
   if vdp_zero!? f then 0 else vdp_getprop(f,'sugar);  % 0 notwendig?

procedure vdp_unit!?(p);
   not vdp_zero!? p and ev_zero!? vdp_evlmon p;

procedure vdp_make(vbc,vev,form);
   {'vdp,vev,vbc,form,nil};

procedure vdp_monp(u);
   dip_monp vdp_poly u;

procedure vdp_tdeg(u);
   dip_tdeg vdp_poly u;

procedure vdp_fdip(u);
   if null u then
      vdp_zero()
   else
      vdp_make(dip_lbc u,dip_evlmon u,u);

procedure vdp_appendmon(vdp,coef,vev);
   % Add a monomial to the end of a vdp (vev remains unchanged).
   if vdp_zero!? vdp then
      vdp_fmon(coef,vev)
   else if bc_zero!? coef then
      vdp
   else
      vdp_make(vdp_lbc vdp,vdp_evlmon vdp,dip_appendmon(vdp_poly vdp,coef,vev));

procedure vdp_nconcmon(vdp,coef,vev);
   % Add a monomial to the end of a vdp (vev remains unchanged).
   if vdp_zero!? vdp then
      vdp_fmon(coef,vev)
   else if bc_zero!? coef then
      vdp
   else
      vdp_make(vdp_lbc vdp,vdp_evlmon vdp,dip_nconcmon(vdp_poly vdp,coef,vev));

procedure vdp_bcquot(p,c);
   begin scalar r;
      r := vdp_fdip dip_bcquot(vdp_poly p,c);
      vdp_setsugar(r,vdp_sugar p);
      return r
   end;

procedure vdp_content(p);
   dip_contenti vdp_poly p;

procedure vdp_content1(d,c);
   dip_contenti1(vdp_poly d,c);

procedure vdp_length(f);
   dip_length vdp_poly f;

procedure vdp_bcprod(p,b);
   begin scalar r;
      r := vdp_fdip dip_bcprod(vdp_poly p,b);
      vdp_setsugar(r,vdp_sugar p);
      return r
   end;

procedure vdp_cancelmev(p,vev);
   begin scalar r;
      r := vdp_fdip dip_cancelmev(vdp_poly p,vev);
      vdp_setsugar(r,vdp_sugar p);
      return r
   end;

procedure vdp_sum(d1,d2);
   begin scalar r;
      r := vdp_fdip dip_sum(vdp_poly d1,vdp_poly d2);
      vdp_setsugar(r,ev_max!#(vdp_sugar d1,vdp_sugar d2));
      return r
   end;

procedure vdp_prod(d1,d2);
   begin scalar r;
      r := vdp_fdip dip_prod(vdp_poly d1,vdp_poly d2);
      vdp_setsugar(r,vdp_sugar d1 #+ vdp_sugar d2);
      return r
   end;

procedure vdp_zero();
   vdp_make('invalid,'invalid,nil);   

procedure vdp_mred(u);
   begin scalar r;
      r := dip_mred vdp_poly u;
      if null r then
 	 return vdp_zero();
      r := vdp_make(dip_lbc r,dip_evlmon r,r);
      vdp_setsugar(r,vdp_sugar(u));
      return r
  end;

procedure vdp_condense(f);
   dip_condense vdp_poly f;

procedure vdp_setsugar(p,s);
   vdp_putprop(p,'sugar,s);

procedure vdp_setnumber(p,n);
   vdp_putprop(p,'number,n);

procedure vdp_putprop(poly,prop,val);
   begin scalar c,p;
      c := cdr cdddr poly;
      p := atsoc(prop,car c);
      if p then
 	 rplacd(p,val)
      else
 	 rplaca(c,(prop . val) . car c);
      return poly
   end;

procedure vdp_getprop(poly,prop);
   (if p then cdr p else nil) where p=atsoc(prop,vdp_plist poly);

procedure vdp_fmon(coef,vev);
   begin scalar r;
      r := vdp_make(coef,vev,dip_fmon(coef,vev));
      vdp_setsugar(r,ev_tdeg vev);
      return r
   end;

procedure vdp_2a(u);
   dip_2a vdp_poly u;

procedure vdp_2f(u);
   dip_2f vdp_poly u;

procedure vdp_init(vars,sm,sx);
   % Initializing vdp-dip polynomial package.
   dip_init(vars,sm,sx);

procedure vdp_cleanup(l);
   dip_cleanup(l);

procedure vdp_f2vdp(u);
   begin scalar dip;
      dip := dip_f2dip u;
      if null dip then
	 return vdp_zero();
      return vdp_make(dip_lbc dip,dip_evlmon dip,dip)
   end;

procedure vdp_enumerate(f);
  % f is a temporary result. Prepare it for medium range storage and
  % assign a number.
  if vdp_zero!? f or vdp_getprop(f,'number) then
     f
  else
     vdp_putprop(f,'number,(vdp_pcount!* := vdp_pcount!* #+ 1));

procedure vdp_simpcont(p);
   begin scalar q;
      q := vdp_poly p;
      if null q then
 	 return p;
      return vdp_fdip dip_simpcont q
   end;

procedure vdp_lsort(pl);
   % Distributive polynomial list sort. pl is a list of distributive
   % polynomials. vdplsort(pl) returns the sorted distributive
   % polynomial list of pl.
   sort(pl,function vdp_evlcomp);

procedure vdp_evlcomp(p1,p2);  % nicht auf das HM?
   dip_evlcomp(vdp_poly p1,vdp_poly p2);

procedure vdp_ilcomb(v1,c1,t1,v2,c2,t2);
   begin scalar r;
      r := vdp_fdip dip_ilcomb(vdp_poly v1,c1,t1,vdp_poly v2,c2,t2);
      vdp_setsugar(r,ev_max!#(
	 vdp_sugar v1 #+ ev_tdeg t1,vdp_sugar v2 #+ ev_tdeg t2));
      return r
   end;

procedure vdp_ilcombr(v1,c1,v2,c2,t2);
   begin scalar r;
      r := vdp_fdip dip_ilcombr(vdp_poly v1,c1,vdp_poly v2,c2,t2);
      vdp_setsugar(r,ev_max!#(vdp_sugar v1,vdp_sugar v2 #+ ev_tdeg t2));
      return r
   end;

endmodule;  % vdp

end;  % of file


REDUCE Historical
REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]