module vdp2dip;
imports 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!* !*gcd);
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);
begin scalar x;
if not vbcsize(u, -100) or not vbcsize(v, -100)
then return '(1 . 1);
x := 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;
return x
end;
symbolic procedure vbcsize(u, n);
if n #> -1 then nil
else if atom u then n
else
begin
n := vbcsize(car u, n #+ 1);
if null n then return nil;
return vbcsize(cdr u, n);
end;
% 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 gcd;
gcd := !*gcd;
return
begin scalar ep1,ep2,sl,res,sum,z1,z2,p1new,p2new,lptr,bptr,c,!*gcd;
!*gcd := if vbcsize(bc1,-100) and vbcsize(bc2,-100)then gcd;
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
<< if !*gcd and not vbcsize(diplbc p1, -100)
then !*gcd := nil;
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
<< if !*gcd and not vbcsize(diplbc p2, -100)
then !*gcd := nil;
c := bcprod(diplbc p2,bc2);
if not bczero!? c then
<<lptr := dipnconcmon (bptr,c,ep2);
bptr := dipmred lptr>>;
p2 := dipmred p2; p2new := t;
>>
else
<< if !*gcd and (not vbcsize(diplbc p1,-100) or
not vbcsize(diplbc p2,-100)) then !*gcd := nil;
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;
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;
end;