Artifact 984e991c9876c48df7614079fd1b7da2ae391d50405fd00624ebe00a4a618314:
- Executable file
r38/packages/groebner/buchbg.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 35648) [annotate] [blame] [check-ins using] [more...]
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;