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;