Artifact dbbe4c435f95deebde1dbde10a773a46de5656b9bddfb90c3c07de898550e2d7:
- Executable file
r37/packages/groebner/groebidq.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: 15397) [annotate] [blame] [check-ins using] [more...]
module groebidq; % calculation of ideal quotient using a modified Buchberger algorithm. % Authors: H. Melenk, H.M. Moeller, W. Neun % July 1988 fluid '( % switches from the user interface !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!* !*groebstat !*groebdivide !*groebprot !*groebidqbasis 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. groetime!* !*gsugar ); switch groebopt,groebfac,groebrm,groebres,trgroeb,trgroebs,trgroeb1, trgroebr,groebstat,groebprot; !*groebidqbasis := t; %default: basis from idq % variables for counting and numbering fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!* basecount!* hzerocount!*); symbolic procedure groebidq2 (p,f); % setup all global variables for the Buchberger algorithm % printing of statistics begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*, pairsdone!*, factorlvevel!*,!*gsugar; factortime!* := 0; groetime!* := time(); vdponepol(); % we construct dynamically hcount!* := 0; pcount!* := 0; mcount!* := 0; fcount!* := 0; bcount!* := 0; b4count!* := 0; hzerocount!* := 0; basecount!* := 0; if !*trgroeb then << prin2 "IDQ Calculation starting "; terprit 2; >>; spac := gctime(); p1:= groebidq3 (p, f); if !*trgroeb or !*trgroebr or !*groebstat then << spac1 := gctime() - spac; terpri(); prin2t "statistics for IDQ 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 if !*groebidqbasis then car groebner2(p1,nil) else p1; end; symbolic procedure groebidq3 (g0,fff); begin scalar result,x,g,d,d1,d2,p,p1,s,h,g99,one,gi; gi := g0; fff := vdpsimpcont fff; vdpputprop(fff,'number,0); % assign number 0 vdpputprop(fff,'cofact,a2vdp 1); % assign cofactor 1 x := for each fj in g0 collect << fj:=vdpenumerate vdpsimpcont fj; vdpputprop(fj,'cofact,a2vdp 0); % assign cofactor 0 fj>>; g0 := list fff; for each fj in x do g0 := vdplsortin(fj,g0); % 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 := delete (p,d); s := idqspolynom (cadr p, caddr p); idqmess3(p,s); h := idqsimpcont idqnormalform(s,g99,'tree); if vdpzero!? h then << !*trgroeb and groebmess4(p,d); x := vdpgetprop(h,'cofact); if not vdpzero!? x then if vevzero!? vdpevlmon x then one:= t else << result := idqtoresult(x , result); idqmess0 x>>; >> >>; if vdpzero!? h then goto bott; if vevzero!? vdpevlmon h then % base 1 found << idqmess5(p,h); result:=gi; d:=g0:=nil; goto bott;>>; s:= nil; % h polynomial is accepted now h := vdpenumerate h; idqmess5(p,h); % construct new critical pairs d1 := nil; for each f in g do <<d1 := groebcplistsortin( list(tt(f,h),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; if groebbuchcrit4t(cadr p1,caddr p1) then d2 := append (d2, list p1) else <<x := idqdirectelement(cadr p1,caddr p1); if not vdpzero!? x then if vevzero!? vdpevlmon x then one:= t else <<idqmess1 (x,cadr p1,caddr p1); result := idqtoresult(x,result); >> >>; d1 := groebinvokecritm (p1,d1); >>; % D := groebInvokeCritB (h,D); d := groebcplistmerge(d,d2); g := h . g; g99 := groebstreeadd(h, g99); !*trgroeb and groebmess8(g,d); bott: end; % ITERATION % now calculate groebner base from quotient base if one then result := list vdpfmon(a2vbc 1,vevzero()); idqmess2 result; return result; end; % MACROLOOP symbolic procedure idqtoresult(x,r); % x is a new element for the quotient r; % is is reduced by r and then added. <<x := groebsimpcontnormalform groebnormalform(x,r,'sort); if vdpzero!? x then r else vdplsortin(x,r)>>; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Reduction of polynomials % symbolic procedure idqnormalform(f,g,type); % 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,done,fold; fold := f; while not vdpzero!? f and g do begin vev:=vdpevlmon f; c:=vdplbc f; if type = 'sort then while g and vevcompless!? (vev,vdpevlmon (car g)) do g := cdr g; divisor := if type = 'tree then groebsearchinstree (vev,g) else groebsearchinlist (vev,g); if divisor then done := t; % true action indicator if divisor and !*trgroebs then << prin2 "//-"; prin2 vdpnumber divisor; >>; if divisor then if !*vdpinteger then f := idqreduceonestepint(f,nil,c,vev,divisor) else f := idqreduceonesteprat(f,nil,c,vev,divisor) else g := nil; end; return if done then f else fold; %in order to preserve history end; % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % special reduction procedures symbolic procedure idqreduceonestepint(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); gcofa := vdpgetprop(g1,'cofact); 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 (vdpred f, a, vevzero(), vdpred g1,vbcneg b, vevlcm); x := vdpilcomb1 (fcofa, a, vevzero(), gcofa ,vbcneg b, vevlcm); vdpputprop(f,'cofact,x); return f; end; symbolic procedure idqreduceonesteprat(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(vdpred f,a2vbc 1,vevzero(), vdpred g1,x,vev); x := vdpilcomb1 (fcofa,a2vbc 1 , vevzero(), gcofa,x,vev); vdpputprop(f,'cofact,x); return f; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % calculation of an S-polynomial and related things % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% symbolic procedure idqspolynom (p1,p2); begin scalar s,ep1,ep2,ep,rp1,rp2,db1,db2,x,cofac1,cofac2; if vdpzero!? p1 then return p1; if vdpzero!? p2 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 := vdpilcomb1 (cofac2,db1,ep2,cofac1,db2,ep1); vdpputprop(s,'cofact,x); return s; end; symbolic procedure idqdirectelement(p1,p2); % the s-Polynomial is reducable to zero because of % buchcrit 4. So we can calculate the corresponing cofactor % directly (if vdpzero!? c1 and vdpzero!? c2 then c1 else vdpdif(vdpprod(p1,c2),vdpprod(p2,c1)) ) where c1 = vdpgetprop(p1,'cofact), c2 = vdpgetprop(p2,'cofact); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Normailsation with cofactors taken into account % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% symbolic procedure idqsimpcont(p); if !*vdpinteger then idqsimpconti p else idqsimpcontr p; % routines for integer coefficient case: % calculation of contents and dividing all coefficients by it symbolic procedure idqsimpconti (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 vdpzero!? cofac then num:=vbcgcd(num,car vdpcontenti cofac); 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); if not vdpzero!? cofac then cofac:=vdpreduceconti(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 idqsimpcontr (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); if not vdpzero!? cofac then cofac:=vdpreduceconti(cofac,vdplbc p,nil); res := vdpputprop(res,'cofact,cofac); return res; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % trace messages % symbolic procedure idqmess0 x; if !*trgroeb then <<prin2t "adding member to intermediate quotient basis:"; vdpprint x; terpri()>>; symbolic procedure idqmess1 (x,p1,p2); if !*trgroeb then <<prin2 "pair ("; prin2 vdpnumber p1; prin2 ","; prin2 vdpnumber p2; prin2t ") adding member to intermediate quotient basis:"; vdpprint x; terpri()>>; symbolic procedure idqmess2 b; if !*trgroeb then <<prin2t "---------------------------------------------------"; prin2 "the full intermediate base of the ideal quotient is:"; for each x in b do vdpprin3t x; prin2t "---------------------------------------------------"; terpri()>>; symbolic procedure idqmess5(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"; vdpprint vdpgetprop(h,'cofact); groetimeprint() >> >> else if !*trgroeb then << % print for input polys prin2t "candidate from input:"; vdpprint h; groetimeprint() >>; symbolic procedure idqmess3(p,s); if !*trgroebs then << prin2 "S-polynomial "; prin2 " from "; groebpairprint p; vdpprint s; prin2t "with cofactor"; vdpprint vdpgetprop(s,'cofact); groetimeprint(); terprit 3 >>; endmodule; end;