Artifact e3cb8b7cfc8580172a29d25b7e1412edc98279d02d136e726f4d4511273617e2:
- File
r33/groebner.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: 69793) [annotate] [blame] [check-ins using] [more...]
module consel; %/*Constructors and selectors for a distributed polynomial form*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ %/*A distributive polynomial has the following informal syntax: % % <dipoly> ::= dipzero % | <exponent vector> . <base coefficient> . <dipoly>*/ %define dipzero = 'nil; fluid '(dipzero); %/*Until we understand how to define something to nil*/ smacro procedure dipzero!? u; null u; smacro procedure diplbc p; % /* Distributive polynomial leading base coefficient. % p is a distributive polynomial. diplbc(p) returns % the leading base coefficient of p. */ cadr p; smacro procedure dipmoncomp (a,e,p); % /* Distributive polynomial monomial composition. a is a base % coefficient, e is an exponent vector and p is a % distributive polynomial. dipmoncomp(a,e,p) returns a dis- % tributive polynomial with p as monomial reductum, e as % exponent vector of the leading monomial and a as leading % base coefficient. */ e . a . p; smacro procedure dipevlmon p; % /* Distributive polynomial exponent vector leading monomial. % p is a distributive polynomial. dipevlmon(p) returns the % exponent vector of the leading monomial of p. */ car p; smacro procedure dipfmon (a,e); % /* Distributive polynomial from monomial. a is a base coefficient % and e is an exponent vector. dipfmon(a,e) returns a % distributive polynomial with e as exponent vector and % a as base coefficient. */ e . a . dipzero; smacro procedure dipnov p; % /* Distributive polynomial number of variables. p is a distributive % polynomial. dipnov(p) returns a digit, the number of variables % of the distributive polynomial p. */ length car p; smacro procedure dipmred p; % /* Distributive polynomial reductum. p is a distributive polynomial % dipmred(p) returns the reductum of the distributive polynomial p, % a distributive polynomial. */ cddr p; endmodule; module bcoeff; % Computation of base coefficients. %/*Definitions of base coefficient operations for distributive % polynomial package. We assume that only field elements are used, but % no check is made for this. In this module, a standard quotient % coefficient is assumed*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ global '(!*nat); expr procedure bcless!? (a1,a2); % /* Base coefficient less. a1 and a2 are base coefficients. % bcless!?(a1,a2) returns a boolean expression, true if % a1 is less than a2 else false. */ minusf numr addsq(a1,negsq a2); smacro procedure bcminus!? u; % /* Boolean function. Returns true if u is a negative base coeff*/ minusf numr u; smacro procedure bczero!? u; % /* Returns a boolean expression, true if the base coefficient u is % zero*/ null numr u; expr procedure bccomp (a1,a2); % /* Base coefficient compare a1 and a2 are base coefficients. % bccomp(a1,a2) compares the base coefficients a1 and a2 and returns % a digit 1 if a1 greater than a2, a digit 0 if a1 equals a2 else a % digit -1. */ (if bczero!? sl then 0 else if bcminus!? sl then -1 else 1) where sl = bcdif(a1, a2); expr procedure bcfi a; % /* Base coefficient from integer. a is an integer. bcfi(a) returns % the base coefficient a. */ mkbc(a,1); expr procedure bclcmd(u,v); % Base coefficient least common multiple of denominators. % u and v are two base coefficients. bclcmd(u,v) calculates the % least common multiple of the denominator of u and the % denominator of v and returns a base coefficient of the form % 1/lcm(denom u,denom v). if bczero!? u then mkbc(1,denr v) else if bczero!? v then mkbc(1,denr u) else mkbc(1,multf(quotf(denr u,gcdf(denr u,denr v)),denr v)); expr procedure bclcmdprod(u,v); % Base coefficient least common multiple denominator product. % u is a basecoefficient of the form 1/integer. v is a base % coefficient. bclcmdprod(u,v) calculates (denom u/denom v)*nom v/1 % and returns a base coefficient. mkbc(multf(quotf(denr u,denr v),numr v),1); expr procedure bcquod(u,v); % Base coefficient quotient. u and v are base coefficients. % bcquod(u,v) calculates u/v and returns a base coefficient. bcprod(u,bcinv v); expr procedure bcone!? u; % /* Base coefficient one. u is a base coefficient. % bcone!?(u) returns a boolean expression, true if the % base coefficient u is equal 1. */ denr u = 1 and numr u = 1; expr procedure bcinv u; % /* Base coefficient inverse. u is a base coefficient. % bcinv(u) calculates 1/u and returns a base coefficient. */ invsq u; expr procedure bcneg u; % /* Base coefficient negative. u is a base coefficient. % bcneg(u) returns the negative of the base coefficient % u, a base coefficient. */ negsq u; expr procedure bcprod (u,v); % /* Base coefficient product. u and v are base coefficients. % bcprod(u,v) calculates u*v and returns a base coefficient. multsq(u,v); expr procedure mkbc(u,v); % /* Convert u and v into u/v in lowest terms*/ if v = 1 then u ./ v else if v<0 then mkbc(negf u,negf v) else quotf(u,m) ./ quotf(v,m) where m = gcdf(u,v); expr procedure bcquot (u,v); % /* Base coefficient quotient. u and v are base coefficients. % bcquot(u,v) calculates u/v and returns a base coefficient. */ quotsq(u,v); expr procedure bcsum (u,v); % /* Base coefficient sum. u and v are base coefficients. % bcsum(u,v) calculates u+v and returns a base coefficient. */ addsq(u,v); expr procedure bcdif(u,v); % /* Base coefficient difference. u and v are base coefficients. % bcdif(u,v) calculates u-v and returns a base coefficient. */ bcsum(u,bcneg v); expr procedure bcpow(u,n); % /*Returns the base coefficient u raised to the nth power, where % n is an integer*/ exptsq(u,n); expr procedure a2bc u; % /*Converts the algebraic (kernel) u into a base coefficient. simp!* u; expr procedure bc2a u; % /* Returns the prefix equivalent of the base coefficient u*/ prepsq u; expr procedure bcprin u; % /* Prints a base coefficient in infix form*/ begin scalar nat; nat := !*nat; !*nat := nil; sqprint u; !*nat := nat end; endmodule; module expvec; % /*Specific support for distributive polynomial exponent vectors*/ % /* Authors: R. Gebauer, A. C. Hearn, H. Kredel */ % We assume here that an exponent vector is a list of integers. This % version uses small integer arithmetic on the individual exponents % and assumes that a compiled function can be dynamically redefined*/ fluid '(dipsortmode!* dipvars!*); expr procedure evperm (e1,n); % /* Exponent vector permutation. e1 is an exponent vector, n is a % index list , a list of digits. evperm(e1,n) returns a list e1 % permuted in respect to n. */ if null n then nil else evnth(e1, car n) . evperm(e1, cdr n); expr procedure evcons (e1,e2); % /* Exponent vector construct. e1 and e2 are exponents. evcons(e1,e2) % constructs an exponent vector. */ e1 . e2; expr procedure evnth (e1,n); % /* Exponent vector n-th element. e1 is an exponent vector, n is a % digit. evnth(e1,n) returns the n-th element of e1, an exponent. */ if n = 1 then evfirst e1 else evnth(evred e1, n - 1); expr procedure evred e1; % /* Exponent vector reductum. e1 is an exponent vector. evred(e1) % returns the reductum of the exponent vector e1. */ cdr e1; expr procedure evfirst e1; % /* Exponent vector first. e1 is an exponent vector. evfirst(e1) % returns the first element of the exponent vector e1, an exponent. */ car e1; expr procedure evsum0(n,p); % exponent vector sum version 0. n is the length of dipvars!*. % p is a distributive polynomial. if dipzero!? p then evzero1 n else evsum(dipevlmon p, evsum0(n,dipmred p)); expr procedure evzero1 n; % Returns the exponent vector power representation % of length n for a zero power. begin scalar x; for i:=1: n do << x := 0 . x >>; return x end; expr procedure indexcpl(ev,n); % returns a list of indixes of non zero exponents. if null ev then ev else ( if car ev = 0 then indexcpl(cdr ev,n + 1) else ( n . indexcpl(cdr ev,n + 1)) ); expr procedure evzer1!? e; % returns a boolean expression. true if e is null else false. null e; expr procedure evzero!? e; % /* Returns a boolean expression. True if all exponents are zero*/ null e or car e = 0 and evzero!? cdr e; expr procedure evzero; % /* Returns the exponent vector representation for a zero power*/ % for i := 1:length dipvars!* collect 0; begin scalar x; for i := 1:length dipvars!* do <<x := 0 . x>>; return x end; expr procedure mkexpvec u; % /* Returns an exponent vector with a 1 in the u place*/ if not(u member dipvars!*) then typerr(u,"dipoly variable") else for each x in dipvars!* collect if x eq u then 1 else 0; expr procedure evcompless!?(e1,e2); % /* Exponent vector compare less. e1, e2 are exponent vectors % in some order. Evcompless? is a boolean function which returns % true if e1 is ordered less than e2. This function is assigned a % value by the ordering mechanism, so is dummy for now*/ apply(get(dipsortmode!*,'evcompless!?),list(e1,e2)); expr procedure evlexcompless!?(e1,e2); % /* Exponent vector lexicographical compare less. e1, e2 are exponent % vectors in lexicographical order. Evlexcompless?(e1,e2) is a % boolean function which returns true if e1 is ordered less than e2*/ if null e1 then nil else if car e1 = car e2 then evlexcompless!?(cdr e1,cdr e2) else car e1 #> car e2; expr procedure evcomp (e1,e2); % /* Exponent vector compare. e1, e2 are exponent vectors in some % order. Evcomp(e1,e2) returns the digit 0 if exponent vector e1 is % equal exponent vector e2, the digit 1 if e1 is greater than e2, % else the digit -1. This function is assigned a value by the % ordering mechanism, so is dummy for now*/ apply(get(dipsortmode!*,'evcomp),list(e1,e2)); expr procedure evilcompless!?(e1,e2); % /* Exponent vector inverse lexicographical compare less. e1, e2 are % exponent vectors in lexicographical order. Evilcompless?(e1,e2) is % a boolean function which returns true if e1 is ordered less than e2*/ if null e1 then nil else if car e1 = car e2 then evilcompless!?(cdr e1,cdr e2) else car e1 #< car e2; expr procedure evlexcomp(e1,e2); % /* Exponent vector lexicographical compare. e1, e2 are exponent % vectors in lexicographical order. Evlexcomp(e1,e2) returns the % digit 0 if exponent vector e1 is equal exponent vector e2, 1 if e1 % is greater than e2, else the digit -1. */ if null e1 then 0 else if car e1 = car e2 then evlexcomp(cdr e1,cdr e2) else if car e1 #< car e2 then 1 else -1; expr procedure evilcomp (e1,e2); % /* Exponent vector inverse lexicographical compare. The % exponent vectors e1 and e2 are in inverse lexicographical % ordering. evilcomp(e1,e2) returns the digit 0 if exponent % vector e1 is equal exponent vector e2, the digit 1 if e1 is % greater than e2, else the digit -1. */ if null e1 then 0 else if car e1 = car e2 then evilcomp(cdr e1,cdr e2) else if car e1 #> car e2 then 1 else -1; expr procedure evitdcompless!?(e1,e2); % /* Exponent vector inverse total degree compare less. % The exponent vectors e1 and e2 are in inverse total degree % ordering. evitdcompless!?(e1,e2) is a boolean function that % returns true if exponent vector e1 is ordered less than e2*/ if null e1 then nil else if car e1 = car e2 then evitdcompless!?(cdr e1, cdr e2) else (if te1 = te2 then car e1 #< car e2 else te1 #< te2) where te1 = evtdeg e1, te2 = evtdeg e2; expr procedure evtdcompless!?(e1,e2); % /*Exponent vector total degree compare less.*/ if null e1 then nil else if car e1 = car e2 then evtdcompless!?(cdr e1,cdr e2) else (if te1 = te2 then car e1 #> car e2 else te1 #< te2) where te1 = evtdeg e1, te2 = evtdeg e2; expr procedure evitdcomp (e1,e2); % /* Exponent vector inverse total degree compare. % The exponent vectors e1 and e2 are in inverse total degree % ordering. evitdcomp(e1,e2) returns the digit 0 if exponent % vector e1 is equal exponent vector e2, the digit 1 if e1 is % greater than e2, else the digit -1. */ if null e1 then 0 else if car e1 = car e2 then evitdcomp(cdr e1, cdr e2) else (if te1 = te2 then if car e1 #> car e2 then 1 else -1 else if te1 #> te2 then 1 else -1) where te1 = evtdeg e1, te2 = evtdeg e2; expr procedure evtdcomp (e1,e2); % /* ... */ if null e1 then 0 else if car e1 = car e2 then evtdcomp(cdr e1,cdr e2) else (if te1 = te2 then if car e1 #< car e2 then 1 else -1 else if te1 #> te2 then 1 else -1) where te1 = evtdeg e1, te2 = evtdeg e2; expr procedure evtdeg e1; % /* Exponent vector total degree. e1 is an exponent vector. % evtdeg(e1) calculates the total degree of the exponent % e1 and returns an integer. */ (<<while e1 do <<x := car e1 #+ x; e1 := cdr e1>>; x>>) where x = 0; expr procedure evlcm (e1,e2); % /* Exponent vector least common multiple. e1 and e2 are % exponent vectors. evlcm(e1,e2) computes the least common % multiple of the exponent vectors e1 and e2, and returns % an exponent vector. */ % for each lpart in e1 each rpart in e2 collect % if lpart #> rpart then lpart else rpart; begin scalar x; while e1 do <<x := (if car e1 #> car e2 then car e1 else car e2) . x; e1 := cdr e1; e2 := cdr e2>>; return reversip x end; expr procedure evmtest!? (e1,e2); % /* Exponent vector multiple test. e1 and e2 are compatible exponent % vectors. evmtest!?(e1,e2) returns a boolean expression. % True if exponent vector e1 is a multiple of exponent % vector e2, else false. */ null e1 or not(car e1 #< car e2) and evmtest!?(cdr e1,cdr e2); expr procedure evsum (e1,e2); % /* Exponent vector sum. e1 and e2 are exponent vectors. % evsum(e1,e2) calculates the sum of the exponent vectors. % e1 and e2 componentwise and returns an exponent vector. */ % for each lpart in e1 each rpart in e2 collect lpart #+ rpart; begin scalar x; while e1 do <<x := (car e1 #+ car e2) . x; e1 := cdr e1; e2 := cdr e2>>; return reversip x end; expr procedure evdif (e1,e2); % /* Exponent vector difference. e1 and e2 are exponent % vectors. evdif(e1,e2) calculates the difference of the % exponent vectors e1 and e2 componentwise and returns an % exponent vector. */ % for each lpart in e1 each rpart in e2 collect lpart #- rpart; begin scalar x; while e1 do <<x := (car e1 #- car e2) . x; e1 := cdr e1; e2 := cdr e2>>; return reversip x end; expr procedure intevprod(n,e); % /* Multiplies each element of the exponent vector u by the integer n*/ for each x in e collect n #* x; expr procedure expvec2a e; % /* Returns list of prefix equivalents of exponent vector e*/ expvec2a1(e,dipvars!*); expr procedure expvec2a1(u,v); % /* Sub function of expvec2a */ if null u then nil else if car u = 0 then expvec2a1(cdr u,cdr v) else if car u = 1 then car v . expvec2a1(cdr u,cdr v) else list('expt,car v,car u) . expvec2a1(cdr u,cdr v); expr procedure dipevlpri(e,v); % /* Print exponent vector e in infix form. V is a boolean variable % which is true if an element in a product has preceded this one*/ dipevlpri1(e,dipvars!*,v); expr procedure dipevlpri1(e,u,v); % /* Sub function of dipevlpri */ if null e then nil else if car e = 0 then dipevlpri1(cdr e,cdr u,v) else <<if v then dipprin2 "*"; dipprin2 car u; if car e #> 1 then <<dipprin2 "**"; dipprin2 car e>>; dipevlpri1(cdr e,cdr u,t)>>; remprop('torder,'stat); expr procedure torder u; % algebraic mode interface to dipsortingmode. dipsortingmode car u; put('torder,'stat,'rlis); expr procedure dipsortingmode u; % /* Sets the exponent vector sorting mode. Returns the previous mode*/ if not idp u or not flagp(u,'dipsortmode) then typerr(u,"term ordering mode") else begin scalar x; x := dipsortmode!*; dipsortmode!* := u; return x end; flag('(lex invlex totaldegree invtotaldegree),'dipsortmode); put('lex,'evcompless!?,'evlexcompless!?); put('lex,'evcomp,'evlexcomp); put('invlex,'evcompless!?,'evilcompless!?); put('invlex,'evcomp,'evilcomp); put('invtotaldegree,'evcompless!?,'evitdcompless!?); put('invtotaldegree,'evcomp,'evitdcomp); put('totaldegree,'evcompless!?,'evtdcompless!?); put('totaldegree,'evcomp,'evtdcomp); dipsortingmode 'invlex; % /*Default value*/ endmodule; module dipoly; % /*Distributive polnomial algorithms*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ fluid '(dipvars!* dipzero); fexpr procedure polyin p; a2dip car p; expr procedure dipconst!? p; not dipzero!? p and dipzero!? dipmred p and evzero!? dipevlmon p; expr procedure dfcprint pl; % h polynomial factor list of distributive polynomials print. for each p in pl do dfcprintin p; expr procedure dfcprintin p; % factor with exponent print. ( if cdr p neq 1 then << prin2 " ( "; dipprint1(p1,nil); prin2 " )** "; prin2 cdr p; terprit 2 >> else << prin2 " "; dipprint p1>> ) where p1:= dipmonic a2dip prepf car p; expr procedure dfcprin p; % print content, factors and exponents of factorized polynomial p. << terpri(); prin2 " content of factorized polynomials = "; prin2 car p; terprit 2; dfcprint cdr p >>; expr procedure diplcm p; % Distributive polynomial least common multiple of denomiators. % p is a distributive rational polynomial. diplcm(p) calculates % the least common multiple of the denominators and returns a % base coefficient of the form 1/lcm(denom bc1,.....,denom bci). if dipzero!? p then mkbc(1,1) else bclcmd(diplbc p, diplcm dipmred p); expr procedure diprectoint(p,u); % Distributive polynomial conversion rational to integral. % p is a distributive rational polynomial, u is a base coefficient % ( 1/lcm denom p ). diprectoint(p,u) returns a primitive % associate pseudo integral ( denominators are 1 ) distributive % polynomial. if bczero!? u then dipzero else if bcone!? u then p else diprectoint1(p,u); expr procedure diprectoint1(p,u); % Distributive polynomial conversion rational to integral internall 1. % diprectoint1 is used in diprectoint. if dipzero!? p then dipzero else dipmoncomp(bclcmdprod(u,diplbc p),dipevlmon p, diprectoint1(dipmred p,u)); expr procedure dipresul(p1,p2); % test for fast downwards calculation % p1 and p2 are two bivariate distributive polynomials. % dipresul(p1,p2) returns the resultant of p1 and p2 with respect % respect to the first variable of the variable list (car dipvars!*). begin scalar q1,q2,q,ct; q1:=dip2a diprectoint(p1,diplcm p1); q2:=dip2a diprectoint(p2,diplcm p2); ct := time(); q:= a2dip prepsq simpresultant list(q1,q2,car dipvars!*); ct := time() - ct; prin2 " resultant : "; dipprint dipmonic q; terpri(); prin2 " time resultant : "; prin2 ct; terpri(); end; expr procedure dipbcprod (p,a); % /* Distributive polynomial base coefficient product. % p is a distributive polynomial, a is a base coefficient. % dipbcprod(p,a) computes p*a, a distributive polynomial. */ if bczero!? a then dipzero else if bcone!? a then p else dipbcprodin(p,a); expr procedure dipbcprodin (p,a); % /* Distributive polynomial base coefficient product internal. % p is a distributive polynomial, a is a base coefficient, % where a is not equal 0 and not equal 1. % dipbcprodin(p,a) computes p*a, a distributive polynomial. */ if dipzero!? p then dipzero else dipmoncomp(bcprod(a, diplbc p), dipevlmon p, dipbcprodin(dipmred p, a)); expr procedure dipdif (p1,p2); % /* Distributive polynomial difference. p1 and p2 are distributive % polynomials. dipdif(p1,p2) calculates the difference of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if dipzero!? p1 then dipneg p2 else if dipzero!? p2 then p1 else ( if sl = 1 then dipmoncomp(diplbc p1, ep1, dipdif(dipmred p1, p2) ) else if sl = -1 then dipmoncomp(bcneg diplbc p2, ep2, dipdif(p1,dipmred p2)) else ( if bczero!? al then dipdif(dipmred p1, dipmred p2) else dipmoncomp(al, ep1, dipdif(dipmred p1, dipmred p2) ) ) where al = bcdif(diplbc p1, diplbc p2) ) where sl = evcomp(ep1, ep2) where ep1 = dipevlmon p1, ep2 = dipevlmon p2; expr procedure diplength p; % /* Distributive polynomial length. p is a distributive % polynomial. diplength(p) returns the number of terms % of the distributive polynomial p, a digit.*/ if dipzero!? p then 0 else 1 + diplength dipmred p; expr procedure diplistsum pl; % /* Distributive polynomial list sum. pl is a list of distributive % polynomials. diplistsum(pl) calculates the sum of all polynomials % and returns a list of one distributive polynomial. */ if null pl or null cdr pl then pl else diplistsum(dipsum(car pl, cadr pl) . diplistsum cddr pl); expr procedure diplmerge (pl1,pl2); % /* Distributive polynomial list merge. pl1 and pl2 are lists % of distributive polynomials where pl1 and pl2 are in non % decreasing order. diplmerge(pl1,pl2) returns the merged % distributive polynomial list of pl1 and pl2. */ if null pl1 then pl2 else if null pl2 then pl1 else ( if sl >= 0 then cpl1 . diplmerge(cdr pl1, pl2) else cpl2 . diplmerge(cdr pl2, pl1) ) where sl = evcomp(ep1, ep2) where ep1 = dipevlmon cpl1, ep2 = dipevlmon cpl2 where cpl1 = car pl1, cpl2 = car pl2; expr procedure diplsort pl; % /* Distributive polynomial list sort. pl is a list of % distributive polynomials. diplsort(pl) returns the % sorted distributive polynomial list of pl. sort(pl, function dipevlcomp); expr procedure dipevlcomp (p1,p2); % /* Distributive polynomial exponent vector leading monomial % compare. p1 and p2 are distributive polynomials. % dipevlcomp(p1,p2) returns a boolean expression true if the % distributive polynomial p1 is smaller or equal the distributive % polynomial p2 else false. */ not evcompless!?(dipevlmon p1, dipevlmon p2); expr procedure dipmonic p; % /* Distributive polynomial monic. p is a distributive % polynomial. dipmonic(p) computes p/lbc(p) if p is % not equal dipzero and returns a distributive % polynomial, else dipmonic(p) returns dipzero. */ if dipzero!? p then p else dipbcprod(p, bcinv diplbc p); expr procedure dipneg p; % /* Distributive polynomial negative. p is a distributive % polynomial. dipneg(p) returns the negative of the distributive % polynomial p, a distributive polynomial. */ if dipzero!? p then p else dipmoncomp ( bcneg diplbc p, dipevlmon p, dipneg dipmred p ); expr procedure dipone!? p; % /* Distributive polynomial one. p is a distributive polynomial. % dipone!?(p) returns a boolean value. If p is the distributive % polynomial one then true else false. */ not dipzero!? p and dipzero!? dipmred p and evzero!? dipevlmon p and bcone!? diplbc p; expr procedure dippairsort pl; % /* Distributive polynomial list pair merge sort. pl is a list % of distributive polynomials. dippairsort(pl) returns the % list of merged and in non decreasing order sorted % distributive polynomials. */ if null pl or null cdr pl then pl else diplmerge(diplmerge( car(pl) . nil, cadr(pl) . nil ), dippairsort cddr pl); expr procedure dipprod (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if diplength p1 <= diplength p2 then dipprodin(p1, p2) else dipprodin(p2, p1); expr procedure dipprodin (p1,p2); % /* Distributive polynomial product internal. p1 and p2 are distrib % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if dipzero!? p1 or dipzero!? p2 then dipzero else ( dipmoncomp(bcprod(bp1, diplbc p2), evsum(ep1, dipevlmon p2), dipsum(dipprodin(dipfmon(bp1, ep1), dipmred p2), dipprodin(dipmred p1, p2) ) ) ) where bp1 = diplbc p1, ep1 = dipevlmon p1; expr procedure dipprodls (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial % using distributive polynomials list sum (diplistsum). */ if dipzero!? p1 or dipzero!? p2 then dipzero else car diplistsum if diplength p1 <= diplength p2 then dipprodlsin(p1, p2) else dipprodlsin(p2, p1); expr procedure dipprodlsin (p1,p2); % /* Distributive polynomial product. p1 and p2 are distributive % polynomials. dipprod(p1,p2) calculates the product of the % two distributive polynomials p1 and p2, a distributive polynomial % using distributive polynomials list sum (diplistsum). */ if dipzero!? p1 or dipzero!? p2 then nil else ( dipmoncomp(bcprod(bp1, diplbc p2), evsum(ep1, dipevlmon p2), car dipprodlsin(dipfmon(bp1, ep1), dipmred p2)) . dipprodlsin(dipmred p1, p2) ) where bp1 = diplbc p1, ep1 = dipevlmon p1; expr procedure dipsum (p1,p2); % /* Distributive polynomial sum. p1 and p2 are distributive % polynomials. dipsum(p1,p2) calculates the sum of the % two distributive polynomials p1 and p2, a distributive polynomial*/ if dipzero!? p1 then p2 else if dipzero!? p2 then p1 else ( if sl = 1 then dipmoncomp(diplbc p1, ep1, dipsum(dipmred p1, p2) ) else if sl = -1 then dipmoncomp(diplbc p2, ep2, dipsum(p1,dipmred p2)) else ( if bczero!? al then dipsum(dipmred p1, dipmred p2) else dipmoncomp(al, ep1, dipsum(dipmred p1, dipmred p2) ) ) where al = bcsum(diplbc p1, diplbc p2) ) where sl = evcomp(ep1, ep2) where ep1 = dipevlmon p1, ep2 = dipevlmon p2; endmodule; module dipvars; %/* Determine distributive polynomial variables in a prefix form*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ expr procedure dipvars u; % /* Returns list of variables in prefix form u*/ dipvars1(u,nil); expr procedure dipvars1(u,v); if atom u then if constantp u or u memq v then v else u . v else if idp car u and get(car u,'dipfn) then dipvarslist(cdr u,v) else if u memq v then v else u . v; expr procedure dipvarslist(u,v); if null u then v else dipvarslist(cdr u,union(dipvars car u,v)); endmodule; module a2dip; %/*Convert an algebraic (prefix) form to distributive polynomial*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ fluid '(dipvars!* dipzero); expr procedure a2dip u; % /*Converts the algebraic (prefix) form u to a distributive poly. % We assume that all variables used have been previously % defined in dipvars!*, but a check is also made for this*/ if atom u then a2dipatom u else if not atom car u or not idp car u then typerr(car u,"dipoly operator") else (if x then apply(x,list for each y in cdr u collect a2dip y) else a2dipatom u) where x = get(car u,'dipfn); expr procedure a2dipatom u; % /*Converts the atom (or kernel) u into a distributive polynomial*/ if u=0 then dipzero else if numberp u or not(u member dipvars!*) then dipfmon(a2bc u,evzero()) else dipfmon(a2bc 1,mkexpvec u); expr procedure dipfnsum u; % /*U is a list of dip expressions. Result is the distributive poly % representation for the sum*/ (<<for each y in cdr u do x := dipsum(x,y); x>>) where x = car u; put('plus,'dipfn,'dipfnsum); put('plus2,'dipfn,'dipfnsum); expr procedure dipfnprod u; % /*U is a list of dip expressions. Result is the distributive poly % representation for the product*/ % /*Maybe we should check for a zero*/ (<<for each y in cdr u do x := dipprod(x,y); x>>) where x = car u; put('times,'dipfn,'dipfnprod); put('times2,'dipfn,'dipfnprod); expr procedure dipfndif u; % /*U is a list of two dip expressions. Result is the distributive % polynomial representation for the difference*/ dipsum(car u,dipneg cadr u); put('difference,'dipfn,'dipfndif); expr procedure dipfnpow u; % /*U is a pair of dip expressions. Result is the distributive poly % representation for the first raised to the second power*/ (if not fixp n or n<0 then typerr(n,"distributive polynomial exponent") else if n=0 then if dipzero!? v then rederr "0**0 invalid" else w else if dipzero!? v or n=1 then v else if dipzero!? dipmred v then dipfmon(bcpow(diplbc v,n),intevprod(n,dipevlmon v)) else <<while n>0 do <<if not evenp n then w := dipprod(w,v); n := n/2; if n>0 then v := dipprod(v,v)>>; w>>) where n := dip2a cadr u, v := car u, w := dipfmon(a2bc 1,evzero()); put('expt,'dipfn,'dipfnpow); expr procedure dipfnneg u; % /*U is a list of one dip expression. Result is the distributive % polynomial representation for the negative*/ (if dipzero!? v then v else dipmoncomp(bcneg diplbc v,dipevlmon v,dipmred v)) where v = car u; put('minus,'dipfn,'dipfnneg); expr procedure dipfnquot u; % /*U is a list of two dip expressions. Result is the distributive % polynomial representation for the quotient*/ if dipzero!? cadr u or not dipzero!? dipmred cadr u or not evzero!? dipevlmon cadr u then typerr(dip2a cadr u,"distributive polynomial denominator") else dipfnquot1(car u,diplbc cadr u); expr procedure dipfnquot1(u,v); if dipzero!? u then u else dipmoncomp(bcquot(diplbc u,v), dipevlmon u, dipfnquot1(dipmred u,v)); put('quotient,'dipfn,'dipfnquot); endmodule; module dip2a; %/* Functions for converting distributive forms into prefix forms*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ expr procedure dip2a u; % /* Returns prefix equivalent of distributive polynomial u*/ if dipzero!? u then 0 else dipreplus dip2a1 u; expr procedure dip2a1 u; if dipzero!? u then nil else ((if bcminus!? x then list('minus,dipretimes(bc2a bcneg x . y)) else dipretimes(bc2a x . y)) where x = diplbc u, y = expvec2a dipevlmon u) . dip2a1 dipmred u; expr procedure dipreplus u; if atom u then u else if null cdr u then car u else 'plus . u; expr procedure dipretimes u; % /* U is a list of prefix expressions the first of which is a number. % Result is prefix representation for their product*/ if car u = 1 then if cdr u then dipretimes cdr u else 1 else if null cdr u then car u else 'times . u; endmodule; module dipprint; %/* printing routines for distributive polynomials*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ fluid '(dipvars!*); expr procedure diplprint u; % /* Prints a list of distributive polynomials using dipprint*/ for each v in u do dipprint v; expr procedure dipprint u; % /* Prints a distributive polynomial in infix form*/ <<terpri(); dipprint1(u,nil); terpri(); terpri()>>; expr procedure dipprint1(u,v); % /* Prints a distributive polynomial in infix form. % U is a distributive form. V is a flag which is true if a term % has preceded current form*/ if dipzero!? u then if null v then dipprin2 0 else nil else begin scalar bool,w; w := diplbc u; if bcminus!? w then <<bool := t; w := bcneg w>>; if bool then dipprin2 " - " else if v then dipprin2 " + "; (if not bcone!? w or evzero!? x then <<bcprin w; dipevlpri(x,t)>> else dipevlpri(x,nil)) where x = dipevlmon u; dipprint1(dipmred u,t) end; expr procedure dipprin2 u; % /* Prints u, preceding by two EOL's if we have reached column 70*/ <<if posn()>69 then <<terpri(); terpri()>>; prin2 u>>; endmodule; module grinterf; % Interface of Groebner package to REDUCE. % /*Authors: R. Gebauer, A. C. Hearn, M. Moeller. fluid '(dipvars!*); expr procedure groebnereval u; begin integer n; n := length u; if n=1 then return groebner(reval car u,nil) else if n neq 2 then rederr "GROEBNER called with wrong number of arguments" else return groebner(reval car u,reval cadr u) end; put('groebner,'psopfn,'groebnereval); expr procedure greduce u; % /* Polynomial reduction modulo a Groebner basis driver. u is an % expression and v a list of expressions. Greduce calculates the % polynomial u reduced wrt the list of expressions v reduced to a % groebner basis modulo using the optional third argument as the % order of variables. begin integer n; scalar dipvars!*,v; n := length u; v := for each j in getrlist reval cadr u collect if eqexpr j then !*eqn2a j else j; if n=2 then dipvars!* := for each j in gvarlis v collect !*a2k j else if n=3 then dipvars!* := getrlist caddr u else rederr "GREDUCE called with wrong number of arguments"; v := groebner2 for each j in v collect a2dip j; if cdr v then errach list("Groebner",u) else if null cdar v and dip2a caar v = 1 then rederr "Inconsistent Basis"; return dip2a dipnorform(car v,a2dip reval car u) end; put('greduce,'psopfn,'greduce); expr procedure preduce(u,v); % /* Polynomial reduction driver. u is an expression and v a list of % expressions. Preduce calculates the polynomial u reduced wrt the list % of expressions v. */ begin scalar dipvars!*; v := for each j in getrlist reval v collect if eqexpr j then !*eqn2a j else j; dipvars!* := for each j in gvarlis v collect !*a2k j; return dip2a dipnorform(for each j in v collect a2dip j, a2dip reval u) end; flag('(preduce),'opfn); endmodule; module groebner; % Basic Groebner base code using Buchberger algorithm. % /*Authors: R. Gebauer, A. C. Hearn, M. Moeller. fluid '(!*groebopt !*groebfac !*hopt !*trgroeb !*trgroebs !*trgroeb0 !*trgroeb1 dipvars!* dipzero); switch groebopt,groebfac,hopt,trgroeb,trgroebs,trgroeb0,trgroeb1; % /* 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 ' trgroeb0' factorizes the S - polynom */ % /* the S - polynom will not be */ % /* replaced by a factor */ % /* option ' trgroebs ' prints S - polynomials */ % /* on the output file */ % /* option ' hopt ' the H- polynomials are */ % /* optimised using resultant */ % /* and factorisation method */ % /* option ' groebfac ' the H - polynomials are */ % /* factorized. If a H - polynom */ % /* could be factorized, new sub- */ % /* problems are generated and */ % /* option ' fac ' is set to off */ % /* NOTE: this option is not complete */ % /* at the moment and has been suppressed */ % expr procedure bas p; diplprint car groebner(p,dipvars!*); expr procedure groebner(u,v); % /* 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. Groebner(p) calculates the Groebner basis of the % expressions wrt the variables. */ begin scalar dipvars!*,w; w := for each j in getrlist u collect if eqexpr j then !*eqn2a j else j; if null w then rederr "Empty list in Groebner" else if null cdr w then return 'list . w; if null v then v := gvarlis w else v := getrlist v; dipvars!* := for each j in v collect !*a2k j; w := groebner2 for each j in w collect a2dip j; if cdr w then errach list("Groebner",u,dipvars!*); return 'list . for each j in car w collect dip2a j end; expr procedure gvarlis u; % Finds variables (kernels) in the list of expressions u. gvarlis1(u,nil); expr procedure gvarlis1(u,v); if null u then v else union(gvar1(car u,v),gvarlis1(cdr u,v)); expr procedure gvar1(u,v); if null u or numberp u then v else if atom u then if u member v then v else u . v else if car u memq '(plus times expt difference minus) then gvarlis1(cdr u,v) else if car u eq 'quotient then gvar1(cadr u,v) else if u member v then v else u . v; expr procedure groebner2 p; begin scalar tim,spac,spac1,p1; tim := time(); % terprit 3; spac := gctime(); p1:= dipgbase p; spac1 := gctime() - spac; % prin2 " garbage collection time for test "; % prin2 spac1; % prin2 "( not yet available )"; if !*trgroeb then <<prin2 "Computing time for test "; prin2(time() - tim - spac1); % prin2(time() - tim ); prin2t " milliseconds ">>; return p1 end; expr procedure dipindexpol(pl,n); % Distributive polynomial index list. pl is a list of distributive % polynomials; n is an index, an integer. dipindexpol(pl,n) % returns a list of distributive polynomials in the form % ( (n,p1) (n+1,p2) ..... (n+k,pk) ). if null pl then pl else list(n,car pl) . dipindexpol(cdr pl, n + 1); expr procedure dipindexpolspec pl; % Distributive polynomial special list. pl is a list produced % by dipindexpol. dipindexpolspec pl constructs a list of lists % of polynomials in the form ( (p1,.....,pl) (p2,.....,pl).... % ..(pl-1,,pl) (pl) ). for each pl0 on pl collect ( for each pl1 in pl0 collect pl1 ); expr procedure dipcpairlistopt pl; % Distributive critical pair list optimise. pl is a special list % ( constructed by dipcpairlist ) of elements used in the % Groebner calculation. dipcpairlistopt(pl) returns a list which % is optimised using Buchberger criterion 4. if pl then ( if buchcrit4(caddr x, cadddr x, cadr x) then x . dipcpairlistopt cdr pl else dipcpairlistopt cdr pl ) where x = car pl else nil; expr procedure dipcpairlistop(d,d0); % Distributive polynomial critical pair list optimise. % dipcpairlistop(d,d0) returns an optimised critical pair % starting list using the new criteria's. begin scalar x; while d do << x:= dipcpairlistopt1(cadar d,d0,d0); d0:= x; d:= cdr d>>; return x end; expr procedure dipcpairlistopt1(h,d,d0); % Distributive polynomial critical pair list optimise version 1. % dipcpairlistopt1(h,d,d0) returns an optimised critical pair % list. if null d then d0 else ( if evmtest!?(cadar d,ev1) then dipcpairlistopt1(h, cdr d,x) else dipcpairlistopt1(h,cdr d,d0) ) where x= dipcpairlistopt1in(ev1,cadar d,car d,d0) where ev1 = dipevlmon h; expr procedure dipcpairlistopt1in(ev1,ev2,id1,d); % Distributive polynomial critical pair list optimise version 1. % internall. dipcpairlistopt1in is used in dipcpairlistopt1. if ev2 neq evlcm(ev1,dipevlmon caddr id1) and ev2 neq evlcm(ev1,dipevlmon cadddr id1) then dipcpairlistopt1in1(id1,d) else d; expr procedure dipcpairlistopt1in1(d1,d); % Distributive polynomial critical pair list optimise version 1 % internall version 1. dipcpairlistopt1in1 is used in % dipcpairlistopt1in. if null d then nil else if d1 eq car d then dipcpairlistopt1in1(d1,cdr d) else car d . dipcpairlistopt1in1(d1,cdr d); expr procedure dipindexpolrec pl; % Distributive index polynom list reconstruct. pl is a list of % polynomials used in the Groebner calculation. dipindexpolrec(pl) % returns a list of distributive polynomials. for each p in pl collect cadr p; expr procedure dipcplist pl; % Distributive polynomial critical pair list construct. % dipcplist returns a list of elements where an element has the % structure ( (ipi,ipj) lcm(epi,epj) pi pj ). % where ipi is the index of polynomial i, epi is the headterm of % the polynomial pi. for each p in pl conc ( dipcplistopt2(nil, dipcplistin(ep, pi1, reverse cdr p)) ) where ep = dipevlmon cadr pi1 where pi1 = car p; expr procedure dipcplistin(ep,p1,pl); % Distributive polynomial critical pair list construct internall. % dipcplistin is used in dipcplist. if null pl then pl else ( list(list(car p1,car p2), evlcm(ep,dipevlmon cadr p2), cadr p1, cadr p2) . dipcplistin(ep, p1, cdr pl) ) where p2 = car pl; expr procedure dipcplistadd(ind,p,pl); % Distributive polynomial critical pair list add. % dipcplistadd returns a new critical pair list where all % combinations of p with pl are added. if null pl then pl else ( list(list(car ps,ind),evlcm(dipevlmon p1, dipevlmon p),p1,p) . dipcplistadd(ind,p,cdr pl) ) where p1 = cadr ps where ps = car pl; expr procedure dipcplistopt2in(p1,pl); % Distributive polynomial critical pair list optimise version 2 % internall use. dipcplistopt2in(pl1,pl) is used in % dipcplistopt2. if null pl then dipzero else ( if evmtest!?(cadr p1, cadr p) then dipcplistopt2in1(p1,p) else dipcplistopt2in(p1,cdr pl) ) where p = car pl; expr procedure dipcplistopt2in1(p1,p2); % Distributive polynomial critical pair list optimise version 2 % internall use version 1. dipcplistopt2in1(pl1,pl) is used in % dipcplistopt2in. if cadr p1 = cadr p2 then ( if evilcompless!?(reverse car p1, reverse car p2) then p1 else p2 ) else p2; expr procedure dipindexpoloptin(p1,pl); % Distributive index polynomial list optimise internall use. % dipindexpoloptin is used in dipindexpolopt. if null pl then dipzero else ( if evmtest!?(dipevlmon cadr p1, dipevlmon cadr p) then dipindexpoloptin1(p1,p) else dipindexpoloptin(p1,cdr pl) ) where p = car pl; expr procedure dipindexpoloptin1(p1,p2); % Distributive index polynomial list optimise internall version 1. % dipindexpoloptin1 is used in dipindexpoloptin. if dipevlmon cadr p1 = dipevlmon cadr p2 then ( if car p1 < car p2 then p1 else p2 ) else p2; expr procedure dipcplistopt2(pl1,pl2); % Distributive polynomial critical pair list optimise version 2. % dipcplistopt2(pl1,pl2) returns the optimised critical pair list. if null pl2 then pl1 else ( if dipzero!? dipcplistopt2in(p,pl1) and dipzero!? dipcplistopt2in(p,pl0) then dipcplistopt2(cons(p,pl1),pl0) else dipcplistopt2(pl1,pl0) ) where p = car pl2, pl0 = cdr pl2; expr procedure dipindexpolopt(pl1,pl2); % Distributive index polynomial list optimise. pl1 and pl2 % are lists of polynomials used in the Groebner calculation. % dipindexpolopt(pl1,pl2) returns an optimised list of polynomials. if null pl2 then pl1 else ( if dipzero!? dipindexpoloptin(p,pl1) and dipzero!? dipindexpoloptin(p,pl0) then dipindexpolopt(cons(p,pl1),pl0) else dipindexpolopt(pl1,pl0) ) where p = car pl2, pl0 = cdr pl2; expr procedure dipcplistsort pl; % Distributive polynomial critical pair list sort. pl is a % special list for Groebner calculation. dipcplistsort(pl) % returns the sorted list pl; begin scalar tree; if null pl then return nil; tree := list(car pl,nil); while pairp(pl:= cdr pl) do dipcplistsortadd(car pl,tree); return tree2list(tree,nil) end; smacro procedure dipcplistevlcomp(p1,p2); % Distributive polynomial critical pair list exponent vector % compare. p1 and p2 are elements of the critical pair list. % dipcplistevlcomp(p1,p2) returns a boolean expression, true % if exponent vector of p1 is smaller or equal exponent vector % of p2 else false. evcompless!?(cadr p1, cadr p2); expr procedure dipcplistsortadd(item,node); % Distributive polynomial critical pair list sort addition. % add item to a node, using dipcplistevlcomp as an order % predicate. if dipcplistevlcomp(item, car node) then if cadr node then dipcplistsortadd(item, cadr node) else rplaca(cdr node,list(item,nil)) else if cddr node then dipcplistsortadd(item,cddr node) else rplacd(cdr node,list(item,nil)); expr procedure dipcplistmerge(pl1,pl2); % Distributive polynomial critical pair list merge. pl1 and pl2 % are critical pair lists used in the Groebner calculation. % dipcplistmerge(pl1,pl2) returns the merged list. if null pl1 then pl2 else if null pl2 then pl1 else ( if sl then cpl1 . dipcplistmerge(cdr pl1,pl2) else cpl2 . dipcplistmerge(pl1,cdr pl2) ) where sl = evcompless!?(cadr cpl1, cadr cpl2) where cpl1 = car pl1, cpl2 = car pl2; expr procedure buchcrit4(p1,p2,e); % Buchberger criterion 4. p1 and p2 are distributive % polynomials. e is the least common multiple of % the leading exponent vectors of the distributive % polynomials p1 and p2. buchcrit4(p1,p2,e) returns a % boolean expression. True if the reduction of the % distributive polynomials p1 and p2 is necessary % else false. e neq evsum( dipevlmon p1, dipevlmon p2); expr procedure dipgbase pl; % /* Distributive polynomial Groebner base. pl is a list of distributiv % polynomials. dipgbase(pl) calculates the Groebner base of the list % of distributive polynomials pl and returns a list of distributive % polynomials. */ if null pl then nil else if null cdr pl then list pl else if !*groebopt then dipgbasein dipvordopt pl else dipgbasein pl; expr procedure gbprint pl; % Groebner basis list of distributive polynomials print. for each p in pl do dipprint dipmonic p; expr procedure rescheck!?(a,h1,vl); length h1 = a and car h1 = vl - 1; expr procedure rescheck1!?(a,h1,vl); length h1 = a and car h1 = vl - 2 and cadr h1 = vl - 1; expr procedure newhpol(p1,p2,x); begin scalar q1,q2,q; q1:=dip2a diprectoint(p1,diplcm p1); q2:=dip2a diprectoint(p2,diplcm p2); q:=a2dip prepsq simpresultant list(q1,q2,x); return q; end; expr procedure sqpol p1; begin scalar q1,q; q1:=dip2a diprectoint(p1,diplcm p1); q:=a2dip caar sqfrf q1; return q; end; expr procedure dipnorfor (pl,p); % /* Distributive polynom normalform. pl is a list of distributive % polynomials, p is a distributive polynomial. dipnorfor(pl,p) % calculates a distributive polynomial such that the powerproduct % of the distributive % polynomial p is reducible to this modulo the distributive % polynomial list pl and is in normalform with respect to the % distributive polynomial p and returns a distributive polynomial. */ if dipzero!? p or null pl then p else ( if dipzero!? q then p else ( if dipzero!? rq then dipnorfor(pl,dipmred p) else dipnorfor(pl, dipdif(dipmred p, dipprod(rq, dipfmon(bcquot(diplbc p, diplbc q), evdif(ep, dipevlmon q) ) ) ) ) ) where rq = dipmred q ) where q = dipnorformsel(ep, pl) where ep = dipevlmon p; expr procedure dipmingbase pl; % Distributive polynomial minimal ordered Groebner base. pl is a % list of distributive polynomials. dipmingbase(pl) calculates % the minimal normed and ordered Groebner base of the distributive % polyomials pl and returns a list of distributive polyomials. if null cdr pl then pl else dipmingbasein2(nil,dipmingbasein1(nil,pl) ); expr procedure dipgbasein ql; % /* Distributive polynomial Groebner base. pl is a list of distributiv % polynomials. dipgbase(pl) calculates the Groebner base of the list % of distributive polynomials pl and returns a list of distributive % polynomials. */ begin scalar ql0,u,ql1,w,d,ql22,lql1,ql11,lv,h1h0,d1,d0,p1, sp0,n,dl,p2,ct1,sp,h,ct11,h1,h10,hs1,h1h1,h0,hs2; u := 1; w := 1; n := 1; ql0 := nil; ql1:= dipindexpol(ql,1); d:= dipcplistsort dipcpairlistopt dipcplist dipindexpolspec ql1; ql22 := ql; lql1:= length ql1; ql11:=dipindexpolopt(nil, ql1); d:=dipcpairlistop(ql11,d); if !*hopt then << lv:=length dipvars!*; h1h0:=nil>>; d1:=list list(lql1,ql1,ql11,ql22,d); if !*trgroeb1 then << prin2 " list d1 = "; prin2 d1; terpri(); prin2 length d1; terpri() >>; while not null d1 do << d0:= car d1; d1:= cdr d1; lql1:= car d0; ql1:= cadr d0; ql11:= caddr d0; ql22:= cadddr d0; d:= cadddr cdr d0; while not null d do << dl:= car d; d := cdr d; p1:= caddr dl; p2:= cadddr dl; if !*trgroeb then << ct1 := time() >>; sp := dipspolynom(p1,p2); if !*trgroebs then << prin2t "S-polynom:"; dipprint sp; terpri() >>; if !*trgroeb0 then << sp0:= dip2a diprectoint(sp,diplcm sp); sp0:= factorf !*q2f simp sp0; dfcprin sp0; terprit 2 >>; h := dipnorform(ql22, sp); if !*trgroeb then << ct11 := time() - ct1 >>; if dipzero!? h then << if !*trgroeb then << terprit 2; printb 57; terpri(); prin2 " / reduction of polynom "; prin2 caar dl; prin2 " and "; prin2 cadar dl; prin2 " leads to 0 "; prin2 " ( "; prin2 ct11; prin2 " ms )"; terpri(); printb 57; terprit 2 >> >>; if not dipzero!? h then if dipconst!? h then << ql11:= list list(lql1,dipmonic h); d:=nil >> else << h1 := dipmonic h; lql1:= lql1 + 1; if !*trgroeb then << prin2 "h-polynom "; prin2 lql1; prin2 " pair"; prin2 " ( "; prin2 caar dl; prin2 ","; prin2 cadar dl; prin2t " ) :"; dipprint h1; terpri(); prin2 " computing time for h-polynom "; prin2 ct11; terprit 3 >>; % The following option has been suppressed since it is not % complete. if nil and !*groebfac and u = 1 then << h10:= h1; h1:= dip2a diprectoint(h1,diplcm h1); h1:= factorf !*q2f simp h1; hs1:= reverse diplsort makdiplist cdr h1; if !*trgroeb then << prin2 "h-polynom factorized: "; terpri(); dfcprin h1; terpri() >>; h1:= dipmonic car hs1; hs1:= reverse cdr hs1; if not dipzero!? (dipdif(h1,h10)) then << u:= 0 >>; if !*trgroeb then << prin2 " new h-polynom "; terprit 3; dipprint h1; terprit 2 >> >>; if !*hopt and w = 1 then << h1h1:= indexcpl(evsum0(lv,h1),1); if !*trgroeb then << prin2 " index: "; prin2 h1h1; terpri(); prin2 " index: "; prin2 h1h0; terprit 3 >>; if h1h1 = h1h0 and rescheck!?(2,h1h0,lv) then << hs2:= reverse diplsort newhpo(h1,h0,cadr reverse dipvars!*); w:= 0>>; if h1h1 = h1h0 and rescheck1!?(2,h1h0,lv) then << hs2:= reverse diplsort newhpo(h1,h0,caddr reverse dipvars!*); w:= 0 >>; if null hs2 then << w:= 1 >> >>; if u = 0 and not null hs1 then << d0:= maklistd1(hs1,lql1,ql1,ql11,ql22,d); u:= 2; d1:=nconc(d0,d1) >>; %%%%%%% u:= 1; d1:=nconc(d0,d1) >>; d:= dipcpairlistopt1(h1,d,d); if !*trgroeb then << terpri(); prin2 "Restpairs: "; prin2t length d; terpri() >>; d:= dipcplistmerge(dipcplistsort dipcpairlistopt dipcplistopt2(nil,dipcplistadd(lql1,h1,ql11)),d); if !*hopt and w = 1 then << h1h0:=indexcpl(evsum0(lv,h1),1); h0:= h1 >>; ql11:= nconc(list list(lql1,h1),ql11); ql22:= nconc(list(h1),ql22); ql11:= dipindexpolopt(nil,ql11); if !*trgroeb1 then << prin2 " *** d = "; prin2 d; terpri(); prin2t " ql11 "; prin2 ql11; terpri() >>; if w = 0 then << h1:= dipmonic car hs2; hs2:= reverse cdr hs2; lql1:= lql1 + 1; if not null hs2 then << d0:= maklistd1(hs2,lql1,ql1,ql11,ql22,d); w:= 2; d1:= nconc(d0,d1) >>; d:= dipcpairlistopt1(h1,d,d); d:= dipcplistmerge(dipcplistsort dipcpairlistopt dipcplistopt2(nil,dipcplistadd(lql1,h1,ql11)),d); ql11:= nconc(list list(lql1,h1),ql11); ql22:= nconc(list(h1),ql22); ql11:= dipindexpolopt(nil,ql11); if !*trgroeb1 then << prin2 " *** d = "; prin2 d; terpri(); prin2t " ql11 "; prin2 ql11; terpri() >> >> >> >>; ql11:=dipindexpolrec ql11; if !*trgroeb then << prin2t " calculation now in final reduction "; terpri(); ct1 := time() >>; ql:=dipmingbase diplsort ql11; if !*trgroeb then << ct11 := time() - ct1; prin2 " computing time for final calculation "; prin2 ct11; prin2 " milliseconds "; terprit 3; prin2 " Number of Groebner Basis Polynomials := "; prin2t length ql; terprit 2; if n = 1 and null d1 then << prin2t " The Groebner Basis Polynomials "; terpri() >> else << prin2 " The Groebner Basis Polynomials ( Factor "; prin2 n; prin2t " )"; terpri(); n:= n + 1 >>; gbprint ql; if not null d1 then << prin2 " Calculation for Factor "; prin2t n; terprit 4 >> >>; ql0:= ql . ql0 >>; return ql0 end; expr procedure makdiplist pl; % Make list of distributive polynomials from list of polynomials pl. for each p in pl collect a2dip prepf car p; expr procedure terprit n; % print blank lines. for i:=1:n do << terpri() >>; expr procedure printb n; % print special sign ( - ). for i:=1:n do << prin2 "-" >>; expr procedure newhpo(h1,h0,x); % new h-polynom calculation. newhpo(h1,h2,x) calculates % the resultant of the two distributive polynomials h1 and h0 % with respect to x. begin scalar ct00,hh,hh1,hs2; if !*trgroeb then << ct00:= time() >>; hh:= dipmonic newhpol(h1,h0,x); if !*trgroeb then << prin2 " resultant "; terprit 2; dipprint hh; terprit 4 >>; hs2:= nil; if not dipzero!? hh then << hh1:= dip2a diprectoint(hh,diplcm hh); hh1:= factorf !*q2f simp hh1; if !*trgroeb then << prin2 " resultant factorized: "; terprit 2; dfcprin hh1; terprit 2; ct00:= time() - ct00; prin2 " special time for h: "; prin2 ct00; terpri() >>; hs2:= makdiplist cdr hh1 >>; return hs2 end; expr procedure maklistd1(x1,x2,x3,x4,x5,x6); % make list d1. save part time problems. begin scalar x,h1; while x1 do << h1:= car x1; x1:= cdr x1; x:= list(x2,x3, (dipindexpolopt(nil,nconc(list list(x2,h1),x4))), (nconc(list h1,x5)), (dipcplistmerge(dipcplistsort dipcpairlistopt dipcplistopt2(nil,dipcplistadd(x2,h1,x4)), dipcpairlistopt1(h1,x6,x6)))) . x >>; return x end; expr procedure dipmingbasein1 (pl1,pl2); % /* Distributive polynomial minimal ordered Groebner base internal1. % pl1 and pl2 are lists of distributive polynomials. % dipmingbasein1(pl1,pl2) is used in dipmingbase and returns a list % of distributive polynomials. */ if null pl2 then pl1 else ( if dipzero!? dipnorformsel(ep, pl1) and dipzero!? dipnorformsel(ep,cpl2) then dipmingbasein1( cons(p, pl1), cpl2) else dipmingbasein1( pl1, cpl2) ) where ep = dipevlmon p, cpl2 = cdr pl2 where p = car pl2; expr procedure dipmingbasein2 (pl1,pl2); % /* Distributive polynomial minimal ordered Groebner base internal2. % pl1 and pl2 are lists of distributive polynomials. % dipmingbasein2(pl1,pl2) is used in dipmingbase and returns a list % of distributive polynomials. */ if null pl2 then pl1 else ( dipmingbasein2(dipnorform(pl1,dipnorform(rp, p)) . pl1, rp) ) where p = car pl2, rp = cdr pl2; expr procedure dipnorform (pl,p); % /* Distributive polynom normalform. pl is a list of distributive % polynomials, p is a distributive polynomial. dipnorform(pl,p) % calculates a distributive polynomial such that the distributive % polynomial p is reducible to this modulo the distributive % polynomial list pl and is in normalform with respect to the % distributive polynomial p and returns a distributive polynomial. */ if dipzero!? p or null pl then p else ( if dipzero!? q then dipmoncomp(diplbc p, ep, dipnorform(pl, dipmred p) ) else ( if dipzero!? rq then dipnorform(pl, dipmred p) else dipnorform(pl, dipdif(dipmred p, dipprod(rq, dipfmon(bcquot(diplbc p, diplbc q), evdif(ep, dipevlmon q) ) ) ) ) ) where rq = dipmred q ) where q = dipnorformsel(ep, pl) where ep = dipevlmon p; expr procedure dipnorformsel (ep,pl); % /* Distributive polynom normalform select. ep is an exponent vector % of a distributive polynomial. pl is a list of distributive % polynomials. dipnorformsel(ep,pl) returns a distributive % polynomial of pl where ep is a multiple of the leading % exponent vector else dipzero. */ if null pl then dipzero else ( if evmtest!?(ep, dipevlmon q) then q else dipnorformsel(ep, cdr pl) ) where q = car pl; expr procedure dipspolynom (p1,p2); % /* Distributive polynom S polynom. p1 and p2 are distributive % polynomials. dipspolynom(p1,p2) calculates the S polynom of the % distributive polynomials p1 and p2 and returns a distributive % polynomial. */ if dipzero!? p1 or dipzero!? p2 then dipzero else ( if dipzero!? rp1 and dipzero!? rp2 then rp1 else ( if dipzero!? rp1 then dipprod(rp2, dipfmon(bcneg diplbc p1, evdif(ep, ep2) ) ) else if dipzero!? rp2 then dipprod(rp1, dipfmon(diplbc p2, evdif(ep, ep1) ) ) else dipdif( dipprod(rp2, dipfmon(diplbc p1, evdif(ep, ep2) ) ), dipprod(rp1, dipfmon(diplbc p2, evdif(ep, ep1) ) ) ) ) where ep = evlcm(ep1, ep2) where ep1 = dipevlmon p1, ep2 = dipevlmon p2 ) where rp1 = dipmred p1, rp2 = dipmred p2; expr procedure delqip1(u,v); if pairp cdr v then if u eq cadr v then rplacd(v,cddr v) else delqip1(u,cdr v); expr procedure delqip(u,v); % /*Destructive delete of first occurrence of u in v*/ if not pairp v then v else if u eq car v then cdr v else <<delqip1(u,v); v>>; endmodule; module dipopt; % /* Authors: R. Gebauer, A. C. Hearn, H. Kredel */ fluid '(!*trbas dipvars!*); %define ezero = 'nil; fluid '(dipzero ezero); %/*Until we understand how to define something to nil*/ expr procedure dipoptmat1 (el,dpl); % /* Distributive optimisation matrix subfunction 1. el is an % exponent vector, dpl is a degree matix. dipoptmat1(el,dpl) % returns the addition of el to dpl. */ if null el then dpl else dipsum ( dipfmon (bcfi 1, evcons(evfirst el, ezero)), car dpl) . dipoptmat1(evred el, cdr dpl); expr procedure dipoptmat2 (p,pl); % /* Distributive optimisation matrix subfunction 2. p is a % distributive polynomial, pl is a list of distributive % polynomials. dipoptmat1 is used. */ if dipzero!? p then pl else dipoptmat2(dipmred p, dipoptmat1(dipevlmon p, pl)); expr procedure dipoptmat3 (p,pl); % /* Distributive optimisation matrix subfunction 3. p is a % distributive polynomial, pl is a list of distributive % polynomials. dipoptmat2 is used. */ if null p then pl else dipoptmat3(cdr p, dipoptmat2(car p, pl)); expr procedure dipoptmat pl; % /* Distributive optimisation matrix. pl is a list of distributive % polynomials. dipoptmat(pl) returns the optimisation matrix % ( a degree matrix ) of pl, a list of univariate distributive % polynomials. */ if null pl then nil else dipoptmat3(pl, for each x in dipvars!* collect dipzero); expr procedure dipless!? (p1,p2); % /* Distributive polynomial less. p1 and p2 are distributive % polynomials. dipless!?(p1,p2) returns a boolean expression, % true if p1 is less than p2 else false. */ if dipzero!? p1 and dipzero!? p2 then nil else if not dipzero!? p1 then if not dipzero!? p2 then ( if sl < 0 then t else if sl > 0 then nil else ( if bl < 0 then t else if bl > 0 then nil else dipless!?(dipmred p1, dipmred p2) ) where bl = bccomp(diplbc p1, diplbc p2) ) where sl = evcomp(dipevlmon p1, dipevlmon p2) else t else nil; expr procedure pvdema pl; % /* Permutation vector degree matrix. pl is a list of univariate % polynomials in distributive representation. pvdema(pl) returns % a list ( indexlist ) where the elements are digits.*/ pvdema2 sort(pvdema1(pl, 1), 'pvdema3); expr procedure pvdema1(pl,n); % /* Permutation vector degree matrix subfunction 1. pl is a list % of univariate distributive polynomials, n is a digit. % pvdema1 changes the internal structure ( add index for each % polynomial ) and is used in pvdema. */ if null pl then pl else list(car pl, n) . pvdema1(cdr pl, n + 1); expr procedure pvdema2(pl); % /* Permutation vector degree matrix subfunction 2. pl is a list of % univariate distributive polynomials. pvdema2(pl) changes the % internal structure ( delete index for each polynomial ) and % is used in pvdema. */ if null pl then pl else nconc(cdar pl, pvdema2(cdr pl)); expr procedure pvdema3 (p1,p2); % /* Permutation vector degree matrix subfunction 3. p1 and p2 are % distributive univariate polynomials. pvdema3(p1,p2) returns % a boolean expression, true if the distributive polynomial p1 % is less than the distributive polynomial p2 else false. */ dipless!?(car p1, car p2); expr procedure listperm (v,n); % /* List permutation. v is a list ( any kind ) and n is an indexlist. % listperm(v,n) permutates v in respect to n and returns a % permutated list v. */ if null n then nil else nth(v, car n) . listperm(v, cdr n); expr procedure dipreorder (p,n); % /* Distributive polynomial reorder. p is a distributive polynomial, % n is an indexlist. dipreorder(p,n) reorders the exponent vectors % of each term of p in respect to the indexlist n and returns a % distributive polynomial. */ if dipzero!? p then nil else dipsum(dipfmon(diplbc p, evperm(dipevlmon p, n)), dipreorder(dipmred p, n)); expr procedure diplreorder (pl,n); % /* Distributive polynomial list reorder. pl is a list of distributive % polynomials and n is an indexlist. diplreorder(pl,n) reorders the % exponent vectors of each term of each polynomial in the list pl in % respect to the indexlist n and returns a list of distributive % polynomials.*/ for each x in pl collect dipreorder(x, n); expr procedure dipvordopt pl; % /* Distributive polynomial variable ordering optimisation. % pl is a list of distributive polynomials. dipvordopt(pl) % calculates the " optimal representation " and returns a list % of distributive polynomials. % NOTE: dipvordopt can change the global variable list dipvars!* */ begin scalar n,olddipvars,pl1; n := pvdema diopmatin pl; if !*trbas then << prin2t " The new index list :"; terprit 2; prin2t n; terprit 2 >>; olddipvars := dipvars!*; dipvars!* := listperm(dipvars!*, n); if !*trbas then << prin2t " The new variable list :"; terprit 2; prin2t dipvars!*; terprit 2 >>; pl1 := diplreorder(pl, n); if !*trbas then << prin2t " The new polynomial list :"; terprit 2; diplprint pl1; terprit 2 >>; % dipvars!* := olddipvars; return pl1 end; expr procedure diopmatin pl; % print univariate polynomials. begin scalar n1; << if !*trbas then << prin2t " The variable list :"; terprit 2; prin2t dipvars!*; terprit 2; prin2t " The univariate polynomials in each variable :"; terprit 2 >>; n1:=dipoptmat pl; if !*trbas then << dioprin(n1,dipvars!*) >> >>; return n1 end; expr procedure dioprin(pl,d); % print variables. begin scalar dipvars!*; for each x in pair(pl,d) do << dipvars!* := list cdr x; dipprint car x >> end; endmodule; end;