module buchbg;% Central Groebner base code: Buchberger algorithm.
% Authors: H. Melenk,H. M. Moeller,W. Neun
% ZIB Berlin,August 2000
flag('(groebrestriction groebresmax gvarslast groebmonfac
groebprotfile glterms),'share);
groebrestriction:=nil;groebresmax:=300;groebmonfac:=1;
groecontcount!*:=10;
!*gsugar:=t;
!*groelterms:=t;
!*groebfullreduction:=t;
!*groebdivide:=t;
switch groebopt,trgroeb,trgroebs,trgroeb1,
trgroebr,groebfullreduction,groebstat,groebprot;
% Variables for counting and numbering.
% Option'groebopt'"optimizes" the given input
% polynomial set(variable ordering).
% option'trgroeb'Prints intermediate
% results on the output file.
% option'trgroeb1'Prints internal representation
% of critical pair list d.
% option'trgroebs'Prints's'- polynomials
% on the output file.
% option'trgroebr'Prints(intermediate)results and
% computation statistics.
% option'groebstat'The statistics are printed.
%
% option'groebrm'Multiplicities of factors in h-polynomials are reduced
% to simple factors .
% option'groebdivide'The algorithm avoids all divisions(only for modular
% calculation), if this switch is set off.
% option'groebprot'Write a protocol to the variable "groebprotfile".
symbolic procedure buchvevdivides!?(vev1,vev2);
% Test : vev1 divides vev2 ? for exponent vectors vev1,vev2.
vevmtest!?(vev2,vev1)and
(null gmodule!* or gevcompatible1(vev1,vev2,gmodule!*));
symbolic procedure gevcompatible1(v1,v2,g);
% Test whether'v1'and'v2'belong to the same vector column.
if null g then t
else if null v1 then(null v2 or gevcompatible1('(0), v2,g))
else if null v2 then gevcompatible1(v1,'(0), g)else
(car g=0 or car v1=car v2)
and gevcompatible1(cdr v1,cdr v2,cdr g);
symbolic procedure gcompatible(f,h);
null gmodule!* or gevcompatible1(vdpevlmon f,vdpevlmon h,gmodule!*);
symbolic procedure groebmakepair(f,h);
% Construct a pair from polynomials'f'and'h'.
begin scalar ttt,sf,sh;
ttt:=tt(f,h);
return if !*gsugar then
<< sf:=gsugar(f)#+ vevtdeg vevdif(ttt,vdpevlmon f);
sh:=gsugar(h)#+ vevtdeg vevdif(ttt,vdpevlmon h);
{ttt,f,h,max(sf,sh)}>>
else{ttt,f,h}end;
% The 1-polynomial will be constructed at run time
% because the length of the vev is not known in advance.
fluid'(vdpone!*);
symbolic procedure vdponepol;
% Construct the polynomial=1.
vdpone!*:=vdpfmon(a2vbc 1,vevzero());
symbolic procedure groebner2(p,r);
% Setup all global variables for the Buchberger algorithm,
% printing of statistics.
begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*,
pairsdone!*,factorlevel!*,groesfactors!*,!*gcd;
factortime!*:=0;groetime!*:=time();
vdponepol();% we construct dynamically
hcount!*:=0;mcount!*:=0;fcount!*:=0;
bcount!*:=0;b4count!*:=0;hzerocount!*:=0;
basecount!*:=0;!*gcd:=t;glterms:={'list};
groecontcount!*:=10;
if !*trgroeb then
<<prin2"Groebner Calculation starting ";terprit 2;
prin2" groebopt: ";print !*groebopt>>;
spac:=gctime();
p1:= if !*groebfac or null !*gsugar
then groebbasein(p,!*groebfac,r)where !*gsugar=nil
else gtraverso(p,nil,nil);
if !*trgroeb or !*trgroebr or !*groebstat then
<<spac1:=gctime()-spac;terpri();
prin2t"statistics for GROEBNER calculation";
prin2t"===================================";
prin2" total computing time(including gc): ";
prin2(( tim1:=time())- groetime!*);
prin2t" milliseconds ";
if factortime!* neq 0 then
<<prin2"(time spent in FACTOR(excl. gc): ";
prin2 factortime!*;prin2t " milliseconds)">>;
prin2"(time spent for garbage collection: ";
prin2 spac1;prin2t " milliseconds)";terprit 1;
prin2"H-polynomials total: ";prin2t hcount!*;
prin2"H-polynomials zero : ";prin2t hzerocount!*;
prin2"Crit M hits: ";prin2t mcount!*;
prin2"Crit F hits: ";prin2t fcount!*;
prin2"Crit B hits: ";prin2t bcount!*;
prin2"Crit B4 hits: ";prin2t b4count!*>>;return p1 end;
smacro procedure testabort h;
vdpmember(h,abort1)or
'cancel=( abort2:=groebtestabort(h,abort2));
symbolic procedure groebenumerate f;
%'f'is a temporary result. Prepare it for medium range storage
% and ssign a number.
if vdpzero!? f then f else
<<vdpcondense f;
if not vdpgetprop(f,'number)then
<<f:=vdpputprop(f,'number,(pcount!*:=pcount!* #+ 1));
if !*groebprot then
<<groebprotsetq(mkid('poly,pcount!*),'candidate);
vdpputprop(f,'groebprot,mkid('poly,pcount!*))
>> >>;f>>;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Buchberger's Algorithm
%
% INPUT : G0={f1,...,fr} set of nonzero polynomials.
% OUTPUT: groebner basis(list of nonzero polynomials).
%
% internal variables:
%
% problems list of problems to be processed. Problems is non nil,
% if the inital problem was split by a successful factoring.
% results Collection of results from problems.
% g Basis under construction.
% g1 Local pointer to g.
% d List of critical pairs during algorithm.
% d1,d2 Local lists of pairs during update of d.
% f,fj Polynomials.
% p,p1,p2 Pairs.
% s,h Polynomials in the central part of the algorithm
% (the "s-poly" and the "h-poly" selon Buchberger).
% g99 Set of polynomials used for reduction.
% abort1 List of polynomials in factorization context.
% S calculation branch can be cancelled if one of
% these polys is detected.
% abort2 List of sets of polynomials. If a new h polynomial
% is calculated,it should be removed from the sets.
% If one set becomes null,the set restriction is
% fulfilled and the branch can be cancelled.
symbolic procedure groebbasein(g0,fact,abort1);
begin scalar abort2,d,d1,d2,g,g1,g99,h,hlist,lasth,lv,p,problems,
p1,results,s,x;integer gvbc,probcount!*;
groebabort!*:=abort1;lv:=length vdpvars!*;
for each p in g0 do if vdpzero!? p then g0:=delete(p,g0);
if !*groebprereduce then g0:=groebprereduce g0;
x:=for each fj in g0 collect
<<groebsavelterm fj;gsetsugar(vdpenumerate vdpsimpcont fj,nil)>>;
if !*groebprot then
for each f in x do
<<groebprotsetq(mkid('poly,h:=vdpnumber f),vdp2a f);
vdpputprop(f,'groebprot,mkid('poly,h))>>;
g0:=x;
% Establish the initial problem
problems:={{nil,nil,nil,g0,abort1,nil,nil,vbccurrentmode!*,nil,nil}};
!*trgroeb and groebmess1(g,d);
go to macroloop;
macroloop:
while problems and gvbc < groebresmax do
begin
% Pick up next problem
x:=car problems;d:=car x;g:=cadr x;
% g99:=groeblistreconstruct caddr x;
g99:=vdplsort caddr x;g0:=cadddr x;abort1:=nth(x,5);
abort2:=nth(x,6);pairsdone!*:=nth(x,7);h:=nth(x,8);
% vbccurrentmode!*
factorlevel!*:=nth(x,9);groesfactors!*:=nth(x,10);
problems:=cdr problems;
g0:=% Sort'g0',but keep factor in first position
if factorlevel!* and g0 and cdr g0 then car g0.vdplsort cdr g0
else vdplsort g0;x:=nil;lasth:=nil;
!*trgroeb and groebmess23(g0,abort1,abort2);
while d or g0 do
begin
if groebfasttest(g0,g,d,g99)then go to stop;
!*trgroeb and groebmess50 g;
if g0 then<<h:=car g0;g0:=cdr g0;gsetsugar(h,nil);
groebsavelterm h;p:={nil,h,h}>>else
<<p:=car d;d:=delete(p,d);
s:=groebspolynom(cadr p,caddr p);
if fact then
pairsdone!*:=(vdpnumber cadr p.vdpnumber caddr p).pairsdone!*;
!*trgroeb and groebmess3(p,s);
h:=groebnormalform0(s,g99,'tree,fact);groebsavelterm h;
h:=groebsimpcontnormalform h;
if vdpzero!? h then !*trgroeb and groebmess4(p,d);
% Test for possible chains
if not vdpzero!? h then % only for real h's
<<s:=groebchain(h,cadr p,g99);
if s=h then h:=groebchain(h,caddr p,g99);
if secondvalue!* then g:=delete(secondvalue!*,g)>> >>;
if vdpzero!? h then go to bott;
if vevzero!? vdpevlmon h then % base 1 found
<<!*trgroeb and groebmess5(p,h);go to stop>>;
if testabort(h)then
<<!*trgroeb and groebmess19(h,abort1,abort2);go to stop>>;
s:= nil;
% Look for implicit or explicit factorization
hlist:=nil;
if groebrestriction!* then hlist:=groebtestrestriction(h,abort1);
if not hlist and fact then hlist:=groebfactorize(h,abort1,g,g99);
if hlist='zero then go to bott;
if groefeedback!* then g0:=append(groefeedback!*,g0);
groefeedback!*:=nil;
% Factorisation found but only one factor survived
if hlist and length hlist=2 then<<h:=car cadr hlist;hlist:=nil>>;
if hlist then
<<if hlist neq'cancel then
problems:= groebencapsulate(hlist,d,g0,g,g99,abort1,abort2,problems,fact);
go to stop>>;
%'h'polynomial is accepted now
h:=groebenumerate h;!*trgroeb and groebmess5(p,h);
% Construct new critical pairs
d1:=nil;
!*trgroeb and groebmess50(g);
for each f in g do
if(car p or % That means "not an input polynomial"
not member(vdpnumber h.vdpnumber f,pairsdone!*)
)and gcompatible(f,h)then
<<d1:=groebcplistsortin(groebmakepair(f,h),d1);
if tt(f,h)=vdpevlmon(f)then
<<g:=delete(f,g);
!*trgroeb and groebmess2 f>> >>;
!*trgroeb and groebmess51 d1;
d2:=nil;
while d1 do
<<d1:=groebinvokecritf d1;p1:=car d1;d1:=cdr d1;
d2:=groebinvokecritbuch4(p1,d2);
d1:=groebinvokecritm(p1,d1)>>;
d:=groebinvokecritb(h,d);
d:=groebcplistmerge(d,d2);
% Monomials and binomials
if vdplength h < 3 and car p then
<<g:=groebsecondaryreduction(h,g,g99,d,nil,t);
if g='abort then go to stop;g99:=secondvalue!*;
d:=thirdvalue!*>>;
g:=h.g;lasth:=h;
g99:=groeblistadd(h,g99);
!*trgroeb and groebmess8(g,d);
go to bott;
stop: d:=g:=g0:=nil;
bott:;
end;
g:=vdplsort g;% Such that g descending
x:=groebbasein2(g,g99,problems,abort1,abort2,fact);
g1:=car x;problems:=cdr x;
if g1 then<<results:=g1.results;gvbc:=gvbc+1>>;
!*trgroeb and groebmess13(g1,problems)end;
if gvbc >= groebresmax then
lpriw("########","warning: GROEBRESMAX limit reached");
return groebbasein3 results end;
symbolic procedure groebfasttest(g0,g,d,g99);
<<g:=g0:=d:=g99:=nil;nil>>;% A hook for special techniques
symbolic procedure groebbasein2(g,g99,problems,abort1,abort2,fact);
% Final reduction for a base'g': reduce each polynomial with the
% other members;here again support of factorization.
begin scalar f,g1,!*groebfullreduction,!*groebheufact,!*gsugar,% Saving value
h,hlist,x;integer cnt;
!*groebfullreduction:=t;g1:=nil;
while g do
<<h:=car g;g:=cdr g;
if !*groebprot then
groebprotsetq('candidate,mkid('poly,vdpnumber h));
h:=groebnormalform0(h,g,'sort,nil);
f:=groebsimpcontnormalform h;
!*trgroeb and groebmess26(h,f);
if !*groebprot then groebprotsetq({'gb,cnt:=cnt+1},'candidate);
if vdpone!? f then<<g1:=g:=nil>>;% Base{1} found
% very late now
if fact and not vdpzero!? f then
<<hlist:=groebfactorize(f,abort1,nil,nil);
if not null hlist then
<< % lift structure
hlist:=for each x in cdr hlist collect car x;
% discard superfluous factors
for each h in hlist do
if vdpmember(h,abort1)then
<<hlist:=delete(h,hlist);!*trgroeb and groebmess19(h,abort1,abort2)>>;
% Generate new subproblems
x:=0;
for each h in hlist do
<<hlist:=delete(h,hlist);
h:=groebenumerate h;
problems:=
{nil, % null D
append(g1,g), % base
g99, % g99
{h}, % g0
append(hlist,abort1),
abort2,
pairsdone!*,
vbccurrentmode!*,
(x:=x+1).factorlevel!*,
groesfactors!*}. problems>>;
g:=g1:=nil;% Cancel actual final reduction
f:=nil >> >>;
if f and vdpevlmon h neq vdpevlmon f then
<<g:=vdplsort append(g,f.g1);g1:=nil>>else
if f and not vdpzero!? f then g1:=append(g1,{f})>>;
return g1.problems end;
symbolic procedure groebbasein3 results;
% Final postprocessing: remove multiple bases from the result.
begin scalar x,g,f,p1,p2;
x:=nil;g:=results;p1:=p2:=0;
while results do
<<if vdpone!? car car results % Exclude multiple{1}
then p1:=p1+1 % count ones
else
<<f:=for each p in car results % Delete props for member
collect vdpremallprops p;
if member(f,x) % each base only once
then p2:=p2+1 % count multiples
else if not groebabortid( f,groebabort!*)
then x:=f.x;
results:=cdr results>> >>;
results:=if null x then{{vdpone!*}}else x;
return results end;
fluid'(!*vbccompress);
symbolic procedure groebchain(h,f,g99);
% Test if a chain from h-plynomials can be computed from the'h'.
begin scalar count,found,h1,h2,h3;
secondvalue!*:=nil;
return h;% Erst einmal.
if not buchvevdivides!?(vdpevlmon h,vdpevlmon f)
then return h;
h2:=h;h1:=f;found:=t;count:=0;
while found do
<<h3:=groebspolynom(h1,h2);
h3:=groebnormalform0(h3,g99,'tree,t);
h3:=vdpsimpcont h3;
if vdpzero!? h3 then
<<found:=nil;
prin2t "chain---------------------------";
vdpprint h1;vdpprint h2;vdpprint h3;
secondvalue!*:=f;count:=9999>>
else if vdpone!? h3 then
<<found:=nil;
prin2t "chain---------------------------";
vdpprint h1;vdpprint h2;vdpprint h3;
h2:=h3;count:=9999>>
else if buchvevdivides!?(vdpevlmon h3,vdpevlmon h2)then
<<found:=t;
prin2t "chain---------------------------";
vdpprint h1;vdpprint h2;vdpprint h3;
h1:=h2;h2:=h3;count:=count+1>>
else found:=nil>>;
return if count > 0 then
<<prin2 "CHAIN :";prin2t count;h2>>else h end;
symbolic procedure groebencapsulate(hlist,d,g0,g,g99,
abort1,abort2,problems,fact);
%'hlist'is a factorized h-poly. This procedure has the job to
% form new problems from hlist and to add them to problems.
% Result is problems.
% Standard procedure: only creation of subproblems.
begin scalar factl, % List of factorizations under way.
u,y,z;integer fc;
if length vdpvars!* > 10 or car hlist neq'factor then
return groebencapsulatehardcase(hlist,d,g0,g,g99,
abort1,abort2,problems,fact);
% Encapsulate for each factor.
factl:=groebrecfactl{hlist};
!*trgroeb and groebmess22(factl,abort1,abort2);
for each x in reverse factl do
<<y:=append(car x,g0);
z:=vdpunion(cadr x,abort1);
u:=append(caddr x,abort2);
problems:={d,g,
g, % future g99
y, % as new problem
z, % abort1
u, % abort2
pairsdone!*, % pairsdone!*
vbccurrentmode!*,
(fc:=fc+1).factorlevel!*,
groesfactors!*
}. problems>>;
return problems end;
symbolic procedure groebencapsulatehardcase(hlist,d,g0,g,g99,
abort1,abort2,problems,fact);
%'hlist'is a factorized h-poly. This procedure has the job to
% form new problems from hlist and to add them to problems.
% Result is problems.
% First the procedure tries to compute new h-polynomials from the
% remaining pairs which are not affected by the factors in hlist.
% Purpose is to find further factorizations and to do calculations
% in common for all factors in order to shorten the separate later
% branches.
begin scalar factl, % List of factorizations under way.
factr, % Variables under factorization.
break,d1,d2,f,fl1,gc,h,p,pd,p1,s,u,y,z;
integer fc;
factl:={hlist};factr:=vdpspace car cadr hlist;
for each x in cdr hlist do
for each p in x do factr:=vevunion(factr,vdpspace p);
% ITER:
% Now process additional pairs.
while d or g0 do
begin
break:=nil;
if g0 then
<< % Next poly from input.
s:=car g0;g0:=cdr g0;p:={nil,s,s}>>
else
<< % Next poly fropm pairs.
p:=car d;d:=delete(p,d);
if not vdporthspacep(car p,factr)then
s:=nil else
<<s:=groebspolynom(cadr p,caddr p);
!*trgroeb and groebmess3(p,s)>> >>;
if null s or not vdporthspacep(vdpevlmon s,factr)then
<< % Throw away s polynomial .
f:=cadr p;
if not vdpmember3(f,g0,g,gc)then gc:=f.gc;
f:=caddr p;
if car p and not vdpmember3(f,g0,g,gc)
then gc:=f.gc;go to bott>>;
h:=groebnormalform(s,g99,'tree);
if vdpzero!? h and car p then !*trgroeb and groebmess4(p,d);
if not vdporthspacep(vdpevlmon h,factr)then
<< % Throw away h-polynomial.
f:=cadr p;
if not vdpmember3(f,g0,g,gc)then gc:=f.gc;
f:=caddr p;
if car p and not vdpmember3(f,g0,g,gc)then gc:=f.gc;
go to bott>>;
%%% if car p then
%%% pairsdone!*:=(vdpnumber cadr p.vdpnumber caddr p).pairsdone!*;
if vdpzero!? h then go to bott;
if vevzero!? vdpevlmon h then % Base 1 found.
go to stop;
h:=groebsimpcontnormalform h;% Coefficients normalized.
if testabort h then
<<!*trgroeb and groebmess19(h,abort1,abort2);
go to stop>>;
s:=nil;hlist:=nil;
if groebrestriction!* then hlist:=groebtestrestriction(h,abort1);
if hlist='cancel then go to stop;
if not hlist and fact then
hlist:=groebfactorize(h,abort1,g,g99);
if groefeedback!* then g0:=append(groefeedback!*,g0);
groefeedback!*:=nil;
if hlist and length hlist=2 then
<<h:=car cadr hlist;hlist:=nil>>;
if hlist then
<<for each x in cdr hlist do
for each h in x do factr:=vevunion(factr,vdpspace h);
factl:=hlist.factl;% Add to factors.
go to bott>>;
h:=groebenumerate h; % Ready now.
!*trgroeb and groebmess5(p,h);
% Construct new critical pairs.
d1:=nil;
for each f in g do
if tt(f,h)=vdpevlmon(f)and gcompatible(f,h)then
<<g:=delete(f,g);
d1:=groebcplistsortin(groebmakepair(f,h),d1);
!*trgroeb and groebmess2 f>>;
!*trgroeb and groebmess51 d1;
d2:=nil;
while d1 do
<<d1:=groebinvokecritf d1;p1:=car d1;d1:=cdr d1;
d2:=groebinvokecritbuch4(p1,d2);
d1:=groebinvokecritm(p1,d1)>>;
d:=groebinvokecritb(h,d);d:=groebcplistmerge(d,d2);
if vdplength h < 3 then
<<g:=groebsecondaryreduction(h,g,g99,d,gc,t);
if g='abort then go to stop;g99:=secondvalue!*;
d:=thirdvalue!*;gc:=fourthvalue!*>>;
g:=h.g;
g99:=groeblistadd(h,g99);
!*trgroeb and groebmess8(g,d);
go to bott;
stop: d:=g:=g0:=gc:=factl:=nil;
bott: end;%ITER
% Now collect all relvevant polys.
g0:=vdpunion(g0,vdpunion(g,gc));
% Encapsulate for each factor.
if factl then
<<factl:=groebrecfactl factl;
!*trgroeb and groebmess22(factl,abort1,abort2)>>;
for each x in reverse factl do
<<fl1:=(fc:=fc+1).factorlevel!*;
break:= nil;y:=append(car x,g0);
z:=vdpunion(cadr x,abort1);
u:=append(caddr x,abort2);
if vdpmember(vdpone!*,y)then break:=vdpone!*;
% Inspect the unreduced list first .
if not break then for each p in z do
if vdpmember(p,y)then break:=p;
% Now prepare the reduced list.
if not break then
<<y:=append(car x,groebreducefromfactors(g0,car x));
pd:=secondvalue!*;
if vdpmember(vdpone!*,y)then break:=vdpone!* else
for each p in z do if vdpmember(p,y)then break:=p;
pd:=subla(pd,pairsdone!*)>>;
if not break then
problems:={
nil, % new d
nil, % new g
nil, % future g99
y, % as new problem
z, % abort1
u, % abort2
nil, % pairsdone!*
vbccurrentmode!*,
fl1, % factorlevel!*,
groesfactors!* % factor db
}.problems else !*trgroeb and groebmess19a(break,fl1)>>;
return problems end;
symbolic procedure groebrecfactl(factl);
% Factl is a list of factorizations:a list of lists of vdps
% generate calculation pathes from them.
begin scalar rf,res,type;
if null factl then return{{nil,nil,nil}};
rf:=groebrecfactl(cdr factl);
factl:=car factl;
type:=car factl;% FACTOR or RESTRICT
factl:=cdr factl;
while factl do
<<for each y in rf do
if vdpdisjoint!?(car factl,cadr y)then
res:={vdpunion(car factl,car y),
(if type='factor then
append(for each x in cdr factl collect car x, cadr y)
else cadr y),
(if type neq'factor then append(cdr factl,caddr y)
else caddr y)}.res;
factl:=cdr factl>>;
return res end;
symbolic procedure groebtestabort(h,abort2);
% Tests if h is member of one of the sets in abort2.
% if yes, it is deleted. If one wet becomes null,the message
% "CANCEL is returned, otherwise the updated abort2.
begin scalar x,break,res;
% First test the occurence.
x:=abort2;
while x and not break do
<<if vdpmember(h,car x)then break:=t;x:=cdr x>>;
if not break then return abort2;% not relvevant
break:=nil;
while abort2 and not break do
<<x:=vdpdeletemember(h,car abort2);
if null x then break:=t;res:=x.res;
abort2:=cdr abort2>>;
!*trgroeb and groebmess25(h,res);
if break then return'cancel;
return res end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Reduction of polynomials.
%
symbolic procedure groebnormalform(f,g,type);
groebnormalform0(f,g,type,nil);
symbolic procedure groebnormalform0(f,g,type,m);
% General procedure for reduction of one polynomial from a set
%'f'is a polynomial,'g'is a set of polynomials either in
% a search tree or in a sorted list.
%'f'has to be reduced modulo'g'.
%'m'is indicator,whether a selection('m'is true)is wanted.
begin scalar a,break,c,divisor,done,f0,f1,f2,fold,gl,vev;
integer n,s1,s2;
scalar zzz;
if !*groebweak and !*vdpinteger
and groebweakzerotest( f,g,type)then return f2vdp nil;
fold:=f;f1:=vdpzero();a:= vbcfi 1;
gsetsugar(f1,gsugar f);
while not vdpzero!? f do
begin
vev:=vdpevlmon f;c:=vdplbc f;
if not !*groebfullreduction and not vdpzero!? f1 then g:=nil;
if null g then
<<f1:=vdpsum(f1,f);f:=vdpzero();break:=t;divisor:=nil;go to ready>>;
divisor:=groebsearchinlist(vev,g);
if divisor then<<done:=t;% true action indicator
if m and vdpsortmode!*='revgradlex
and vdpzero!? f1 then gl:=f.gl;
if !*trgroebs then
<<prin2"//-";prin2 vdpnumber divisor>> >>;
if divisor then
if vdplength divisor=1 then
f:=vdpcancelmvev(f,vdpevlmon divisor)else
if !*vdpinteger or not !*groebdivide then
<<f:=groebreduceonestepint(f,f1,c,vev,divisor);
f1:=secondvalue!*;n:=n+1;
if not vdpzero!? f and n #> groecontcount!* then
<<f0:=f;
f:=groebsimpcont2(f,f1);
groecontentcontrol(f neq f0);
f1:=secondvalue!*;n:=0>> >>
else
f:=groebreduceonesteprat(f,nil,c,vev,divisor)
else
<<!*gsugar and<<s1:=gsugar(f);s2:=gsugar(f1)>>;
f1:=vdpappendmon(f1,vdplbc f,vdpevlmon f);
f:=vdpred f;
!*gsugar and<<gsetsugar(f,s1);
gsetsugar(f1,max(s2,s1))>> >>;
ready:
end;
if !*groebprot then groebreductionprotocolborder();
if vdpzero!? f1 then go to ret;
zzz:=f1;
if not done then f1:=fold else
if m and vdpsortmode!*='revgradlex then
<<if not vdpzero!? f1 then gl:=f1.gl;
while gl do
<<f2:=groebnormalformselect car gl;
if f2 then<<f1:=f2;gl:=nil>>else gl:=cdr gl>> >>;
ret: return f1 end;
symbolic procedure groecontentcontrol u;
%'u'indicates,that a substantial content reduction was done;
% update content reduction limit from'u'.
groecontcount!*:=if not numberp groecontcount!* then 10 else
if u then max(0,groecontcount!*-1)
else min(10,groecontcount!*+1);
symbolic procedure groebvbcbig!? a;
% Test if'a'is a "big" coefficient.
(if numberp x then(x > 1000000000000 or x <-1000000000000)
else t)where x=vbcnumber a;
symbolic procedure groebnormalformselect v;
% Select the vdp'v',if the'vdplastvar*'- variable occurs in all
% terms(then return it)or don't select it(then return'nil').
if countlastvar(v,t)#> 0 then v;
symbolic procedure groebsimpcontnormalform h;
% SimpCont version preserving the property SUGAR.
if vdpzero!? h then h else
begin scalar sugar,c;
sugar:=gsugar h;c:=vdplbc h;
h:=vdpsimpcont h;gsetsugar(h,sugar);
if !*groebprot and not(c=vdplbc h)then groebreductionprotocol2
reval{'quotient,vbc2a vdplbc h,vbc2a c};
return h end;
symbolic procedure groebsimpcont2(f,f1);
% Simplify two polynomials with the gcd of their contents.
begin scalar c,s1,s2;
s1:=gsugar f;s2:=gsugar f1;
c:=vdpcontent f;
if vbcone!? vbcabs c then go to ready;
if not vdpzero!? f1 then
<<c:=vdpcontent1(f1,c);
if vbcone!? vbcabs c then go to ready;
f1:= vdpdivmon(f1,c,nil)>>;
f:=vdpdivmon(f,c,nil);
!*trgroeb and groebmess28 c;
groebsaveltermbc c;
gsetsugar(f,s1);gsetsugar(f1,s2);
ready:secondvalue!*:=f1;return f end;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
% Special case reductions.
%
symbolic procedure groebprereduce g;
% Reduce the polynomials in g with themselves.
% The reduction is continued until headterms are stable is possible.
begin scalar res,work,oldvev,f,oldf,!*groebweak,
!*groebfullreduction;integer count;
if !*trgroebs then
<<g:=for each p in g collect vdpenumerate p;
for each p in g do vdpprint p>>;
res:=nil;% Delete zero polynomials from'g'.
for each f in g do if not vdpzero!? f then res:=f.res;
work:=g:=res:=reversip res;
while work do
<<g:=vdplsort res;% Sort prvevious result.
if !*trgroebs then prin2t "Starting cycle in prereduction.";
res:=nil;count:=count+1;work:=nil;
while g do
<<oldf:=f:= car g;g:=cdr g;
oldvev:=vdpevlmon f;
f:=vdpsimpcont groebnormalform(f,g,'sort);
if(!*trgroebs or !*groebprot)and not vdpequal(f,oldf)then
<<f:=vdpenumerate f;
if !*groebprot then
if not vdpzero!? f then
groebprotsetq(mkid('poly,vdpnumber f),vdp2a f)
else groebprotval 0;
if !*trgroebs then
<<prin2t "reducing";vdpprint oldf;prin2t "to";vdpprint f>> >>;
if not vdpzero!? f then
<<if oldvev neq vdpevlmon f then work:=t;
res:=f.res>> >> >>;
return for each f in res collect vdpsimpcont f end;
symbolic procedure groebreducefromfactors(g,facts);
% Reduce the polynomials in G from those in facts.
begin scalar new,gnew,f,nold,nnew,numbers;
if !*trgroebs then
<<prin2t "replacing from factors:";
for each x in facts do vdpprin2t x>>;
while g do
<<f:=car g;g:=cdr g;nold:=vdpnumber(f);
new:= groebnormalform(f,facts,'list);
if vdpzero!? new then
<<if !*trgroebs then<<prin2 "vanishes ";
prin2 vdpnumber f>> >>
else if vevzero!? vdpevlmon new then
<<if !*trgroebs then<<prin2 "ONEPOL ";prin2 vdpnumber f>>;
g:=nil;
gnew:={vdpone!*}>>
else<<if new neq f then
<<new:=vdpenumerate vdpsimpcont new;
nnew:=vdpnumber new;
numbers:=(nold.nnew).numbers;
if !*trgroebs then<<prin2 "replacing ";
prin2 vdpnumber f;
prin2 " by ";
prin2t vdpnumber new>> >>;
gnew:=new.gnew>> >>;
secondvalue!*:=numbers;
return gnew end;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
% Support for reduction by "simple" polynomials.
symbolic procedure groebnormalform1(f,p);
% Short version;reduce f by p;
% special case: p is a monomial.
if vdplength p=1 then vdpcancelmvev(f,vdpevlmon p)
else groebnormalform(f,{p},nil);
symbolic procedure groebprofitsfromvev(p,vev);
% Tests,if at least one monomial from p would be reduced by vev.
if vdpzero!? p then nil
else if buchvevdivides!?(vev,vdpevlmon p)then t
else groebprofitsfromvev(vdpred p,vev);
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
% Special reduction procedures.
symbolic procedure groebreduceonestepint(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 vevlcm,a,b,cg,x,rg1;
% Trivial case: g1 single monomial.
if vdpzero!?(rg1:=vdpred g1)
then return<<f:=vdpred f;secondvalue!*:=f1;f>>;
vevlcm:=vevdif(vev,vdpevlmon g1);
cg:=vdplbc g1;
% Calculate coefficient factors .
x:=if not !*groebdivide then vbcfi 1 else vbcgcd(c,cg);
a:=vbcquot(cg,x);
b:=vbcquot(c,x);
% Multiply relvevant parts from f and f1 by a(vbc).
if f1 and not vdpzero!? f1 then f1:=vdpvbcprod(f1,a);
if !*groebprot then groebreductionprotocol(a,vbcneg b,vevlcm,g1);
f:= vdpilcomb1(vdpred f,a,vevzero(),
rg1,vbcneg b,vevlcm);
% Return with f and f1.
secondvalue!*:= f1;thirdvalue!*:=a;return f end;
symbolic procedure groebreduceonesteprat(f,dummy,c,vev,g1);
% Reduction step for rational case:
% calculate f= f-g/vdpLbc(f).
begin scalar x,rg1,vevlcm;
% Trivial case: g1 single monomial.
dummy:=nil;
if vdpzero!?(rg1:=vdpred g1)then return vdpred f;
% Calculate coefficient factors.
x:=vbcneg vbcquot(c,vdplbc g1);
vevlcm:=vevdif(vev,vdpevlmon g1);
if !*groebprot then
groebreductionprotocol( a2vbc 1,x,vevlcm,g1);
return vdpilcomb1(vdpred f,a2vbc 1,vevzero(),
rg1,x,vevlcm)end;
symbolic procedure groebreductionprotocol(a,b,vevlcm,g1);
if !*groebprot then
groebprotfile:=
if not vbcone!? a then
append(groebprotfile,
{{'equal,'candidate,
{'times,'candidate,vbc2a a}},
{'equal,'candidate,
{'plus,'candidate,
{'times,vdp2a vdpfmon(b,vevlcm),
mkid('poly,vdpnumber g1)}}}
})
else
append(groebprotfile,
{{'equal,'candidate,
{'plus,'candidate,
{'times,vdp2a vdpfmon(b,vevlcm),
mkid('poly,vdpnumber g1)}}}
});
symbolic procedure groebreductionprotocol2 a;
if !*groebprot then
groebprotfile:=
if not(a=1)then
append(groebprotfile,
{{'equal,'candidate,{'times,'candidate,a}}});
symbolic procedure groebreductionprotocolborder();
append(groebprotfile,'!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+!+.nil);
symbolic procedure groebprotsetq(a,b);
groebprotfile:=append(groebprotfile,{{'equal,a,b}});
symbolic procedure groebprotval a;
groebprotfile:=
append(groebprotfile,{{'equal,'intermediateresult,a}});
symbolic procedure subset!?(s1,s2);
% Tests,if s1 is a subset of s2.
if null s1 then t else
if member(car s1,s2)then subset!?(cdr s1,s2)
else nil;
symbolic procedure vevsplit vev;
% Split vev such that each exponent vector has only one 1.
begin scalar e,vp;integer n;
for each x in vev do
<<n:=n+1;
if x neq 0 then
<<e:=append(vdpevlmon vdpone!*,nil);
rplaca(pnth(e,n),1);
vp:=e.vp>> >>;return vp end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Calculation of an S-polynomial.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General strategy:
%
% groebspolynom4 calculates the traditional s-polynomial from p1,p2
%(linear combination such that the highest term vanishes).
% groebspolynom2 subtracts multiples of p2 from the s-polynomial such
% that head terms are eliminated early.
symbolic procedure groebspolynom(p1,p2);
groebspolynom2(p1,p2);
symbolic procedure groebspolynom2(p1,p2);
if vdpzero!? p1 then p2 else if vdpzero!? p2 then p1 else
begin scalar cand,s,tp1,tp2,ts;
s:=groebspolynom3(p1,p2);
if vdpzero!? s or vdpone!? s or !*groebprot then return s;
tp1:=vdpevlmon p1;tp2:=vdpevlmon p2;
while not vdpzero!? s
and(( buchvevdivides!?(tp2,(ts:=vdpevlmon s)) and(cand:=p2))
or(buchvevdivides!?(tp1,(ts:=vdpevlmon s))
and(cand:=p1)))
do<<if !*vdpinteger then
s:=% vdpsimpcont
groebreduceonestepint(s,nil,vdplbc s,ts,cand)
else
% Rational, float and modular case .
s:=groebreduceonesteprat(s,nil,vdplbc s,ts,cand)>>;
return s end;
symbolic procedure groebspolynom3(p,q);
begin scalar r;r:=groebspolynom4(p,q);
groebsavelterm r;return r end;
symbolic procedure groebspolynom4(p1,p2);
begin scalar db1,db2,ep1,ep2,ep,r,rp1,rp2,x;
ep1:=vdpevlmon p1;ep2:=vdpevlmon p2;
ep:=vevlcm(ep1,ep2);
rp1:=vdpred p1;rp2:=vdpred p2;
gsetsugar(rp1,gsugar p1);gsetsugar(rp2,gsugar p2);
r:=(if vdpzero!? rp1 and vdpzero!? rp2 then rp1
else(if vdpzero!? rp1 then
<<db2:=a2vbc 0;
vdpprod( rp2,vdpfmon(db1:=a2vbc 1,vevdif(ep,ep2)))>>
else if vdpzero!? rp2 then
<<db1:=a2vbc 0;
vdpprod(rp1,vdpfmon(db2:=a2vbc 1,
vevdif(ep,ep1)))>>
else
<<db1:=vdplbc p1;
db2:=vdplbc p2;
if !*vdpinteger then
<<x:= vbcgcd(db1,db2);
if not vbcone!? x then
<<db1:=vbcquot(db1,x);db2:=vbcquot(db2,x)>> >>;
vdpilcomb1(rp2,db1,vevdif(ep,ep2),
rp1,vbcneg db2,vevdif(ep,ep1)) >>
));
if !*groebprot then
groebprotsetq('candidate,
{'difference,
{'times,vdp2a vdpfmon(db2,vevdif(ep,ep2)),
mkid('poly,vdpnumber p1)},
{'times,vdp2a vdpfmon(db1,vevdif(ep,ep1)),
mkid('poly,vdpnumber p2)}});
return r end;
symbolic procedure groebsavelterm r;
if !*groelterms and not vdpzero!? r then groebsaveltermbc vdplbc r;
symbolic procedure groebsaveltermbc r;
<<r:=vbc2a r;
if not numberp r and not constant_exprp r then
for each p in cdr fctrf numr simp r do
<<p:=prepf car p;
if not member(p,glterms)then nconc(glterms,{p})>> >>;
symbolic procedure sfcont f;
% Calculate the integer content of standard form f.
if domainp f then f else gcdf(sfcont lc f,sfcont red f);
symbolic procedure vdplmon u;vdpfmon(vdplbc u,vdplbc u);
symbolic procedure vdpmember3(p,g1,g2,g3);
% Test membership of p in one of then lists g1,g2,g3.
vdpmember(p,g1)or vdpmember(p,g2)or vdpmember(p,g3);
symbolic procedure groebabortid(base,abort1);
% Test whether one of the elements in abort1 is
% member of the ideal described by base. Definite
% test here.
if null abort1 then nil else
vdpzero!?(groebnormalform(car abort1,base,'list))
or groebabortid(base,cdr abort1);
endmodule;;end;