reduce-historical/r38/packages/cgb/cgb.red

% ----------------------------------------------------------------------
% $Id: cgb.red,v 1.30 2004/05/03 16:38:25 sturm Exp $
% ----------------------------------------------------------------------
% Copyright (c) 1999-2003 Andreas Dolzmann and Thomas Sturm
% ----------------------------------------------------------------------
% $Log: cgb.red,v $
% Revision 1.30  2004/05/03 16:38:25  sturm
% Hopefully clean solution for deadlock with CGB/REDLOG compilation.
%
% Revision 1.29  2003/10/21 10:24:36  gilch
% Incorporated new module gbsc.
%
% Revision 1.28  2003/10/12 14:50:30  sturm
% The bootstrapping technique via remflag('(load!-package),'eval); does
% not work for CSL. Added corresponding preprocessing directive for now.
% As a consequence, under CSL "redlog" has to be loaded explicitly when
% using CGB.
%
% Revision 1.27  2003/07/17 06:30:38  dolzmann
% Added new argument xvarl to Groebner system computation. xvarl is a list
% of variables. If cgbgen is on gsys makes no assumptions on variaboles
% in xvarl.
%
% Revision 1.26  2003/05/20 08:17:38  dolzmann
% Moved cd_init to the beginning of cgb_interface!$.
% This may be neccessary for the a2s procedures.
%
% Revision 1.25  2003/05/20 07:38:26  dolzmann
% Do not load modules belongig to this package.
%
% Revision 1.24  2003/05/20 07:24:46  dolzmann
% Moved macro *_mkinterface to the right place.
%
% Revision 1.23  2003/05/19 10:27:18  dolzmann
% Adapted to the Groebner simplifier using gb.
% Added interface generator as in gb.
% Used interface generator for creating all interfaces.
% Added first version of a wrapper for non-parametric input.
% Removed old interface code.
% Introduced initial theory for Gröbner and green Gröbner system
% computation.
%
% Revision 1.22  2003/04/17 16:14:45  dolzmann
% Added AM interface ggsys for computing green Groebner systems.
% Added switch cgbsgreen. If it is on then the green Groebner system is
% computed by computing a regular Groebner system and finally colorying
% it green. If off (derfault) the green groebner system is computed by a
% modified Groebner system computation.
% Renamed cgb_greengsysf(u) to cgb_ggsysf(u);
%
% Revision 1.21  2003/04/16 09:43:18  dolzmann
% Added (inefficient) procedure for computing green Groebner systems.
% Added procedure for computing groebner systems for SFs.
% Added and corrected some comments.
%
% Revision 1.20  1999/04/13 20:56:04  dolzmann
% Added default setting for switches.
%
% Revision 1.19  1999/04/13 18:41:21  dolzmann
% Dropped zeroes and duplicates in the input system.
% Sort Groebner systems, conditions, and (partial) bases.
% Removed switch gsugar.
% Removed main variable list and parameter list arguments from the
% entire call tree.
% Renamed all gb switches and fluids to cgb.
%
% Revision 1.18  1999/04/11 11:31:37  dolzmann
% Introduced wrappers for using the gb package in case of non-parametric
% problems.
%
% Revision 1.17  1999/04/11 09:49:05  dolzmann
% Completely rewritten the interface code for the AM and the standard
% form interface.
%
% Revision 1.16  1999/04/07 15:54:16  dolzmann
% Fixed a bug in cgb_gsys2cgb: Rewritten procedure cgb_rtgsys to handle
% the case of no main variable.
%
% Revision 1.15  1999/04/07 12:37:00  dolzmann
% Fixed a bug in cgp_monp.
% Added comments to all procedures.
%
% Revision 1.14  1999/04/07 09:27:08  dolzmann
% Added switch cgbgen and related code for computing only the generic branch.
%
% Revision 1.13  1999/04/06 12:13:59  dolzmann
% Moved procedures dip_append, dip_cp, dip_dcont, and dip_dcont1 from
% module dipto into module dip.
% Moved procedures bc_mkat, bc_dcont, and bc_2d from module bcto into the
% bc modules of the dip package.
%
% Revision 1.12  1999/04/05 09:16:46  sturm
% Do not load Redlog during complilation.
%
% Revision 1.11  1999/04/05 09:06:09  sturm
% Locally bind !*rlgsvb for calls to rl_gsd.
%
% Revision 1.10  1999/04/04 18:30:52  sturm
% Provide a standard form interface cgb_cgbf to cgb's.
%
% Revision 1.9  1999/04/04 16:46:07  sturm
% Changed cgb_groebnereval into cgb_gsys.
% Added copyright and CVS fluids.
% Added create!-package.
%
% Revision 1.8  1999/04/04 14:50:37  sturm
% Implemented switch tdusetorder.
%
% Revision 1.7  1999/04/04 14:09:31  sturm
% Moved dip_ilcomb and dip_ilcombr from cgb.red to dp.red.
% Created vdp_ilcomb and vdp_ilcombr for gb.red.
%
% Revision 1.6  1999/04/04 12:20:00  dolzmann
% The counter gb_hzerocount!* works now.
% Fixed a bug in cgp_2scpl: It was possible that the condition becomes
% inconsistent.
%
% Revision 1.5  1999/04/03 13:37:21  sturm
% cgb_groebner1a runs under errorset.
% Adapted to new dip_init/dip_cleanup.
% Bind !*msg during rl_set.
% Replaced cgb_surep and cgb_gsd by correct versions with non-renamed dipoly
% fluids.
%
% Revision 1.4  1999/04/03 11:07:29  dolzmann
% Fixed some bugs.
% The test file runs without Reduce errors.
%
% Revision 1.3  1999/04/03 10:16:16  dolzmann
% Code completely rewritten:
% Introduced splitted polynomials, data types for the Groebner system,
% for branches, and for critical pairs.
% Procedure cgb_groebner1 sets the Redlog context for the condition
% handling.
%
% Revision 1.2  1999/03/31 14:05:22  sturm
% Simple examples run.
% cgb_spolsc is mathematically not correct.
%
% Revision 1.1  1999/03/24 15:10:23  sturm
% Initial check-in. Copy of gb.red 1.16.
%
% ----------------------------------------------------------------------
lisp <<
   fluid '(cgb_rcsid!* cgb_copyright!*);
   cgb_rcsid!* := "$Id: cgb.red,v 1.30 2004/05/03 16:38:25 sturm Exp $";
   cgb_copyright!* := "Copyright (c) 1999-2003 by A. Dolzmann and T. Sturm"
>>;

% TODO:
% - Normalize green groebner systems: Detect branches containing a unit
% - Detect green monomials in RP
% - Final simplification with groebner simplifier
% - Computing reduced or pseudo reduced groeber systems.
% - Computing relatively generic and local groebner systems.

module cgb;

create!-package('(cgb gb dp gbsc),nil);

load!-package 'ezgcd;

if 'psl member lispsystem!* then
   if filestatus("$reduce/lisp/psl/$MACHINE/red/redlog.b",nil) then
      load!-package 'redlog;

if 'csl member lispsystem!* then
   if modulep 'redlog then
      load!-package 'redlog;
   
switch cgbstat,cgbfullred,cgbverbose,cgbcontred,cgbgs,cgbreal,cgbgen,
   cgbsgreen;

fluid '(!*cgbstat !*cgbfullred !*cgbverbose !*cgbcontred !*cgbgs !*cgbreal
   !*cgbgen !*cgbsloppy !*cgbcdsimpl);

off1 'cgbstat;
on1 'cgbfullred;
off1 'cgbverbose;
off1 'cgbcontred;
off1 'cgbgs;
off1 'cgbreal;
off1 'cgbgen;
off1 'cgbsgreen;  % Simulate green. Compute Gsys and colore it green
!*cgbsloppy := T;
!*cgbcdsimpl := T;

fluid '(!*cgbgreen);  % pseudo switch for computing green Gsys'

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

fluid '(cgp_pcount!* cgb_hashsize!*);

cgb_hashsize!* := 65521;  % The size of the hash table for BETA (in gbsc).

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

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

