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