Artifact 1480712104316344a66030779bf8ec01bbfccbfadd72796f9551f9aec214be3f:
- Executable file
r37/packages/dipoly/vdp2dip.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 34151) [annotate] [blame] [check-ins using] [more...]
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;