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

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;