reduce-historical/r37/packages/dipoly/vdpcom.red

module vdpcom;
 
% common routines to all vdp mappings
 

fluid '(intvdpvars!* vdpvars!* secondvalue!* vdpsortmode!* !*groebrm
	!*vdpinteger !*trgroeb !*groebdivide pcount!* !*gsugar
         vdpsortextension!* );


global '(vdpprintmax);
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 e1 and e2 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 ssign 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>>;

%smacro procedure vdpNumber f;
%     vdpGetProp(f,'NUMBER) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SUGAR of polynomials
%

symbolic procedure gsugar p;
   if !*gsugar then
    ((s or 
        << print list("*** 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 therefor 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 (po,pl);
%    po is a polynomial, pl is a list of polynomials.
%    po is inserted into pl at its place determined by vevlcompless?.
%    the result is the updated pl;
   if null pl then list po
   else if vevcompless!? (vdpevlmon po, vdpevlmon car pl)
       then  car pl . vdplsortin (po, cdr pl)
     else po . pl;
 
symbolic procedure  vdplsortinreplacing (po,pl);
%    po is a polynomial, pl is a list of polynomials.
%    po is inserted into pl at its place determined by vevlcompless?.
%    if there is a multiple of the exponent of pl, this is deleted
%    the result is the updated pl;
   if null pl then list po
   else if vevdivides!? (vdpevlmon po, vdpevlmon car pl) 
       then vdplsortinreplacing (po, cdr pl)
   else if vevcompless!? (vdpevlmon po, vdpevlmon car pl)
       then  car pl . vdplsortinreplacing (po, cdr pl)
     else po . pl;





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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,
                  list("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,
                  list("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 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 return u;
      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;