Artifact 3903ce97b9518ae791df2079d21c5a600f04e1b2274db8b2ac28cbe8c23d0776:
- Executable file
r38/packages/dipoly/vdpcom.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: 11056) [annotate] [blame] [check-ins using] [more...]
module vdpcom; % Common routines to all vdp mappings. flag('(vdpprintmax),'share); vdpprintmax:=5; % Repeat of smacros defined in vdp2dip. smacro procedure vdppoly u;cadr cddr u; smacro procedure vdpzero!? u;null u or null vdppoly u; smacro procedure vdpevlmon u;cadr u; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % manipulating of exponent vectors % symbolic procedure vevnth(a,n); % Extract nth element from 'a'. if null a then 0 else if n=1 then car a else vevnth(cdr a,n #- 1); % Unrolled code for zero test(very often called). smacro procedure vevzero!? u; null u or(car u=0 and vevzero!?1 cdr u); symbolic procedure vevzero!?1 u; null u or(car u=0 and vevzero!? cdr u); symbolic procedure veveq(vev1,vev2); if null vev1 then vevzero!? vev2 else if null vev2 then vevzero!? vev1 else(car vev1=car vev2 and vevequal(cdr vev1,vev2)); symbolic procedure vevmaptozero e; % Generate an exponent vector with same length as e and zeros only. vevmaptozero1(e,nil); symbolic procedure vevmaptozero1(e,vev); if null e then vev else vevmaptozero1(cdr e,0 .vev); symbolic procedure vevmtest!?(e1,e2); % Exponent vector multiple test.'e1' and 'e2' are compatible exponent % vectors.vevmtest?(e1,e2) returns a boolean expression. % True if exponent vector 'e1' is a multiple of exponent % vector 'e2', else false. if null e2 then t else if null e1 then if vevzero!? e2 then t else nil else not(car e1 #< car e2)and vevmtest!?(cdr e1,cdr e2); symbolic procedure vevlcm(e1,e2); % Exponent vector least common multiple.'e1' and 'e2' are % exponent vectors.'vevlcm(e1,e2)' computes the least common % multiple of the exponent vectors 'e1' and 'e2', and returns % an exponent vector. begin scalar x; while e1 and e2 do <<x:=(if car e1 #> car e2 then car e1 else car e2).x; e1:=cdr e1;e2:=cdr e2>>; x:=reversip x; if e1 then x:=nconc(x,e1)else if e2 then x:=nconc(x,e2); return x end; symbolic procedure vevmin(e1,e2); % Exponent vector minima. begin scalar x; while e1 and e2 do <<x:=(if car e1 #< car e2 then car e1 else car e2).x; e1:=cdr e1;e2:=cdr e2>>; while x and 0=car x do x:=cdr x;% Cut trailing zeros. return reversip x end; symbolic procedure vevsum(e1,e2); % Exponent vector sum.'e1' and 'e2' are exponent vectors. % 'vevsum(e1,e2)' calculates the sum of the exponent vectors % 'e1' and 'e2' componentwise and returns an exponent vector. begin scalar x; while e1 and e2 do <<x:=(car e1 #+ car e2).x;e1:=cdr e1;e2:=cdr e2>>; x:=reversip x; if e1 then x:=nconc(x,e1)else if e2 then x:=nconc(x,e2); return x end; symbolic procedure vevtdeg u; % Calculate the total degree of u. if null u then 0 else car u #+ vevtdeg cdr u; symbolic procedure vdptdeg u; if vdpzero!? u then 0 else max(vevtdeg vdpevlmon u,vdptdeg vdpred u); symbolic procedure vevdif(ee1,ee2); % Exponent vector difference.'e1' and 'e2' are exponent % vectors.'vevdif(e1,e2)' calculates the difference of the % exponent vectors componentwise and returns an exponent vector. begin scalar x,y,break,e1,e2; e1:=ee1;e2:=ee2; while e1 and e2 and not break do <<y:=(car e1 #- car e2);x:=y.x;break:=y #< 0; e1:=cdr e1;e2:=cdr e2>>; if break or(e2 and not vevzero!? e2)then <<print ee1;print ee2;if getd 'backtrace then backtrace(); return rerror(dipoly,5,"Vevdif, difference would be < 0")>>; return nconc(reversip x,e1)end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Numbering of polynomials. % symbolic procedure vdpenumerate f; % 'f' is a temporary result.Prepare it for medium range storage % and sign a number. if vdpzero!? f then f else <<f:=vdpsave f; if not vdpgetprop( f,'number)then f:=vdpputprop(f,'number,(pcount!*:=pcount!* #+ 1));f>>; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % SUGAR of polynomials. % symbolic procedure gsugar p; if !*gsugar then (( s or <<print{"*** missing sugar",p};backtrace(); s:=vdptdeg p;gsetsugar(p,s);s>> )where s= if vdpzero!? p then 0 else vdpgetprop(p,'sugar)); symbolic procedure gsetsugar(p,s); !*gsugar and vdpputprop(p,'sugar,s or vdptdeg p)or p; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Operations on sets of polynomials. % symbolic procedure vdpmember(p1,l); % Test membership of a polynomial in a list of polys. if null l then nil else if vdpequal(p1,car l)then l else vdpmember(p1,cdr l); symbolic procedure vdpunion(s1,s2); % 's1' and 's2' are two sets of polynomials; % union of the sets using vdpMember as crit. if null s1 then s2 else if vdpmember(car s1,s2)then vdpunion(cdr s1,s2) else car s1.vdpunion(cdr s1,s2); symbolic procedure vdpintersection(s1,s2); % 's1' and 's2' are two sets of polynomials; % intersection of the sets using vdpmember as crit. if null s1 then nil else if vdpmember(car s1,s2)then car s1.vdpunion(cdr s1,s2) else vdpunion(cdr s1,s2); symbolic procedure vdpsetequal!?(s1,s2); % Tests if 's1' and 's2' have the same polynomials as members. if not(length s1=length s2)then nil else vdpsetequal!?1(s1,append(s2,nil)); symbolic procedure vdpsetequal!?1(s1,s2); % Destroys its second parameter(is therefore copied when called). if null s1 and null s2 then t else if null s1 or null s2 then nil else (if hugo then vdpsetequal!?1(cdr s1,groedeletip(car hugo,s2)) else nil)where hugo=vdpmember(car s1,s2); symbolic procedure vdpsortedsetequal!?(s1,s2); % Tests if 's1' and 's2' have the same polynomials as members % here assuming, that both sets are sorted by the same principles. if null s1 and null s2 then t else if null s1 or null s2 then nil else if vdpequal(car s1,car s2)then vdpsortedsetequal!?(cdr s1,cdr s2)else nil; symbolic procedure vdpdisjoint!?(s1,s2); % 's1' and 's2' are two sets of polynomials; % test that there are no common members. if null s1 then t else if vdpmember(car s1,s2)then nil else vdpdisjoint!?(cdr s1,s2); symbolic procedure vdpsubset!?(s1,s2); not(length s1 > length s2)and vdpsubset!?1(s1,s2); symbolic procedure vdpsubset!?1(s1,s2); % 's1' and 's2' are two sets of polynomials. % Test if 's1' is subset of 's2'. if null s1 then t else if vdpmember(car s1,s2)then vdpsubset!?1(cdr s1,s2)else nil; symbolic procedure vdpdeletemember(p,l); % Delete polynomial 'p' from list 'l'. if null l then nil else if vdpequal(p,car l)then vdpdeletemember(p,cdr l) else car l.vdpdeletemember(p,cdr l); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Sorting of polynomials. % symbolic procedure vdplsort pl; % Distributive polynomial list sort.'pl' is a list of % distributive polynomials.'vdplsort(pl)' returns the % sorted distributive polynomial list of pl. sort(pl,function vdpvevlcomp); symbolic procedure vdplsortin(p,pl); % 'p' is a polynomial, 'pl' is a list of polynomials. % 'p' is inserted into 'pl' at its place determined by vevlcompless?. % The result is the updated pl. if null pl then{p}else<<vdplsortin1(p,pl,pl);pl>>; symbolic procedure vdplsortin1(p,pl,oldpl); if null pl then cdr oldpl:= p.nil else if vevcompless!?(vdpevlmon p,vdpevlmon car pl) then vdplsortin1(p,cdr pl,pl) else <<cdr pl:=car pl.cdr pl;car pl:=p>>; symbolic procedure vdplsortinreplacing(po,pl); % 'po' is a polynomial, 'pl' is a linear list of polynomials(sorted). % 'po' is inserted into 'pl' at its place determined by 'vevlcompless?'. % If there is a multiple of the first exponent of a polynomial in 'pl', % this one is deleted from 'pl'.The result is the updated 'pl'. % 'opl' is the initial value of 'pl',(initial multiples of 'po' are % removed);'oopl' a working version of 'opl'. begin scalar oopl,opl;if pl then go to bb; aa : return po.nil; bb : if vevdivides!?(vdpevlmon po,vdpevlmon car pl)then <<pl:=cdr pl;if null pl then go to aa;go to bb>>; opl:=pl; cc : if null pl then <<oopl:=lastpair opl;cdr oopl:=po.nil;return opl>>; if not(pl eq opl)and vevdivides!?(vdpevlmon po,vdpevlmon car pl)then <<if null cdr pl then <<oopl:=lastpair1 opl;cdr oopl:=nil;pl:=nil>> else <<car pl:=cadr pl;cdr pl:=cddr pl>>; go to cc>>; if vevcompless!?(vdpevlmon po,vdpevlmon car pl)then <<pl:=cdr pl;go to cc>>; cdr pl:=car pl.cdr pl;car pl:=po;% Insert 'po'. return opl end; symbolic procedure lastpair1 opl; % Determine the last full pair(the 'cdr' non-nil)of the linear list % 'opl';if the routine is called with 'nil' or cdr='nil', % return 't'. if null opl or null cdr opl then t else <<while cddr opl do opl:=cdr opl;opl>>; symbolic procedure countlastvar(a,m); % Count the monomials with the last variable of a vdp-polynomial % 'a';'m' determines, whether the first('m' is true)non-factor % of the last variable leads to the result '0';if the polynomial has more % than 25 elements, return '0' if 'm' is false or if the polynomial has % more than 25 terms(divisible by the last variable). begin integer n,nn;a:=vdppoly a; aa: if atom a then return n; nn:=nn #+ 1;if nn #> 25 then return 0;n:=n #+ 1; if countlastvar1 dipevlmon a #< 1 then(if m then return 0 else n:=0); a:=dipmred a;go to aa end; symbolic procedure countlastvar1 b; begin scalar n;n:=1; aa : if atom b then return 0 else if n=vdplastvar!* then return car b; b:=cdr b;n:=n #+ 1;go to aa end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Property lists for polynomials. % symbolic procedure vdpputprop(poly,prop,val); begin scalar c,p; if not pairp poly or not pairp(c:=cdr poly)or not pairp(c:=cdr c) or not pairp(c:=cdr c)or not pairp(c:=cdr c) then rerror(dipoly,6, {"vdpputprop given a non-vdp as 1st parameter",poly,prop,val}); p:=assoc(prop,car c); if p then rplacd(p,val)else rplaca(c,(prop.val).car c); return poly end; symbolic procedure vdpgetprop(poly,prop); if null poly then nil % nil is a legal variant of vdp=0 else if not eqcar(poly,'vdp) then rerror( dipoly,7, {"vdpgetprop given a non-vdp as 1st parameter", poly,prop}) else(if p then cdr p else nil) where p=assoc(prop,cadr cdddr poly); symbolic procedure vdpremallprops u; begin scalar c; if not(not pairp u or not pairp(c:=cdr u)or not pairp(c:=cdr c) or not pairp(c:=cdr c)or not pairp(c:=cdr c)) then rplaca(c,nil);return u end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Groebner interface to power substitution. % fluid'(!*sub2); symbolic procedure groebsubs2 q;(subs2 q)where !*sub2=t; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % And a special print. % symbolic procedure vdpprintshort u; begin scalar m; m:=vdpprintmax;vdpprintmax:= 2;vdpprint u;vdpprintmax:=m end; endmodule;;end;