Artifact 1276f1e4b066bfe9bae281809104055782a6e10264128b053a308fa32ecabeb1:
- Executable file
r36/src/dipoly.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: 95301) [annotate] [blame] [check-ins using] [more...]
module dipoly; % Header module for dipoly package. %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ % % Significant modifications: H. Melenk. % %% Modifications: % % 14-Dec-1994 (HM): Term order GRADED added. % % 17-Sep-1994 (HM): The ideal variables are now declared in the TORDER % statement. The calling conventions can be still % used, but are removed from the documents. % % 12-Sep-1994 (HM): Make the base coefficient arithmatic call subs2 if % the switch *bcsub2 is on. This is turned on if % there are roots in the coefficient domain. Without % subs2 the zero detection would be incomplete in % such cases. % Term order MATRIX added. % % 5-Jun-1994 (HM): Introduced zero divisor list for the base % coefficients. These are polynomial variants of let % rules which Groebner has found for the parameters. % For the time being, this contains the smacros that used to be in % consel, and repeats those in bcoeff. %---------------------------------------------------------------- % For compatibility with REDUCE 3.5: macro procedure dipoly!-compat(u); if version!*="REDUCE 3.5" then <<fluid '(!*arbvars !*varopt); {'progn,'(setq !*arbvars t), '(setq !*varopt t), if 'psl memq lispsystem!* then << newtok'((!# !=) eq); put('iequal,'quotenewnam,'eq); nil >> else << newtok'((!# !=) eqn); put('eqn,'op,get('eq,'op)); put('eqn,'infix,get('eq,'infix)); put('iequal,'quotenewnam,'eqn); nil >> } >>; dipoly!-compat(); %---------------------------------------------------------------- %/*Constructors and selectors for a distributed polynomial form*/ %/*A distributive polynomial has the following informal syntax: % % <dipoly> ::= dipzero % | <exponent vector> . <base coefficient> . <dipoly>*/ % Vdp2dip modules included. They could be in a separate package. create!-package('(dipoly a2dip bcoeff dip2a dipoly1 dipvars expvec torder vdp2dip vdpcom condense), '(contrib dipoly)); put('dipoly,'version,3.7); %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 a2dip; %/*Convert an algebraic (prefix) form to distributive polynomial*/ %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ % Modified by: H. Melenk. fluid '(dipvars!* dipzero !*vdpinteger); symbolic 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") % Handling expt separately because the exponents should % not be simplified as domain elements. else if car u = 'expt then dipfnpow(a2dip cadr u,caddr u) 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); symbolic procedure dipfnpow(v,n); % V is a dip. Result is the distributive poly v**n. (if not fixp n or n<0 then typerr(n,"distributive polynomial exponent") else if n=0 then if dipzero!? v then rerror(dipoly,1,"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 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); symbolic 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 or (!*vdpinteger and not bcone!? diplbc 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 bcoeff; % Computation of base coefficients. % Definitions of base coefficient operations for distributive % polynomial package. Fields and rings are supported as coefficient % domains. Side relations (computing modulo an ideal) are supported % if the list bczerodivl is non-zero. % % In this module, a standard quotient coefficient is assumed, unless % !*grmod!* is true, in which case it is a small modular number. %/*Authors: R. Gebauer, A. C. Hearn, H. Kredel*/ % H. Melenk: added routines for faster computation with % quotients representing integers. % June 94: Added side relations (zero divisor list) bczerodivl!*; % Nov 94: Using quotfx and quotientx for divison with known zero % remainder. fluid '(dmode!* bczerodivl!* % List of expressions (sf) equivalent to zero. !*bcsubs2 % check bc for substitutions. !*grmod!* % indicating modular coefficients !*vdpinteger % if t: arithmetic restricted to ring mode. ); symbolic procedure bcint2op(a1,a2,op); null dmode!* and 1=denr a1 and numberp (a1:=numr a1) and 1=denr a2 and numberp (a2:=numr a2) and % (a1 := idapply(op,list(a1,a2))) and (a1 := apply2(op,a1,a2)) and ((if a1=0 then nil else a1) ./ 1); fluid '(!*nat); % symbolic 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); % The following two could be smacros. However, they would then need to % be included in dipoly, thus destroying the modularity of the base % coefficient code. symbolic procedure bcminus!? u; % /* Boolean function. Returns true if u is a negative base coeff*/ null !*grmod!* and minusf numr u; symbolic procedure bczero!? u; % /* Returns a boolean expression, true if the base coefficient u is % zero*/ if !*grmod!* then eqn(u,0) else null numr u; %symbolic procedure bcxzero!? u; % extended zero test: check for zero divisors. % begin scalar l,d; % l:=bczerodivl!*; % u := numr u; % while not d and l do <<d:=quotf(u,car l);l:=cdr l>>; % return d; % end; % symbolic 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); symbolic procedure bcfd a; % /* Base coefficient from domain. a is a domain element. bcfd(a) % returns the base coefficient a. */ if null !*grmod!* then mkbc(a,1) else if fixp a then bcfi a else if not(car a eq '!:mod!:) then rederr list("Invalid modular coefficient",a) else bcfi cdr a; symbolic procedure bcfi a; % /* Base coefficient from integer. a is an integer. bcfi(a) % returns the base coefficient a. */ (if u<0 then current!-modulus-u else u) where u=remainder(a,current!-modulus); symbolic procedure bcdomain!? u; % True if base coefficient u is a domain element. !*grmod!* or (denr u =1 and domainp numr u); symbolic 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(quotfx(denr u,gcdf(denr u,denr v)),denr v)); symbolic 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(quotfx(denr u,denr v),numr v),1); % symbolic 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); symbolic 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. */ if !*grmod!* then eqn(u,1) else denr u = 1 and numr u = 1; symbolic procedure bcinv u; % /* Base coefficient inverse. u is a base coefficient. % bcinv(u) calculates 1/u and returns a base coefficient. */ if !*grmod!* then modular!-reciprocal u else invsq u; symbolic procedure bcneg u; % /* Base coefficient negative. u is a base coefficient. % bcneg(u) returns the negative of the base coefficient % u, a base coefficient. */ if !*grmod!* then if eqn(u,0) then u else current!-modulus #- u else negsq u; symbolic procedure bcprod (u,v); % /* Base coefficient product. u and v are base coefficients. % bcprod(u,v) calculates u*v and returns a base coefficient. if !*grmod!* then remainder(u*v,current!-modulus) else bcint2op(u,v,'times2) or bccheckz multsq(u,v); symbolic procedure mkbc(u,v); % /* Convert u and v into u/v in lowest terms*/ if !*grmod!* then remainder(u*general!-modular!-reciprocal v,current!-modulus) else if v = 1 then u ./ v else if minusf v then mkbc(negf u,negf v) else quotfx(u,m) ./ quotfx(v,m) where m = gcdf(u,v); if null getd 'quotientx then copyd('quotientx,'quotient); symbolic procedure bcquot(u,v); % /* Base coefficient quotient. u and v are base coefficients. % bcquot(u,v) calculates u/v and returns a base coefficient. */ if !*grmod!* then remainder(u*general!-modular!-reciprocal v,current!-modulus) else if !*vdpinteger then (bcint2op(u,v,'quotientx) or !*f2q quotfx(numr u,numr v)) else quotsq(u,v); symbolic procedure bcsum(u,v); % /* Base coefficient sum. u and v are base coefficients. % bcsum(u,v) calculates u+v and returns a base coefficient. */ if !*grmod!* then ((if w#<current!-modulus then w else w #- current!-modulus) where w=u #+ v) else bcint2op(u,v,'plus2) or bccheckz addsq(u,v); symbolic procedure bccheckz u; % Reduce a sum/difference result by members of bczerodivl!*. if null numr u then u else if !*bcsubs2 then subs2 u else <<while l and n do <<n:=cdr qremf(n,car l);l := cdr l>>; n./d>> where l=bczerodivl!*,n=numr u,d=denr u; symbolic procedure bcdif(u,v); % /* Base coefficient difference. u and v are base coefficients. % bcdif(u,v) calculates u-v and returns a base coefficient. */ if !*grmod!* then if u#>v then u #- v else u #- v #+ current!-modulus else bcint2op(u,v,'difference) or bcsum(u,bcneg v); symbolic procedure bcpow(u,n); % /*Returns the base coefficient u raised to the nth power, where % n is an integer*/ if !*grmod!* then modular!-expt(u,n) else exptsq(u,n); symbolic procedure a2bc u; % /*Converts the algebraic (kernel) u into a base coefficient. if !*grmod!* then if not domainp u then rederr list("Invalid coefficient",u) else bcfd u else simp!* u; symbolic procedure bc2a u; % /* Returns the prefix equivalent of the base coefficient u*/ if !*grmod!* then u else prepsq u; fluid '(!*groebigpos !*groebigneg !*groescale); !*groescale := 20; !*groebigpos:= 10** !*groescale; !*groebigneg := - 10** !*groescale; symbolic procedure bcprin u; % /* Prints a base coefficient in infix form*/ if !*grmod!* then prin2 u else begin scalar nat; nat := !*nat; !*nat := nil; if cdr u = 1 and numberp car u and (car u>!*groebigpos or car u<!*groebigneg) then bcprin2big car u else if cdr u neq 1 or not numberp car u then <<prin2!* "[";sqprint u; prin2!* "]">> else sqprint u; !*nat := nat end; symbolic procedure bcprin2big u; << if u<0 then<< prin2 "-"; u:= -u>>; bcprin2big1(u,0)>>; symbolic procedure bcprin2big1 (u,n); if u>!*groebigpos then bcprin2big1 (u/!*groebigpos,n#+!*groescale) else <<prin2 u; prin2 "e"; prin2 n>>; 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 dipoly1; % Distributive polynomial algorithms. % Authors: R. Gebauer, A. C. Hearn, H. Kredel. % Modification for REDUCE > 3.3: H. Melenk. fluid '(dipvars!* dipzero); symbolic procedure dipconst!? p; not dipzero!? p and dipzero!? dipmred p and evzero!? dipevlmon p; symbolic procedure terprit n; for i:=1:n do terpri(); symbolic procedure dfcprint pl; % h polynomial factor list of distributive polynomials print. for each p in pl do dfcprintin p; symbolic 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; symbolic 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 >>; symbolic procedure diplcm p; % Distributive polynomial least common multiple of denominators. % 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); symbolic 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); symbolic procedure diprectoint1(p,u); % Distributive polynomial conversion rational to integral internal 1. % diprectoint1 is used in diprectoint. if dipzero!? p then dipzero else dipmoncomp(bclcmdprod(u,diplbc p),dipevlmon p, diprectoint1(dipmred p,u)); symbolic 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; symbolic 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); symbolic 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)); symbolic 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; symbolic 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; symbolic 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); symbolic 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; symbolic 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); symbolic 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); symbolic 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); symbolic 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 ); symbolic 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; symbolic 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); symbolic 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); symbolic 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; symbolic 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); symbolic 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; %symbolic 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; symbolic 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. % Iterative version, better suited for very long polynomials. % Warning: this routine uses "dipmred" == "cdr cdr" for a % destrcutive concatenation. begin scalar w,rw,sl,ep1,ep2,nt,al,done; while not done do <<if dipzero!? p1 then <<nt:=p2; done:=t>> else if dipzero!? p2 then <<nt:=p1; done:=t>> else <<ep1 := dipevlmon p1; ep2 := dipevlmon p2; sl := evcomp(ep1, ep2); % Compute the next term. if sl #= 1 then <<nt := dipmoncomp(diplbc p1, ep1, nil); p1 := dipmred p1>> else if sl #= -1 then <<nt := dipmoncomp(diplbc p2, ep2, nil); p2 := dipmred p2>> else <<al := bcsum(diplbc p1, diplbc p2); nt := if not bczero!? al then dipmoncomp(al,ep1,nil); p1 := dipmred p1; p2 := dipmred p2>>; >>; % Append the term to the sum polynomial. if nt then if null w then w:=rw:=nt else <<cdr cdr rw := nt; rw := nt>>; >>; return w; end; 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 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*/ % Modification H. Melenk (August 1988) % 1. Most ev-routines handle exponent vectors with variable length: % the convention is, that trailing zeros may be omitted. % 2. evcompless!? is mapped to evcomp such that each term order mode % is supported by exactly one procedure entry. % 3. complete exponent vector compare collected in separate module % TORDER (TORD33) fluid '(dipsortmode!* dipvars!*); symbolic 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); symbolic procedure evcons (e1,e2); % /* Exponent vector construct. e1 and e2 are exponents. evcons(e1,e2) % constructs an exponent vector. */ e1 . e2; symbolic 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 null e1 then 0 else if n = 1 then evfirst e1 else evnth(evred e1, n - 1); symbolic procedure evred e1; % /* Exponent vector reductum. e1 is an exponent vector. evred(e1) % returns the reductum of the exponent vector e1. */ if e1 then cdr e1 else nil; symbolic procedure evfirst e1; % /* Exponent vector first. e1 is an exponent vector. evfirst(e1) % returns the first element of the exponent vector e1, an exponent. */ if e1 then car e1 else 0; symbolic 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)); symbolic 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; symbolic procedure indexcpl(ev,n); % returns a list of indexes 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)) ); symbolic procedure evzer1!? e; % returns a boolean expression. true if e is null else false. null e; symbolic procedure evzero!? e; % /* Returns a boolean expression. True if all exponents are zero*/ null e or car e = 0 and evzero!? cdr e; symbolic 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; symbolic 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; symbolic 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 and e2 do <<x := (if car e1 #> car e2 then car e1 else car e2) . x; e1 := cdr e1; e2 := cdr e2>>; return reversip x end; symbolic 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. */ if e1 and e2 then not(car e1 #< car e2) and evmtest!?(cdr e1,cdr e2) else evzero!? e2 ; symbolic 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 and e2 do <<x := (car e1 #+ car e2) . x; e1 := cdr e1; e2 := cdr e2>>; x := reversip x; return if e1 then nconc(x,e1) else if e2 then nconc(x,e2) else x; end; symbolic 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 e2 do <<if null e1 then e1 := '(0); x := (car e1 #- car e2) . x; e1 := cdr e1; e2 := cdr e2>>; return nconc (reversip x,e1); end; symbolic procedure intevprod(n,e); % /* Multiplies each element of the exponent vector u by the integer n*/ for each x in e collect n #* x; symbolic procedure expvec2a e; % /* Returns list of prefix equivalents of exponent vector e*/ expvec2a1(e,dipvars!*); symbolic 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); symbolic 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); symbolic 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 "*"; if atom car u or null get(caar u,'dipprifn) then dipprin2 car u else apply1(get(caar u,'dipprifn),car u); if car e #> 1 then <<dipprin2 "**"; dipprin2 car e>>; dipevlpri1(cdr e,cdr u,t)>>; endmodule; module torder; % Term order modes for distributive polynomials. % H. Melenk, ZIB Berlin. % The routines of this module should be coded as efficiently as % possible. fluid '(dipsortmode!* dipsortevcomp!* olddipsortmode!* dipvars!*); fluid '(vdpsortmode!* vdpsortextension!* vdpmatrix!* global!-dipvars!* compiled!-orders!*); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % switching between term order modes: TORDER statement. % global!-dipvars!*:='(list); symbolic procedure torder u; begin scalar oldmode,oldex,oldvars,w; oldmode := vdpsortmode!*; oldex := vdpsortextension!*; oldvars := global!-dipvars!*; global!-dipvars!* := '(list); a: w:=reval car u; u:=cdr u; if eqcar(w,'list) and null u then<<u:=cdr w; goto a>>; if eqcar(w,'list) then <<global!-dipvars!*:=w; w:=reval car u; u:=cdr u>>; vdpsortmode!* := w; vdpsortextension!* := for each x in u join (if eqcar(x:=reval x,'list) then cdr x else {x}); if flagp(vdpsortmode!*,'dipsortextension) and null vdpsortextension!* then rederr "term order needs additional parameter(s)"; return 'list . oldvars . oldmode . oldex; end; remprop('torder,'number!-of!-args); put('torder,'psopfn,'torder); symbolic procedure dipsortingmode u; % /* Sets the exponent vector sorting mode. Returns the previous mode*/ begin scalar x,z; if not idp u or not flagp(u,'dipsortmode) then return typerr(u,"term ordering mode"); x := dipsortmode!*; dipsortmode!* := u; % saves thousands of calls to GET; dipsortevcomp!* := get(dipsortmode!*,'evcomp); if not getd dipsortevcomp!* then rerror(dipoly,2, "No compare routine for term order mode found"); if (z:=get(dipsortmode!*,'evcompinit)) then apply(z,nil); if (z:=get(dipsortmode!*,'evlength)) and z neq length dipvars!* then rederr "wrong variable number for fixed length term order"; return x end; flag('(lex gradlex revgradlex),'dipsortmode); put('lex,'evcomp,'evlexcomp); put('gradlex,'evcomp,'evgradlexcomp); put('revgradlex,'evcomp,'evrevgradlexcomp); symbolic 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. % Mapped to evcomp 1 = evcomp(e2,e1); symbolic 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. % IDapply would be better here, but is not within standard LISP! apply(dipsortevcomp!*,list(e1,e2)); symbolic procedure evlexcomp (e1,e2); % /* Exponent vector lexicographical compare. The % exponent vectors e1 and e2 are in lexicographical % ordering. evLexComp(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 null e2 then evlexcomp(e1,'(0)) else if car e1 #= car e2 then evlexcomp(cdr e1,cdr e2) else if car e1 #> car e2 then 1 else -1; symbolic procedure evinvlexcomp (e1,e2); % Exponent vector inverse lexicographical compare. if null e1 then if null e2 then 0 else evinvlexcomp('(0),e2) else if null e2 then evlexcomp(e1,'(0)) else if car e1 #= car e2 then evinvlexcomp(cdr e1,cdr e2) else (if not(n#=0) then n else if car e2 eq car e1 then 0 else if car e2 #> car e1 then 1 else -1) where n = evinvlexcomp(cdr e1,cdr e2); symbolic procedure evgradlexcomp (e1,e2); % /* Exponent vector graduated lex compare. % The exponent vectors e1 and e2 are in graduated lex % ordering. evGradLexComp(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 null e2 then evgradlexcomp(e1,'(0)) else if car e1 #= car e2 then evgradlexcomp(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; symbolic procedure evrevgradlexcomp (e1,e2); % /* Exponent vector reverse graduated lex compare. % The exponent vectors e1 and e2 are in reverse graduated lex % ordering. evRevGradLexcomp(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 null e2 then evrevgradlexcomp(e1,'(0)) else if car e1 #= car e2 then evrevgradlexcomp(cdr e1, cdr e2) else (if te1 #= te2 then evinvlexcomp(e1,e2) else if te1 #> te2 then 1 else -1) where te1 = evtdeg e1, te2 = evtdeg e2; symbolic 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; % The following secion contains additional term order modes. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % gradlexgradlex % % this order can have several steps % torder gradlexgradlex,3,2,4; % flag ('(gradlexgradlex),'dipsortmode); flag ('(gradlexgradlex),'dipsortextension); put('gradlexgradlex,'evcomp,'evgradgradcomp); symbolic procedure evgradgradcomp (e1,e2); evgradgradcomp1 (e1,e2,car vdpsortextension!*, cdr vdpsortextension!*); symbolic procedure evgradgradcomp1 (e1,e2,n,nl); if null e1 then 0 else if null e2 then evgradgradcomp1(e1,'(0),n,nl) else if n#=0 then if null nl then evgradlexcomp(e1,e2) else evgradgradcomp1 (e1,e2,car nl,cdr nl) else if car e1 #= car e2 then evgradgradcomp1(cdr e1,cdr e2,n#-1,nl) else (if te1 #= te2 then if car e1 #> car e2 then 1 else -1 else if te1 #> te2 then 1 else -1) where te1 = evpartdeg(e1,n), te2 = evpartdeg(e2,n); symbolic procedure evpartdeg(e1,n); evpartdeg1(e1,n,0); symbolic procedure evpartdeg1(e1,n,sum); if n #= 0 or null e1 then sum else evpartdeg1(cdr e1,n #-1, car e1 #+ sum); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % gradlexrevgradlex % % flag ('(gradlexrevgradlex),'dipsortmode); flag ('(gradlexrevgradlex),'dipsortextension); put('gradlexrevgradlex,'evcomp,'evgradrevgradcomp); symbolic procedure evgradrevgradcomp (e1,e2); evgradrevgradcomp1 (e1,e2,car vdpsortextension!*); symbolic procedure evgradrevgradcomp1 (e1,e2,n); if null e1 then 0 else if null e2 then evgradrevgradcomp1(e1,'(0),n) else if n#=0 then evrevgradlexcomp(e1,e2) else if car e1 #= car e2 then evgradrevgradcomp1(cdr e1,cdr e2,n#-1) else (if te1 #= te2 then if car e1 #< car e2 then 1 else -1 else if te1 #> te2 then 1 else -1) where te1 = evpartdeg(e1,n), te2 = evpartdeg(e2,n); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % LEXGRADLEX % % flag ('(lexgradlex),'dipsortmode); flag ('(lexgradlex),'dipsortextension); put('lexgradlex,'evcomp,'evlexgradlexcomp); symbolic procedure evlexgradlexcomp (e1,e2); evlexgradlexcomp1 (e1,e2,car vdpsortextension!*); symbolic procedure evlexgradlexcomp1 (e1,e2,n); if null e1 then (if evzero!? e2 then 0 else -1) else if null e2 then evlexgradlexcomp1(e1,'(0),n) else if n#=0 then evgradlexcomp(e1,e2) else if car e1 #= car e2 then evlexgradlexcomp1(cdr e1,cdr e2,n#-1) else if car e1 #> car e2 then 1 else -1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % LEXREVGRADLEX % % flag ('(lexrevgradlex),'dipsortmode); flag ('(lexrevgradlex),'dipsortextension); put('lexrevgradlex,'evcomp,'evlexrevgradlexcomp); symbolic procedure evlexrevgradlexcomp (e1,e2); evlexrevgradlexcomp1 (e1,e2,car vdpsortextension!*); symbolic procedure evlexrevgradlexcomp1 (e1,e2,n); if null e1 then (if evzero!? e2 then 0 else -1) else if null e2 then evlexrevgradlexcomp1(e1,'(0),n) else if n#=0 then evrevgradlexcomp(e1,e2) else if car e1 #= car e2 then evlexrevgradlexcomp1(cdr e1,cdr e2,n#-1) else if car e1 #> car e2 then 1 else -1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % WEIGHTED % % flag ('(weighted),'dipsortmode); flag ('(weighted),'dipsortextension); put('weighted,'evcomp,'evweightedcomp); symbolic procedure evweightedcomp (e1,e2); (if dg1 #= dg2 then evlexcomp(e1,e2) else if dg1 #> dg2 then 1 else -1 ) where dg1=evweightedcomp1(e1,vdpsortextension!*), dg2=evweightedcomp1(e2,vdpsortextension!*); symbolic procedure evweightedcomp1 (e,w); % scalar product of exponent and weight vector if null e then 0 else if null w then evweightedcomp1 (e,'(1 1 1 1 1)) else car e #* car w #+ evweightedcomp1(cdr e,cdr w); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % GRADED term order % cascading a graded sorting with another term order. % % The grade of a term is defined as a scalar product of the exponent % vector and a grade vector which contains non-negative integers. % In contrast to a weight vector the grade vector may contain also % zeros. A vector of ones is used if no vector is given explicitly. % fluid '(gradedrec!*); flag ('(graded),'dipsortmode); flag ('(graded),'dipsortextension); put('graded,'evcomp,'evgradedcomp); put('graded,'evcompinit,'evgradedinit); symbolic procedure evgradedinit(); begin scalar w,gvect,vse; vse:=vdpsortextension!*; while pairp vdpsortextension!* and numberp car vdpsortextension!* do <<gvect:=car vdpsortextension!* . gvect; vdpsortextension!* := cdr vdpsortextension!*>>; if vdpsortextension!* then <<w:=car vdpsortextension!*; vdpsortextension!* := cdr vdpsortextension!*>> else w:='lex; dipsortingmode w; gradedrec!*:={reversip gvect,dipsortevcomp!*,vdpsortextension!*}; dipsortevcomp!* := 'evgradedcomp; dipsortmode!* := 'graded; vdpsortextension!* := vse; end; symbolic procedure evgradedcomp (e1,e2); (if dg1 #= dg2 then apply(cadr gradedrec!*,{e1,e2}) where vdpsortextension!*=caddr gradedrec!* else if dg1 #> dg2 then 1 else -1 ) where dg1=ev!-gamma e1, dg2=ev!-gamma e2; symbolic procedure ev!-gamma(ev); % compute the grade of an exponent vector; evweightedcomp1 (ev, if dipsortmode!* = 'graded then car gradedrec!* else if dipsortmode!* = 'weighted then vdpsortextension!* else nil); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % MATRIX % % % In the following routines I assume that 99 percent of the matrix % entries will be 0 or 1 such that the special branches for these % numbers makes sense. We save lots of memory read and % multiplication is needed only entries other than 0 and 1. % % I could do the same optimization step for -1, but I don't % expect that many people will use term orders with negative % numbers. % % This package includes a compilation mode for matrix term orders % for fixed length variable lists. Compilation is done implicilty % when *comp is on, or explicitly by callint torder_compile. flag ('(matrix),'dipsortmode); flag ('(matrix),'dipsortextension); put('matrix,'evcomp,'evmatrixcomp); put('matrix,'evcompinit,'evmatrixinit); symbolic procedure evmatrixcomp(e1,e2); evmatrixcomp1(e1,e2,vdpmatrix!*); symbolic procedure evmatrixcomp1(e1,e2,m); if null m then 0 else (if w1 #= w2 then evmatrixcomp1(e1,e2,cdr m) else % #= if w1 #> w2 then 1 else -1) where w1= evmatrixcomp2 (e1,car m,0), w2= evmatrixcomp2 (e2,car m,0); symbolic procedure evmatrixcomp2(e,l,w); if null e or null l then w else (if l1 #= 0 then evmatrixcomp2(cdr e,cdr l,w) else if l1 #= 1 then evmatrixcomp2(cdr e,cdr l,w #+ car e) else evmatrixcomp3(e,l1,l,w)) where l1=car l; symbolic procedure evmatrixcomp3(e,l1,l,w); evmatrixcomp2(cdr e,cdr l,w #+ car e #* l1); symbolic procedure evmatrixinit1(w,mode); begin scalar m,mm; if not eqcar(w,'mat) or mode and length cadr w neq length dipvars!* then typerr(w,"term order matrix for". dipvars!*); for each row in cdr w do <<row:=for each c in row collect ieval c; mm:=row . mm; row:=reverse row; while eqcar(row,0) do row := cdr row; m:=reversip row . m>>; m:=reversip m; mm:=reversip mm; if m neq vdpmatrix!* then <<if length cadr w > length cdr w then lprim "Warning: non-square matrix used in torder" else if 0=reval{'det,w} then typerr(w,"term order (singular matrix)"); if not evmatrixcheck mm then typerr(w,"term order (non admissible)") >>; return m end; symbolic procedure evmatrixinit(); begin scalar c,m,w; w:=reval car vdpsortextension!*; m:=evmatrixinit1(w,t); if (c:=assoc(m,compiled!-orders!*)) then dipsortevcomp!* := cdr c else if !*comp then dipsortevcomp!* := evmatrixcompile m; vdpmatrix!*:=m; end; symbolic procedure evmatrixcheck m; % Check the usability of the term order matrix: the % top elements of each column must be positive. This % approach goes back to a recommendation of J. Apel. begin scalar bad,c,w; integer i,j,r; r:=length m; for i:=1:length car m do <<c:=nil; for j:=1:r do if (w:=nth(nth(m,j),i)) neq 0 and null c then <<c:=t; bad:=w < 0>> >>; return not bad; end; symbolic procedure evmatrixcompile m; begin scalar w; w:= evmatrixcompile1 m; putd(car w,'expr,caddr w); compiled!-orders!* := (m.car w).compiled!-orders!*; return car w; end; symbolic procedure evmatrixcompile1 m; begin scalar c,n,x,w,lvars,code; integer ld,p,k; for each row in m do k:=max(k,length row); lvars := for i:=1:k collect gensym(); code := {{'setq,car lvars, '(idifference (car e1) (car e2))}}; ld:=1; for each row in m do <<p:=0; w:=nil; while row do <<c:=car row; row := cdr row; p:=p+1; if c neq 0 then << % load the differences up to the current point. for i:=ld+1:p do << code:= append(code,'((setq e1(cdr e1))(setq e2(cdr e2)))); code := append(code, {{'setq,nth(lvars,i), '(idifference (car e1) (car e2))}}); ld:=i; >>; % collect the terms of the row sum x:=nth(lvars,p); if c = -1 then x := {'iminus,x} else if c neq 1 then x:={'itimes2,c,x}; w:=if w then {'iplus2,x,w} else x; >>; >>; if not atom w then <<code:=append(code,{{'setq,'w,w}});w:='w>>; code:=append(code, {{'cond,{{'iequal,w,0},t}, {{'igreaterp,w,0},'(return 1)}, '(t (return -1))}}); >>; % common trailor code:=append(code,'((return 0))); n:=gensym(); return {n,k,evform {'lambda,'(e1 e2), 'prog.('w.lvars). code}}; end; symbolic procedure evform(u); % Let form play on the generated code; form1(u,nil,'symbolic); symbolic procedure torder_compile_form(w,c,m); begin scalar n; if length w < 3 then rederr "illegal arguments"; m:=evmatrixinit1(eval form caddr w,nil); c:=evmatrixcompile1 m; n:=eval form cadr w; return {'progn, {'putd,mkquote n,mkquote 'expr,mkquote caddr c}, {'setq,'compiled!-orders!*, {'cons,{'cons,mkquote m,mkquote n}, 'compiled!-orders!*}}, {'put,mkquote n,''evcomp,mkquote n}, {'put,mkquote n,''evlength,cadr c}, {'flag,mkquote{n},''dipsortmode}, mkquote n}; end; put('torder_compile,'formfn,'torder_compile_form); endmodule; module vdp2dip; imports dipoly; % create!-package('(vdp2dip vdpcom vdp2dip1),'(contrib dipoly)); % load!-package 'dipoly; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % interface for Virtual Distributive Polynomials (VDP) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % "Distributive representation" with respect to a given set of % variables ("vdpvars") means for a polynomial, that the polynomial % is regarded as a sequence of monomials, each of which is a % product of a "coefficient" and of some powers of the variables. % This internal representation is very closely connected to the % standard external (printed) representation of a polynomial in % REDUCE if nothing is factored out. The monomials are logically % ordered by a term order mode based on the ordering which is % given bye the sequence "vdpvars"; with respect to this ordering % the representation of a polynomial is unique. The "highest" term % is the car one. Monomials are represented by their coefficient % ("vbc") and by a vector of the exponents("vev") (in the order % corresponding to the vector vars). The distributive representation % is good for those algorithms, which base their decisions on the % complete ledading monomial: this representation guarantees a % fast and uniform access to the car monomial and to the reductum % (the cdr of the polynomial beginning with the cadr monomial). % The algorithms of the Groebner package are of this type. The % interface defines the distributive polynomials as abstract data % objects via their acess functions. These functions map the % distributive operations to an arbitrary real data structure % ("virtual"). The mapping of the access functions to an actual % data structure is cdrricted only by the demand, that the typical % "distributive operations" be efficient. Additionally to the % algebraic value a VDP object has a property list. So the algorithms % using the VDP interface can assign name-value-pairs to individual % polynomials. The interface is defined by a set of routines which % create and handle the distributive polynomials. In general the % first letters of the routine name classifies the data its works on: % vdp... complete virt. polynomial objects % vbc... virt. base coefficients % vev... virt. exponent vectors % 0. general control % % vdpinit(dv) initialises the vdp package for the variables % given in the list "dv". vdpinit modifies the % torder and returns the prvevious torder as its % result. vdpinit sets the global variable % vdpvars!*; % 1. conversion % % a2vdp algebraic (prefix) to vdp % f2vdp standard form to vdp % a2vbc algebraic (prefix) to vbc % vdp2a vdp to algebraic (prefix) % vdp2f vdp to standard form % vbc2a vbc to algebraic (prefix) % 2. composing/decomposing % % vdpfmon make a vdp from a vbc and an vev % vdpMonComp add a monomial (vbc and vev) to the front of a vdp % vdpAppendMon add a monomial (vbc and vev) to the bottom of a vdp % vdpMonAdd add a monomial (vbc and vev) to a vdp, not yet % knowing the place of the insertion % vdpAppendVdp concat two vdps % % vdpLbc extract leading vbc % vdpevlmon extract leading vev % vdpred reductum of vdp % vdplastmon last monomial of polynomial % vevnth nth element from exponent vector % 3. testing % % vdpZero? test vdp = 0 % vdpredZero!? test rductum of vdp = 0 % vdpOne? test vdp = 1 % vevZero? test vev = (0 0 ... 0) % vbczero? test vbc = 0 % vbcminus? test vbc <= 0 (not decidable for algebraic vbcs) % vbcplus? test vbc >= 0 (not decidable for algebraic vbcs) % vbcone!? test vbc = 1 % vbcnumberp test vbc is a numeric value % vevdivides? test if vev1 < vev2 elementwise % vevlcompless? test ordering vev1 < vev2 % vdpvevlcomp calculate ordering vev1 / vev1: -1, 0 or +1 % vdpEqual test vdp1 = vdp2 % vdpMember member based on "vdpEqual" % vevequal test vev1 = vev2 % 4. arithmetic % % 4.1 vdp arithmetic % % vdpsum vdp + vdp % special routines for monomials: see above (2.) % vdpdif vdp - vdp % vdpprod vdp * vdp % vdpvbcprod vbc * vdp % vdpDivMon vdp / (vbc,vev) divisability presumed % vdpCancelvev substitute all multiples of monomial (1,vev) in vdp by 0 % vdpLcomb1 vdp1*(vbc1,vev1) + vdp2*(vbc2,vev2) % vdpContent calculate gcd over all vbcs % 4.2 vbc arithmetic % % vbcsum vbc1 + vbc2 % vbcdif vbc1 - vbc2 % vbcneg - vbc % vbcprod vbc1 * vbc2 % vbcquot vbc1 / vbc2 divisability assumed if domain = ring % vbcinv 1 / vbc only usable in field % vbcgcd gcd(vbc1,vbc2) only usable in Euclidean field % 4.2 vev arithmetic % % vevsum vev1 + vev2 elementwise % vevdif vev1 - vev2 elementwise % vevtdeg sum over all exponents % vevzero generate a zero vev % 5. auxiliary % % vdpPutProp assign indicator-value-pair to vdp % the property "number" is used for printing. % vdpGetProp read value of indicator from vdp % vdplSort sort list of polynomials with respect to ordering % vdplSortIn sort a vdp into a sorted list of vdps % vdpprint print a vdp together with its number % vdpprin2t print a vdp "naked" % vdpprin3t print a vdp with closing ";" % vdpcondense replace exponent vectors by equal objects from % global list dipevlist!* in order to save memory %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % RECORD STRUCTURE % % a virtual polynomial here is a record (list) with the entries % ('vdp <vdpevlmon> <vdplbc> <form> <plist>) % % 'vdp a type tag % <vdpevlmon> the exponents of the variables in the leading % leading monomial; the positions correspond to % the positions in vdpvars!*. Trailing zeroes % can be omitted. % % <lcoeff> the "coefficient" of the leading monomial, which % in general is a standard form. % % <form> the complete polynomial, e.g. as REDUCE standard form. % % <plist> an asso list for the properties of the polynomial % % The components should not be manipulated only via the interface % functions and macros, so that application programs remain % independent from the internal representation. % The only general assumption made on <form> is, that the zero % polynomial is represented as NIL. That is the case e.g. for both, % REDUCE standard forms and DIPOLYs. % Conventions for the usage: % ------------------------- % % vdpint has to be called prveviously to all vdp calls. The list of % vdp paraemters is passed to vdpinit. The value of vdpvars!* % and the current torder must remain unmodfied afterwards. % usual are simple id's, e.g. % % % Modifications to vdpvars!* during calculations % ---------------------------------------------- % % This mapping of vdp operations to standard forms offers the % ability to enlarge vdpvars during the calculation in order % to add new (intermediate) variables. Basis is the convention, % that exponent vectors logically have an arbitrary number % of trailing zeros. All routines processing exponent vectors % are able to handle varying length of exponent vectors. % A new call to vdpinit is necessary. % % During calculation vdpvars may be enlarged (new variables % suffixed) without needs to modify existing polynomials; only % korder has to be set to the new variable sequence. % modifications to the sequence in vdpvars requires a % new call to vdpinit and a reordering of exisiting % polynomials, e.g. by % vdpint newvdpvars; % f2vdp vdp2f p1; f2vdp vdp2f p2; ..... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DECLARATION SECTION % % this module must be present during code generation for modules % using the vdp - sf interface fluid '(vdpvars!* intvdpvars!* secondvalue!* vdpsortmode!* !*groebrm !*vdpinteger !*trgroeb !*trgroebs !*groebdivide pcount!* !*gsugar dipevlist!*); global '(vdpprintmax groebmonfac); flag('(vdpprintmax),'share); % basic internal constructor of vdp-record smacro procedure makevdp (vbc,vev,form); list('vdp,vev,vbc,form,nil); % basic selectors (conversions) smacro procedure vdppoly u; cadr cddr u; smacro procedure vdplbc u; caddr u; smacro procedure vdpevlmon u; cadr u; % basic tests smacro procedure vdpzero!? u; null u or null vdppoly u; smacro procedure vevzero!? u; null u or (car u = 0 and vevzero!?1 cdr u); smacro procedure vdpone!? p; not vdpzero!? p and vevzero!? vdpevlmon p; % base coefficients % manipulating of exponent vectors smacro procedure vevdivides!? (vev1,vev2); vevmtest!?(vev2,vev1); smacro procedure vevzero(); vevmaptozero1(vdpvars!*,nil); smacro procedure vdpnumber f; vdpgetprop(f,'number) ; % the code for checkpointing is factored out % This version: NO CHECKPOINT FACILITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % interface for DIPOLY polynomials as records (objects). % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % fluid '(intvdpvars!* vdpvars!* secondvalue!* vdpsfsortmode!* !*groebrm !*vdpinteger !*trgroeb !*trgroebs !*groebdivide pcount!* !*groebsubs); fluid '(vdpsortmode!*); global '(vdpprintmax groebmonfac); flag('(vdpprintmax),'share); fluid '(dipvars!* !*vdpinteger); symbolic procedure dip2vdp u; % Is used when u can be empty. (if dipzero!? uu then makevdp(a2bc 0,nil,nil) else makevdp(diplbc uu,dipevlmon uu,uu)) where uu = if !*groebsubs then dipsubs2 u else u; % some simple mappings smacro procedure makedipzero(); nil; symbolic procedure vdpredzero!? u; dipzero!? dipmred vdppoly u; symbolic procedure vdplastmon u; % Return bc . ev of last monomial of u. begin u:=vdppoly u; if dipzero!? u then return nil; while not dipzero!? u and not dipzero!? dipmred u do u:=dipmred u; return diplbc u . dipevlmon u; end; symbolic procedure vbczero!? u; bczero!? u; symbolic procedure vbcnumber u; if pairp u and numberp car u and 1=cdr u then cdr u else nil; symbolic procedure vbcfi u; bcfd u; symbolic procedure a2vbc u; a2bc u; symbolic procedure vbcquot(u,v); bcquot(u,v); symbolic procedure vbcneg u; bcneg u; symbolic procedure vbcabs u; if vbcminus!? u then bcneg u else u; symbolic procedure vbcone!? u; bcone!? u; symbolic procedure vbcprod (u,v); bcprod(u,v); % initializing vdp-dip polynomial package symbolic procedure vdpinit2(vars); begin scalar oldorder; vdpcleanup(); oldorder := kord!*; if null vars then rerror(dipoly,8,"Vdpinit: vdpvars not set"); vdpvars!* := dipvars!* := vars; torder2 vdpsortmode!*; return oldorder; end; symbolic procedure vdpcleanup(); dipevlist!*:={nil}; symbolic procedure vdpred u; begin scalar r,s; r:=dipmred vdppoly u; if dipzero!? r then return makevdp(nil ./ nil,nil,makedipzero()); r:=makevdp(diplbc r,dipevlmon r,r); if !*gsugar and (s:=vdpgetprop(u,'sugar)) then gsetsugar(r,s); return r; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % coefficient handling; here we assume that coefficients are % standard quotients; % symbolic procedure vbcgcd (u,v); if denr u = 1 and denr v = 1 then if fixp numr u and fixp numr v then gcdn(numr u,numr v) ./ 1 else gcdf!*(numr u,numr v) ./ 1 else 1 ./ 1; % Cofactors: compute (q,v) such that q*a=v*b. symbolic procedure vbc!-cofac(bc1,bc2); % Compute base coefficient cofactors. <<if vbcminus!? bc1 and vbcminus!? bc2 then gcd:=vbcneg gcd; vbcquot(bc2,gcd).vbcquot(bc1,gcd) >> where gcd=vbcgcd(bc1,bc2); symbolic procedure vev!-cofac(ev1,ev2); % Compute exponent vector cofactors. (vevdif(lcm,ev1) . vevdif(lcm,ev2)) where lcm =vevlcm(ev1,ev2); % the following functions must be redefinable symbolic procedure vbcplus!? u; (numberp v and v>0) where v = numr u; symbolic procedure bcplus!? u; (numberp v and v>0) where v = numr u; symbolic procedure vbcminus!? u; (numberp v and v<0) where v = numr u; symbolic procedure vbcinv u; bcinv u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % conversion between forms, vdps and prefix expressions % % prefix to vdp symbolic procedure a2vdp u; if u=0 or null u then makevdp(nil ./ nil,nil,makedipzero()) else (makevdp(diplbc r,dipevlmon r,r) where r = a2dip u); % vdp to prefix symbolic procedure vdp2a u; dip2a vdppoly u; symbolic procedure vbc2a u; bc2a u; % form to vdp symbolic procedure f2vdp(u); if u=0 or null u then makevdp(nil ./ nil,nil,makedipzero()) else (makevdp(diplbc r,dipevlmon r,r) where r = f2dip u); % vdp to form symbolic procedure vdp2f u; dip2f vdppoly u; % vdp from monomial symbolic procedure vdpfmon (coef,vev); begin scalar r; r:= makevdp(coef,vev,dipfmon(coef,vev)); if !*gsugar then gsetsugar(r,vevtdeg vev); return r; end; % add a monomial to a vdp in front (new vev and coeff) symbolic procedure vdpmoncomp(coef,vev,vdp); if vdpzero!? vdp then vdpfmon(coef,vev) else if vbczero!? coef then vdp else makevdp(coef,vev,dipmoncomp(coef,vev,vdppoly vdp)); %add a monomial to the end of a vdp (vev remains unchanged) symbolic procedure vdpappendmon(vdp,coef,vev); if vdpzero!? vdp then vdpfmon(coef,vev) else if vbczero!? coef then vdp else makevdp(vdplbc vdp,vdpevlmon vdp, dipsum(vdppoly vdp,dipfmon(coef,vev))); % add monomial to vdp, place of new monomial still unknown symbolic procedure vdpmonadd(coef,vev,vdp); if vdpzero!? vdp then vdpfmon(coef,vev) else (if c = 1 then vdpmoncomp(coef,vev,vdp) else if c = -1 then makevdp (vdplbc vdp,vdpevlmon vdp, dipsum(vdppoly vdp,dipfmon(coef,vev))) else vdpsum(vdp,vdpfmon(coef,vev)) ) where c = vevcomp(vev,vdpevlmon vdp); symbolic procedure vdpzero(); a2vdp 0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % comparing of exponent vectors % % symbolic procedure vdpvevlcomp (p1,p2); dipevlcomp (vdppoly p1,vdppoly p2); symbolic procedure vevilcompless!?(e1,e2); 1 = evilcomp(e2,e1); symbolic procedure vevilcomp (e1,e2); evilcomp (e1,e2); symbolic procedure vevcompless!?(e1,e2); 1 = evcomp(e2,e1); symbolic procedure vevcomp (e1,e2); evcomp (e1,e2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % routines traversing the "coefficients" % % CONTENT of a vdp % The content is the gcd of all coefficients. symbolic procedure vdpcontent d; if vdpzero!? d then a2bc 0 else <<d := vdppoly d; dipnumcontent(dipmred d,diplbc d)>>; symbolic procedure vdpcontent1(d,c); dipnumcontent(vdppoly d,c); symbolic procedure dipnumcontent(d,c); if bcone!? c or dipzero!? d then c else dipnumcontent(dipmred d,vbcgcd(c,diplbc d)); symbolic procedure dipcontenti p; % the content is a pair of the lcm of the coefficients and the % exponent list of the common monomial factor. if dipzero!? p then 1 else (if dipzero!? rp then diplbc p . (if !*groebrm then dipevlmon p else nil) else dipcontenti1(diplbc p, if !*groebrm then dipevlmon p else nil,rp) ) where rp=dipmred p; symbolic procedure dipcontenti1 (n,ev,p1); if dipzero!? p1 then n . ev else begin scalar nn; nn := vbcgcd (n,diplbc p1); if ev then ev := dipcontevmin(dipevlmon p1,ev); if bcone!? nn and null ev then return nn . nil else return dipcontenti1 (nn,ev,dipmred p1) end; % CONTENT and MONFAC (if groebrm on) symbolic procedure vdpcontenti d; vdpcontent d . if !*groebrm then vdpmonfac d else nil; symbolic procedure vdpmonfac d; dipmonfac vdppoly d; symbolic procedure dipmonfac p; % exponent list of the common monomial factor. if dipzero!? p or not !*groebrm then evzero() else (if dipzero!? rp then dipevlmon p else dipmonfac1(dipevlmon p,rp) ) where rp=dipmred p; symbolic procedure dipmonfac1(ev,p1); if dipzero!? p1 or evzero!? ev then ev else dipmonfac1(dipcontevmin(ev,dipevlmon p1),dipmred p1); % vdpCoeffcientsFromDomain!? symbolic procedure vdpcoeffcientsfromdomain!? w; dipcoeffcientsfromdomain!? vdppoly w; symbolic procedure dipcoeffcientsfromdomain!? w; if dipzero!? w then t else (if bcdomain!? v then dipcoeffcientsfromdomain!? dipmred w else nil) where v =diplbc w; symbolic procedure vdplength f; diplength vdppoly f; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % polynomial operations: % coefficient normalization and reduction of monomial % factors % symbolic procedure vdpequal(p1,p2); p1 eq p2 or (n1 and n1 = n2 % number comparison is faster most times or dipequal(vdppoly p1,vdppoly p2) where n1 = vdpgetprop(p1,'number), n2 = vdpgetprop(p2,'number)); symbolic procedure dipequal(p1,p2); if dipzero!? p1 then dipzero!? p2 else if dipzero!? p2 then nil else diplbc p1 = diplbc p2 and evequal(dipevlmon p1,dipevlmon p2) and dipequal(dipmred p1,dipmred p2); symbolic procedure evequal(e1,e2); % test equality with variable length exponent vectors if null e1 and null e2 then t else if null e1 then evequal('(0),e2) else if null e2 then evequal(e1,'(0)) else 0=(car e1 #- car e2) and evequal(cdr e1,cdr e2); symbolic procedure vdplcm p; diplcm vdppoly p; symbolic procedure vdprectoint(p,q); dip2vdp diprectoint(vdppoly p,q); symbolic procedure vdpsimpcont(p); begin scalar r,q; q := vdppoly p; if dipzero!? q then return p; r := dipsimpcont q; p := dip2vdp cdr r; % the polynomial r := car r; % the monomial factor if any if not evzero!? r and (dipmred q or evtdeg r>1) then vdpputprop(p,'monfac,r); return p; end; symbolic procedure dipsimpcont (p); if !*vdpinteger or not !*groebdivide then dipsimpconti p else dipsimpcontr p; % routines for integer coefficient case: % calculation of contents and dividing all coefficients by it symbolic procedure dipsimpconti (p); % calculate the contents of p and divide all coefficients by it begin scalar co,lco,res,num; if dipzero!? p then return nil . p; co := bcfd 1; co := if !*groebdivide then dipcontenti p else if !*groebrm then co . dipmonfac p else co . nil; num := car co; if not bcplus!? num then num := bcneg num; if not bcplus!? diplbc p then num := bcneg num; if bcone!? num and cdr co = nil then return nil . p; lco := cdr co; if groebmonfac neq 0 then lco := dipcontlowerev cdr co; res := p; if not(bcone!? num and lco = nil) then res := dipreduceconti (p,num,lco); if null cdr co then return nil . res; lco := evdif(cdr co,lco); return(if lco and not evzero!? evdif(dipevlmon res,lco) then lco else nil).res; end; symbolic procedure vdpreduceconti (p,co,vev); % divide polynomial p by monomial from co and vev vdpdivmon(p,co,vev); % divide all coefficients of p by cont symbolic procedure dipreduceconti (p,co,ev); if dipzero!? p then makedipzero() else dipmoncomp ( bcquot (diplbc p,co), if ev then evdif(dipevlmon p,ev) else dipevlmon p, dipreduceconti (dipmred p,co,ev)); % routines for rational coefficient case: % calculation of contents and dividing all coefficients by it symbolic procedure dipsimpcontr (p); % calculate the contents of p and divide all coefficients by it begin scalar co,lco,res; if dipzero!? p then return nil . p; co := dipcontentr p; if bcone!? diplbc p and co = nil then return nil . p; lco := dipcontlowerev co; res := p; if not(bcone!? diplbc p and lco = nil) then res := dipreducecontr (p,bcinv diplbc p,lco); return (if co then evdif(co,lco) else nil) . res; end; symbolic procedure dipcontentr p; % the content is the exponent list of the common monomial factor. (if dipzero!? rp then (if !*groebrm then dipevlmon p else nil) else dipcontentr1(if !*groebrm then dipevlmon p else nil,rp) ) where rp=dipmred p; symbolic procedure dipcontentr1 (ev,p1); if dipzero!? p1 then ev else begin if ev then ev := dipcontevmin(dipevlmon p1,ev); if null ev then return nil else return dipcontentr1 (ev,dipmred p1) end; % divide all coefficients of p by cont symbolic procedure dipreducecontr (p,co,ev); if dipzero!? p then makedipzero() else dipmoncomp ( bcprod (diplbc p,co), if ev then evdif(dipevlmon p,ev) else dipevlmon p, dipreducecontr (dipmred p,co,ev)); symbolic procedure dipcontevmin (e1,e2); % calculates the minimum of two exponents; if one is shorter, trailing % zeroes are assumed. % e1 is an exponent vector. e2 is a list of exponents begin scalar res; while e1 and e2 do <<res := (if ilessp(car e1,car e2) then car e1 else car e2) . res; e1 := cdr e1; e2 := cdr e2>>; while res and 0=car res do res := cdr res; return reversip res; end; symbolic procedure dipcontlowerev (e1); % subtract a 1 from those elements of an exponent vector which % are greater than 1. % e1 is a list of exponents, the result is an exponent vector. begin scalar res; while e1 do <<res := (if igreaterp(car e1,0) then car e1 - 1 else 0) . res; e1 := cdr e1>>; while res and 0 = car res do res := cdr res; if res and !*trgroebs then <<prin2 "***** exponent reduction:"; prin2t reverse res>>; return reversip res; end; symbolic procedure dipappendmon(dip,bc,ev); append(dip,dipfmon(bc,ev)); smacro procedure dipnconcmon(dip,bc,ev); nconc(dip,dipfmon(bc,ev)); smacro procedure dipappenddip(dip1,dip2); append(dip1,dip2); smacro procedure dipnconcdip(dip1,dip2); nconc(dip1,dip2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % basic polynomial arithmetic: % symbolic procedure vdpsum(d1,d2); begin scalar r; r:=dip2vdp dipsum(vdppoly d1,vdppoly d2); if !*gsugar then gsetsugar(r,max(gsugar d1,gsugar d2)); return r; end; symbolic procedure vdpdif(d1,d2); begin scalar r; r:= dip2vdp dipdif(vdppoly d1,vdppoly d2); if !*gsugar then gsetsugar(r,max(gsugar d1,gsugar d2)); return r; end; symbolic procedure vdpprod(d1,d2); begin scalar r; r:= dip2vdp dipprod(vdppoly d1,vdppoly d2); if !*gsugar then gsetsugar(r,gsugar d1 + gsugar d2); return r; end; % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % linear combination: the Buchberger Workhorse % % LCOMB1: calculate mon1 * vdp1 + mon2 * vdp2 symbolic procedure vdpilcomb1(d1,vbc1,vev1,d2,vbc2,vev2); begin scalar r; r:= dip2vdp dipilcomb1 (vdppoly d1,vbc1,vev1,vdppoly d2,vbc2,vev2); if !*gsugar then gsetsugar(r,max(gsugar d1+vevtdeg vev1, gsugar d2+vevtdeg vev2)); return r; end; symbolic procedure dipilcomb1 (p1,bc1,ev1,p2,bc2,ev2); % same asl dipILcomb, exponent vectors multiplied in already begin scalar ep1,ep2,sl,res,sum,z1,z2,p1new,p2new,lptr,bptr,c; z1 := not evzero!? ev1; z2 := not evzero!? ev2; p1new := p2new := t; lptr := bptr := res := makedipzero(); loop: if p1new then << if dipzero!? p1 then return if dipzero!? p2 then res else dipnconcdip(res, dipprod(p2,dipfmon(bc2,ev2))); ep1 := dipevlmon p1; if z1 then ep1 := evsum(ep1,ev1); p1new := nil;>>; if p2new then << if dipzero!? p2 then return dipnconcdip(res, dipprod(p1,dipfmon(bc1,ev1))); ep2 := dipevlmon p2; if z2 then ep2 := evsum(ep2,ev2); p2new := nil; >>; sl := evcomp(ep1, ep2); if sl = 1 then << c := bcprod(diplbc p1,bc1); if not bczero!? c then <<lptr := dipnconcmon (bptr,c,ep1); bptr := dipmred lptr>>; p1 := dipmred p1; p1new := t; >> else if sl = -1 then << c := bcprod(diplbc p2,bc2); if not bczero!? c then <<lptr := dipnconcmon (bptr,c,ep2); bptr := dipmred lptr>>; p2 := dipmred p2; p2new := t; >> else << sum := bcsum (bcprod(diplbc p1,bc1), bcprod(diplbc p2,bc2)); if not bczero!? sum then << lptr := dipnconcmon(bptr,sum,ep1); bptr := dipmred lptr>>; p1 := dipmred p1; p2 := dipmred p2; p1new := p2new := t; >>; if dipzero!? res then <<res := bptr := lptr>>; % initial goto loop; end; symbolic procedure vdpvbcprod(p,a); (if !*gsugar then gsetsugar(q,gsugar p) else q) where q=dip2vdp dipbcprod(vdppoly p,a); symbolic procedure vdpdivmon(p,c,vev); (if !*gsugar then gsetsugar(q,gsugar p) else q) where q=dip2vdp dipdivmon(vdppoly p,c,vev); symbolic procedure dipdivmon(p,bc,ev); % divides a polynomial by a monomial % we are sure that the monomial ev is a factor of p if dipzero!? p then makedipzero() else dipmoncomp ( bcquot(diplbc p,bc), evdif(dipevlmon p,ev), dipdivmon (dipmred p,bc,ev)); symbolic procedure vdpcancelmvev(p,vev); (if !*gsugar then gsetsugar(q,gsugar p) else q) where q=dip2vdp dipcancelmev(vdppoly p,vev); symbolic procedure dipcancelmev(f,ev); % cancels all monomials in f which are multiples of ev dipcancelmev1(f,ev,makedipzero()); symbolic procedure dipcancelmev1(f,ev,res); if dipzero!? f then res else if evmtest!?(dipevlmon f,ev) then dipcancelmev1(dipmred f,ev,res) else dipcancelmev1(dipmred f,ev, % dipAppendMon(res,diplbc f,dipevlmon f)); dipnconcmon(res,diplbc f,dipevlmon f)); % some prehistoric routines needed in resultant operation symbolic procedure vevsum0(n,p); % exponent vector sum version 0. n is the length of vdpvars!*. % p is a distributive polynomial. if vdpzero!? p then vevzero1 n else vevsum(vdpevlmon p, vevsum0(n,vdpred p)); symbolic procedure vevzero1 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; symbolic procedure vdpresimp u; % fi domain changes, the coefficients have to be resimped dip2vdp dipresimp vdppoly u; symbolic procedure dipresimp u; if null u then nil else (for each x in u collect <<toggle := not toggle; if toggle then simp prepsq x else x>> ) where toggle = t; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % printing of polynomials % symbolic procedure vdpprin2t u; << vdpprint1(u,nil,9999); terpri()>>; symbolic procedure vdpprin3t u; << vdpprint1(u,nil,9999); prin2t ";">>; symbolic procedure vdpprint u; <<vdpprin2 u; terpri()>>; symbolic procedure vdpprin2 u; <<(if x then <<prin2 "P("; prin2 x; if s then<<prin2 "/",prin2 s>>; prin2 "): ">>) where x=vdpgetprop(u,'number), s=vdpgetprop(u,'sugar); vdpprint1(u,nil,vdpprintmax)>>; symbolic procedure vdpprint1(u,v,max); vdpprint1x(vdppoly u,v,max); symbolic procedure vdpprint1x(u,v,max); % /* 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 % max limits the number of terms to be printed if dipzero!? u then if null v then dipprin2 0 else nil else if max = 0 then % maximum of terms reached << terpri(); prin2 " ### etc ("; prin2 diplength u; prin2 " terms) ###"; terpri();>> 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; vdpprint1x(dipmred u,t, max - 1) end; symbolic procedure dipprin2 u; <<if posn()>69 then terprit 2 ; prin2 u>>; symbolic procedure vdpsave u; u; % switching between term order modes symbolic procedure torder2 u; dipsortingmode u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % additional conversion utilities % conversion dip to standard form / standard quotient symbolic procedure dip2f u; (if denr v neq 1 then <<print u; rerror(dipoly,9, "Distrib. poly. with rat coeff cannot be converted")>> else numr v) where v = dip2sq u; symbolic procedure dip2sq u; % convert a dip into a standard quotient. if dipzero!? u then nil ./ 1 else addsq(diplmon2sq(diplbc u,dipevlmon u),dip2sq dipmred u); symbolic procedure diplmon2sq(bc,ev); %convert a monomial into a standard quotient. multsq(bc,dipev2f(ev,dipvars!*) ./ 1); symbolic procedure dipev2f(ev,vars); if null ev then 1 else if car ev = 0 then dipev2f(cdr ev,cdr vars) else multf(car vars .** car ev .* 1 .+ nil, dipev2f(cdr ev,cdr vars)); % evaluate SUBS2 for the coefficients of a dip symbolic procedure dipsubs2 u; begin scalar v,secondvalue!*; secondvalue!* := 1 ./ 1; v := dipsubs21 u; return diprectoint(v,secondvalue!*); end; symbolic procedure dipsubs21 u; begin scalar c; if dipzero!? u then return u; c := groebsubs2 diplbc u; if null numr c then return dipsubs21 dipmred u; if not(denr c = 1) then secondvalue!* := bclcmd(c,secondvalue!*); return dipmoncomp(c,dipevlmon u,dipsubs21 dipmred u); end; % conversion standard form to dip symbolic procedure f2dip u; f2dip1(u,evzero(),bcfd 1); symbolic procedure f2dip1 (u,ev,bc); % f to dip conversion: scan the standard form. ev % and bc are the exponent and coefficient parts collected % so far from higher parts. if null u then nil else if domainp u then dipfmon(bcprod(bc,bcfd u),ev) else dipsum(f2dip2(mvar u,ldeg u,lc u,ev,bc), f2dip1(red u,ev,bc)); symbolic procedure f2dip2(var,dg,c,ev,bc); % f to dip conversion: % multiply leading power either into exponent vector % or into the base coefficient. <<if ev1 then ev := ev1 else bc := multsq(bc,var.**dg.*1 .+nil./1); f2dip1(c,ev,bc)>> where ev1=if memq(var,dipvars!*) then evinsert(ev,var,dg,dipvars!*) else nil; symbolic procedure evinsert(ev,v,dg,vars); % f to dip conversion: % Insert the "dg" into the ev in the place of variable v. if null ev or null vars then nil else if car vars eq v then dg . cdr ev else car ev . evinsert(cdr ev,v,dg,cdr vars); symbolic procedure vdpcondense f; dipcondense car cdddr f; endmodule; module vdpcom; % common routines to all vdp mappings fluid '(intvdpvars!* vdpvars!* secondvalue!* vdpsortmode!* !*groebrm !*vdpinteger !*trgroeb !*groebdivide pcount!* !*gsugar vdpsortextension!* ); global '(vdpprintmax); flag('(vdpprintmax),'share); vdpprintmax := 5; % Repeat of smacros defined in vdp2dip. smacro procedure vdppoly u; cadr cddr u; smacro procedure vdpzero!? u; null u or null vdppoly u; smacro procedure vdpevlmon u; cadr u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % manipulating of exponent vectors % symbolic procedure vevnth (a,n); % extract nth element from a if null a then 0 else if n=1 then car a else vevnth(cdr a,n #- 1); % unrolled code for zero test (very often called) smacro procedure vevzero!? u; null u or (car u = 0 and vevzero!?1 cdr u); symbolic procedure vevzero!?1 u; null u or (car u = 0 and vevzero!? cdr u); symbolic procedure veveq(vev1,vev2); if null vev1 then vevzero!? vev2 else if null vev2 then vevzero!? vev1 else (car vev1 = car vev2 and vevequal(cdr vev1,vev2)); symbolic procedure vevmaptozero e; % generate an exponent vector with same length as e and zeros only vevmaptozero1(e,nil); symbolic procedure vevmaptozero1(e,vev); if null e then vev else vevmaptozero1(cdr e, 0 . vev); symbolic procedure vevmtest!? (e1,e2); % /* Exponent vector multiple test. e1 and e2 are compatible exponent % vectors. vevmtest?(e1,e2) returns a boolean expression. % True if exponent vector e1 is a multiple of exponent % vector e2, else false. */ if null e2 then t else if null e1 then if vevzero!? e2 then t else nil else not(car e1 #<car e2)and vevmtest!?(cdr e1,cdr e2); symbolic procedure vevlcm (e1,e2); % /* Exponent vector least common multiple. e1 and e2 are % exponent vectors. vevlcm(e1,e2) computes the least common % multiple of the exponent vectors e1 and e2, and returns % an exponent vector. */ begin scalar x; while e1 and e2 do <<x := (if car e1 #> car e2 then car e1 else car e2) . x; e1 := cdr e1; e2 := cdr e2>>; x := reversip x; if e1 then x := nconc(x,e1) else if e2 then x := nconc(x,e2); return x; end; symbolic procedure vevmin (e1,e2); % Exponent vector minima begin scalar x; while e1 and e2 do <<x := (if car e1 #< car e2 then car e1 else car e2) . x; e1 := cdr e1; e2 := cdr e2>>; while x and 0=car x do x := cdr x; % cut trailing zeros return reversip x; end; symbolic procedure vevsum (e1,e2); % /* Exponent vector sum. e1 and e2 are exponent vectors. % vevsum(e1,e2) calculates the sum of the exponent vectors. % e1 and e2 componentwise and returns an exponent vector. */ begin scalar x; while e1 and e2 do <<x := (car e1 #+ car e2) . x;e1 := cdr e1; e2 := cdr e2>>; x := reversip x; if e1 then x := nconc(x,e1) else if e2 then x := nconc(x,e2); return x; end; symbolic procedure vevtdeg u; % calculate the total degree of u if null u then 0 else car u #+ vevtdeg cdr u; symbolic procedure vdptdeg u; if vdpzero!? u then 0 else max(vevtdeg vdpevlmon u,vdptdeg vdpred u); symbolic procedure vevdif (ee1,ee2); % Exponent vector difference. e1 and e2 are exponent % vectors. vevdif(e1,e2) calculates the difference of the % exponent vectors e1 and e2 componentwise and returns an % exponent vector. begin scalar x,y,break,e1,e2; e1 := ee1; e2 := ee2; while e1 and e2 and not break do <<y := (car e1 #- car e2); x := y . x; break := y #< 0; e1 := cdr e1; e2 := cdr e2>>; if break or (e2 and not vevzero!? e2) then <<print ee1; print ee2; if getd 'backtrace then backtrace(); return rerror(dipoly,5,"Vevdif, difference would be < 0")>>; return nconc(reversip x,e1); end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % numbering of polynomials % symbolic procedure vdpenumerate f; % f is a temporary result. Prepare it for medium range storage % and ssign a number if vdpzero!? f then f else << f := vdpsave f; if not vdpgetprop(f,'number) then f := vdpputprop(f,'number,(pcount!* := pcount!* #+ 1)); f>>; %smacro procedure vdpNumber f; % vdpGetProp(f,'NUMBER) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % SUGAR of polynomials % symbolic procedure gsugar p; if !*gsugar then ((s or << print list("*** missing sugar",p); backtrace(); s:=vdptdeg p; gsetsugar(p,s); s>> ) where s= if vdpzero!? p then 0 else vdpgetprop(p,'sugar)); symbolic procedure gsetsugar(p,s); !*gsugar and vdpputprop(p,'sugar,s or vdptdeg p) or p; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % operations on sets of polynomials % symbolic procedure vdpmember(p1,l); % test membership of a polynomial in a list of polys if null l then nil else if vdpequal(p1,car l) then l else vdpmember(p1,cdr l); symbolic procedure vdpunion (s1,s2); % s1 and s2 are two sets of polynomials. % union of the sets using vdpMember as crit if null s1 then s2 else if vdpmember(car s1,s2) then vdpunion(cdr s1,s2) else car s1 . vdpunion(cdr s1,s2); symbolic procedure vdpintersection (s1,s2); % s1 and s2 are two sets of polynomials. % intersection of the sets using vdpMember as crit if null s1 then nil else if vdpmember(car s1,s2) then car s1 . vdpunion(cdr s1,s2) else vdpunion(cdr s1,s2); symbolic procedure vdpsetequal!?(s1,s2); % tests if s1 and s2 have the same polynomials as members if not (length s1 = length s2) then nil else vdpsetequal!?1(s1,append(s2,nil)); symbolic procedure vdpsetequal!?1(s1,s2); % destroys its second parameter (is therefor copied when called) if null s1 and null s2 then t else if null s1 or null s2 then nil else (if hugo then vdpsetequal!?1(cdr s1,groedeletip(car hugo,s2)) else nil) where hugo = vdpmember(car s1,s2); symbolic procedure vdpsortedsetequal!?(s1,s2); % tests if s1 and s2 have the same polynomials as members % here assuming, that both sets are sorted by the same % principles if null s1 and null s2 then t else if null s1 or null s2 then nil else if vdpequal(car s1,car s2) then vdpsortedsetequal!?(cdr s1,cdr s2) else nil; symbolic procedure vdpdisjoint!? (s1,s2); % s1 and s2 are two sets of polynomials. % test that there are no common members if null s1 then t else if vdpmember(car s1,s2) then nil else vdpdisjoint!?(cdr s1,s2); symbolic procedure vdpsubset!? (s1,s2); not(length s1 > length s2) and vdpsubset!?1(s1,s2); symbolic procedure vdpsubset!?1 (s1,s2); % s1 and s2 are two sets of polynomials. % test if s1 is subset of s2 if null s1 then t else if vdpmember(car s1,s2) then vdpsubset!?1 (cdr s1,s2) else nil; symbolic procedure vdpdeletemember(p,l); % delete polynomial p from list l if null l then nil else if vdpequal(p,car l) then vdpdeletemember(p,cdr l) else car l . vdpdeletemember(p,cdr l); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % sorting of polynomials % symbolic procedure vdplsort pl; % /* Distributive polynomial list sort. pl is a list of % distributive polynomials. vdplsort(pl) returns the % sorted distributive polynomial list of pl. sort(pl, function vdpvevlcomp); symbolic procedure vdplsortin (po,pl); % po is a polynomial, pl is a list of polynomials. % po is inserted into pl at its place determined by vevlcompless?. % the result is the updated pl; if null pl then list po else if vevcompless!? (vdpevlmon po, vdpevlmon car pl) then car pl . vdplsortin (po, cdr pl) else po . pl; symbolic procedure vdplsortinreplacing (po,pl); % po is a polynomial, pl is a list of polynomials. % po is inserted into pl at its place determined by vevlcompless?. % if there is a multiple of the exponent of pl, this is deleted % the result is the updated pl; if null pl then list po else if vevdivides!? (vdpevlmon po, vdpevlmon car pl) then vdplsortinreplacing (po, cdr pl) else if vevcompless!? (vdpevlmon po, vdpevlmon car pl) then car pl . vdplsortinreplacing (po, cdr pl) else po . pl; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % property lists for polynomials % symbolic procedure vdpputprop (poly,prop,val); begin scalar c,p; if not pairp poly or not pairp (c := cdr poly) or not pairp (c := cdr c) or not pairp (c := cdr c) or not pairp (c := cdr c ) then rerror(dipoly,6, list("vdpPutProp given a non-vdp as 1st parameter", poly,prop,val)); p := assoc(prop,car c); if p then rplacd(p,val) else rplaca(c,(prop . val) . car c); return poly; end; symbolic procedure vdpgetprop (poly,prop); if null poly then nil % nil is a legal variant of vdp=0 else if not eqcar(poly,'vdp) then rerror(dipoly,7, list("vdpGetProp given a non-vdp as 1st parameter", poly,prop)) else (if p then cdr p else nil) where p=assoc(prop,cadr cdddr poly); symbolic procedure vdpremallprops u; begin scalar c; if not pairp u or not pairp (c := cdr u) or not pairp (c := cdr c) or not pairp (c := cdr c) or not pairp (c := cdr c) then return u; rplaca(c,nil); return u; end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Groebner interface to power substitution fluid '(!*sub2); symbolic procedure groebsubs2 q; (subs2 q) where !*sub2=t; % and a special print symbolic procedure vdpprintshort u; begin scalar m; m := vdpprintmax; vdpprintmax:= 2; vdpprint u; vdpprintmax:=m; end; endmodule; module condense; % unify exponent vectors for lower memory consumption. % Author: Herbert Melenk fluid '(dipevlist!*); dipevlist!*:={nil}; symbolic procedure dipcondense f; begin scalar dl,ev; dl:=dipevlist!*; while f do <<ev := dipevlmon f; while cdr dl and evcompless!?(dipevlmon f,cadr dl) do dl:=cdr dl; if cdr dl and ev=cadr dl then car f := cadr dl else cdr dl:= ev.cdr dl; f:=dipmred f; >>; end; endmodule; end;