reduce-historical/r38/packages/crack/conlaw0.red

%            CONLAW file with subroutines for CONLAW1/2/3/4

%                   by Thomas Wolf, September 1997

%----------------------------------------------------------------------

symbolic fluid '(reducefunctions_ print_)$ 

algebraic procedure print_claw(eqlist,qlist,plist,xlist)$
begin scalar n$
 n:=length eqlist$
 while qlist neq {} do <<
  if length qlist < n then write "+"$
  write"( ",first qlist," ) * ( ",first eqlist," )"$
  qlist:=rest qlist; eqlist:=rest eqlist
 >>$
 write" = "$
 n:=length xlist$
 while plist neq {} do <<
  if length plist < n then write "+"$
  write"df( ",first plist,", ",first xlist," )"$
  plist:=rest plist;
  xlist:=rest xlist
 >>
end$

symbolic operator lhsli$
symbolic procedure lhsli(eqlist)$
% lhslist1 will be a list of all those lhs's which are a derivative or
%          a power of a derivative which is used to fix dependencies
%          of q_i or p_j
% lhslist2 will be a list of all lhs's of all equations in their
%          order with those lhs's set to 0 which can not be used 
%          for substitutions
begin scalar lhslist1,lhslist2,h1,flg1,flg2$
  for each h1 in cdr eqlist do <<
    flg1:=nil$    % no assignment to lhslist1 done yet
    if (pairp h1) and (car h1 = 'EQUAL) then <<
      h1:=reval cadr h1;
      if (pairp h1) and 
         (car h1='EXPT) and 
         (numberp caddr h1) then <<flg2:=nil;h1:=cadr h1>>
                            else   flg2:=t;
      if (not numberp h1) and
         ((atom h1) or ((car h1='DF) and (atom cadr h1) )) then 
      <<lhslist1:=cons(h1,lhslist1)$
        if flg2 then <<lhslist2:=cons(h1,lhslist2)$
                       flg1:=t>>
      >>
    >>;
    if null flg1 then lhslist2:=cons(0,lhslist2);
  >>$
  return list('LIST,cons('LIST,lhslist1),cons('LIST,lhslist2))
end$

symbolic operator chksub$
symbolic procedure chksub(eqlist,ulist)$
% eqlist is a list of equations   df(f,x,2,y) = ...
% this procedure tests whether 
% - for any equation a derivative on the rhs is equal or a derivative of
%   the lhs?
% - any lhs is equal or the derivative of any other lhs
begin scalar h1,h2,derili,complaint$
 derili:=for each e1 in cdr eqlist collect 
 cons( all_deriv_search(cadr  e1,cdr ulist),      % lhs
       all_deriv_search(caddr e1,cdr ulist) );    % rhs
 %--- Is for any equation a derivative on the rhs equal or a derivative of
 %--- the lhs?
 for each e1 in derili do 
 if car e1 then <<
  h1:=caaar e1;                  % e.g. h1 = (f x 2 y)
  for each h2 in cdr e1 do
  if (car h1 = caar h2) and null which_deriv(cdar h2,cdr h1) then <<
   complaint:=t;
   write "The left hand side ",
         if length h1 = 1 then car h1
                          else cons('DF,h1)$terpri()$
   write " is not a leading derivative in its equation!"$ terpri()
  >>
 >>$
 %--- Is any lhs equal or the derivative of any other lhs? 
 if derili then
 while cdr derili do <<
  if caar derili then <<
   h1:=caaaar derili$
   for each h2 in cdr derili do
   if (car h2)           and
      (car h1=caaaar h2) and
      ((null which_deriv(cdr h1,cdaaar h2)) or
       (null which_deriv(cdaaar h2,cdr h1))    ) then <<
    complaint:=t;
    write"--> One left hand side (lhs) contains a derivative which"$
    terpri()$
    write"is equal or a derivative of a derivative on another lhs!"$
    terpri()$
   >>$
  >>$
  derili:=cdr derili
 >>;
 if complaint then terpri()$
end$

%==== Procedures as in LIEPDE:

symbolic procedure comparedif1(u1l,u2l)$
% checks whether u2l has more or at least equally many 1's, 2's, ...
% contained as u1l.
% returns a list of 1's, 2's, ... which are in excess in u2l
% compared with u1l. The returned value is 0 if both are identical
begin
 scalar ul;
 if u2l=nil then if u1l neq nil then return nil 
				else return 0
	    else if u1l=nil then return u2l
			    else % both are non-nil
 if car u1l < car u2l then return nil else 
 if car u1l = car u2l then return comparedif1(cdr u1l,cdr u2l) else <<
  ul:=comparedif1(u1l,cdr u2l);
  return if not   ul then nil          else 
	 if zerop ul then list car u2l else
			  cons(car u2l,ul)
 >>
end$ % of comparedif1

%-------------

symbolic procedure comparedif2(u1,u1list,du2)$
% checks whether du2 is a derivative of u1 differentiated
% wrt. u1list
begin
 scalar u2l;
 u2l:=combidif(du2)$ % u2l=(u2, 1, 1, ..)
 if car u2l neq u1 then return nil else
 return comparedif1(u1list, cdr u2l)
end$ % of comparedif2

%-------------

symbolic procedure listdifdif1(du1,deplist)$
% lists all elements of deplist which are *not* derivatives
% of du1
begin
 scalar u1,u1list,res,h$
 h:=combidif(du1);
 u1:=car h;
 u1list:=cdr h;
 for each h in deplist do
 if not comparedif2(u1,u1list,h) then res:=cons(h,res);
 return res
end$ % of listdifdif1

%-------------

symbolic operator listdifdif2$
symbolic procedure listdifdif2(lhslist,deplist)$
% lists all elements of deplist which are *not* derivatives
% of any element of lhslist
begin
  scalar h;
  deplist:=cdr reval deplist;
  lhslist:=cdr reval lhslist;
  for each h in lhslist do
  deplist:=listdifdif1(h,deplist);
  return cons('LIST,deplist)
end$ % of listdifdif2

%-------------

symbolic operator totdeg$
symbolic procedure totdeg(p,f)$
%   Ordnung (total) der hoechsten Ableitung von f im Ausdruck p
eval(cons('PLUS,ldiff1(car ldifftot(reval p,reval f),fctargs reval f)))$

%-------------

% symbolic operator totordpot$
% symbolic procedure totordpot(p,f)$
% %   Ordnung (total) der hoechsten Ableitung von f im Ausdruck p
% %   und hoechste Potenz der hoechsten Ableitung
% %   currently not used
% begin scalar a;
%  a:=ldifftot(reval p,reval f);
%  return
%  cons(eval(cons('PLUS,ldiff1(car a,fctargs reval f))), cdr a)
% end$

%-------------

symbolic procedure diffdeg(p,v)$
%   liefert Ordnung der Ableitung von p nach v$
%   p Liste Varible+Ordnung der Ableitung, v Variable (Atom)
if null p then 0                        %  alle Variable bearbeitet ?
else if car p=v then                    %  v naechste Variable ?
     if cdr p then
	if numberp(cadr p) then cadr p  %  folgt eine Zahl ?
				else 1
	else 1
     else diffdeg(cdr p,v)$             %  weitersuchen

%-------------

symbolic procedure ldiff1(l,v)$
%  liefert Liste der Ordnungen der Ableitungen nach den Variablen aus v
%  l Liste (Variable + Ordnung)$ v Liste der Variablen
if null v then nil                      %  alle Variable abgearbeitet ?
else cons(diffdeg(l,car v),ldiff1(l,cdr v))$
					%  Ordnung der Ableitung nach
					%  erster Variable anhaengen

%-------------

symbolic procedure ldifftot(p,f)$
%  leading derivative total degree ordering
%  liefert Liste der Variablen + Ordnungen mit Potenz
%  p Ausdruck in LISP - Notation, f Funktion
ldifftot1(p,f,fctargs f)$

%-------------

symbolic procedure ldifftot1(p,f,vl)$
%  liefert Liste der Variablen + Ordnungen mit Potenz
%  p Ausdruck in LISP - Notation, f Funktion, lv Variablenliste
begin scalar a$
  a:=cons(nil,0)$
  if not atom p then
%    if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,
%	                 'QUOTIENT,'DF,'EQUAL)) then
    if member(car p,REDUCEFUNCTIONS_) then
					%  erlaubte Funktionen
    <<if (car p='PLUS) or (car p='TIMES) or 
         (car p='QUOTIENT) or (car p='EQUAL) then
      <<p:=cdr p$
        while p do 
        <<a:=diffreltot(ldifftot1(car p,f,vl),a,vl)$
          p:=cdr p
        >> 
      >>                      else 
      if car p='MINUS then a:=ldifftot1(cadr p,f,vl) else 
      if car p='EXPT then               %  Exponent
%      if numberp caddr p then
%      <<a:=ldifftot1(cadr p,f,vl)$
%        a:=cons(car a,times(caddr p,cdr a))
%      >>                 else a:=cons(nil,0)
      <<a:=ldifftot1(cadr p,f,vl)$
        if (numberp caddr p) and
           (numberp cdr a) then a:=cons(car a,times(caddr p,cdr a))
                           else a:=cons(car a,10000)
      >>
                                        %  Potenz aus Basis wird mit
					%  Potenz multipliziert
                     else 
      if car p='DF then                 %  Ableitung
      if cadr p=f then a:=cons(cddr p,1)
					%  f wird differenziert?
	          else a:=cons(nil,0)
                   else                 %  any other non-linear function
      <<p:=cdr p$
        while p do 
        <<a:=diffreltot(ldifftot1(car p,f,vl),a,vl)$
          p:=cdr p
        >>;
        a:=cons(car a,10000)
      >>  
    >> else                             %  sonst Konstante bzgl. f
               
    if p=f then a:=cons(nil,1)	        %  Funktion selbst
           else a:=cons(nil,0)          %  alle uebrigen Fkt. werden
  else if p=f then a:=cons(nil,1)$        %  wie Konstante behandelt
  return a
end$

%-------------

symbolic procedure diffreltot(p,q,v)$
%   liefert komplizierteren Differentialausdruck$
if diffreltotp(p,q,v) then q
		   else p$

%-------------

symbolic procedure diffreltotp(p,q,v)$
%   liefert t, falls p einfacherer Differentialausdruck, sonst nil
%   p, q Paare (liste.power), v Liste der Variablen
%   liste Liste aus Var. und Ordn. der Ableit. in Diff.ausdr.,
%   power Potenz des Differentialausdrucks
begin scalar n,m$
m:=eval(cons('PLUS,ldiff1(car p,v)))$
n:=eval(cons('PLUS,ldiff1(car q,v)))$
return
 if m<n then t
 else if n<m then nil
      else diffrelp(p,q,v)$
end$

%-------------

algebraic procedure subdif1(xlist,ylist,ordr)$
% A list of lists of derivatives of one order for all functions
begin
 scalar allsub,revx,i,el,oldsub,newsub;
 revx:=reverse xlist;
 allsub:={};
 oldsub:= for each y in ylist collect y=y;
 for i:=1:ordr do      %  i is the order of next substitutions
 <<oldsub:=for each el in oldsub join nextdy(revx,xlist,el);
   allsub:=cons(oldsub,allsub)
 >>;
 return allsub
end$

%-------------

algebraic procedure nextdy(revx,xlist,dy)$
% generates all first order derivatives of lhs dy
% revx = reverse xlist; xlist is the list of variables;
%                       dy the old derivative
begin
  scalar x,n,ldy,rdy,ldyx,sublist;
  x:=first revx; revx:=rest revx;
  sublist:={};
  ldy:=lhs dy;
  rdy:=rhs dy;
 
  while lisp(not member(prepsq simp!* algebraic x,
	     prepsq simp!* algebraic ldy))
	and (revx neq {}) do 
  <<x:=first revx; revx:=rest revx>>;
 
  n:=length xlist;
  if revx neq {} then                % dy is not the function itself
  while first xlist neq x do xlist:=rest xlist;
  xlist:=reverse xlist;

  % New higher derivatives
  while xlist neq {} do
  <<x:=first xlist;
    ldyx:=df(ldy,x);
    sublist:=cons((lisp reval algebraic ldyx)=
                  mkid(mkid(rdy,!`),n), sublist);
    n:=n-1;
    xlist:=rest xlist
  >>;
  return sublist
end$

%-------------

symbolic procedure combidif(s)$
% extracts the list of derivatives from s: % u`1`1`2 --> (u, 1, 1, 2)
begin scalar temp,ans,no,n1;
  s:=reval s; % to guarantee s is in true prefix form
  temp:=reverse explode s;
  
  while not null temp do
  <<n1:=<<no:=nil;
	  while (not null temp) and (not eqcar(temp,'!`)) do
	  <<no:=car temp . no;temp:=cdr temp>>; 
	  compress no
	>>;
    if (not fixp n1) then n1:=intern n1;
    ans:=n1 . ans;
    if eqcar(temp,'!`) then <<temp:=cdr temp; temp:=cdr temp>>;
  >>;
  return ans
end$

%-------------

symbolic operator dif$
symbolic procedure dif(s,n)$
% e.g.:   dif(fnc!`1!`3!`3!`4, 3) --> fnc!`1!`3!`3!`3!`4
begin scalar temp,ans,no,n1,n2,done;
  s:=reval s; % to guarantee s is in true prefix form
  temp:=reverse explode s;
  n2:=reval n;
  n2:=explode n2;

  while (not null temp) and (not done) do
  <<n1:=<<no:=nil;
	  while (not null temp) and (not eqcar(temp,'!`)) do
	  <<no:=car temp . no;temp:=cdr temp>>; 
	  compress no
	>>;
    if (not fixp n1) or ((fixp n1) and (n1 leq n)) then
    <<ans:=nconc(n2,ans); ans:='!` . ans; ans:='!! . ans; done:=t>>;
    ans:=nconc(no,ans);
    if eqcar(temp,'!`) then <<ans:='!` . ans; ans:='!! . ans; 
			      temp:=cdr temp; temp:=cdr temp>>; 
  >>;
  return intern compress nconc(reverse temp,ans);
end$

%-------------

algebraic procedure depnd(y,xlist)$
for each xx in xlist do
for each x  in xx    do depend y,x$

%==== Other procedures:

symbolic operator totdif$
symbolic procedure totdif(s,x,n,dylist)$
% total derivative of s(x,dylist) w.r.t. x which is the n'th variable
begin
  scalar tdf,el1,el2;
  tdf:=simpdf {s,x};
  <<dylist:=cdr dylist;
    while dylist do
    <<el1:=cdar dylist;dylist:=cdr dylist;
      while el1 do
      <<el2:=car el1;el1:=cdr el1;
	tdf:=addsq(tdf ,multsq( simp!* dif(el2,n), simpdf {s,el2}))
      >>
    >>
  >>;
  return prepsq tdf
end$

%-------------

algebraic procedure simppl(pllist,ulist,tt,xx)$
begin
 scalar pl,hh,td,xd,lulist,ltt,lxx,ltd,dv,newtd,e1,deno,ok,
        newpllist,contrace; 
% contrace:=t;
 lisp <<
  lulist:=cdr reval algebraic ulist;
  lxx:=reval algebraic xx;
  ltt:=reval algebraic tt;
 >>;
 newpllist:={};
 for each pl in pllist do <<
  td:=first  pl;
  xd:=second pl;
  repeat <<
   lisp <<
    ltd:=reval algebraic td;
    if contrace then <<write"ltd1=",ltd;terpri()>>$
    dv:=nil;
    newtd:=nil;
    deno:=nil;
    if (pairp ltd) and (car ltd='QUOTIENT)   and
      my_freeof(caddr ltd,ltt) and
      my_freeof(caddr ltd,lxx) 
    then <<deno:=caddr ltd;ltd:=cadr ltd>>;
    ok:=t;

    if (pairp ltd) and (car ltd = 'PLUS) then ltd:= cdr ltd else
    if (pairp ltd) and (car ltd neq 'TIMES) then ok:=nil
                                            else ltd:=list ltd;   
    if contrace then <<write"ltd2=",ltd;terpri()>>$
    if ok then <<
     for each e1 in ltd do <<
      hh:=intpde(e1, lulist, list(lxx,ltt), lxx, t);
      if null hh then hh:=list(nil,e1);
      dv   :=cons(car hh,dv);
      newtd:=cons(cadr hh,newtd);
     >>;
     dv   :=reval cons('PLUS,dv);
     newtd:=reval cons('PLUS,newtd);
     if deno then <<newtd:=list('QUOTIENT,newtd,deno);
                    dv   :=list('QUOTIENT,dv   ,deno) >>;
     if contrace then <<write"newtd=",newtd;terpri();
                        write"dv=",dv      ;terpri() >>$

     td:=newtd;
     if contrace then <<write"td=",td;terpri()>>$
     if (dv neq 0) and (dv neq nil) then <<
      xd:=reval(list('PLUS,xd,list('DF,dv,tt)));
      if contrace then <<write"xd=",xd;terpri()>>$
      %algebraic mode:
      %hh:=lisp gensym()$
      %sbb:=absorbconst({td*hh,xd*hh},{hh})$
      %if (sbb neq nil) and (sbb neq 0) then
      %<<td:=sub(sbb,td*hh)/hh;  xd:=sub(sbb,xd*hh)/hh>>;
      % cllist would have to be scaled as well
     >>
    >>
   >>
  >>
  until lisp(dv)=0;
  newpllist:=cons({td,xd}, newpllist);
 >>;
 return reverse newpllist
end$ % simppl

%-------------

symbolic operator fdepterms$
symbolic procedure fdepterms(td,f)$
% fdepterms regards td as a fraction where f occurs only in the
% numerator. It determines all terms of the numerator in
% which f occurs divided through the denominator.
begin
  scalar nu,de,e1,sm;
  td:=reval td;
  if pairp td then
  if car td='QUOTIENT then <<nu:=cadr td;de:=caddr td>>;
  if null nu then nu:=td;
  if not pairp nu then if freeof(nu,f) then sm:=0
				       else sm:=nu
		  else <<
    if car nu = 'PLUS then nu:=cdr nu
		      else nu:=list nu;
    for each e1 in nu do
    if not freeof(e1,f) then sm:=cons(e1,sm);
    if null sm then sm:=0 else
    if length sm = 1 then sm:=car sm
		     else sm:=cons('PLUS,sm)
  >>;
  if de then sm:=list('QUOTIENT,sm,de);
  return sm
end$ % of fdepterms             

%-------------

symbolic procedure subtract_diff(d1,d2)$
% assumes d1,d2 to be equally long lists of numbers (at least one)
% that are orders of derivatives (which may be 0),
% These lists ca be produced using the procedure maxderivs(),
% returns nil if any number in d2 is bigger than the corresponding
% number in d1, returns list of differences otherwise
begin scalar d;
 return
 if car d2 > car d1 then nil else
 if null cdr d1 then {car d1 - car d2} else 
 if d:=subtract_diff(cdr d1,cdr d2) then cons(car d1 - car d2,d)
                                    else nil
end$

%-------------

symbolic procedure transfer_fctrs(h,flist)$
begin scalar fctrs;
%algebraic write"begin: caar h=",lisp caar h," cdar h =",lisp cdar h;
 if (pairp cdar h) and (cadar h='MINUS) then 
 rplaca(h,cons(reval {'MINUS,caar h},cadr cdar h));

 if (pairp cdar h) and (cadar h='TIMES) then 
 for each fc in cddar h do
 if freeoflist(fc,flist) then fctrs:=cons(fc,fctrs);
 if fctrs then <<
  if cdr fctrs then fctrs:=cons('TIMES,fctrs)
               else fctrs:=car fctrs;
  rplaca(h,cons(reval {'TIMES   ,caar h,fctrs},
                reval {'QUOTIENT,cdar h,fctrs} ))
 >>
%;algebraic write"end:   caar h=",lisp caar h," cdar h =",lisp cdar h;
end$

%-------------

symbolic operator partintdf$
symbolic procedure partintdf(eqlist,qlist,plist,xlist,flist,jlist,sb)$
% eqlist ... list of equations
% qlist  ... list of characteristic functions
% plist  ... list of components of conserved current
% xlist  ... list of independent variables
% flist  ... list of the arbitrary function occuring in this conservation law
% jlist  ... list of all jet-variables
% eqlist and qlist are in order.
% plist and xlist are in order.
% The aim is to remove all derivatives of f in the conservation law
% At first terms with derivatives of f in qlist are partially integrated.
% Then terms with derivatives of f in plist are partially integrated.
begin scalar f,n,d,deltali,subli,lhs,rhs,cof,x,y,cpy,newpl,lowd,su,vle,
             idty,idtysep,sbrev,dno,lsb,h0,h1,h2,h3,h4,h5,h6,h7,ldh1,ldh2,
             reductions_to_do,ld1,ld2,h0_changed;

 % 0. check that plist is homogeneous in flist
 algebraic <<
  cpy:=plist$
  for each f in flist do cpy:=sub(f=0,cpy)$
  while (cpy neq {}) and (first cpy = 0) do cpy:=rest cpy$ 
 >>$
 if cpy neq {'LIST} then return nil$

 eqlist:=cdr eqlist$
 qlist :=cdr  qlist$
 plist :=cdr  plist$
 xlist :=cdr  xlist$
 flist :=cdr  flist$
 jlist :=cdr  jlist$

 % 0. check that flist functions do only depend on xlist variables
 d:=t;
 for each f in flist do 
 if not_included(fctargs f,xlist) then d:=nil$
 if null d then return nil$

 terpri()$
 write"An attempt to factor out linear differential operators:"$terpri()$
 n:=0;
 while eqlist do <<
  n:=add1 n;
  su:=print_;print_:=nil;
  d:=newfct('eq_,xlist,n);
  print_:=su;
  deltali:=cons(d,deltali);
  algebraic write d,":=",lisp car eqlist$
  subli:=cons({'EQUAL,d,car eqlist},subli)$
  lhs:=cons({'TIMES,car qlist,d % car eqlist
     },lhs);
  eqlist:=cdr eqlist;
  qlist:=cdr qlist
 >>;
 lhs:=reval cons('PLUS,lhs)$
 subli:=cons('LIST,subli)$
 for each f in flist do << 
  f:=reval f$

  % removing f-derivatives from the lhs
  repeat <<
   d:=car ldiffp(lhs,f)$ %  liefert Liste der Variablen + Ordnungen mit Potenz
   if d then <<
    % correcting plist
    cpy:=d;
    while cpy and ((numberp car cpy) or freeof(xlist,car cpy)) do cpy:=cdr cpy;
    if null cpy then d:=nil
		else <<
     cof:=coeffn(lhs,cons('DF,cons(f,d)),1);
     lhs:=reval {'DIFFERENCE,lhs,cons('DF,cons({'TIMES,cof,f},d))}$
     x:=car cpy;
     lowd:=lower_deg(d,x)$ % the derivative d reduced by one
     su:=if lowd then cons('DF,cons({'TIMES,cof,f},lowd))
		 else               {'TIMES,cof,f}$

     cpy:=xlist;
     newpl:=nil;
     while cpy and (x neq car cpy) do <<
      newpl:=cons(car plist,newpl);
      plist:=cdr plist;
      cpy:=cdr cpy
     >>;
     plist:=cons({'DIFFERENCE,car plist,su},cdr plist);
     while newpl do <<
      plist:=cons(car newpl,plist)$
      newpl:=cdr newpl
     >>
    >>
   >>
  >> until null d;     % until no derivative of f occurs
  plist:=cdr algebraic(sub(subli,lisp cons('LIST,plist)))$

  % Now we add trivial conservation laws in order to get rid of
  % derivatives of f in the conserved current
  repeat <<
   newpl:=nil;
   cpy:=xlist;
   while plist and null(d:=car ldiffp(car plist,f)) do <<
    newpl:=cons(car plist,newpl);
    plist:=cdr plist;
    cpy:=cdr cpy
   >>;
   if d and (car d neq car cpy) then <<   % otherwise infinte loop
    cof:=coeffn(car plist,cons('DF,cons(f,d)),1);
    x:=car d;
    lowd:=lower_deg(d,x)$ % the derivative d reduced by one
    su:=if lowd then {'TIMES,cof,cons('DF,cons(f,lowd))}
		else {'TIMES,cof,              f       }$

    plist:=cons(reval reval {'DIFFERENCE,car plist,{'DF,su,x}},cdr plist);
    while newpl do <<
     plist:=cons(car newpl,plist)$
     newpl:=cdr newpl
    >>$

    % adding the correction to the other component of plist
    y:=car cpy;
    cpy:=xlist;
    while x neq car cpy do <<
     newpl:=cons(car plist,newpl);
     plist:=cdr plist;
     cpy:=cdr cpy
    >>$
    plist:=cons(reval reval {'PLUS,car plist,{'DF,su,y}},cdr plist);
    while newpl do <<
     plist:=cons(car newpl,plist)$
     newpl:=cdr newpl
    >>
   >> else <<d:=nil;plist:=append(reverse newpl,plist)>>
  >> until null d;
 >>;

 vle:=length xlist;

 newpl:=algebraic absorbconst(lisp cons('LIST,append(qlist,plist)),
                                   cons('LIST,flist))$
 if newpl then newpl:=cdadr newpl;

 % Now factorizing out a linear differential operator
 % 2. extend dependencies of functions from flist and add extra conditions
 for each f in flist do <<
  depl!*:=delete(assoc(f,depl!*),depl!*);
  depl!*:=cons(cons(f,xlist),depl!*);
 >>$
 % 3. compute coefficients of the conditions in the identity 
 idty:=algebraic(sub(subli,lhs))$
 for n:=1:vle do 
 if not zerop nth(plist,n) then 
 idty:={'DIFFERENCE,idty,{'DF,nth(plist,n),nth(xlist,n)}}$
 % 4. separate idty into conditions with multiplicities
 sbrev:=cons('LIST,for each d in cdr sb collect {'EQUAL,caddr d,cadr d})$
 idty:=reval reval idty$
 dno:=algebraic den idty;
 if dno neq 1 then idty:=algebraic num idty$

 idty:=algebraic(sub(sbrev,idty))$
 su:=print_;print_:=nil;
 idtysep:=separ(reval idty,flist,jlist,nil)$
 print_:=su;
 idtysep:=for each d in idtysep collect 
 cons(algebraic(sub(sb,lisp car d)),cdr d);

 % 5. integrations of cdr of the elements of idty have to be done:
 %    - sufficiently often so that there are not more conditions
 %      than functions in flist
 %    - as few as possible to have factored out afterall an as
 %      high as possible operator

 reductions_to_do:=length idtysep - length flist;
 if reductions_to_do>0 then <<
  h0:=idtysep;
  while h0 do <<
   rplaca(h0,cons(reval caar h0, reval cdar h0));
   transfer_fctrs(h0,flist); h0:=cdr h0
  >>$

%write"Separation gives:"$terpri()$
%for each d in idtysep do 
%algebraic write "0 = (",lisp car d,") * (",lisp cdr d,")"$

  h0:=idtysep;
  repeat <<  % check whether cdar h0 is a derivative of another condition
   h0_changed:=nil;
   h1:=cdar h0;
%algebraic write"caar h0=",lisp caar h0," cdar h0 =",lisp cdar h0;
   % find a function appearing in h1 and its leading derivative
   cpy:=flist;
   while cpy and freeof(h1,car cpy) do cpy:=cdr cpy;
   % if null cpy then error!
     
   ld1:=car ldiffp(h1,car cpy)$
   ldh1:=maxderivs(nil,ld1,xlist)$
   ld1:=if null ld1 then car cpy
                    else cons('DF,cons(car cpy,ld1))$
     
   h2:=idtysep;
   while h2 do 
   % is h1 a derivative of car h2 or car h2 a derivative of h1?
   if (h2 eq h0) or freeof(cdar h2,car cpy) then h2:=cdr h2
                                            else <<

%algebraic write"caar h2=",lisp caar h2," cdar h2 =",lisp cdar h2;
    ld2:=car ldiffp(cdar h2,car cpy)$
    ldh2:=maxderivs(nil,ld2,xlist)$
    ld2:=if null ld2 then car cpy
                     else cons('DF,cons(car cpy,ld2))$

    % is h1 a derivative of car h2?
    h3:=subtract_diff(ldh1,ldh2);
    if null h3 then h2:=cdr h2
               else << 
     % the leading derivative in h1 is a derivative of 
     % the leading derivative in cdar h2
     h4:=cdar h2;
%write"h4=",h4;terpri()$
     if pairp h4 and (car h4 = 'PLUS) then <<
      for n:=1:vle do if not zerop nth(h3,n) then
      h4:={'DF,h4,nth(xlist,n),nth(h3,n)};
      if null freeoflist(h5:=algebraic(h1/h4),flist) then h2:=cdr h2
                                                     else <<
       % h1 = h5 * derivative of (cdar h2)
       h6:={'TIMES,caar h0,reval h5};
       for n:=1:vle do <<
        h7:=nth(h3,n);
        if not zerop h7 then
        h6:={'TIMES,{'EXPT,-1,h7},{'DF,h6,nth(xlist,n),h7}};
       >>;
       rplaca(h2,cons(reval {'PLUS,caar h2,h6},cdar h2));
       rplaca(h0,cons(0,0));
%algebraic write"Change(1):"$
%algebraic write"caar h2=",lisp caar h2," cdar h2 =",lisp cdar h2;
%algebraic write"caar h0=",lisp caar h0," cdar h0 =",lisp cdar h0;
       reductions_to_do:=sub1 reductions_to_do;
       h2:=nil
      >>
     >>                               else <<
      % Update of car h2
      h6:=algebraic(lisp(caar h0)*coeffn(h1,ld1,1));
      for n:=1:vle do <<
       h7:=nth(h3,n);
       if not zerop h7 then
       h6:={'TIMES,{'EXPT,-1,h7},{'DF,h6,nth(xlist,n),h7}};
      >>;
      rplaca(h2,cons(reval {'PLUS,caar h2,h6},cdar h2));
%;algebraic write"Change(2):"$
%algebraic write"caar h2=",lisp caar h2," cdar h2 =",lisp cdar h2;
      % Update of car h0
      h1:=reval {'DIFFERENCE,h1,{'TIMES,coeffn(h1,ld1,1),ld1}}$

      if zerop h1 then <<rplaca(h0,cons(0,0));h2:=nil;
                         reductions_to_do:=sub1 reductions_to_do;>>
                  else <<rplaca(h0,cons(caar h0,h1));
                         transfer_fctrs(h0,flist);
                         h1:=cdar h0;
                         cpy:=flist;
                         while cpy and freeof(h1,car cpy) do cpy:=cdr cpy;
                         ld1:=car ldiffp(h1,car cpy)$
                         ldh1:=maxderivs(nil,ld1,xlist)$
                         ld1:=if null ld1 then car cpy
                                          else cons('DF,cons(car cpy,ld1))$
                         h2:=cdr h2;h0_changed:=t>>
%;algebraic write"caar h0=",lisp caar h0," cdar h0 =",lisp cdar h0;
     >>

    >>
   >>;
   if (null h0_changed) or (zerop caar h0) then h0:=cdr h0
  >> until (reductions_to_do=0) or (null h0);

%write"After correction the separation gives:"$terpri()$
%for each d in idtysep do 
%if not zerop car d then
%algebraic write "0 = (",lisp car d,") * (",lisp cdr d,")"$

 >>$

 % Now the number of f in flist should be equal the number of conditions
 % or as low as possible
 n:=0;
 rhs:=nil;
 for each d in idtysep do 
 if not zerop car d then << % for each condition
  n:=add1 n;
  su:=print_;print_:=nil;
  x:=newfct('l_,xlist,n);
  print_:=su;
  su:=if dno=1 then car d 
               else reval {'QUOTIENT,car d,dno}$
  algebraic write x,":=",su$
  lsb:=cons({'EQUAL,x,su},lsb);
  % 5. for each condition integrate all terms
  y:=cdr d;
  cpy:=flist;
  while y and not zerop y do <<
   repeat <<
    d:=ldiffp(y,car cpy)$
    if zerop cdr d then
    if null cpy then <<write"The backintegration is faulty."$terpri()>>
                else cpy:=cdr cpy
   >> until not zerop cdr d;
   if car d = nil then <<
    cof:=coeffn(y,car cpy,1);
    rhs:={'PLUS,{'TIMES,x,cof,car cpy},rhs};
    y:=reval reval {'DIFFERENCE,y,{'TIMES,cof,car cpy}}
   >>             else <<
    cof:=coeffn(y,cons('DF,cons(car cpy,car d)),1);
    rhs:=reval {'PLUS,rhs,{'TIMES,cons('DF,cons({'TIMES,x,cof},car d)),
                                  car cpy,{'EXPT,{'MINUS,1},absdeg(car d)}}};
    y:=reval reval {'DIFFERENCE,y,{'TIMES,cof,cons('DF,cons(car cpy,car d))}}
   >>
  >>
 >>$
 lsb:=cons('LIST,lsb)$
 flist:=cons('LIST,flist)$
 algebraic <<
  d:=gcd(den lhs,den rhs);
  lhs:=lhs*d; rhs:=rhs*d;

  %--- Correctness test
  d:=sub(subli,lhs)-sub(lsb,rhs);
  if d neq 0 then write "Not identically zero : ",d$

  for each f in flist do algebraic <<
   x:=coeffn(num lhs,f,1);  y:=coeffn(num rhs,f,1);
   d:=gcd(x,y);
   algebraic write x/d/den lhs," = ",y/d/den rhs$
  >>  
 >>$

end$

%-------------

end$