module groebtra;
% calculation of a Groebner base with the Buchberger algorithm
% including the backtracking information which denotes the
% dependency between base and input polynomials
% Authors: H. Melenk, H.M. Moeller, W. Neun
% November 1988
fluid '( % switches from the user interface
!*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
!*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
!*groebstat !*groebdivide !*groebprot
vdpvars!* % external names of the variables
!*vdpinteger !*vdpmodular % indicating type of algorithm
vdpSortMode!* % term ordering mode
secondvalue!* thirdvalue!* % auxiliary: multiple return values
fourthvalue!*
factortime!* % computing time spent in factoring
factorlvevel!* % bookkeeping of factor tree
pairsdone!* % list of pairs already calculated
probcount!* % counting subproblems
vbcCurrentMode!* % current domain for base coeffs.
groetags!* % starting point of tag variables
groetime!*
);
global '(gvarslast);
switch groebopt,groebfac,groebres,trgroeb,trgroebs,trgroeb1,
trgroebr,groebstat,groebprot;
% variables for counting and numbering
fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
basecount!* hzerocount!*);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interface
symbolic procedure groebnertraeval u;
% backtracking Groebner calculation
begin integer n; scalar !*groebfac,!*groebrm,!*groebprot,!*gsugar;
n := length u;
if n=1 then return groebnertra1(reval car u,nil,nil)
else if n neq 2
then rerror(groebnr2,10,
"GROEBNERT called with wrong number of arguments")
else return groebnertra1(reval car u,reval cadr u,nil)
end;
put('groebnert,'psopfn,'groebnertraeval);
smacro procedure vdpnumber f;
vdpgetprop(f,'number) ;
symbolic procedure groebnertra1(u,v,mod1);
% Buchberger algorithm system driver. u is a list of expressions
% and v a list of variables or NIL in which case the variables in u
% are used.
begin scalar vars,w,y,z,x,np,oldorder,groetags!*,tagvars;
integer pcount!*,nmod;
u := for each j in getrlist u collect
<<if not eqcar(j,'equal)
or not (idp (w:=cadr j) or (pairp w and
w = reval w and
get(car w,'simpfn)='simpiden))
then rerror(groebnr2,11,
"groebnert parameter not {...,name = polynomial,...}")
else cadr j . caddr j
>>;
if null u then rerror(groebnr2,12,"Empty list in Groebner")
else if null cdr u then
return 'list . list('equal,cdar u,caar u);;
w := for each x in u collect cdr x;
if mod1 then
<<z := nil;
for each vect in w do
<<if not eqcar(vect,'list) then
rerror(groebnr2,13,"Non list given to groebner");
if nmod=0 then nmod:= length cdr vect else
if nmod<(x:=length cdr vect) then nmod=x;
z := append(cdr vect,z);
>>;
if not eqcar(mod1,'list) then
rerror(groebnr2,14,"Illegal column weights specified");
vars := groebnervars(z,v);
tagvars := for i:=1:nmod collect mkid('! col,i);
w := for each vect in w collect
<<z:=tagvars; y := cdr mod1;
'plus . for each p in cdr vect collect
'times . list('expt,car z,car y) .
<<z := cdr z; y := cdr y;
if null y then y := '(1); list p>>
>>;
groetags!* := length vars + 1;
vars := append(vars,tagvars);
>>
else vars := groebnervars(w,v);
groedomainmode();
if null vars then vdperr 'groebner;
oldorder := vdpinit vars;
% cancel common denominators
u := pair(for each x in u collect car x,w);
u := for each x in u collect
<<z := simp cdr x;
multsq(simp car x,denr z ./ 1) . reorder numr z>>;
w := for each j in u collect cdr j;
% optimize varable sequence if desired
if !*groebopt then<< w:=vdpvordopt (w,vars); vars := cdr w;
w := car w; vdpinit vars>>;
w := pair (for each x in u collect car x,w);
w := for each j in w collect
<<z := f2vdp cdr j; vdpPutProp(z,'cofact,car j)>>;
if not !*vdpInteger then
<<np := t;
for each p in w do
np := if np then vdpcoeffcientsfromdomain!? p else nil;
if not np then <<!*vdpModular:= nil; !*vdpinteger := T>>;
>>;
w := groebtra2 w;
w := if mod1 then groebnermodres(w,nmod,tagvars) else
groebnertrares w;
setkorder oldorder;
gvarslast := 'list . vars; return w;
end;
symbolic procedure groebnertrares w;
begin scalar c,u;
return 'list . for each j in w collect
<<c := prepsq vdpgetprop(j,'cofact);
u := vdp2a j;
if c=0 then u else list('equal,u,c) >>;
end;
symbolic procedure groebnermodres(w,n,tagvars);
begin scalar x,c,oldorder;
c := for each u in w collect prepsq vdpgetprop(u,'cofact);
oldorder := setkorder tagvars;
w := for each u in w collect
'list . <<u := numr simp vdp2a u;
for i := 1:n collect
prepf if not domainp u and mvar u = nth(tagvars,i)
then <<x := lc u; u := red u; x>> else nil >>;
setkorder oldorder;
% reestablish term order for output
w := for each u in w collect vdp2a a2vdp u;
w := pair(w,c);
return 'list . for each p in w collect
if cdr p=0 then car p else
list('equal,car p,cdr p);
end;
symbolic procedure preduceteval pars;
% trace version of PREDUCE
% parameters:
% 1 expression to be reduced
% formula or equation
% 2 polynomials or equations; base for reduction
% must be equations with atomic lhs
% 3 optional: list of variables
begin scalar vars,x,y,u,v,w,z,oldorder,!*factor,!*exp,!*gsugar;
integer pcount!*; !*exp := t;
pars := groeparams(pars,2,3);
y := car pars; u := cadr pars; v:= caddr pars;
u := for each j in getrlist u
collect
<<if not eqcar(j,'equal) then
j . j else cadr j . caddr j>>;
if null u then rerror(groebnr2,15,"Empty list in preducet");
w := for each p in u collect cdr p; % the polynomials
groedomainmode();
vars := if null v then
for each j in gvarlis w collect !*a2k j
else getrlist v;
if not vars then vdperr 'preducet;
oldorder := vdpinit vars;
u := for each x in u collect
<<z := simp cdr x;
multsq(simp car x,denr z ./ 1) . reorder numr z>>;
w := for each j in u collect
<<z := f2vdp cdr j; vdpputprop(z,'cofact,car j)>>;
if not eqcar(y,'equal) then y := list('equal,y,y);
x := a2vdp caddr y; % the expression
vdpputprop(x,'cofact,simp cadr y); % the lhs (name etc.)
w := tranormalform(x,w, 'sort,'f);
u := list('equal,vdp2a w,prepsq vdpgetprop(w,'cofact));
setkorder oldorder;
return u;
end;
put('preducet,'psopfn,'preduceteval);
symbolic procedure groebnermodeval u;
% Groebner for moduli calculation
(if n=0 or n>3 then
rerror(groebnr2,16,
"groebnerm called with wrong number of arguments")
else
groebnertra1(reval car u,
if n>=2 then reval cadr u else nil,
if n>=3 then reval caddr u else '(list 1))
) where n = length u;
put('groebnerm,'psopfn,'groebnermodeval);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some macros for local usage only
smacro procedure tt(s1,s2);
% lcm of leading terms of s1 and s2
vevlcm(vdpevlmon s1,vdpevlmon s2);
symbolic procedure groebtra2 (p);
% setup all global variables for the Buchberger algorithm
% printing of statistics
begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*,
pairsdone!*,factorlvevel!*;
factortime!* := 0;
groetime!* := time();
vdponepol(); % we construct dynamically
hcount!* := pcount!* := mcount!* := fcount!* :=
bcount!* := b4count!* := hzerocount!* := basecount!* := 0;
if !*trgroeb then
<< prin2 "Groebner Calculation with Backtracking starts ";
terprit 2 >>;
spac := gctime();
p1:= groebTra3 (p);
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 ";
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;
symbolic procedure groebTra3 (g0);
begin scalar x,g,d,d1,d2,p,p1,s,h,g99,one;
x := for each fj in g0 collect vdpenumerate trasimpcont fj;
for each fj in x do g := vdplsortin(fj,g0);
g0 := g; g := nil;
% iteration :
while (d or g0) and not one do
begin if g0 then
<< % take next poly from input
h := car g0; g0 := cdr g0;
p := list(nil,h,h) >>
else
<< % take next poly from pairs
p := car d; d := cdr d;
s := traspolynom (cadr p, caddr p); tramess3(p,s);
h := groebnormalform(s,g99,'tree); %piloting wo cofact
if vdpzero!? h then groebmess4(p,d)
else h := trasimpcont tranormalform(s,g99,'tree,'h) >>;
if vdpzero!? h then goto bott;
if vevzero!? vdpevlmon h then % base 1 found
<< tramess5(p,h);
g0 := d := nil; g:= list h;
goto bott>>;
s:= nil;
% h polynomial is accepted now
h := vdpenumerate h;
tramess5(p,h);
% construct new critical pairs
d1 := nil;
for each f in g do
if groebmoducrit(f,h) then
<<d1 := groebcplistsortin( list(tt(f,h),f,h) , d1);
if tt(f,h) = vdpevlmon(f) then
<<g := delete (f,G); groebmess2 f>> >>;
groebmess51(d1);
d2 := nil;
while d1 do
<<d1 := groebinvokecritf d1;
p1 := car d1;
d1 := cdr d1;
if groebbuchcrit4(cadr p1,caddr p1,car p1)
then d2 := append (d2, list p1);
d1 := groebinvokecritm (p1,d1) >>;
d := groebinvokecritb (h,d);
d := groebcplistmerge(d,d2);
g := h . g;
g99 := groebstreeadd(h, g99);
groebmess8(g,d);
bott:
end;
return groebtra3post g;
end;
symbolic procedure groebtra3post (g);
% final reduction
begin scalar r,p;
g := vdplsort g;
while g do
<<p := tranormalform(car G,cdr G,'sort,'f);
if not vdpzero!? p then r := trasimpcont p . r;
g := cdr g>>;
return reversip r;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Reduction of polynomials
%
symbolic procedure tranormalform(f,g,type,mode);
% 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
% type describes the ordering of the set G:
% 'TREE G is a search tree
% 'SORT G is a sorted list
% 'LIST G is a list, but not sorted
% f has to be reduced modulo G
% version for idealQuotient: doing side effect calculations for
% the cofactors; only headterm reduction
begin scalar c,vev,divisor,break;
while not vdpzero!? f and not break do
begin
vev:=vdpevlmon f; c:=vdpLbc f;
divisor :=
if type = 'tree then groebsearchinstree (vev,g)
else groebsearchinlist (vev,g);
if divisor and !*trgroebs then
<< prin2 "//-";
prin2 vdpnumber divisor >>;
if divisor then
if !*vdpInteger then
f := trareduceonestepint(f,nil,c,vev,divisor)
else
f := trareduceonesteprat(f,nil,c,vev,divisor)
else
break := t;
end;
if mode = 'f then f := tranormalform1(f,g,type,mode);
return f
end;
symbolic procedure tranormalform1(f,g,type,mode);
% reduction of subsequent terms
begin scalar c,vev,divisor,break,f1;
mode := nil;
f1 := f;
while not vdpzero!? f and not vdpzero!? f1 do
<<f1 := f; break := nil;
while not vdpzero!? f1 and not break do
<<vev:=vdpevlmon f1; c:=vdpLbc f1;
f1 := vdpred f1;
divisor :=
if type = 'tree then groebsearchinstree (vev,g)
else groebsearchinlist (vev,g);
if divisor and !*trgroebs then
<< prin2 "//-";
prin2 vdpnumber divisor >>;
if divisor then
<< if !*vdpInteger then
f := trareduceonestepint(f,nil,c,vev,divisor)
else
f := trareduceonesteprat(f,nil,c,vev,divisor);
break := t>> >> >>;
return f;
end;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
%
% special reduction procedures
symbolic procedure trareduceonestepint(f,dummy,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,fcofa,gcofa;
dummy := nil;
fcofa := vdpgetprop(f,'cofact);
if null fcofa then fcofa := nil ./ 1;
gcofa := vdpgetprop(g1,'cofact);
if null gcofa then gcofa := nil ./ 1;
vevlcm := vevdif(vev,vdpevlmon g1);
cg := vdpLbc g1;
% calculate coefficient factors
x := vbcgcd(c,cg);
a := vbcquot(cg,x);
b := vbcquot(c,x);
f:= vdpilcomb1 (f, a, vevzero(), g1,vbcneg b, vevlcm);
x := vdpilcomb1tra (fcofa, a, vevzero(),
gcofa ,vbcneg b, vevlcm);
vdpputprop(f,'cofact,x);
return f;
end;
symbolic procedure trareduceonesteprat(f,dummy,c,vev,g1);
% reduction step for rational case:
% calculate f= f - g/vdpLbc(f)
%
begin scalar x,fcofa,gcofa,vev;
dummy := nil;
fcofa := vdpgetprop(f,'cofact);
gcofa := vdpgetprop(g1,'cofact);
vev := vevdif(vev,vdpevlmon g1);
x := vbcneg vbcquot(c,vdplbc g1);
f := vdpilcomb1(f,a2vbc 1,vevzero(), g1,x,vev);
x := vdpilcomb1tra (fcofa,a2vbc 1 , vevzero(),
gcofa,x,vev);
vdpputprop(f,'cofact,x);
return f;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% calculation of an S-polynomial
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure traspolynom (p1,p2);
begin scalar s,ep1,ep2,ep,rp1,rp2,db1,db2,x,cofac1,cofac2;
if vdpzero!? p1 then return p1; if vdpzero!? p1 then return p2;
cofac1 := vdpgetprop(p1,'cofact); cofac2 := vdpGetProp(p2,'cofact);
ep1 := vdpevlmon p1; ep2 := vdpevlmon p2;
ep := vevlcm(ep1, ep2);
rp1 := vdpred p1; rp2 := vdpred p2;
db1 := vdpLbc p1; db2 := vdpLbc p2;
if !*vdpinteger then
<<x:=vbcgcd(db1,db2); db1:=vbcquot(db1,x); db2:=vbcquot(db2,x)>>;
ep1 := vevdif(ep,ep1); ep2 := vevdif(ep,ep2); db2 := vbcneg db2;
s := vdpilcomb1 (rp2,db1,ep2,rp1,db2,ep1);
x := vdpilcomb1tra (cofac2,db1,ep2,cofac1,db2,ep1);
vdpputprop(s,'cofact,x);
return s;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Normalisation with cofactors taken into account
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
symbolic procedure trasimpcont(p);
if !*vdpInteger then trasimpconti p else trasimpcontr p;
% routines for integer coefficient case:
% calculation of contents and dividing all coefficients by it
symbolic procedure trasimpconti (p);
% calculate the contents of p and divide all coefficients by it
begin scalar res,num,cofac;
if vdpzero!? p then return p;
cofac := vdpgetprop(p,'cofact);
num := car vdpcontenti p;
if not vbcplus!? num then num := vbcneg num;
if not vbcplus!? vdpLbc p then num:=vbcneg num;
if vbcone!? num then return p;
res := vdpreduceconti (p,num,nil);
cofac:=vdpreducecontitra(cofac,num,nil);
res := vdpputprop(res,'cofact,cofac);
return res;
end;
% routines for rational coefficient case:
% calculation of contents and dividing all coefficients by it
symbolic procedure trasimpcontr (p);
% calculate the contents of p and divide all coefficients by it
begin scalar res,cofac;
cofac := vdpgetprop(p,'cofact);
if vdpzero!? p then return p;
if vbcone!? vdpLbc p then return p;
res := vdpreduceconti (p,vdplbc p,nil);
cofac:=vdpreducecontitra(cofac,vdplbc p,nil);
res := vdpputprop(res,'cofact,cofac);
return res;
end;
symbolic procedure vdpilcomb1tra (cofac1,db1,ep1,cofac2,db2,ep2);
% the linear combination, here done for the cofactors
% (standard quotients)
addsq(multsq(cofac1,vdp2f vdpfmon(db1,ep1) ./ 1),
multsq(cofac2,vdp2f vdpfmon(db2,ep2) ./ 1));
symbolic procedure vdpreducecontitra(cofac,num,dummy);
% divide the cofactor by a number
<<dummy := nil; quotsq(cofac,simp vbc2a num)>>;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% special handling of moduli
%
symbolic procedure groebmoducrit(p1,p2);
null groetags!*
or pnth(vdpevlmon p1,groetags!*) = pnth(vdpevlmon p2,groetags!*);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% trace messages
%
symbolic procedure tramess0 x;
if !*trgroeb then
<<prin2t "adding member to intermediate quotient base:";
vdpprint x;
terpri()>>;
symbolic procedure tramess1 (x,p1,p2);
if !*trgroeb then
<<prin2 "pair ("; prin2 vdpnumber p1; prin2 ",";
prin2 vdpnumber p2;
prin2t ") adding member to intermediate quotient base:";
vdpprint x;
terpri()>>;
symbolic procedure tramess5(p,h);
if car p then % print for true h-Polys
<< hcount!* := hcount!* + 1;
if !*trgroeb then << terpri();
prin2 "H-polynomial ";
prin2 pcount!*;
groebmessff(" from pair (",cadr p,nil);
groebmessff(",", caddr p,")");
vdpprint h;
prin2t "with cofactor";
writepri(mkquote prepsq vdpgetprop(h,'cofact),'only);
groetimeprint() >> >>
else
if !*trgroeb then << % print for input polys
prin2t "candidate from input:";
vdpprint h;
groetimeprint() >>;
symbolic procedure tramess3(p,s);
if !*trgroebs then <<
prin2 "S-polynomial ";
prin2 " from ";
groebpairprint p;
vdpprint s;
prin2t "with cofactor";
writepri(mkquote prepsq vdpGetProp(s,'cofact),'only);
groetimeprint();
terprit 3 >>;
endmodule;
end;