1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936 |
- 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;
|