fluid '(!*rlgsvb !*rlspgs !*rlsithok);

%DS
% <AMCGB> ::= <AMPSYS>
% <AMPSYS> ::= ('list,...,<Lisp-prefix-form>,...)
% <AMGSYS> ::= ('list,...,<AMBRANCH>,...)
% <AMBRANCH> ::= ('list,<RL-Formula>,<AMPSYS>)
% <FPSYS> ::= (...,<SF>,...)
% <FGSYS> ::= (...,<FBRANCH>,...)
% <FBRANCH> ::= (<CDLIST>,<FPSYS>)
% <CDLIST> ::= (...,<RL_Formula>,...)

macro procedure cgb_mkinterface(argl);
   begin
      scalar a2sl1,a2sl2,defl,xvfn,s2a,s2s,s,
	 args,bname,len,sm,prgn,ami,smi,psval,postfix,modes;
      bname := eval nth(argl,2);
      a2sl1 := eval nth(argl,3);
      a2sl2 := eval nth(argl,4);
      defl := eval nth(argl,5);
      xvfn := eval nth(argl,6);
      s2a := eval nth(argl,7);
      s2s := eval nth(argl,8);
      s := eval nth(argl,9);
      postfix := eval nth(argl,10);
      modes := eval nth(argl,11);
      len := length a2sl1;
      args := for i := 1:len+3 collect mkid('a,i);
      if (null modes or modes eq 'sm) then <<	 
      	 sm := intern compress append('(!c !g !b !_),explode bname);     
      	 % Define the symbolic mode interface
      	 smi := intern compress nconc(explode sm,explode postfix);      
      	 prgn := {'put,mkquote smi,''number!-of!-args,len+3} . prgn;
      	 prgn := {'de,smi,args,{'cgb_interface!$,mkquote sm, mkquote a2sl1,
	    mkquote a2sl2,mkquote defl,mkquote xvfn,mkquote
	       s2a,mkquote s2s,mkquote s,T,'list . args}} . prgn
      >>;
      if (null modes or modes eq 'am) then <<
      	 % Define the algebraic mode interface      
      	 ami := bname;
      	 %      ami := intern compress append('(!g !b),explode bname);
      	 psval := intern compress nconc(explode ami,'(!! !$));
      	 prgn := {'put,mkquote ami,''psopfn,mkquote psval} . prgn;
      	 prgn := {'put,mkquote psval,''number!-of!-args,1} . prgn;
      	 prgn := {'put,mkquote psval,''cleanupfn,''cgb_cleanup} . prgn;
      	 prgn := {'de,psval,'(argl),{'cgb_interface!$,mkquote sm,
	    mkquote a2sl1,mkquote a2sl2,mkquote defl,mkquote
	       xvfn,mkquote s2a,mkquote s2s,mkquote s,nil,'argl}} . prgn;
      >>;
      return 'progn . prgn
   end;

cgb_mkinterface('cgb,'(cgb_a2s!-psys),'(cgb_a2s2!-psys),
   nil,'cgb_xvars!-psys,'cgb_s2a!-cgb,'cgb_s2s!-cgb,T,'f,nil);
cgb_mkinterface('gsys,'(cgb_a2s!-psys cgb_a2s!-cd cgb_a2s!-varl),
   '(cgb_a2s2!-psys cgb_a2s2!-cd cgb_a2s2!-varl),
   {'true,'(list)},'cgb_xvars!-psys3,'cgb_s2a!-gsys,'cgb_s2s!-gsys,T,'f,nil);
%cgb_mkinterface('gsys,'(cgb_a2s!-psys cgb_a2s!-cd),
%   '(cgb_a2s2!-psys cgb_a2s2!-cd),
%   {'true},'cgb_xvars!-psys2,'cgb_s2a!-gsys,'cgb_s2s!-gsys,T,'f,nil);
cgb_mkinterface('ggsys,'(cgb_a2s!-psys cgb_a2s!-cd cgb_a2s!-varl),
   '(cgb_a2s2!-psys cgb_a2s2!-cd cgb_a2s2!-varl),
   {'true,'(list)},'cgb_xvars!-psys3,'cgb_s2a!-gsys,'cgb_s2s!-gsys,T,'f,nil);
cgb_mkinterface('gsys2cgb,'(cgb_a2s!-gsys),'(cgb_a2s2!-gsys),
   nil,'cgb_xvars!-gsys,'cgb_s2a!-cgb,'cgb_s2s!-cgb,T,'f,nil);

put('cgb_cgb,'gb_wrapper,{'gb_gb,'(gb_a2s!-psys),'(gb_a2s2!-psys),
   nil,'gb_xvars!-psys,'gb_s2a!-gbx,'gb_s2s!-gb,T});
put('cgb_gsys,'gb_wrapper,{'gb_gbgsys,'(gb_a2s!-psys),'(gb_a2s2!-psys),
   nil,'gb_xvars!-psys,'gb_s2a!-gsys,'gb_s2s!-gsys,T});

procedure cgb_a2s!-psys(l);
   % Comprehensive 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 reval l do <<
      	 w := numr simp j;
	 if w and not(w member resl) then
	    resl := w . resl
      >>;
      return sort(resl,'ordp)
   end;

procedure cgb_a2s2!-psys(fl);
   for each x in fl collect cgp_f2cgp x;

procedure cgb_xvars!-psys(l,vl);
   cgb_vars(l,vl);

procedure cgb_xvars!-psys2(l,cd,vl);
   cgb_vars(l,vl);

procedure cgb_xvars!-psys3(l,cd,xvl,vl);
   cgb_vars(l,vl);

procedure cgb_s2a!-cgb(u);
   % Comprehensive Groebner bases symbolic mode to algebraic mode CGB.
   % [u] is a list of CGP's. Returns an AMPSYS.
   'list . for each x in u collect cgp_2a x;

procedure cgb_s2s!-cgb(l);
   cgb_cgb!-sfl l;

procedure cgb_s2a!-gsys(u);
   % Comprehensive Groebner bases symbolic mode to algebraic mode
   % Groebner system. [u] is a GSY. Returns an AMGSYS.
   'list . for each bra in u collect cgb_s2a!-bra bra;

procedure cgb_s2a!-bra(bra);
   % Comprehensive Groebner bases symbolic mode to algebraic mode
   % branch. [u] is a BRA. Returns an AMBRANCH.
   {'list,rl_mk!*fof rl_smkn('and,bra_cd bra),
      'list . for each x in bra_system bra collect cgp_2a x};

procedure cgb_s2s!-gsys(u);
   for each bra in u collect cgb_s2s!-bra bra;

procedure cgb_s2s!-bra(bra);
   {bra_cd bra,cgb_s2s!-cgb bra_system bra};

procedure cgb_a2s!-gsys(u);
   % Comprehensive Groebner bases algebraic mode to symbolic mode
   % Groebner system. [u] is AMGSYS. Returns an FGSYS.
   begin scalar sys,w;
      sys := getrlist reval u;
      return for each bra in sys collect <<
	 w := getrlist bra;
	 bra_mk(cd_for2cd rl_simp car w,cgb_a2s!-psys cadr w,nil)
      >>
   end;

procedure cgb_a2s2!-gsys(sys);
   for each bra in sys collect
      bra_mk(car bra,cgb_a2s2!-psys cadr bra,nil);

procedure cgb_xvars!-gsys(sys,vl);
   begin scalar w;
      w := for each bra in sys join
	 bra_system bra;
      return cgb_vars(w,vl)
   end;

procedure cgb_a2s!-cd(cd);
   cd_for2cd rl_simp reval cd;

procedure cgb_a2s2!-cd(cd);
   cd;

procedure cgb_a2s!-varl(varl);
   cdr varl;

procedure cgb_a2s2!-varl(varl);
   varl;

procedure cgb_cleanup(u,v);  % Do not use reval.
   u;

procedure cgb_interface!$(fname,a2sl1,a2sl2,defl,xvfn,s2a,s2s,s,smp,argl);
   % fname is a function, the name of the procedure to be called;
   % [a2sl1] and [as2sl2] are a list of functions, called to be
   % transform algebraic arguments to symbolic arguments; [defl] is a
   % list of algebraic defualt arguments; xvfn is a procedure for
   % extracting the variables from all arguments; [s2a] is procedure
   % for transforming the symbolic return value to an algebraic mode
   % return value; [argl] is the list of arguments; [s] is a flag;
   % [smp] is a flag. Return an S-expr. If [s] is on then second stage
   % of argument processing is done with the results of the first one.
   begin scalar w,vl,nargl,oenv,ocdenv,m,c,x;
      ocdenv := cd_init();  % early setup for a2s procedures...
      if not smp then <<
	 nargl := cgb_am!-pargl(fname,a2sl1,argl,defl);
	 vl := apply(xvfn,append(nargl,{td_vars()}));
	 if null cdr vl and (w:=get(fname,'gb_wrapper)) then <<
	    cd_cleanup ocdenv;
	    return apply('gb_interface!$,append(w,{smp,argl}))
	 >>;
	 oenv := cgp_init(car vl,td_sortmode(),td_sortextension());
      >> else <<
	 w := cgb_sm!-pargl argl;
	 nargl := car w;
	 m := cadr w;
	 c := caddr w;
	 x := cadddr w;
	 vl := apply(xvfn,append(nargl,{m}));
	 if null cdr vl and (w:=get(fname,'gb_wrapper)) then <<
	    cd_cleanup ocdenv;
	    return apply('gb_interface!$,append(w,{smp,argl}))
	 >>;
	 oenv := cgp_init(car vl,c,x);
      >>;	    
      w := errorset({'cgb_interface1!$,
	 mkquote fname,mkquote a2sl2,mkquote s2a,mkquote s2s,mkquote s,
	 mkquote smp,mkquote argl, mkquote nargl,mkquote car vl,
	 mkquote cdr vl},T,!*backtrace);
      cd_cleanup ocdenv;
      cgp_cleanup oenv;
      if errorp w then
      	 rederr {"Error during ",fname};
      return car w
   end;

procedure cgb_sm!-pargl(argl);
   begin scalar nargl,m,c,x;
      nargl := reverse argl;
      x := car nargl;
      nargl := cdr nargl;
      c := car nargl;
      nargl := cdr nargl;
      m := car nargl;
      nargl := cdr nargl;
      nargl := reversip nargl;
      return {nargl,m,c,x}
   end;

procedure cgb_am!-pargl(fname,a2sl1,argl,defl);
   % process argument list for algebraic mode.
   begin integer l1,l2,l3,noa,da; scalar w,nargl,scargl,scdefl;
      l1 := length argl;
      l2 := length a2sl1;
      l3 := l2 - length defl;
      if l1 < l3 or l1 > l2 then
	 rederr {fname,"called with",l1,"arguments instead of",l3,"-",l2};
      scargl := argl;
      scdefl := defl;
      da := l2 - length defl;
      noa := 1;
      nargl := for each x in a2sl1 collect <<
	 if scargl then <<
	    w := car scargl;
	    scargl := cdr scargl
	 >> else <<
	    w := car scdefl;
	 >>;
	 if noa>da then
	    scdefl := cdr scdefl;
	 noa := noa+1;
	 apply(x,{w})
      >>;
      return nargl
   end;

procedure cgb_interface1!$(fname,a2sl2,s2a,s2s,s,smp,argl,nargl,m,p);
   begin scalar w,pl;
      pl := if s then nargl else argl;
      argl := for each x in a2sl2 collect <<
	 w := car pl;
	 pl := cdr pl;
	 apply(x,{w})
      >>;
%      w := apply(fname,nconc(argl,{m,p}));
      w := apply(fname,argl);
      w := if smp then
	 apply(s2s,{w})
      else
      	 apply(s2a,{w});
      return w
   end;

procedure cgb_greengsysf(u,m,sm,sx,theo,xvarl);
   cgb_ggsysf(u,m,sm,sx,theo,xvarl);
         
procedure cgb_gsys2green(u,theo);
   % Comprehensive Groebner bases Groebner system to gree Groebner
   % system. [u] is a GSY; [theo] is a CD. Returns a GSY, in which
   % all polynomials are colored green, i.e., the green colore head
   % part is deleted.
   for each bra in u collect
      bra_mk(bra_cd bra,cgb_cgpl2green(bra_system bra,append(theo,bra_cd bra)),
	 bra_cprl bra);

procedure cgb_cgpl2green(l,theo);  % TODO: delete green monomials in RP.
   % Comprehensive Groebner bases CGP list 2 green CGP list. [l] is a
   % list of CGP's; [theo] is a CD. Returns a list of CGP's. All CGP's
   % in the returned list are colred green, i.e., the green colored
   % head part is deleted.
   for each cgp in l collect
      cgp_green cgp;

procedure cgb_domainchk();
   % Comprehensive Groebner bases domain check. No argument. Return
   % value not defined. Raises an error if the current domain is not
   % valid for CGB computations.
   if not memq(dmode!*,'(nil)) then
      rederr bldmsg("cgb does not support domain: %w",get(dmode!*,'dname));
      
procedure cgb_vars(l,vl);
   % Comprehensive 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.
   begin scalar w,m,p;
      for each f in l do
	 w := union(w,kernels f);
      if vl then <<
      	 m := cgb_intersection(vl,w);
      	 p := setdiff(w,vl)
      >> else
	 m := w;
      return m . p
   end;

procedure cgb_varsgsys(gsys,vl);
   % Comprehensive Groebner bases variables in a Groebner system.
   % [gsys] is FGSYS; [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.
   begin scalar w,m,p;
      for each bra in gsys do
      	 for each f in bra_system bra do
	    w := union(w,kernels f);
      m := cgb_intersection(vl,w);
      p := setdiff(w,vl);
      return m . p
   end;

procedure cgb_intersection(a,b);
   % Comprehensive 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 cgb_cgb(u);
   % Comprehensive Groebner bases CGB computation. [u] is a list of
   % CGP's. Returns a list of CGP's.
   cgb_gsys2cgb cgb_gsys(u,nil,nil);

procedure cgb_gsys2cgb(u);
   % Comprehensive Groebner bases CGB to Groebner system conversion.
   % [u] is a GSY. Returns a list of CGP's.
   begin scalar cgbase;
      for each bra in u do
	 for each p in bra_system bra do
	    if not (p member cgbase) then  % TODO: cgp_member?
	       cgbase := p . cgbase;
      return cgp_lsort cgbase
   end;
   
procedure cgb_cgb!-sfl(u);
   % Comprehensive Groebner bases CGB to SF list. [u] is a list of
   % CGP's. Returns a list of SF's.
   for each p in u collect cgp_2f p;
   
smacro procedure cgb_tt(s1,s2);
   % Comprehensive Groebner bases tt. [s1] and [s2] are CGP's. Returns
   % an EV, the lcm of the leading terms of [s1] and [s2].
   ev_lcm(cgp_evlmon s1,cgp_evlmon s2);

procedure cgb_gsys(u,theo,xvarl);
   % Comprehensive Groebner bases Groebner system computation. [u] is
   % a list of CGP's; [theo] is the inital theory. Returns a GSY, the
   % Groebner system of [u].
   gsy_normalize cgb_gsys1(cgp_lsort u,theo,xvarl);

procedure cgb_ggsys(u,theo,xvarl);
   % Comprehensive Groebner bases green Groebner system computation.
   % [u] is a list of CGP's; [theo] is the initial theory. Returns a
   % GSY, the green Groebner system of [u].
   begin scalar w,!*cgbgreen,sgreen;
      if !*cgbsgreen then
	 return gsy_normalize
	    cgb_gsys2green(cgb_gsys1(cgp_lsort u,theo,xvarl),theo);
      sgreen := !*cgbgreen;
      !*cgbgreen := T;
      w := cgb_gsys(u,theo,xvarl);
      !*cgbgreen := sgreen;
      return w
   end;

procedure cgb_gsys1(u,theo,xvarl);
   % Comprehensive Groebner bases Groebner system computation
   % subroutine. [u] is a list of CGP's; [theo] is the initaila
   % theory. Returns a GSY, the Groebner system of [u].
   begin
      scalar spac,stime,p1,!*factor,!*exp,!*gcd,!*ezgcd,cgb_cd!*,!*cgbverbose;
      integer cgp_pcount!*,cgb_contcount!*,cgb_hcount!*,cgb_hzerocount!*,
	 cgb_tr1count!*,cgb_tr2count!*,cgb_tr3count!*,cgb_b4count!*,
	 cgb_strangecount!*,cgb_paircount!*,cgb_gbcount!*,cgb_contcount!*;
      !*exp := !*gcd := !*ezgcd := t;
      cgb_contcount!* := cgb_mincontred!*;
      if !*cgbstat then <<
      	 spac := gctime();
	 stime := time()
      >>;
      p1 := cgb_traverso(u,theo,xvarl);
      if !*cgbstat then <<
	 ioto_tprin2t "Statistics for GB computation:";
	 ioto_prin2t {"Time: ",time() - stime," ms plus GC time: ",
	    gctime() - spac," ms"};
	 ioto_prin2t {"H-polynomials total: ",cgb_hcount!*};
	 ioto_prin2t {"H-polynomials zero: ",cgb_hzerocount!*};
	 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!*};
%	 ioto_prin2t {"Strange reductions: ",cgb_strangecount!*}
      >>;
      return p1
   end;

procedure cgb_traverso(g0,theo,xvars);
   % Comprehensive Groebner bases Traverso. [g0] is a list of CGP's;
   % [theo] is a initial theory. Returns a GSY of [g0].  
   begin scalar bra,gsys,resl,bral;
      g0 := for each fj in g0 collect
	 cgp_simpdcont fj;
      gsys := gsy_init(g0,theo,xvars);
      while gsys do <<
	 bra := car gsys;	 
	 gsys := cdr gsys;
	 if bra_cprl bra eq 'final or null bra_cprl bra then
	    resl := bra . resl
	 else <<
	    bral := cgb_traverso1(bra,xvars);
	    gsys := nconc(bral,gsys)
	 >>
      >>;
      return resl  % TODO: reduction
   end;
	 
procedure cgb_traverso1(bra,xvars);
   % Comprehensive Groebner bases Traverso subroutine. [bra] is a BRA.
   % Returns a GSY. Performs one step in the computation of a GSY.
   begin scalar g,d,s,h,p;
      cgb_cd!* := bra_cd bra;
      g := bra_system bra;
      d := bra_cprl bra;
      if !*cgbverbose then <<
      	 ioto_prin2 {"[",cgb_paircount!*,"] "};
	 cgb_paircount!* := cgb_paircount!* #- 1
      >>;
      p := car d;
      d := cdr d;
      s := cgb_spolynomial p;
      h := cgb_normalform(s,g,xvars);
      h := cgp_simpdcont h;
      if !*cgbstat then
	 cgb_hcount!* := cgb_hcount!* #+ 1;
      if cgp_zerop h then
 	 cgb_hzerocount!* := cgb_hzerocount!* #+ 1;
      return bra_split(bra_mk(cgb_cd!*,g,d),h,xvars)
   end;

procedure cgb_spolynomial(pr);
   % Comprehensive Groebner bases S-polynomial. [pr] is a CPR. Returns
   % a CGP the S-polynomial of [pr] possibly reduced wrt. the
   % polynomials in [pr].
   begin scalar s;
      s := cgb_spolynomial1 pr;  % TODO: updcondition
      return s;  % TODO: Strange reduction 
   end;

procedure cgb_spolynomial1(pr);
   % Comprehensive Groebner bases S-polynomial subroutine. [pr] is a
   % CPR. Returns a CGP. the S-polynomial of [pr].
   begin scalar p1,p2,ep,ep1,ep2,rp1,rp2,db1,db2,x,spol;
      p1 := cpr_p1 pr;
      p2 := cpr_p2 pr;
      ep := cpr_lcm pr;
      ep1 := cgp_evlmon p1;
      ep2 := cgp_evlmon p2;
      rp1 := cgp_mred p1;
      rp2 := cgp_mred p2;
      if cgp_greenp rp1 and cgp_greenp rp2 then
 	 return cgp_zero();
      db1 := cgp_lbc p1;
      db2 := cgp_lbc p2;
      x := bc_gcd(db1,db2);
      db1 := bc_quot(db1,x);
      db2 := bc_quot(db2,x);
      spol := cgp_ilcomb(rp1,db2,ev_dif(ep,ep1),rp2,bc_neg db1,ev_dif(ep,ep2));
      if cgp_greenp spol then
	 return cgp_zero();
      return spol
   end;

procedure cgb_normalform(f,g,xvars);
   % Comprehensive Groebner bases normal form computation. [f] is a
   % CGP; [g] is a list of CGP's with red HT's. Returns a CGP $p$.
   % Depends on switch [!*cgbfullred]. $p$ is computed by
   % reducing [f] with polynomials in [g].
   begin scalar fold,c,tai,divisor;
      if null g then
	 return f;
      if cgp_greenp f then
	 return cgp_zero();
      fold := f;
      f := cgp_hpcp f;
      f := cgp_shift(f,xvars);
      c := T; while c and cgp_rp f do <<
	 divisor := cgb_searchinlist(cgp_evlmon f,g);
	 if divisor then <<
	    tai := T;
	    f := cgb_reduce(f,divisor)
	 >> else if !*cgbfullred then
	    f := cgp_shiftwhite f
	 else
	    c := nil;
	 if c then
	    f := cgp_shift(f,xvars)
      >>;
      if not tai then
	 return fold;
      return cgp_backshift f  % TODO: updccondition
   end;

procedure cgb_searchinlist(vev,g);
   % Comprehensive Groebner bases search for a polynomial in a list.
   % [vev] is a EV; [g] is a CGP. Returns a CGP $p$, such that the RP
   % of [g] is reducible wrt. $p$.
   if null g then
      nil
   else if cgb_buch!-ev_divides!?(cgp_evlmon car g,vev) then
      car g
   else
      cgb_searchinlist(vev,cdr g);

procedure cgb_buch!-ev_divides!?(vev1,vev2);
   % Comprehensive Groebner bases Buchberger exponent vector divides.
   % [vev1] and [vev2] are EV's. Returns non-[nil] if [vev1] divides
   % [vev2].
   ev_mtest!?(vev2,vev1);

procedure cgb_reduce(f,g1);
   % Comprehensive Groebner bases reduce. [f] is a CGP; [g1] is a CGP,
   % such that the RP of [f] is reducible wrt. [g1]. Returns a CGP
   % $p$. $p$ is computed by reducing [f] with [g1].
   if cgp_monp g1 then
      cgp_cancelmev(cgp_bcprod(f,cgp_lbc g1),cgp_evlmon g1)  % TODO: numberp
   else
      cgb_reduceonestep(f,g1);  % TODO: Content reduction

procedure cgb_reduceonestep(f,g);
   % Comprehensive Groebner bases reduce one step. [f] is a CGP; [g]
   % is a CGP, such that the RP of [f] is top-reducible wrt. [g].
   % Returns a CGP $p$. $p$ is computed by performing one
   % top-reduction.
   begin scalar cot,hcf,hcg,x,a,b;
      cot := ev_dif(cgp_evlmon f,cgp_evlmon g);
      hcf := cgp_lbc f;
      hcg := cgp_lbc g;
      x := bc_gcd(hcf,hcg);
      a := bc_quot(hcg,x);
      b := bc_quot(hcf,x);
      return cgp_setci(cgp_ilcombr(f,a,g,bc_neg b,cot),cgp_ci f) 
   end;  % TODO: updccondition

endmodule;   % cgb;


module cd;
% Conditions.

% DS
% <CD> ::= (...,<Atomic Formula>,...)

procedure cd_init();
   % Condition init. No argument. Return value describes the current
   % context. Depends on switch [!*cgbreal]. Sets up the environment
   % for handling conditions in the choosen context.
   (if !*cgbreal then rl_set '(ofsf) else rl_set '(acfsf)) where !*msg=nil;

procedure cd_cleanup(oc);
   % Condition clean-up. [oc] decsribes the context wich should be
   % selected. Return value unspecified.
   rl_set oc where !*msg=nil;

procedure cd_falsep(cd);
   % Condion false predicate. [cd] is a CD. Returns bool. If [t] is
   % retunred then the condion [cd] is inconsistent.
   eqcar(cd,'false);

procedure cd_siadd(atl,sicd);
   % Condion simplify add. [atl] is a list of atomic formulas; [sicd]
   % is a CD. Returns a CD, the union of [cd] and [atl].
   begin scalar w;
      if not !*cgbcdsimpl then
      	 return nconc(atl,sicd);
      w := if !*cgbgs then
	 cd_gsd(rl_smkn('and,nconc(atl,sicd)),nil)
      else
	 rl_siaddatl(atl,rl_smkn('and,sicd));
      return cd_for2cd w
   end;

procedure cd_for2cd(f);
   % Condition formula to condition. [f] is either ['false] , ['true],
   % or a conjunction of atomic formulas. Returns a CD equivalent to
   % [f]. Formula to condition.
   if f eq 'true then
      nil
   else if f eq 'false then
      '(false)
   else if cl_cxfp f then
      rl_argn f
   else
      {f};

procedure cd_surep(f,cd);
   % Condition sure predicate. [f] is an atomic formuls; [cd] is a CD.
   % If [T] is returned, then [cd] implies [f].
   begin scalar !*rlgsvb;
      return rl_surep(f,cd) where !*rlspgs=!*cgbgs,!*rlsithok=T;
   end;

procedure cd_gsd(f,cd);
   % Condition Groebner simplifier. [f] is a formula; [cd] is a
   % condition. Simplies [f] wrt. the theory [cd].
   begin scalar !*rlgsvb;
      return rl_gsd(f,cd)
   end;

procedure cd_ordp(cd1,cd2);
   % Condition order predicate. [cd1] and [cd2] are conditions sorted
   % wrt. ['cd_ordatp]. Returns bool.
   if null cd1 then
      T
   else if null cd2 then
      nil
   else if car cd1 neq car cd2 then
      cd_ordatp(car cd1,car cd2)
   else
      cd_ordp(cdr cd1,cdr cd2);

procedure cd_ordatp(a1,a2);
   % Condition order atomic formula predicate. [a1] and [a2] are
   % atomic formulas. Returns bool.
   if car a1 eq 'neq and car a2 eq 'equal then
      T
   else if car a1 eq 'equal and car a2 eq 'neq then
      nil
   else
      ordp(cadr a1,cadr a2);

endmodule;  % cd


module cpr;
% Critical pairs.

%DS
% <CPRL> ::= (...,<CPR>,...)
% <CPR> ::= (<LCM>,<P1>,<P2>,<SUGAR>);

procedure cpr_mk(f,h);
   % Critical pair make. [f], and [h] are CGP's. Returns a CPR.
   % Construct a pair from polynomials [f] and [h].
   begin scalar ttt,sf,sh;
      ttt := cgb_tt(f,h);
      sf := cgp_sugar(f) #+ ev_tdeg ev_dif(ttt,cgp_evlmon f);
      sh := cgp_sugar(h) #+ ev_tdeg ev_dif(ttt,cgp_evlmon h);
      return cpr_mk1(ttt,f,h,ev_max!#(sf,sh))
   end;

procedure cpr_mk1(lcm,p1,p2,sugar);
   % Critical pair make subroutine. [lcm] is an EV, the lcm of [evlmon
   % p1] and [evlmon p2]; [p1] and [p2] are CGP's with red HC; [sugar]
   % is a machine integer, the sugar of the S-polynomials of [p1] and
   % [p2]. Returns a CPR.
   {lcm,p1,p2,sugar};

procedure cpr_lcm(cpr);
   % Critical pair lcm. [cpr] is a critical pair. Returns the lcm part
   % of [cpr].
   car cpr;

procedure cpr_p1(cpr);
   % Critical pair p1. [cpr] is a critical pair. Returns the p1 part
   % of [cpr].
   cadr cpr;

procedure cpr_p2(cpr);
   % Critical pair p2. [cpr] is a critical pair. Returns the p2 part
   % of [cpr].
   caddr cpr;

procedure cpr_sugar(cpr);
   % Critical pair suger. [cpr] is a critical pair. Returns the sugar
   % part of [cpr].
   cadddr cpr;

procedure cpr_traverso!-pairlist(gk,g,d);
   % Critical pair Travero pair list. [gk] is a CGP with red HT; [g]
   % is a list of CGP's with red HT's; [d] is a sorted list of CPR's.
   % Returns a sorted list of CPR's the result of updating [w] with
   % critical pairs construction by combining [gk] with polynomials in
   % [g].
   begin scalar ev,r,n;
      d := cpr_traverso!-pairs!-discard1(gk,d);
      % build new pair list:
      ev := cgp_evlmon gk;
      for each p in g do
     	 if not cpr_buchcrit4t(ev,cgp_evlmon p) then <<
	    if !*cgbstat then
	       cgb_b4count!* := cgb_b4count!* #+ 1;
 	    r := ev_lcm(ev,cgp_evlmon p) . r
       	 >> else
 	    n := cpr_mk(p,gk) . n;
      n := cpr_tr2crit(n,r);
      n := cpr_listsort(n,!*cgbsloppy);
      n := cpr_tr3crit n;
      if !*cgbverbose and n then <<
	 cgb_paircount!* := cgb_paircount!* #+ length n;
	 ioto_cterpri();
	 ioto_prin2 {"(",cgb_gbcount!*,") "}
      >>;
      return cpr_listmerge(d,reversip n)
   end;

procedure cpr_tr2crit(n,r);
   % Critical pair Travero 2 criterion. [n] is a list of CPR's; [r] is
   % a list of EV's. Returns a list of CPR's. Delete equivalents to
   % coprime lcm
   for each p in n join 
      if ev_member(cpr_lcm p,r) then <<
	 if !*cgbstat then
	    cgb_tr2count!* := cgb_tr2count!* #+ 1;
	 nil
      >> else
	 {p};

procedure cpr_tr3crit(n);
   % Critical pair Travero 3 criterion. [n] is a sorted list of CPR's;
   % [r] is a list of EV's. Returns a sorted list of CPR's.
   begin scalar newn,scannewn,q;
      for each p in n do <<
	 scannewn := newn;
	 q := nil;
 	 while scannewn do
	    if ev_divides!?(cpr_lcm car scannewn,cpr_lcm p) then <<
	       q := t;
	       scannewn := nil;
	       if !*cgbstat then
	       	  cgb_tr3count!* := cgb_tr3count!* #+ 1
	    >> else
	       scannewn := cdr scannewn;
	 if not q then
	    newn := cpr_listsortin(p,newn,nil)
      >>;
      return newn
   end;

procedure cpr_traverso!-pairs!-discard1(gk,d);
   % Critical pairs Traverso pairs discard 1. [gk] is a CGP with red
   % HT; [d] is a sorted list of CPR's. Returns a list of [cpr]'s.
   % Criterion B. Delete triange relations.
   for each pij in d join
      if cpr_traverso!-trianglep(cpr_p1 pij,cpr_p2 pij,gk,cpr_lcm pij) then <<
	 if !*cgbstat then
	    cgb_tr1count!* := cgb_tr1count!* #+ 1;
	 if !*cgbverbose then
	    cgb_paircount!* := cgb_paircount!* #- 1;
	 nil
      >> else
	 {pij};

procedure cpr_traverso!-trianglep(gi,gj,gk,tij);
   % Critical pairs Traverso triangle predicate. [gi], [gj], and [gk]
   % are CGP's with red HT; [tij] is an EV.
   ev_sdivp(cgb_tt(gi,gk),tij) and ev_sdivp(cgb_tt(gj,gk),tij);

procedure cpr_buchcrit4t(e1,e2);
   % Critical pair Buchbergers criterion 4. [e1], [e2] are EV's.
   % Returns [T] if [e1] and [e2] are disjoint.
   not ev_disjointp(e1,e2);
   
procedure cpr_listsort(g,sloppy);
   % Critical pair list sort. [g] is a list of CPR's, [sloppy] is
   % bool. Returns a list of CPR'S. Destructively sorts [g]
   begin scalar gg;
      for each p in g do
 	 gg := cpr_listsortin(p,gg,sloppy);
      return gg
   end;

procedure cpr_listsortin(p,pl,sloppy);
   % Critical pair list sort into. [p] is a CPR; [pl] is a sorted list
   % of CPR's, [sloppy] is bool. Destructively sorts [p] into [pl].
   if null pl then
      {p}
   else <<
      cpr_listsortin1(p,pl,sloppy);
      pl
   >>;

procedure cpr_listsortin1(p,pl,sloppy);
   % Critical pair list sort into. [p] is a CPR; [pl] is a non-empty,
   % sorted list of CPR's; [sloppy] is bool. Destructively sorts [p]
   % into [pl].
   if not cpr_lessp(car pl,p,sloppy) then <<
      rplacd(pl,car pl . cdr pl);
      rplaca(pl,p)
   >> else if null cdr pl then
      rplacd(pl,{p})
   else
      cpr_listsortin1(p,cdr pl,sloppy);

procedure cpr_lessp(pr1,pr2,sloppy);
   % Critical pair less predicate. [p1] and [p2] are CPR's; [sloppy]
   % is bool. Returns [T] is [p1] is less than [p2]. Compare 2 pairs
   % wrt. their sugar or their lcm.
   if sloppy then
      ev_compless!?(cpr_lcm pr1,cpr_lcm pr2)
   else
      cpr_lessp1(pr1,pr2,cpr_sugar pr1 #- cpr_sugar pr2,
	 ev_comp(cpr_lcm pr1,cpr_lcm pr2));

procedure cpr_lessp1(pr1,pr2,d,q);
   % Critical pair less predicate subroutine. [p1] and [p2] are CPR's.
   % Returns [T] is [p1] is less than [p2]. Compare 2 pairs wrt. their
   % sugar or their lcm.
   if not(d #= 0) then
      d #< 0
   else if not(q #= 0) then
      q #< 0
   else
      cgp_number cpr_p2 pr1 #< cgp_number cpr_p2 pr2;

procedure cpr_listmerge(pl1,pl2);  % TODO: Rekursiv, konstruktiv !!!
   % Critical pair list merge. [pl1] and [pl2] are sorted list of
   % CPR's. Returns a sorted list of CPR's the restult of merging the
   % lists [pl1] and [pl2].
   begin scalar cpl1,cpl2;
      if null pl1 then
 	 return pl2;
      if null pl2 then
 	 return pl1;
      cpl1 := car pl1;
      cpl2 := car pl2;
      return if cpr_lessp(cpl1,cpl2,nil) then
 	 cpl1 . cpr_listmerge(cdr pl1,pl2)
      else
 	 cpl2 . cpr_listmerge(pl1,cdr pl2)
   end;

endmodule;  % cpr


module bra;

%DS
% <BRA> ::= (<CD>,<SYSTEM>,<CPRL>)

procedure bra_cd(br);
   % Branch condition. [br] is a BRA. Returns a CD, the condition part
   % of [br].
   car br;

procedure bra_system(br);
   % Branch system. [br] is a BRA. Returns a list of CGP's, the
   % system part of [br].
   cadr br;

procedure bra_cprl(br);
   % Branch critical pair list. [br] is a BRA. Returns a list of
   % CPR's, the pairs part of [br].
   caddr br;

procedure bra_mk(cd,system,cprl);
   % Branch make. [cd] is a CD; [system] is a list of CGP's with red
   % HT's; [cprl] is a list of CPR's. Returns a BRA.
   {cd,system,cprl};

procedure bra_split(bra,p,xvars);
   % Branch split. [bra] is a BRA; [p] is a CGP. Returns a GSY.
   if cgp_greenp p then
      {bra}
   else if bra_cprl bra eq 'final then
      {bra}
   else
      bra_split1(bra,cgp_enumerate cgp_condense p,xvars);

procedure bra_split1(bra,p,xvars);
   % Branch split subroutine. [bra] is a BRA; [p] is a CGP. Returns a GSY.
   for each pr in cgp_2scpl(p,bra_cd bra,xvars) collect
      bra_ext(bra,car pr,cdr pr);

procedure bra_ext(bra,cd,scp);
   % Branch extend. [bra] is a BRA; [cd] is a CD; [scp] is CGP with
   % red HT. Returns a BRA.
   begin scalar sy,d;
      if cgp_unitp scp then
	 return bra_mk(cd,{scp},'final);
      sy := for each p in bra_system bra collect cgp_cp p;  % TODO: Copy?
      d := for each pr in bra_cprl bra collect pr;  % TODO: Copy?
      if cgp_greenp scp then
	 return bra_mk(cd,sy,d);
      d := cpr_traverso!-pairlist(scp,sy,d);
      return bra_mk(cd,nconc(sy,{scp}),d)
   end;

procedure bra_ordp(b1,b2);
   % Branch order predicate. [b1] and [b2] are branches. Returns bool.
   cd_ordp(bra_cd b1,bra_cd b2);

endmodule;  % bra


module gsy;
% Groebner system.

%DS
% <GSY> ::= (...,<BRA>,...)

procedure gsy_init(l,theo,xvars);
   % Groebner system initialize. [l] is a list of CGP's. Returns a
   % GSY. We construct a case distinction wrt. to the parametric
   % coefficients in the elements of [l].
   begin scalar s;
      s := {bra_mk(theo,nil,nil)};
      for each x in l do
	 s := for each y in s join
	    bra_split(y,x,xvars);
      return s
   end;

procedure gsy_normalize(l);
   % Groebner system normalize. [l] is a GSY. Returns a GSY.
   sort(gsy_normalize1 l,'bra_ordp);

procedure gsy_normalize1(l);
   % Groebner system normalize subroutine. [l] is a GSY. Returns a GSY.
   for each bra in l collect
      bra_mk(sort(bra_cd bra,'cd_ordatp),
	 cgp_lsort for each x in bra_system bra collect cgp_normalize x,
	 bra_cprl bra);
	 
endmodule;  % gsy


module cgp;
% Comprehensive Groebner basis polynomial.

%DS
% <CGP> ::= ('cgp,<HP>,<RP>,<SUGAR>,<NUMBER>,<CI>)
% <HP> ::= <DIP>
% <RP> ::= <DIP>
% <SUGAR> ::= <Machine Integer> | nil
% <NUMBER> ::= <Machine Integer> | nil
% <CI> ::= 'unknown | 'red | 'green | 'zero | ('mixed . <WTL>) | green_colored
% <WTL> ::= (...,<EV>,...)

procedure cgp_mk(hp,rp,sugar,number,ci);
   % CGP make. [hp] and [rp] are DIP's; [sugar] and [number] are
   % machine numbers; [ci] is an S-expr.
   {'cgp,hp,rp,sugar,number,ci};

procedure cgp_hp(cgp);
   % CGP head polynomial. [cgp] is a CGP. Returns a DIP, the head
   % polynomial part of [cgp]. 
   cadr cgp;

procedure cgp_rp(cgp);
   % CGP rest polynomial. [cgp] is a CGP. Returns a DIP, the rest
   % polynomial part of [cgp].
   caddr cgp;

procedure cgp_sugar(cgp);
   % CGP sugar. [cgp] is a CGP. Returns a machine number, the sugar
   % part of [cgp].
   cadddr cgp;

procedure cgp_number(cgp);
   % CGP number. [cgp] is a CGP. Returns a machine number, the number
   % part of [cgp].
   nth(cgp,5);

procedure cgp_ci(cgp);
   % CGP number. [cgp] is a CGP. Returns an S-expr, the coloring %
   %  information of [cgp].
   nth(cgp,6);

procedure cgp_init(vars,sm,sx);
   % CGP init. [vars] is a list of variables. Returns an S-expr.
   % Initializing the DIP package.
   dip_init(vars,sm,sx);

procedure cgp_cleanup(l);
   % CGP clean-up. [l] is an S-expr returned by calling [cgp_init].
   dip_cleanup(l);

procedure cgp_lbc(u);
   % CGP leading base coefficient. [u] is a CGP. Returns the HC of the
   % rest part of [u].
   dip_lbc cgp_rp u;

procedure cgp_evlmon(u);
   % CGP exponent vector of leading monomial. [u] is a CGP. Returns
   % the HT of the rest part of [u].
   dip_evlmon cgp_rp u;

procedure cgp_zerop(u);
   % CGP zero predicate. [u] is a CGP. Returns [T] if [u] is the zero
   % polynomial.
   null cgp_hp u and null cgp_rp u;

procedure cgp_greenp(u);
   % CGP green predicate. [u] is a CGP. Returns [T] if [u] is
   % completely green colored.
   null cgp_rp u;

procedure cgp_monp(u);
   % CGP monomial predicate. [u] is a CGP. Returns [T] if [u] is a monomial.
   null cgp_hp u and dip_monp cgp_rp u;

procedure cgp_zero();
   % CGP zero. No argument. Returns the zero polynomial.
   cgp_mk(nil,nil,nil,nil,'zero); 

procedure cgp_one();
   % CGP one. No argument. Returns a CGP, the polynomial one in CGP
   % representation.
   cgp_mk(nil,dip_one(),0,nil,'red);
   
procedure cgp_tdeg(u);
   % CGP total degree. [u] is a CGP. Returns the total degree of the
   % rest polynomial of [u].
   dip_tdeg cgp_rp u;

procedure cgp_mred(cgp);
   % CGP monomial reductum. [cgp] is a CGP. Returns a CGP $p$. $p$ is
   % computed from [cgp] by deleting the HM of the rest part of [cgp].
   cgp_mk(cgp_hp cgp,dip_mred cgp_rp cgp,cgp_sugar cgp,nil,'unknown);

procedure cgp_cp(cgp);
   % CGP copy. [cgp] is a CGP. Returns a CGP, the top-level copy of
   % [cgpl
   cgp_mk(cgp_hp cgp,cgp_rp cgp,cgp_sugar cgp,cgp_number cgp,cgp_ci cgp);
      
procedure cgp_f2cgp(u);
   % CGP form to cgp. [u] is a SF. Returns a CGP. 
   cgp_mk(nil,dip_f2dip u,nil,nil,'unknown);

procedure cgp_2a(u);
   % CGP to algebraic. [u] is a CGP. Returns the AM representation of
   % [u].
   dip_2a dip_append(cgp_hp u,cgp_rp u);

procedure cgp_2f(u);
   % CGP to algebraic. [u] is a CGP. Returns the AM representation of
   % [u].
   dip_2f dip_append(cgp_hp u,cgp_rp u);

procedure cgp_enumerate(p);
   % CGP enumerate. [p] is a CGP. Returns a CGP. Sets the number of
   % [p] destructively to the next free number.
   cgp_setnumber(p,cgp_pcount!* := cgp_pcount!* #+ 1);

procedure cgp_unitp(p);
   % CGP unit predicate. [p] is a CGP with red HT. Returns [T] if [p]
   % is a unit.
   cgp_rp p and ev_zero!? cgp_evlmon p;

procedure cgp_setnumber(p,n);
   % CGP set number. [p] is a CGP; [n] is a machine number. Returns a
   % CGP. Sets the number of [p] destructively to [n].
   <<
      nth(p,5) := n;
      p
   >>;

procedure cgp_setsugar(p,s);
   % CGP set sugar. [p] is a CGP; [s] is a machine number. Returns a
   % CGP. Sets the sugar of [p] destructively to [s].
   <<
      nth(p,4) := s;
      p
   >>;

procedure cgp_setci(p,tg);
   % CGP set coloring information. [p] is a CGP; [tg] is an S-expr.
   % Returns a CGP. Sets the coloring information of [p] destructively
   % to [s].
   <<
      nth(p,6) := tg;
      p
   >>;

procedure cgp_condense(p);
   % CGP condense. [p] is a CGP. Returns a CGP. Condenses both the
   % head and the rest polynomial of [p].
   <<
      dip_condense cgp_hp p;
      dip_condense cgp_rp p;
      p
   >>;

procedure cgp_2scpl(p,cd,xvars);
   % CGP to strong cpl. [p] is a CGP; [cd] is a CD. Returns a list of
   % pairs $(...,(\gamma . p'),...)$, where $\gamma$ is a condition
   % and $p'$ is a CGP with red HC.
   if !*cgbgen and null xvars then
      cgp_2scpl!-gen(p,cd)
   else
      cgp_2scpl1(p,cd,xvars);

procedure cgp_2scpl1(p,cd,xvars);
   % CGP to strong cpl subroutine. [p] is a CGP; [cd] is a CD. Returns
   % a list of pairs $(...,(\gamma . p'),...)$, where $\gamma$ is a
   % condition and $p'$ is a CGP with red HC.
   begin scalar hp,rp,s,n,hc,ht,l,ncdeq,ncdneq;
      hp := cgp_hp p;
      if !*cgbgreen and hp then
	 rederr {"cgp_2scpl1: Non empty hp",p};
      rp := cgp_rp p;      
      s := cgp_sugar p;
      n := cgp_number p;
      while rp do <<
	 hc := dip_lbc rp;
	 ht := dip_evlmon rp;
	 ncdeq := ncdneq := nil;
	 if cd_surep(bc_mkat('neq,hc),cd) or
	    eqcar(ncdeq := cd_siadd({bc_mkat('equal,hc)},cd),'false)
	 then <<
	    l := (cd . cgp_mk(hp,rp,s or dip_tdeg rp,n,'red)) . l;
	    hc := 'break;
	    rp := nil
	 >> else if !*cgbgen and null intersection(xvars,bc_vars hc) then <<
	    ncdneq := cd_siadd({bc_mkat('neq,hc)},cd);
	    l := (ncdneq . cgp_mk(hp,rp,s or dip_tdeg rp,n,'red)) . l;
	    hc := 'break;
	    rp := nil	    
	 >> else <<
	    if not (cd_surep(bc_mkat('equal,hc),cd) or
	       eqcar(ncdneq := cd_siadd({bc_mkat('neq,hc)},cd),'false))
	    then <<
	       ncdneq := ncdneq or cd_siadd({bc_mkat('neq,hc)},cd);
	       ncdeq := ncdeq or cd_siadd({bc_mkat('equal,hc)},cd);
	       l := (ncdneq . cgp_mk(hp,rp,s or dip_tdeg rp,n,'red)) . l;
	       cd := ncdeq;
	    >>;	 
	    rp := dip_mred rp;
	    if not(!*cgbgreen) then
	       hp := dip_appendmon(hp,hc,ht);
	 >>
      >>;	 
      if hc neq 'break then
 	 l := (cd . cgp_zero()) . l;
      return reversip l
   end;

procedure cgp_2scpl!-gen(p,cd);
   % CGP to strong cpl generic case. [p] is a CGP; [cd] is a CD. Returns
   % a list of one pair $((\gamma . p'))$, where $\gamma$ is a
   % condition and $p'$ is a CGP with red HC.
   begin scalar hp,rp;
      hp := cgp_hp p;
      rp := cgp_rp p;      
      if null rp then
	 return {cd . cgp_zero()};
      cd := cd_siadd({bc_mkat('neq,dip_lbc rp)},cd);
      return {cd . cgp_mk(hp,rp,cgp_sugar p or dip_tdeg rp,cgp_number p,'red)}
   end;

procedure cgp_ilcomb(p1,c1,t1,p2,c2,t2);
   % CGP integer linear combination. [p1], [p2] are CGP's; [c1], [c2]
   % are BC's; [t1], [t2] are EV's. Returns a CGP. Computes
   % $p1*c1^t1+p2*c2^t2$.
   begin scalar hp,rp,s;
      hp := dip_ilcomb(cgp_hp p1,c1,t1,cgp_hp p2,c2,t2);
      rp := dip_ilcomb(cgp_rp p1,c1,t1,cgp_rp p2,c2,t2);
      s := ev_max!#(cgp_sugar p1 #+ ev_tdeg t1,cgp_sugar p2 #+ ev_tdeg t2);
      return cgp_mk(hp,rp,s,nil,'unknown)  % TODO: Summe ?????
   end;

procedure cgp_ilcombr(p1,c1,p2,c2,t2);
   % CGP integer linear combination for reduction. [p1], [p2] are
   % CGP's; [c1], [c2] are BC's; [t2] is a EV's. Returns a CGP.
   % Computes $p1*c1+p2*c2^t2$.
   begin scalar hp,rp,s;
      hp := dip_ilcombr(cgp_hp p1,c1,cgp_hp p2,c2,t2);
      rp := dip_ilcombr(cgp_rp p1,c1,cgp_rp p2,c2,t2);
      s := ev_max!#(cgp_sugar p1,cgp_sugar p2 #+ ev_tdeg t2);
      return cgp_mk(hp,rp,s,nil,'unknown)
   end;

procedure cgp_hpcp(cgp);
   % CGP head polynomial copy. [cgp] is a CGP. Returns a CGP, in which
   % the head polynomial is copied.
   cgp_mk(dip_cp cgp_hp cgp,cgp_rp cgp,cgp_sugar cgp,
      cgp_number cgp,cgp_ci cgp);

procedure cgp_shift(p,xvars);
   % CGP shift. [p] is a CGP, which is neither zero nor green. Returns
   % a [CGP]. Shifts all leading green monomials from the rest part
   % into the head part.
   if !*cgbgen and null xvars then
      cgp_shift!-gen p
   else
      cgp_shift1(p,xvars);

procedure cgp_shift1(p,xvars);
   % CGP shift subroutine. [p] is a CGP, which is neither zero nor
   % green. Returns a [CGP]. Shifts all leading green monomials from
   % the rest part into the head part.
   begin scalar hp,rp,ht,hc,c;
      hp := cgp_hp p;
      rp := cgp_rp p;
      c := T;
      while c and rp do <<
	 ht := dip_evlmon rp;
	 hc := dip_lbc rp;
	 if cd_surep(bc_mkat('equal,hc),cgb_cd!*) then <<
	    if not(!*cgbgreen) then
	       hp := dip_nconcmon(hp,hc,ht);
	    rp := dip_mred rp
	 >> else
	    c := nil
      >>;
      if null rp and idp cgp_ci p then
	 return cgp_zero();
      return cgp_mk(hp,rp,cgp_sugar p,cgp_number p,cgp_ci p)
   end;

procedure cgp_shift!-gen(p);
   % CGP shift generic case. [p] is a CGP, which is neither zero nor
   % green. Returns a [CGP]. Shifts all leading green monomials from
   % the rest part into the head part, i.e. we do nothing because
   % there are no green BC's.
   p;

procedure cgp_shiftwhite(p);
   % CGP shift white. [p] is a CGP, which is neither zero nor green.
   % Returns a [CGP]. Shifts the leading white monomials from the rest
   % part into the head part and set the wtl accordingly.
   begin scalar nhp,nci;
      nhp := dip_nconcmon(cgp_hp p,cgp_lbc p,cgp_evlmon p);
      nci := cgp_ci p;
      nci := 'mixed . (cgp_evlmon p . if idp nci then nil else cdr nci);      
      return cgp_mk(nhp,dip_mred cgp_rp p,cgp_sugar p,cgp_number p,nci)
   end;

procedure cgp_backshift(p);
   % CGP back shift. [p] is a CGP. Returns a CGP. Shifts all white
   % monomials from the head part into the rest part using the wtl.
   begin scalar ci;
      ci := cgp_ci p;
      if not pairp ci or pairp ci and null cdr ci then
	 return p;
      if cgp_rp p then
	 rederr "cgp_backshift: Rest polynomial must be zero";
      return cgp_backshift1 p
   end;

procedure cgp_backshift1(p);
   % CGP back shift subroutine. [p] is a CGP. Returns a CGP. Shifts
   % all white monomials from the head part into the rest part using
   % the wtl.
   begin scalar hp,wtl,nhp;
      hp := cgp_hp p;
      wtl := cdr cgp_ci p;
      % TODO: Update condition
      while hp and not ev_member(dip_evlmon hp,wtl) do <<  % TODO: Destructive?
	 nhp := dip_nconcmon(nhp,dip_lbc hp,dip_evlmon hp);
	 hp := dip_mred hp
      >>;
      if hp then
	 return cgp_mk(nhp,hp,cgp_sugar p,cgp_number p,'unknown);      
      return cgp_zero()
   end;

procedure cgp_cancelmev(p,ev);
   % CGP cancel monomoial ev's. [p] is a CGP; [ev] is an EV. Returns a
   % CGP. Cancels all monomials in f which are multiples of [ev].
   cgp_mk(cgp_hp p,dip_cancelmev(cgp_rp p,ev),
      cgp_sugar p,cgp_number p,cgp_ci p);

procedure cgp_bcquot(p,c);
   % CGP base coefficient procuct. [p] is a CGP; [c] is a BC. Returns
   % a CGP. Computes $(1/[c])[p]$.
   cgp_mk(dip_bcquot(cgp_hp p,c),dip_bcquot(cgp_rp p,c),
      cgp_sugar p,cgp_number p,cgp_ci p);

procedure cgp_bcprod(p,c);
   % CGP base coefficient procuct. [p] is a CGP; [c] is a BC. Returns
   % a CGP. Computes $[c][p]$.
   cgp_mk(dip_bcprod(cgp_hp p,c),dip_bcprod(cgp_rp p,c),
      cgp_sugar p,cgp_number p,cgp_ci p);

procedure cgp_simpdcont(p);
   % CGP simplify domain content. [p] is a CGP. Returns a CGP $p'$
   % such that $p'$ is primitive as a multivariate polynomial over Z
   % and there is an integer $c$ such that $[p]=cp'$.
   begin scalar c;
      if cgp_zerop p then
 	 return p;
      c := cgp_dcont p;
      if bc_minus!? cgp_rlbc p then
	 c := bc_neg c;
      return cgp_mk(dip_bcquot(cgp_hp p,c),dip_bcquot(cgp_rp p,c),
      	 cgp_sugar p,cgp_number p,cgp_ci p)
   end;

procedure cgp_rlbc(p);
   % CGP real leading base coefficient. [p] is a CGP. Returns a BC,
   % the coefficient of the largest term in both the head polynomial
   % and the rest polynomial part.
   if cgp_zerop p then
      bc_fd 0
   else if cgp_hp p then 
      dip_lbc cgp_hp p
   else
      cgp_lbc p;

procedure cgp_dcont(p);
   % CGP domain content. [p] is a CGP. Returns a BC, the domain
   % content of [p], i.e. the content of [p] considered as an
   % multivariate polynomial over Z.
   begin scalar c;
      c := dip_dcont cgp_hp p;
      if bc_one!? c then
	 return c;
      return dip_dcont1(cgp_rp p,c)
   end;

procedure cgp_normalize(u);
   % CGP normalize. [u] is a CGP. Returns a unique representation of
   % [u] as a CGP.
   cgp_mk(nil,dip_append(cgp_hp u,cgp_rp u),nil,nil,'unknown);

procedure cgp_green(u);
   % CGP green. [u] is A CGP. Returns a green CGP, i.e. a CGP in which
   % the green head part is cancelled.
   cgp_mk(nil,cgp_rp u,nil,nil,'green_colored);
   
procedure cgp_lsort(pl);
   % CGP list sort. pl is a list of CGP's. Returns a list of CGP's.
   sort(pl,function cgp_comp);

procedure cgp_comp(p1,p2);
   dip_comp(cgp_rp p1,cgp_rp p2);

endmodule;  % cgp

end;  % of file