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

%*********************************************************************
module gensep_lin$
%*********************************************************************
%  Routines for generalized separation of de's
%  Author: Andreas Brand, Thomas Wolf 1990 1994 1997
%          Thomas Wolf since 1997

symbolic procedure quick_gen_separation(arglist)$
% Indirect separation of a pde
if vl_ then % otherwise not possible --> save time
begin scalar p,l,l1,pdes,stp$
 % pdes:=clean_up(car arglist)$
 % if pdes then l:=list(pdes,cadr arglist)$ 
 % because the bookeeping of to_drop is not complete instead:
 pdes:=car arglist$
 % to return the new list of pdes in case gensep is not successful
 if expert_mode then <<
   l1:=selectpdes(pdes,1);
   if get(car l1,'starde) then flag(l1,'to_gensep)
 >>             else l1:=cadddr arglist$
 if (p:=get_gen_separ_pde(l1,t,t)) then % high priority
    <<l:=gensep(p,pdes)$
      pdes:=delete(p,pdes)$
      for each a in l do <<
        pdes:=eqinsert(a,pdes)$
        if member(a,pdes) and (stp:=get(a,'starde)) then
        to_do_list:=cons(list(if cdr stp=0 then 'separation
                                           else 'gen_separation,
                              %pdes,cadr arglist,caddr arglist,
                              list a),
                         to_do_list)
      >>$
      l:=list(pdes,cadr arglist)>>$
 return l$
end$

symbolic procedure gen_separation(arglist)$
% Indirect separation of a pde
if vl_ then % otherwise not possible --> save time
begin scalar p,l,l1,pdes,stp$
 % pdes:=clean_up(car arglist)$
 % if pdes then l:=list(pdes,cadr arglist)$ 
 % because the bookeeping of to_drop is not complete instead:
 pdes:=car arglist$
 % to return the new list of pdes in case gensep is not successful
 if expert_mode then <<
   l1:=selectpdes(pdes,1);
   if get(car l1,'starde) then flag(l1,'to_gensep)
 >>             else l1:=cadddr arglist$
 if (p:=get_gen_separ_pde(l1,nil,t)) then % low priority
    <<l:=gensep(p,pdes)$
      pdes:=delete(p,pdes)$
      for each a in l do <<
        pdes:=eqinsert(a,pdes)$
        if member(a,pdes) and (stp:=get(a,'starde)) then
        to_do_list:=cons(list(if cdr stp=0 then 'separation
                                           else 'gen_separation,
                              %pdes,cadr arglist,caddr arglist,
                              list a),
                         to_do_list)
      >>$
      l:=list(pdes,cadr arglist)>>$
 return l$
end$

symbolic procedure maxnoargs(fl,v)$
% determines the maximal number of arguments of any of the 
% functions of fl
begin scalar f,n,m;
  n:=0;
  for each f in fl do
  <<m:=fctargs f;
    m:=if m and not_included(v,m) then length m
            else 0;
    if n<m then n:=m
  >>;
  return n
end$

symbolic procedure get_gen_separ_pde(pdes,high_priority,lin)$
% looking for a pde in pdes which can be indirectly separated
% if lin then only a liner PDE
% p ...the next equation that will be chosen
% dw...whether p was already delt with
% na...number of variables in the equation
% nv...maximal number of arguments of any of the functions of p
% nf...min number of functions to be dropped before direct sep.
% leng...length of p
begin scalar p,nv,nf,dw,len,h1,h2,h3,h4,nvmax$ %na,h5
%  ncmax:=nvmax$
  if high_priority then <<
    nvmax:=0;
    for each p in pdes do if (h1:=get(p,'nvars))>nvmax then nvmax:=h1;
    p:=nil
  >>$
  while pdes do <<
    if flagp(car pdes,'to_gensep)                  and
       (null lin or get(car pdes,'linear_))        and
         % not too many terms or enough terms
      <<h1:=get(car pdes,'length)$
        (null high_priority) or
        (get(car pdes,'nvars)=nvmax) or
        ( low_gensep>h1)     or      
        (high_gensep<h1)  
      >>                                           and 
         % no single function depending on all variables:
      (h3:=get(car pdes,'starde)                 ) and
         % not directly separable: 
      (cdr h3 neq 0                              ) and
         % Each time the equation is investigated and differentiated
         % wrt the first element of car h3, this element is dropped.
         % --> The equation should not already have been differentiated
         % wrt all variables:
      (not null car h3                           ) and     
         % If equations have been investigated by generalized
         % separation or if equations resulted from generalized
         % separation then they get the flag used_ to be solved
         % first, not to have too many unevaluated new functions
         % at a time
      ((h4:=flagp(car pdes,'used_)          ) or 
       (null dw)                                 ) and
         % The variables in h3 are the ones wrt which direct separation
         % shall be achieved after differentiation, therefore functions
         % of these variables have to be thrown out. The remaining
         % functions shall be of as many as possible arguments to
         % make quick progress:
      ((null p                              ) or
       (len > h1                            ) or   % neu
       ((len = h1) and (                           % neu
       (nv < (h2:=maxnoargs(
                  get(car pdes,'fcts),
                  car h3   ))               ) or
       ((nv = h2) and (
%       (na < (h5:=get(car pdes,'nvars))     ) or
%       ((na = h5) and (
       ((null dw) and flagp(car pdes,'used_)) or
       (( nf > cdr h3       ) or
        ((nf = cdr h3 ) and          
         (len > h1    )     )    )         )    ))))      then
    <<p:=car pdes;
      nv:=if (null nv) or (null h2) then maxnoargs(get(p,'fcts),car h3)
                                    else h2;
%      na:=if (null na) or (null h5) then get(car pdes,'nvars)
%                                    else h5;
      if h4 then dw:=h4;
      nf:=cdr h3;
      len:=h1>>$
    pdes:=cdr pdes$
  >>;
  return p
end$

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

symbolic procedure gensep(p,pdes)$
%  generalized separation of pde p
if zerop cdr get(p,'starde) then separate(p,pdes)   % be dropped?
else                                                % e.g. a=((x y z).2)
% POSSIBLE REASONS FOR FAILURE: 
% - After the sequence of divisions and differentiations in the direct 
%   separation, if there explicit v-dependent coefficients are taken 
%   out which also contain later integration variables, then the integrands 
%   are not total derivatives anymore --> no backintegration is possible.
% - This corresponds to the case when all variables occur explicitly but
%   in a non-product expression, like sin(x*y*z)
begin    scalar a,pl$
  if print_ then <<terpri()$write "generalized separation of ",p>>$
  if tr_gensep then
  <<a:=get(p,'starde)$
    terpri()$write "de to be separated : "$
    typeeqlist list p$
    terpri()$write "variable list for separation : ",car a$
    terpri()$write "for each of these variables ",cdr a$
    if (cdr a)=1 then write" function does"
                 else write" functions"$
    write" depend on it "
  >>$

  %--- so far only one DE p in the pool starlist of equations
  pl:=partitn(get(p,'val),    % expression
              nil,            % history of divisions/diff so far
              get(p,'fcts),   % functions
              get(p,'vars),   % variables
              car get(p,'starde)  % separation charac.
             );

  if pl then <<
    pl:=append(for each a in car pl collect cdr a,cadr pl);
    pl:=mkeqlist(pl,fctsort union(ftem_,get(p,'fcts)),get(p,'vars),
                 cons('to_drop,allflags_),t,get(p,'orderings),pdes)$
    drop_pde(p,nil,nil);
    flag(pl,'used_);
    if print_ then <<terpri()$
      a:=length pl$
      write "separation yields ",a," new equation"$
      if a > 1 then write"s"$ write" : "$
      if tr_gensep then typeeqlist pl
                   else listprint pl$
      terpri()
    >> 
  >> else <<
    remflag1(p,'to_gensep)$
    pl:=list p
  >>$
  return pl$
end$

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

symbolic procedure partitn(q,old_histy,ftem,vl,a)$
%       This procedure calls itself recursively!
% q:    a **-expression to be separated
% old_histy: a list of elements {denom,v,{(divisor . expr_before),..}}
%       where a sequence of divisions through factors from the
%       list of divisors and differentiations wrt. v and 
%       afterwards multiplication with denom resulted in q
% ftem: functions in the expression
% vl:   variables in the expression
% a:    the variables for direct separation=car starp()
%
% RETURNS {{(c1.q1),(c2.q2),(c3.q3),..},{s1,s2,s3,..},{r1,r2,..},{f1,f2,..}}
% where qi=0 are necessary consequences,
% qi are not *-expressions,
% sum_i ci*qi = q
% si=0 are consistency conditions determining constants/functions
% of integration 
% ri=0 are other cases to be checked --> case distinctions
begin scalar histy,l1,l4,nv,vl1,nv1,h,x,f,ft,aa,bb,cc,y,
      ruli,extra_cond,par,cases,newf$

  %--- ft: the list of functions to drop from q by differentiation
  %--- to do a direct separation wrt x:
  % x = any one variable in a on which a function with as 
  % many as possible variables does not depend on
  % Find such a function and variable x
  ft:=ftem;    
  nv:=0;
  while ft do <<
    vl1:=fctargs car ft;
    nv1:=if vl1 then length vl1
                else 0;
    if nv1 > nv then <<
      h:=setdiff(a,vl1);
      if h then <<
        x:=car h;
        % if possible find an x occuring explicitly in q
        while h and freeof(q,car h) do h:=cdr h;
        if h then x:=car h;
        f:=car ft;
        nv:=nv1
      >>
    >>;
    ft:=cdr ft
  >>; 
  if nv=0 then x:=car a;  % no x was found

  if tr_gensep then
  <<terpri()$write "--- The aim is to separate directly w.r.t. ",x$
    write " the expression : "$deprint list q  >>$

  % Find all functions ft in q depending on x
  ft:=nil$
  for each f in ftem do
  if member(x,fctargs f) and not freeof(q,f) then ft:=cons(f,ft)$
  ft:=fctsort reverse ft$      % sorting w.r.t. number of args

  ruli:=start_let_rules()$

  %--- throwing out functions ft until ft=nil
  %--- or until the expression lost the *-property


  while ft do % for each function to get rid of 
              % (possibly each time a different diff variable)   
  if null (l1:=felim(q,car ft,ftem,vl)) then ft:=nil  % to stop
                                        else <<
%prettyprint l1;
    for each h in cdadr l1 do               % special extra cases
    if not freeoflist(car h,ftem) then cases:=cons(car h,cases);

%write"cadr l1=",cadr l1$terpri()$
    if zerop car l1 then <<
      q:=reval {'QUOTIENT,cdr cadadr l1,car cadadr l1}; % new expression
      cc:=cons(car cadr l1,cddadr l1);
    >>              else <<
      q:=car l1$                            % new expression
      cc:=cadr l1;
    >>$
    if (pairp q) and (car q='QUOTIENT) then <<
      bb:=caddr q;                          % we take off the denimonator
      q:=cadr q
    >>                                 else bb:=1$
    histy:=cons(cons(bb,cc),histy)$         % extended history
%terpri()$write"q=",q$terpri()$
%write"histy=",histy$terpri()$

    ftem:=smemberl(ftem,q)$                 % functions still in q
    aa:=stardep(ftem,argset(ftem))$         % still *-expression?
    if not aa or 
       zerop cdr aa then ft:=nil                   % to stop
                    else ft:=smemberl(cdr ft,ftem) % remain. fcts of x
  >>$

  stop_let_rules(ruli)$
  if null l1 then if tr_gensep then <<terpri()$
                    write"felim or newgensep gave nil!!"$terpri()$
                    write"q=",q;terpri()
                  >>           else
             else
  RETURN <<
    if tr_gensep then <<
      terpri()$
      write"Now ready for direct separation."
    >>;

    %--- prepare list of variables vl1 for direct separation
    vl1:=nil$
    for each y in vl do if my_freeof(ftem,y) then vl1:=cons(y,vl1);

    %--- direct separation if useful (i.e. if aa(=stardep) neq nil)
    if vl1 and not (q=0) then
    <<if tr_gensep then
      <<terpri()$write "trying direct separation of "$
        deprint list q$
        write "Remaining variables: ",vl1>>$
      l1:=separ(q,ftem,vl1,nil)$   % direct separation of the numerator
      if tr_gensep then
      <<terpri()$
        write "The result of direct separation: "$
        deprint for each bb in l1 collect cdr bb>>$
    >>             else l1:=list cons(1,q)$

    if tr_gensep then <<
      terpri()$
      write"Separation gave ",length l1," condition(s)"
    >>;

    % Although the vaiable x does not occur anymore
    % (each felim call removed one function of x
    %  and direct separation removed explicit occurences of x)
    % the remaining expression may still be indirectly separable
    % --> recursive call of partitn
    % l4 becomes a list of pairs (sep_coefficient . sep_remainding_factor) 

    for each h in l1 do <<
      ft:=smemberl(ftem,cdr h);
      vl1:=argset(ft)$
      if null (aa:=stardep(ft,vl1)) then l4:=cons(h,l4)
                                    else <<
        par:=partitn(cdr h,   % expression
                     append(histy,      % history so far, 
                            old_histy), % needed to add new functions
                              % of integration properly differentiated to be
                              % able to integrate below
                     ft,      % functions
                     vl1,     % variables
                     car aa   % separation charac.
                    );
	% RETURNS {{(c1.q1),(c2.q2),(c3.q3),..},{s1,s2,s3,..},
	%          {r1,r2,..},{f1,f2,..}                      }
	% where qi=0 are necessary consequences,
	% qi are not *-expressions,
	% sum_i ci*qi = q
	% si=0 are consistency conditions determining constants/functions
	% of integration 
	% ri=0 are other cases to be checked --> case distinctions

        l4:=append(l4,for each f in car par collect 
                        ({'TIMES,car h,car f} . cdr f));
        extra_cond:=append(extra_cond,cadr par); % collecting any
                                                 % appearing conditions
        cases:=append(cases,caddr par);
        newf:=cadddr par;
        ftem:=append(ftem,newf);
      >>
    >>$

    %--- backintegration
    par:=backint(l4,old_histy,histy,ftem,vl)$
    extra_cond:=append(extra_cond,cadr par); % collecting any conditions

    {car par,extra_cond,cases,append(newf,caddr par)}
  >>
end$

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

symbolic procedure felim(q,f,ftem,vl)$
  % returns:  nil if not successful (quotient) otherwise
  %  {the expression after all the divisions and differentiations,
  %   {diff variable, sequence of (factor . expression before)}   }
begin scalar a,b,l,l1,ft1,v,prflag$

  %--- getting rid of f through diff. wrt. v
  v:=car setdiff(vl,fctargs f)$  

  %--- ft1 are all v-independent functions
  %--- In the call to separ one has to disregard v-dep. functions
  ft1:=nil$
  for each f in ftem do if my_freeof(f,v) then ft1:=cons(f,ft1)$

  %--- To run separ, functions ft1 should not be in the denominator
  %--- ?????? What if nonl. Problems?
  if not (pairp q and (car q='QUOTIENT) and smemberl(ft1,caddr q)) then
  % This exceptional case should not occure anymore
  <<prflag:=print_$print_:=nil$
    l:=separ(q,ft1,list v,nil)$ % det. all lin. ind. factors with v
    for each a in reverse idx_sort for each b in l 
                                   collect cons(delength car b,b)
    collect cdr a$

    if tr_gensep then
    <<terpri()$write "To get rid of ",f,
		     " we differentiate w.r.t. : ",v>>$
    print_:=prflag$

    %--- l is a list of dotted pairs a each representing a term in q 
    % where car(a) is the product of v-dep. factors and cdr(a) the
    % product of v-independent factors
    %--- A list l1 of car(a) is generated for which cdr(a) depends
    % on f. l1 is the list of divisions to be done before differen.
    l1:=nil$
    while l do
    <<a:=car l$
      b:=nil$
      if not freeof(cdr a,f) and (not zerop car a) then
      l1:=cons(car a,l1)$
      l:=cdr l
    >>$
    if tr_gensep then
    <<terpri()$
      write v," - depending coefficients of terms containing ",f," : "$
      for each ss in l1 do eqprint ss>>$

    %--- Now the divisions and differentiations are done
    while l1 do
    <<b:=reval aeval car l1$ %--- b is the v-dep. coefficient
      l1:=cdr l1$
      %--- ????? If for non-linear Problems b includes ftem functions
      % then the extra case 0 = b has to be considered
      if not zerop b then
      <<a:=reval aeval list('QUOTIENT,q,b)$
        %--- for later backward integrations: extension of the history
        l:=cons(b . q ,l)$  %--- new: q is the equ. before division & diff.
                            %    formerly: l:=cons(b ,l)$      
        % l will be returned later by felim 
        %--- l1 has to be updated as the coefficients
        % change through division and differentiation
	l1:=for each c in l1 collect
	    reval list('DF,list('QUOTIENT,c,b),v)$
        %--- the differentiation
	q:=reval list('DF,a,v)$
	if tr_gensep then
	<<write "The new equation: "$
	  eqprint q$
	  write "The remaining factors:"$
	  for each ss in l1 do eqprint ss
	>>
      >>
    >>$
    %if l then part_histy:=cons(v,l)$

    %--- output
    if tr_gensep then
    <<terpri()$write "new expression (should not depend on ",f,") : "$
      eqprint q$>>$
    if tr_gensep and l then
    <<write"To invert the last steps one has to integr. wrt. ",v$ 
      terpri()$
      write "each time before multiplying with "$
      for each aa in l do eqprint car aa
    >>$

    l1:=list(q,cons(v,l))
  >>$
  return l1
end$

symbolic procedure backint(l4,old_histy,histy,ftem,vl)$
% l4 is a list of pairs (sep_coefficient . 
%                        sep_remainding_factor_to_be_integrated) 
% old_histy, histy are lists of elements 
%       {denom,v,{(divisor . expr_before),..}}
%       where a sequence of divisions through factors from the
%       list of divisors and differentiations wrt. v and 
%       afterwards multiplication with denom resulted in q
%       Integrations should only be done inverting histy, but
%       in choosing functions of integration, both should be used
%
% returns {- integrated equivalent of l4 where the cdr of each element
%            is the integrated expression,
%          - a list of check_sum conditions,
%          - a list of new functions}
begin scalar succ,ft,q,l,v,v1,vf,s1,s2,p,f1,f2,fctr,check_sum,
             allfnew,new_cond,denomi$
  % start of the backintegration
  succ:=t$
  while histy and succ do
  <<l:=car histy$     % l={diff variable, sequence of (factor . exp before)}
    histy:=cdr histy$
    denomi:=car l$
    v:=cadr l$
    l:=cddr l$

    % At first putting the denominator back in
    % l4 is a list of pairs (sep_coefficient . sep_remainding_factor) 
    if denomi neq 1 then
    l4:=for each h in l4 collect (car h . {'QUOTIENT,cdr h,denomi});

    if tr_gensep then <<terpri()$ write "backward integration w.r.t. ",v>>$

    % Now the sequence of integrations wrt v
    % l is the list of (factor . expression before division & diff)
    while l do <<         % while l and q do
      fctr:=caar l;
      check_sum:=cdar l;
      l:=cdr l;
      if tr_gensep then
      <<terpri()$
        write "The integrals of the following partitioned subexpressions"$
        terpri()$
        write "added up should be equal the original expression: "$
        terpri()$
        eqprint check_sum
      >>$
%write"l4="$
%prettyprint l4; 
      % l4 is a list of pairs (sep_coefficient . sep_remainding_factor) 
      l4:=for each h in l4 collect if null car h then h % one integration
                                                        % was not succ.ful
                                                 else <<
        ft:=smemberl(ftem,cdr h)$
        fnew_:=nil$
        if tr_gensep then
        <<terpri()$
          write "Backintegration of: "$eqprint cdr h
        >>$
        q:=integratepde(cdr h,ft,v,nil,nil)$ % genflag:=nil, potflag=nil
        if null q then <<
          succ:=nil$
          if print_ then <<
            terpri()$
            write "#### Back integration of "$
            eqprint cdr h$
            write " wrt ",v," during generalized ",
                  "separation was not successful ####."$
            terpri()$
            write "The coeff. dropped in direct separation was "$
            mathprint car h
          >>
        >>        else <<
          q:=reval list('TIMES,fctr,car q)$  
          % fctr is the next integrating factor
          
          % Neccessary: Substituting the new functions of integration by
          % derivatives of them such that back-integration can be made
          % hat fnew_ nur ein element, d.h. wird nur eine Integration gemacht
          % oder mehrere?
          for each f1 in fnew_ do
          <<f2:=f1$

            vf:=setdiff(vl,fctargs f1)$
            for each s1 in reverse append(histy,old_histy) do
            <<v1:=cadr s1$
              % The following also if not my_freeof(f1,v1)
              % The reason is that divisors may contain variables which
              % are later integration variables
              s2:=reverse cddr s1$
              while s2 do
              <<% only divisions through factors that can be swallowed by f1
                if not smemberl(vf,caar s2) then 
                f2:=list('QUOTIENT,f2,caar s2)$
                if not my_freeof(f1,v1) then f2:=reval list('DF,f2,v1)$
                % actually only dividing through those factors of (caar s2)
                % which depend on v1 and which do not contain variables
                % which f2 does not depent on.
                s2:=cdr s2
              >>$
              if not smemberl(vf,car s1) then f2:=list('TIMES,f2,car s1)$
            >>$
            % the remaining integrations in the current element of histy
            if histy then <<
              s2:=reverse l$ 
              while s2 do
              <<if not smemberl(vf,caar s2) then f2:=list('QUOTIENT,f2,caar s2)$
                if not my_freeof(f1,v1) then f2:=reval list('DF,f2,v1)$
                s2:=cdr s2
              >>;
            >>;

            if f1 neq f2 then
            <<if tr_gensep then <<terpri()$
                                  write f1," is replaced by "$
                                  eqprint f2>>$
              q:=subst(f2,f1,q)$
            >>
          >>$
          allfnew:=append(fnew_,allfnew)$
          ftem:=union(fnew_,ftem);
          if succ then check_sum:={'DIFFERENCE,check_sum,{'TIMES,q,car h}};
          % car h is the coefficient dropped through direct separation
        >>$
        (car h . q) % the value to be collected to give the new l4
      >>;
      check_sum:=reval check_sum$
      if succ then new_cond:=cons(check_sum,new_cond)$
      if succ and tr_gensep then
      <<terpri()$
        write "Consistency condition: "$eqprint check_sum
      >>$
    >>
  >>$
  for each f in allfnew do ftem_:=fctinsert(f,ftem_)$
  if tr_gensep then if succ then <<terpri()$write "yields : "$
			 	   eqprint p$>>
			    else <<terpri()$
				   write "was not successful.">>$
  fnew_:=nil$
  return {l4,new_cond,allfnew}
end$

endmodule$

%*********************************************************************
module gensep_non_lin$
%*********************************************************************
%  Routines for generalized separation of de's
%  Author: Thomas Wolf since 1997

symbolic procedure my_smemberl(p,vl)$
begin scalar l,v;
 for each v in vl do
 if not my_freeof(p,v) then l:=cons(v,l);
 return reverse l
end$

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

symbolic procedure stripcond(conds)$
begin scalar newconds,condi;
 newconds:=nil;
 while conds do <<
  condi:=cdar conds;
  conds:=cdr conds;
  if length condi=1 then condi:=car condi
		    else condi:=cons('PLUS,condi);
  newconds:=cons(condi,newconds)
 >>;
 return newconds
end$

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

symbolic procedure checkli(exlist,condi)$
begin
 scalar ok,isincondi,isinexli,excopy,n;
 ok:=t;
 while condi and ok do <<
  % all i in the condition car condi
  isincondi:=car condi; %smemberl(ilist,car condi);
  n:=length isincondi;
  % are all isincondi contained in one of the elements of exlist?
  excopy:=exlist;
  while excopy and ok do <<
   isinexli:=smemberl(isincondi,car excopy);
   if isinexli then
   if length(isinexli) = n then ok:=nil;
   excopy:=cdr excopy
  >>;
  condi:=cdr condi
 >>;
 return ok
end$

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

symbolic procedure longst(exlist)$
% returns the element of exlist which (is a list and)
% has the most elements
begin
 scalar lo;
 while exlist do <<
  if not lo then lo:=car exlist else
  if length(lo)<length(car exlist) then lo:=car exlist;
  exlist:=cdr exlist
 >>;
 return lo
end$

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

symbolic procedure starequ(n,alindep,blindep)$
% alindep is a list of lists of factors ai which are all non-zero and
%         are all linear independent from each other within such a list
% blindep like alindep
% generates all cases each with all conditions with _i representing
% ai or bi, equations and new functions are not generated
begin
 comment
  The equation to separate has the form 0 = sum_i ai*bi
  where the bi do not depend on some variable x. The
  procedure starequ generates cases:
  cases          ... ( all cases )
  each case      ... ( list of all a-conditions,
		       list of all b-conditions)
  each condition ... ( the ai,bi contained in the condition
		       with _i representing ai and bi )
 ;
 scalar i,j,cases,oldcases,case,ai,bi,ci,oldaconds,oldbconds,
	newaconds,newbcond,newbconds,newacond,
	ilist,cona,conb,unin,el,pri;  % ,oldpri

 % Determine the longest union of two list, one, ai, being element of
 % alindep and one, bi, being from blindep

 %pri:=t;
 i:=0;
 if alindep then for each cona in alindep do
 if blindep then for each conb in blindep do
 if (j:=length union(cona,conb)) > i then <<ai:=cona;bi:=conb;i:=j>>
				     else
	    else % no blindep given
 if (j:=length cona) > i then <<ai:=cona;i:=j>>
			 else
	    else % no alindep given
 if blindep then for each conb in blindep do
 if (j:=length conb) > i then <<bi:=conb;i:=j>>;

 % ai, bi will now be determined 

 % preparation of the sequence ilist of extensions
 ilist:=for i:=1:n collect i;
 if pri then <<write"222";terpri()>>$
 if i neq 0 then <<
  if ai then i:=length ai else i:=0;
  if bi then j:=length bi else j:=0;

  unin:=union(ai,bi);
  % extensions through bch should be done when elements from
  % bi are treated. This is coded in ilist through negative numbers
  ilist:=setdiff(ilist,unin);

  if i>j then <<
   for each el in setdiff(unin,ai) do ilist:=cons(-el,ilist);
   for each el in ai do ilist:=cons( el,ilist)
  >>     else <<
   for each el in setdiff(unin,bi) do ilist:=cons( el,ilist);
   for each el in bi do ilist:=cons(-el,ilist)
  >>;
  ilist:=reverse ilist
 >>;
 % ilist is prepared now
 if pri then <<write"333 ilist=",ilist;terpri()>>$

 % `cases' is a list of cases, each is a dotted pair with
 % the car being the a-conditions and cdr being the b-conditions
 % The first two cases
 i:=car ilist;
 if i<0 then i:=-i;
 ci:=mkid('_,i);
 cases:=list(list(list(list(ci)),       nil     ),
	     list(     nil      , list(list(ci))) );

% oldpri:=print_;
% print_:=nil;
 ilist:=cdr ilist;
 if pri then <<write"555";terpri()>>$
 while ilist do <<
  i:=car ilist;ilist:=cdr ilist;
  if pri then <<write"iii=",i;terpri()>>$

  if i>0 then ci:=mkid('_, i)
	 else ci:=mkid('_,-i);
  if pri then <<
   write"666 car ilist=",i;
   terpri()
  >>$

  % if i>0 then the list of cases is extended with ai else with bi

  oldcases:=cases;
  cases:=nil;
  while oldcases do <<  % for each old case do:
   case:=car oldcases;
   if pri then <<write"old case:",case;terpri()>>$
   oldcases:=cdr oldcases;
   if i>0 then <<
    oldaconds:=car case;
    if pri then <<write"888 oldaconds=",oldaconds;terpri()>>$
    % at first add condition i=0 to the case
    cases:=cons(cons(cons(list(ci),oldaconds), cdr case) ,
		cases);
    if pri then <<write"999 new case=",car cases;terpri()>>$
    % then: - add to each a-condition i
    %       - add one new b-condition with the first element `i'
    %         and furtherelements which are the first elements of
    %         the a-lists which have been extended
    newaconds:=nil;
    newbcond:=list(ci);
    while oldaconds do <<
     j:=caar oldaconds;
     newaconds:=cons(cons(j,cons(ci,cdar oldaconds)),
		     newaconds);
     newbcond:=cons(j,newbcond);
     oldaconds:=cdr oldaconds
    >>;
    if pri then <<write"newaconds=",newaconds,
		       " rev newbcond=",reverse newbcond;
		  terpri()>>$
    cases:=cons(list(newaconds, cons(reverse newbcond,cadr case)) ,
		cases);
    if pri then <<write"000 new case=",car cases;terpri()>>$
   >>            else <<
    oldbconds:=cadr case;
    if pri then <<write"888 oldbconds=",oldbconds;terpri()>>$
    % at first add condition bi=0 to the case
    cases:=cons(list(car case, cons(list(ci),oldbconds)),
		cases);
    if pri then <<write"999 new case=",car cases;terpri()>>$
    % then: - add to each b-condition i
    %       - add one new a-condition with the first element `i'
    %         and further elements which are the first elements of
    %         the b-lists which have been extended
    newbconds:=nil;
    newacond:=list(ci);
    while oldbconds do <<
     j:=caar oldbconds;
     newbconds:=cons(cons(j,cons(ci,cdar oldbconds)),
		     newbconds);
     newacond:=cons(j,newacond);
     oldbconds:=cdr oldbconds
    >>;
    cases:=cons(list(cons(reverse newacond,car case), newbconds),
		cases);
    if pri then <<write"000 new case=",car cases;terpri()>>$
   >>

  >>;
 >>;

 % Throwing out cases which are forbidden by alindep and blindep
 alindep:=for each ci in alindep collect
          for each i in ci collect mkid('_,i);
 blindep:=for each ci in blindep collect
          for each i in ci collect mkid('_,i);
 oldcases:=nil;
 % ilist:=for i:=1:n collect mkid('_,i);
 while cases do <<
  if checkli(alindep,caar cases%,ilist
) then
  if checkli(blindep,cadar cases%,ilist
) then
  oldcases:=cons(car cases,oldcases);
  cases:=cdr cases
 >>;

% print_:=oldpri;
 return oldcases
end$ % of starequ

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

symbolic procedure pickfac(ex,indx)$
% returns the n'th element of ex where indx has the form _n
nth(ex,compress cdr explode indx)$

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

symbolic procedure find_cond(bcons,ai)$
% find the element of bcons with ai as (automatically first) element
% (there must be an b-condition with ai as first element
% if that has not already been dropped
% because ai is not the first element of the a-condition)
begin
 while (pairp bcons) and (caar bcons neq ai) do bcons:=cdr bcons;
 return if pairp bcons then car bcons
		       else nil
end$

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

symbolic procedure starsep(exx,ex,ftem,vl,x)$
% exx is the original expression to be separated
% ex is a *-expression returned from SEPAR
% vl are all variables which really occur in ex or 
%    on which ex depends on
% x is a variable on which the bi do not depend on
%
% RETURNS a list of cases, each case having the form
% {{new constants/functions of integration},
%  eq1,eq2,eq3,...} 
% where eqi are a set of necc and suff conditions
begin
 scalar cases,newcases,acons,bcons,acond,newca,alindep,blindep,aco,bco,
	ai,bi,ci,a1,avars,bvars,i,ili,cilist,ali,n,addex,bcopy,cntr,
	pri;

% pri:=t;

 ili:=for i:=1:length ex collect mkid('_,i);
 n:=length vl;

 % at first determining lists of non-vanishing and linear independent
 % a-factors alindep and b-factors blindep
 % does so far only the obvious test which is useful essentially for 
 % linear problems
 cntr:=0;
 for each ci in ex do <<
   cntr:=add1 cntr;
   if pri then <<
     write"a",cntr," = "$mathprint car ci$
     write"b",cntr," = "$mathprint cdr ci$
   >>$
   if null smemberl(ftem,car ci) then alindep:=cons(cntr,alindep)$
   if null smemberl(ftem,cdr ci) then blindep:=cons(cntr,blindep)$
 >>;
 if alindep then alindep:=list alindep$
 if blindep then blindep:=list blindep$

 % generation of all logical cases with the factors ai,bi
 % substituted by _i
 cases:=starequ(length ex,alindep,blindep);
 if pri then <<terpri()$write"Returned from STAREQU: ",cases>>;

 % newcases will be the new list of all cases
 newcases:=nil;

 while cases do << 
  % car cases is one case with alltogether n conditions which 
  % The conditions for the a-factors are called below acons
  % and for the b-factors bcons.

  acons:= caar cases;
  bcons:=cadar cases;
  cases:= cdr cases;
  if pri then <<write"acons=",acons,"  bcons=",bcons;terpri()>>$

  % The case will now be formulated with the terms of the expression ex
  newca:=nil;
  cilist:=nil;
  addex:=nil;
  bcopy:=nil;

  % extract at first all b-conditions with only one term as they do not
  % need constants of integration and are used for any grade of ex
  while bcons do <<
   if length car bcons=1 then
   newca:=cons(cdr pickfac(ex,caar bcons),newca)
			 else bcopy:=cons(car bcons,bcopy);
   bcons:=cdr bcons
  >>;
  bcons:=bcopy;

  % The a-factors may depend on all variables vl whereas the
  % b-factors do at least not depend on x.
  while acons do <<   % formulating all a-conditions
   aco:=car acons;    % aco encodes one condition for a-factors
   acons:=cdr acons;
   if pri then <<write"aco=",aco;terpri()>>$

   a1:=car aco;       % e.g. a1 = _7
   acond:=list car pickfac(ex,a1);
   if pri then <<write"acond=",acond;terpri()>>$
   % acond becomes a list of all a-factors encoded by aco

   % adding all a-conditions to the new case newca
   if (length aco)=1 then newca:=cons(car acond,newca)
		     else <<

    ali:=for each i in aco collect car pickfac(ex,i);
    avars:=my_smemberl(ali,vl);
    if length avars neq n then <<
     % The progress in this case is that all new equations will
     % have less variables.
     % Now determination on which variables the constants of back-
     % integration would depend on. This is the intersection of all
     % variables avars in the a-condition and the variables bvars in
     % the b-condition in which the constant of integration occurs. The
     % a-condition is known, it will be made from acond. The relevant
     % b-condition is determined through the current index of aco
     % from which the coefficient is to be determined (which is not the
     % first index of aco)

     aco:=cdr aco;    % a1 already in acond
     while aco do <<
      ai:=car aco;aco:=cdr aco;

      % find the list of variables bvers of the b-condition bco
      % which includes the b-factor corresponding to ai=car aco
      % disadvantage of this way: if bco has m elements then all
      % variables of bco are determined m-1 times.

      % determining bco as the b-condition which contains ai
      bco:=find_cond(bcopy,ai);
      % bcopy is used instead of bcons
      % the condition with ai may already have been dropped from
      % bcons because of ai depending on all variables, i.e. the
      % new functions in the b-conditions would depend on all
      % variables and be no help.

      if pri then <<write"bco=",bco;terpri()>>$
      % bcoli:=smemberl(ili,bco); % in case bco has already had subst.

      % determining bvars
      bvars:=nil;
      for each bi in bco do
      bvars:=union(my_smemberl(cdr pickfac(ex,bi),vl),bvars);
      if pri then <<write"bvars=",bvars$terpri()>>$

      % introducing new constants of integration
      ci:=newfct(fname_,intersection(avars,bvars),nfct_)$
      cilist:=cons(ci,cilist);
      nfct_:=nfct_+1;

      acond:=cons(list('MINUS,list('TIMES,ci,car pickfac(ex,ai))),acond);
      if pri then <<write"acond=",acond;terpri()>>$

      if bco:=find_cond(bcons,ai) then <<
       bcons:=subst(subst(list('TIMES,ci,a1),a1,bco),bco,bcons);
       if pri then <<write"bcons=",bcons;terpri()>>$
      >>

     >>;
     acond:=cons('PLUS,acond)

    >>
	       else << % The condition aco is a *-condition
     % progress in this case is that new a-conditions
     % have less functions
     addex:=t;
     ali:=reverse ali;
     aco:=reverse aco;
     while length ali > 1 do <<
      if pri then <<write"ali1=",ali;terpri()>>$
      % Generate the a-condition
      if pri then <<write"###";terpri()>>$

      ali:=
      if not zerop car ali then
      for each i in cdr ali collect
      reval list('DF,list('QUOTIENT,i,car ali),x)
			   else cdr ali;
      if pri then <<write"ali2=",ali;terpri()>>$

      % Drop that element from bcons which includes
      % car ali (as first element)
      if bco:=find_cond(bcons,car aco) then 
      bcons:=setdiff(bcons,list bco);
      aco:=cdr aco

     >>;
     acond:=car ali;
     if (pairp acond) and (car acond = 'QUOTIENT) then
     acond:=cadr acond;

    >>;

    newca:=cons(acond,newca)

   >>;
   if pri then <<write"newca1=",newca;terpri()>>;
  >>;  % all a-conditions have been dealt with
  if pri then <<write"newca2=",newca;terpri()>>;

  % completing all b-conditions
  for each bi in ili do bcons:=subst(cdr pickfac(ex,bi),bi,bcons);

  % adding all b-conditions to the new case newca
  while bcons do <<
   if length car bcons = 1 then newca:=cons(caar bcons,newca)
			   else newca:=cons(cons('PLUS,car bcons),
						 newca);
   bcons:=cdr bcons
  >>;
  % if ex is an *-expression with grade>1 then possibly b-conditions
  % had to be dropped, so ex must be added
  if addex then newca:=cons(exx,newca);

  if pri then <<write"newca3=",newca;terpri()>>;

  % adding the list of constants of integeration
  newca:=cons(cilist,newca);
  if pri then <<write"cilist=",cilist;terpri()>>;

  newcases:=cons(newca,newcases)
 >>;
 return newcases
end$ % of starsep

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

symbolic procedure separizable(p,ftem,vl)$
begin scalar x,ft,f,ex,v,a,b,vlcp,allvarcaara,print_bak$
  vlcp:=vl;
  repeat <<
    x:=car vl;  vl:=cdr vl;
     % Determining all x-dependent functions ft
    ft:=nil;
    for each f in ftem do
    if member(x,fctargs f) and
       not my_freeof(p,f) then ft:=cons(f,ft)$
    f:=car reverse fctsort ft$     % sorting w.r.t. number of args
    v:=car setdiff(vlcp,fctargs f)$  % getting rid of f by v-differen.

    % preparation of the separ-call, ft are now v-indep. functions
    ft:=nil$
    for each f in ftem do if my_freeof(f,v) then ft:=cons(f,ft)$
 %  ex:=separ(p,ft,list v,nil)$ % det. all lin. ind. factors
    print_bak:=print_;  print_:=nil;
    ex:=separ2(p,ft,list v)$ % det. all lin. ind. factors
    print_:=print_bak;
    a:=ex;

    while a and <<
      b:=vlcp;
      while b and not my_freeof(caar a,car b) do b:=cdr b;
      b
    >> do a:=cdr a;
    if a then allvarcaara:=cons(caar a,allvarcaara);

  >> until (null a) or (null vl);
  % if a then null vl then whatever x was, there is allways one
  % element (car a) of ex such that car of this element (caar a)
  % does depend on all variables --> no separability possible,
  % new functions would depend on all variables
  
  % if a then test whether separation would be possible by getting
  % rid of functions through differentiation 
  % (this would not be the case if e.g. sin(all variables) would occur)
  % --> use of smemberl

  vl:=vlcp;
  while allvarcaara and not not_included(vlcp,smemberl(vlcp,car allvarcaara)) do
  <<allvarcaara:=cdr allvarcaara; vl:=cdr vl>>$

  return if a and null allvarcaara then nil          % no chance
         else if a then {nil,car allvarcaara,car vl} % deleting functions first
                   else <<                           % separation now possible
    if tr_gensep then
    <<terpri()$write "To separate directly wrt. ",x$
      write " the expression : "$deprint list p$
      write "will be differentiated wrt. ",v," to get rid of ",ft," ">>$
    {ex,v}
  >>
end$

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

symbolic procedure newgensep(p,starpro,ftem,vl)$
% ftem, vl should be correct:
% ftem:=smemberl(ftem_,p)$
% vl:=varslist(p,ftem,vl)$
% starpro:=stardep(ftem,vl)$
%    returns what starsep returns
begin 
scalar pl,v,ex,a,b$
%            ,gense,el1,el2,conds,newcali,l,clist$
%  if pairp p and (car p = 'QUOTIENT) then <<casecheck(caddr p,vl)$
%                                            p:=cadr p>>$
%  ftem:=smemberl(ftem,p)$
%  vl:=varslist(p,ftem,vl)$
%  if not (starpro:=stardep(ftem,vl)) then % then no *-equation
%  pl:=list list(nil,p)              % one case, no new functions
%			       else % e.g. starpro=((x y z).2)
%  if zerop cdr starpro then pl:=nil% ##############################
%                  %list cons(nil,separate(p,ftem,vl)) % direct sep
%		 else                 
%  if delength(p) leq gensep_ then   % generalized separation
%  <<
    if print_ then <<terpri()$write "generalized separation ">>$
    if tr_gensep then
    <<terpri()$write "de to be separated : "$
      deprint list p$
      terpri()$write "variable list for separation : ",car starpro$
      terpri()$write "for each of these variables ",cdr starpro,
		     " function(s) depend(s) on it.">>$

    for each v in car starpro do vl:=delete(v,vl);
    vl:=append(car starpro,vl);

    a:=separizable(p,ftem,vl)$
    if null a then return nil else 
    if null car a then return <<
      % functions to be deleted before separation are those in cadr a
      % ft:=smemberl(ftem,cadr a);
      if print_ then <<terpri()$
        write"In order to be separable with this procedure at first"$terpri()$
        write"one or more functions have to be eliminated through"$terpri()$
        write"differentiation and algebraic elimination, for example,"$terpri()$
        write"the functions: ",smemberl(ftem,cadr a)$terpri()$
      >>;
      nil
    >>            else <<ex:=car a;v:=cadr a>>$

    for each a in 
      reverse idx_sort for each b in ex collect cons(delength car b,b)
    collect cdr a$

    if tr_gensep then
    <<terpri()$write "Return from SEPAR: "$terpri()$prettyprint ex>>$

    % with v and v-dep. functions as first factors in the terms in ex
    pl:=starsep(p,ex,ftem,vl,v);
    if tr_gensep then
    <<terpri()$write "Return from STARSEP: "$terpri()$prettyprint pl>>$
%    print_:=oldpri$

%%############################################################
%    % l is a list of cases each (list of new fncts, cond1, cond2, ...)
%    % for each condition (neq p) in all cases calling gensep again
%    % if needed
%    pl:=nil;           % the final list of cases of only non-*-equ.
%    while l do         % checking all cases
%    <<clist:=caar l;   % list of new functions
%      conds:=cdar l;   % list of conditions
%      l:=cdr l;
%      newcali:=list list clist;
%      % newcali will be the list of new cases which starts as
%      % only one case and from this only the list of new functions
%      % but no conditions
%      while conds do % checking all conditions of one case
%      <<
%%        if car conds = ex then
%%        % ex aready investigated, not again
%%        % append ex to the conditions of all new cases
%%        newcali:=for each el1 in newcali collect
%%                 cons(car el1,cons(ex,cdr el1))
%%                          else <<
%          gense:=newgensep(car conds,append(ftem,clist),vl);
%          newcali:=for each el1 in gense join
%                   for each el2 in newcali collect
%                   cons(append(car el1,car el2),
%                        append(cdr el1,cdr el2));
%%        >>; 
%        conds:=cdr conds
%      >>;
%      pl:=append(newcali,pl)
%    >>
%  >>;
  return pl
end$ % of newgensep

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

symbolic procedure gen_separation2(arglist)$
% Indirect separation of a pde, 2nd version for non-linear PDEs
begin scalar p,h,fl,l,l1,pdes,forg,n,result,d,contrad,newpdes$%,newfdep,bak,sol
  pdes:=car arglist$
  forg:=cadr arglist$
  if expert_mode then <<
    l1:=selectpdes(pdes,1);
    if get(car l1,'starde) then flag(l1,'to_gensep)
  >>             else l1:=pdes$
  if (p:=get_gen_separ_pde(l1,nil,nil)) then
  if l:=newgensep(get(p,'val),
                  get(p,'starde),
                  get(p,'fcts),
                  get(p,'vars)) then
  if cdr l then <<
   if print_ then <<
    terpri()$
    write"The indirect separation leads to ",length l," cases."$
    %terpri()$
   >>$
   contrad:=t$
   n:=0;
   remflag1(p,'to_gensep)$
%   bak:=backup_pdes(pdes,forg)$ % must come before drop_pde(...
   backup_to_file(pdes,forg,nil)$
%   newfdep:=nil$
   while l do <<
    d:=car l; l:=cdr l;
    if not memberl(cdr d,ineq_) then << % non of the equations is an inequality
     if n neq 0 then <<

      h:=restore_and_merge(l1,pdes,forg)$
      pdes:= car h; 
      forg:=cadr h; % was not assigned above as it has not changed probably

%      h:=restore_pdes(bak);
%      pdes:=car h; forg:=cadr h
     >>;
     n:=n+1$
     level_:=cons(n,level_)$
     if print_ then <<
      print_level(t)$
      terpri()$
      write "CRACK is now called with the assumption : "$
      deprint(cdr d)
     >>$
     % formulation of new equations
     for each h in car d do ftem_:=fctinsert(h,ftem_);
     fl:=append(get(p,'fcts),car d);
     newpdes:=pdes$
     for each h in cdr d do 
     newpdes:=eqinsert(mkeq(h,fl,vl_,allflags_,t,list(0),nil,newpdes),newpdes);
     % further necessary step to call crackmain():
     recycle_fcts:=nil$  % such that functions generated in the sub-call 
                         % will not clash with existing functions
     l1:=crackmain(newpdes,forg)$
%     for each sol in l1 do
%     if sol then <<
%      for each f in caddr sol do 
%      if h:=assoc(f,depl!*) then newfdep:=cons(h,newfdep);
%     >>;
     if not contradiction_ then contrad:=nil$ 
     if l1 and not contradiction_ then
        result:=union(l1,result);
     contradiction_:=nil$                   

    >>
   >>;
   delete_backup()$
%   % newfdep are additional dependencies of the new functions in l1
%   depl!*:=append(depl!*,newfdep);

   contradiction_:=contrad$ 
   if contradiction_ then result:=nil$
   if print_ then <<
    terpri()$
    write"This completes the investigation of all cases of an ",
         "indirect separation."$
    terpri()$
   >>$
   result:=list result % to tell crackmain that computation is completed
  >>       else <<     % only one case
   l:=car l;
   for each h in car l do ftem_:=fctinsert(h,ftem_);
   fl:=append(get(p,'fcts),car l);
   pdes:=drop_pde(p,pdes,nil)$
   for each h in cdr l do 
   pdes:=eqinsert(mkeq(h,fl,vl_,allflags_,t,list(0),nil,pdes),pdes);
   result:=list(pdes,forg)
  >>$
  return result$
end$

endmodule$

end$