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

%********************************************************************
module utilities$
%********************************************************************
%  Routines for finding leading derivatives and others
%  Author: Andreas Brand 1990 1994
%          Thomas Wolf since 1994

%%%%%%%%%%%%%%%%%%%%%%%%%
%  properties of pde's  %
%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure drop_dec_with(de1,de2,rl)$
% drop de1 from the 'dec_with or 'dec_with_rl list of de2
% currently for all orderings
begin scalar a,b,c$
  a:=if rl then get(de2,'dec_with_rl)
           else get(de2,'dec_with)$
  for each b in a do << % for each ordering b
    b:=delete(de1,b);
    if length b>1 then c:=cons(b,c);
  >>; 
  if rl then put(de2,'dec_with_rl,c)
        else put(de2,'dec_with   ,c)
end$

symbolic procedure add_dec_with(ordering,de1,de2,rl)$
% add (ordering de1) to 'dec_with or 'dec_with_rl of de2 
begin scalar a,b$
  a:=if rl then get(de2,'dec_with_rl)
           else get(de2,'dec_with)$
  b:=assoc(ordering,a)$
  a:=delete(b,a)$
  if b then b:=cons(ordering,cons(de1,cdr b))
       else b:=list(ordering,de1)$
  if rl then put(de2,'dec_with_rl,cons(b,a))
        else put(de2,'dec_with   ,cons(b,a))$
end$

symbolic procedure add_both_dec_with(ordering,de1,de2,rl)$
% add (ordering de1) to 'dec_with or 'dec_with_rl of de2  and 
% add (ordering de2) to 'dec_with or 'dec_with_rl of de1 
begin
  add_dec_with(ordering,de1,de2,rl)$
  add_dec_with(ordering,de2,de1,rl)$
end$

symbolic procedure drop_rl_with(de1,de2)$
% drop de1 from the 'rl_with list of de2
put(de2,'rl_with,delete(de1,get(de2,'rl_with)))$

symbolic procedure add_rl_with(de1,de2)$
% add de1 to 'rl_with of de2 and vice versa
<<put(de2,'rl_with,cons(de1,get(de2,'rl_with)))$
  put(de1,'rl_with,cons(de2,get(de1,'rl_with)))>>$

symbolic procedure prevent_simp(v,de1,de2)$
% it is df(de1,v) = de2
% add dec_with such that de2
% will not be simplified to 0=0
begin scalar a,b$
  % a:=get(de1,'fcts)$
  a:=list(0);   % all orderings for which de1 is used (-->ord)
  for each b in a do if member(v,fctargs(b)) then 
  <<add_dec_with(b,de2,de1,nil);add_dec_with(b,de2,de1,t)>>;
  % a:=get(de2,'fcts)$
  a:=list(0);   % all orderings for which de2 is used (-->ord)
  for each b in a do if member(v,fctargs(b)) then 
  <<add_dec_with(b,de1,de2,nil);add_dec_with(b,de1,de2,t)>>;
end$

symbolic procedure termread$
begin scalar val, !*echo;  % Don't re-echo tty input
 if not null old_history then <<
  val:=car old_history$
  if print_ then <<write"old input: ",val$terpri()>>$
  old_history:=cdr old_history
 >>                      else <<
  rds nil; wrs nil$         % Switch I/O to terminal
  val := read()$
  if ifl!* then rds cadr ifl!*$  %  Resets I/O streams
  if ofl!* then wrs cdr ofl!*$
 >>$
 history_:=cons(val,history_)$
 return val
end$

symbolic procedure termxread$
begin scalar val, !*echo;  % Don't re-echo tty input
 if not null old_history then <<
  val:=car old_history$
  if print_ then <<write"old input: ",val$terpri()>>$
  old_history:=cdr old_history
 >>                      else <<
  rds nil; wrs nil$         % Switch I/O to terminal
  val := xread(nil)$
  if ifl!* then rds cadr ifl!*$  %  Resets I/O streams
  if ofl!* then wrs cdr ofl!*$
 >>$
% history_:=cons(compress(append(explode val,list('$))),history_)$
 history_:=cons(val,history_)$
 return val
end$

symbolic procedure termlistread()$
begin scalar l;
  l:=termxread()$
  if (not null l) and
     ((atom l) or
      (pairp l and (car l neq '!*comma!*)))
  then l:=list('!*comma!*,l);
  if l and ((not pairp l) or (car l neq '!*comma!*)) then
  <<terpri()$write"Error: not a legal list of elements.";terpri()$
    l:=nil>>
  else if pairp l then l:=cdr l; % dropping '!*comma!*
  return l
end$

symbolic procedure mkeqlist(vallist,ftem,vl,flaglist,simp_flag,orderl,pdes)$
%  make a list of equations
%    vallist:  list of expressions
%    ftem:     list of functions
%    vl:       list of variables
%    flaglist: list of flags
%    orderl:   list of orderings where the equations are valid
%    pdes:     list of all equations by name to update inequalities
%              within update()
begin scalar l1$
 for each a in vallist do
     l1:=eqinsert(mkeq(a,ftem,vl,flaglist,simp_flag,orderl,
                       nil,append(l1,pdes)),l1)$
return l1 
end$

symbolic procedure mkeq(val,ftem,vl,flaglist,simp_flag,orderl,hist,pdes)$
%  make a new equation
%    val:      expression
%    ftem:     list of functions
%    vl:       list of variables
%    flaglist: list of flags
%    orderl:   list of orderings where the equation is valid
%    hist:     the history of val
%    pdes:     list of all equations by name to update inequalities
%              within update()
%  If the new equation to be made is only to exist temporarily then
%  call mkeq with pdes=nil to avoid lasting effects of the temprary pde.
% 
begin scalar s$
 s:=new_pde()$
 if record_hist and hist then put(s,'histry_,reval hist)$
 for each a in flaglist do flag1(s,a)$
 if not update(s,val,ftem,vl,simp_flag,orderl,pdes) then 
   <<drop_pde(s,nil,nil)$
   s:=nil>>$
 if record_hist and null hist and s then put(s,'histry_,s)$
 return s
end$

symbolic procedure no_of_derivs(equ)$
begin scalar h,dl;
 h:=0;
 dl:=get(equ,'derivs);
 while dl do <<
  if (pairp caar dl) and (cdaar dl) then h:=add1 h;
  dl:=cdr dl
 >>;
 return h
end$

symbolic procedure update(equ,val,ftem,vl,simp_flag,orderl,pdes)$
%  update the properties of a pde
%    equ:      pde
%    val:      expression
%    ftem:     list of functions
%    vl:       list of variables
%    orderl:   list of orderings where the equation is valid
%    pdes:     to call addineq at end
%  *** important ***:
%  call afterwards also drop_pde_from_idties(p,pdes,pval) and
%                       drop_pde_from_properties() 
%  if this is now a new equation
begin scalar l$
  if val and not zerop val then <<
    %ftem:=reverse union(smemberl(ftem,val),nil)$
    ftem:=sort_according_to(smemberl(ftem,val),ftem_)$
    put(equ,'terms,no_of_terms(val))$    
    if simp_flag then <<
      % the following test to avoid factorizing big equations
      val:=% if get(equ,'terms)>max_factor then simplifypde(val,ftem,nil,equ)
           %                               else 
           simplifypde(val,ftem,t,equ)$  
      if val and not zerop val then <<
	ftem:=reverse union(smemberl(ftem,val),nil)$
        put(equ,'terms,no_of_terms(val))$    
      >>
    >>$
  >>$
  depl!*:=delete(assoc(reval equ,depl!*),depl!*)$
  if val and not zerop val then <<
    put(equ,'val,val)$
    put(equ,'fcts,ftem)$
    for each v in vl do
      if not my_freeof(val,v) then l:=cons(v,l)$
    vl:=sort_according_to(l,vl_);
    put(equ,'vars,vl)$  
    if vl then
    depl!*:=cons(cons(equ,vl),depl!*)$ % needed in expressions in idnties_
    put(equ,'nvars,length vl)$
    put(equ,'level,level_)$
    put(equ,'derivs,sort_derivs(all_deriv_search(val,ftem),ftem,vl))$
    if struc_eqn then put(equ,'no_derivs,no_of_derivs(equ));
    put(equ,'fcteval_lin,nil)$
    put(equ,'fcteval_nca,nil)$
    put(equ,'fcteval_nli,nil)$   
    put(equ,'fct_nli_lin,nil)$   
    put(equ,'fct_nli_nca,nil)$   
    put(equ,'fct_nli_nli,nil)$   
    put(equ,'fct_nli_nus,nil)$   
%    put(equ,'terms,no_of_terms(val))$    
    put(equ,'length,pdeweight(val,ftem))$    
    put(equ,'printlength,delength val)$
    put(equ,'rational,nil)$
    put(equ,'nonrational,nil)$
    put(equ,'allvarfcts,nil)$
    put(equ,'orderings,orderl)$	% Orderings !
    for each f in reverse ftem do
      if rationalp(val,f) then 
	 <<put(equ,'rational,cons(f,get(equ,'rational)))$
	   if fctlength f=get(equ,'nvars) then
	      put(equ,'allvarfcts,cons(f,get(equ,'allvarfcts)))>>
      else put(equ,'nonrational,cons(f,get(equ,'nonrational)))$
%    put(equ,'degrees,          % too expensive
%     if linear_pr then cons(1,for each l in get(equ,'rational) 
%                              collect (l . 1))
%                  else fct_degrees(val,get(equ,'rational))    )$
    put(equ,'partitioned,nil)$
    put(equ,'starde,stardep(ftem,vl))$
    flag1(equ,'to_eval)$
    if (l:=get(equ,'starde)) then
       <<%remflag1(equ,'to_eval)$
       remflag1(equ,'to_int)$
       remflag1(equ,'to_fullint)$
       if simp_flag and (zerop cdr l) then flag1(equ,'to_sep)$
       % remflag1(equ,'to_diff)
       >>
    else remflag1(equ,'to_gensep)$
    if get(equ,'starde) and
       (zerop cdr get(equ,'starde)      ) % or (get(equ,'length)<=gensep_)) 
       then 
       else remflag1(equ,'to_sep)$
    if get(equ,'nonrational) then 
       <<%remflag1(equ,'to_decoup)$
       if not freeoflist(get(equ,'allvarfcts),get(equ,'nonrational)) then
	  remflag1(equ,'to_eval)>>$
%    if not get(equ,'allvarfcts) then remflag1(equ,'to_eval)$
    if (not get(equ,'rational)) or
       ((l:=get(equ,'starde)) and (cdr l = 0)) then remflag1(equ,'to_eval)$
    if homogen_ then <<
      l:=cdr algebraic find_hom_deg(lisp val);
      put(equ,'hom_deg,l)$
%      if car l=1 then << % i.e. linear in flin_
%        l:=get(equ,'derivs);  
%        while l and (null linf or (length linf < 3)) do <<
%          if not freeoflist(car l,flin_) then <<
%            linf:=cons(car l,linf);
%            if member(car l,ineq_) then fd1:=car l
%          >>;
%          l:=cdr l
%        >>;
%        if linf and (length linf = 2) and fd1 then <<<<
%          if NON-ZERO(coeffn(get(equ,'val),fd1,1)) then <<
%            fd2:=car delete(fd1,linf);
%  braucht pdes, was nicht vorhanden ist
%            addineq(pdes,fd2);
%            addineq(pdes,coeffn(get(equ,'val),fd2,1))
%          >>
%        >>
%      >>      
    >>$
    put(equ,'split_test,nil)$
    put(equ,'linear_,if lin_problem then t
                                    else lin_check(val,ftem))$
    new_ineq_from_pde(equ,pdes);
    return equ
  >>
end$

symbolic procedure new_ineq_from_pde(equ,pdes)$
% currently only effective for equations with 2 terms
% If one term of the equation is non-zero then the sum of the 
% remaining terms has to be non-zero too
if pdes and null lin_problem and (get(equ,'terms)=2) % >1) 
then begin scalar valu;
 valu:=get(equ,'val);
 if not (pairp valu and (car valu='PLUS)) then valu:=reval valu; 
 if not (pairp valu and (car valu='PLUS)) then write"err in update"
                                          else 
%    for each l in cdr valu do
%    if null may_vanish l then addineq(pdes,reval{'DIFFERENCE,valu,l})
 if null may_vanish cadr  valu then addineq(pdes,caddr valu) else 
 if null may_vanish caddr valu then addineq(pdes,cadr  valu)
end$

symbolic procedure fct_degrees(pv,ftem)$
% ftem are to be the rational functions
begin
 scalar f,l,ll,h,degs$
 if den pv then pv:=num pv$
 for each f in ftem do <<
  l:=gensym()$
  ll:=cons((f . l),ll)$
  pv:=subst({'TIMES,l,f},f,pv)$
 >>$
 pv:=reval pv$
 for each l in ll do <<
  degs:=cons((car l . deg(pv,cdr l)),degs)$ 
 >>;
 h:=cdar ll$
 for each l in cdr ll do pv:=subst(h,cdr l,pv)$
 pv:=reval pv$
 return cons(deg(pv,h),degs)
end$

symbolic procedure pde_degree(pv,ftem)$
% ftem are to be the rational functions
begin
 scalar f,h$
 if den pv neq 1 then pv:=num pv$
 h:=gensym()$
 for each f in ftem do pv:=subst({'TIMES,h,f},f,pv)$
 pv:=reval pv$
 return deg(pv,h)
end$

symbolic procedure dfsubst_update(f,der,equ)$
%  miniml update of some properties of a pde
%    equ:      pde
%    der:      derivative
%    f:        f new function
begin scalar l$
  for each d in get(equ,'derivs) do
    if not member(cadr der,car d) then l:=cons(d,l)
    else
      <<l:=cons(cons(cons(f,df_int(cdar d,cddr der)),cdr d),l)$
      put(equ,'val,
          subst(reval cons('DF,caar l),reval cons('DF,car d),
                get(equ,'val)))>>$
  put(equ,'fcts,subst(f,cadr der,get(equ,'fcts)))$
  put(equ,'allvarfcts,subst(f,cadr der,get(equ,'allvarfcts)))$
  if get(equ,'allvarfcts) then flag(list equ,'to_eval)$
%  This would reactivate equations which resulted due to
%  substitution of derivative by a function.
%  8.March 98: change again: the line 3 lines above has been reactivated
  put(equ,'rational,subst(f,cadr der,get(equ,'rational)))$
  put(equ,'nonrational,subst(f,cadr der,get(equ,'nonrational)))$
  put(equ,'derivs,sort_derivs(l,get(equ,'fcts),get(equ,'vars)))$
  return equ
end$

symbolic procedure eqinsert(s,l)$
% l is a sorted list
if not (s or get(s,'val)) or zerop get(s,'length) or member(s,l) then l
else if not l then list s
else begin scalar l1,n$
 l1:=proddel(s,l)$
 if car l1 then
  <<n:=get(s,'length)$
    l:=cadr l1$
    l1:=nil$
    while l and (n>get(car l,'length)) do
     <<l1:=cons(car l,l1)$
       l:=cdr l>>$
    l1:=append(reverse l1,cons(s,l))$
  >>
 else if l1 then l1:=cadr l1  % or reverse of it
            else l1:=l$
 return l1$
end$

symbolic procedure not_included(a,b)$
% meaning: not_all_a_in_b = setdiff(a,b)
% Are all elements of a also in b? If yes then return nil else t
% This could be done with setdiff(a,b), only setdiff
% copies expressions and needs extra memory whereas here we only
% want to know one bit (included or not)
begin scalar c$
 c:=t;
 while a and c do <<
  c:=b;
  while c and ((car a) neq (car c)) do c:=cdr c;
  % if c=nil then car a is not in b
  a:=cdr a;
 >>;
 return if c then nil
             else t
end$

symbolic procedure proddel(s,l)$
% delete all pdes from l with s as factor
% delete s if it is a consequence of any known pde from l
begin scalar l1,l2,l3,n,lnew,pdes,s_hist$
 if pairp(lnew:=get(s,'val)) and (car lnew='TIMES) then lnew:=cdr lnew
                                                   else lnew:=list lnew$
 n:=length lnew$
 pdes:=l$
 while l do << 
  if pairp(l1:=get(car l,'val)) and (car l1='TIMES) then l1:=cdr  l1
                                                    else l1:=list l1$
  if n<length l1 then     % s has less factors than car l
    if not_included(lnew,l1) then 
    l2:=cons(car l,l2)    % car l is not a consequ. of s
                             else 
    <<l3:=cons(car l,l3); % car l is a consequ. of s
      drop_pde(car l,nil,{'QUOTIENT,{'TIMES,s,get(car l,'val)},get(s,'val)})
    >>
  else <<
    if null not_included(l1,lnew) then % s is a consequence of car l
    <<if print_ then <<terpri()$write s," is a consequence of ",car l,".">>$
      % one could stop here but continuation can still be useful
      if null s_hist then s_hist:={'QUOTIENT,
                                   {'TIMES,car l,get(s,'val)},
                                   get(car l,'val)            }$
    >>$
    % else 
    if null l3 or (car l3 neq car l) then l2:=cons(car l,l2)$
  >>;
  l:=cdr l
 >>$
 if print_ and l3 then <<
  listprint l3$
  if cdr l3 then write " are consequences of ",s
            else write " is a consequence of ",s;
  terpri()$
 >>$
 if s_hist then <<drop_pde(s,nil,s_hist);s:=nil>>$
 return list(s,reverse l2)$
end$


symbolic procedure myprin2l(l,trenn)$
if l then 
   <<if pairp l then
        while l do
          <<write car l$
          l:=cdr l$
          if l then write trenn>> 
   else write l>>$

symbolic procedure print_stars(s)$
begin scalar b,star,pv$
 pv:=get(s,'val)$
 if (pairp pv) and (car pv='TIMES) then pv:=t else pv:=nil$
 star:=get(s,'starde)$
 if star or pv then <<
  write "("$
  if pv then write"#"$
  if star then for b:=1:(1+cdr star) do write"*"$
  write")"$
 >>$
end$

symbolic procedure typeeq(s)$
%  print equation 
 if (null print_) or (get(s,'printlength)>print_) then begin scalar a,b$
   print_stars(s);
   write " ",(a:=get(s,'terms))," terms"$
   if a neq (b:=get(s,'length)) then write", ",b," factors"$
   write", with derivatives"$
   if get(s,'starde) then <<
    write": "$ terpri()$
    print_derivs(s,nil)$
   >>   else <<
    write" of functions of all variables: "$ terpri()$
    print_derivs(s,t)$
   >>
 end                     else
 mathprint list('EQUAL,0,get(s,'val))$

symbolic procedure print_derivs(p,allvarf)$
begin scalar a,d,dl,avf;
 dl:=get(p,'derivs)$
 if allvarf then <<
  avf:=get(p,'allvarfcts);
  for each d in dl do
  if not freeoflist(d,avf) then a:=cons(d,a);
  dl:=reverse a
 >>$
 dl:=for each d in dl collect <<
  a:=if null cdar d then caar d
                    else cons('DF,car d);
  if cdr d=1 then a else {'EXPT,a,cdr d}
 >>$
 mathprint cons('! ,dl)
end$

symbolic procedure typeeqlist(l)$
%  print equations and their property lists
<<terpri()$
for each s in l do
 <<terpri()$
 write s," : "$ 
 if not print_all then typeeq(s)
                  else
 if (null print_) or (get(s,'printlength)>print_) then 
 <<write get(s,'terms)," terms"$terpri()>>        else
 mathprint list('EQUAL,0,get(s,'val))$
 if print_all then
    <<        write "   derivs     : "$    
    terpri()$print_derivs(s,nil)$
%   if struc_eqn then <<
%    terpri()$write "   no_derivs  : ",get(s,'no_derivs)$
%   >>$
             write "   terms      : ",get(s,'terms)$
    terpri()$write "   fcts       : ",get(s,'fcts)$
    terpri()$write "   vars       : ",get(s,'vars)$
    terpri()$write "   nvars      : ",get(s,'nvars)$
    terpri()$write "   length     : ",get(s,'length)$
    terpri()$write "   printlength: ",get(s,'printlength)$
    terpri()$write "   level      : ",get(s,'level)$
    terpri()$write "   allvarfcts : ",get(s,'allvarfcts)$
    terpri()$write "   rational   : ",get(s,'rational)$ 
    terpri()$write "   nonrational: ",get(s,'nonrational)$ 
    terpri()$write "   degrees    : ",get(s,'degrees)$
    terpri()$write "   starde     : ",get(s,'starde)$ 
    terpri()$write "   fcteval_lin: ",get(s,'fcteval_lin)$
    terpri()$write "   fcteval_nca: ",get(s,'fcteval_nca)$
    terpri()$write "   fcteval_nli: ",get(s,'fcteval_nli)$
    terpri()$write "   fct_nli_lin: ",get(s,'fct_nli_lin)$
    terpri()$write "   fct_nli_nca: ",get(s,'fct_nli_nca)$
    terpri()$write "   fct_nli_nli: ",get(s,'fct_nli_nli)$
    terpri()$write "   fct_nli_nus: ",get(s,'fct_nli_nus)$
    terpri()$write "   rl_with    : ",get(s,'rl_with)$ 
    terpri()$write "   dec_with   : ",get(s,'dec_with)$ 
    terpri()$write "   dec_with_rl: ",get(s,'dec_with_rl)$
%    terpri()$write "   dec_info   : ",get(s,'dec_info)$ 
    terpri()$write "   to_int     : ",flagp(s,'to_int)$ 
    terpri()$write "   to_fullint : ",flagp(s,'to_fullint)$ 
    terpri()$write "   to_sep     : ",flagp(s,'to_sep)$ 
    terpri()$write "   to_gensep  : ",flagp(s,'to_gensep)$ 
    terpri()$write "   to_decoup  : ",flagp(s,'to_decoup)$
    terpri()$write "   to_drop    : ",flagp(s,'to_drop)$
    terpri()$write "   to_eval    : ",flagp(s,'to_eval)$
    terpri()$write "   to_diff    : ",flagp(s,'to_diff)$
    terpri()$write "   to_under   : ",flagp(s,'to_under)$
    terpri()$write "   to_symbol  : ",flagp(s,'to_symbol)$
    terpri()$write "   not_to_eval: ",get(s,'not_to_eval)$
    terpri()$write "   used_      : ",flagp(s,'used_)$
    terpri()$write "   orderings  : ",get(s,'orderings)$
    terpri()$write "   histry_    : ",get(s,'histry_)$
    terpri()$write "   partitioned: ",if get(s,'partitioned) then 
                                      "not nil"              else
                                      "nil"$
    terpri()$write "   split_test : ",get(s,'split_test)$
    terpri()$write "   linear_    : ",get(s,'linear_)$
    if homogen_ then <<
      terpri()$write "   hom_deg    : ",get(s,'hom_deg)
    >>$
    terpri()>>
 >> >>$

symbolic procedure rationalp(p,f)$
%  tests if p is rational in f and its derivatives
not pairp p
or
((car p='QUOTIENT) and 
 polyp(cadr p,f) and polyp(caddr p,f))
or 
((car p='EQUAL) and
 rationalp(cadr p,f) and rationalp(caddr p,f)) 
or 
polyp(p,f)$

symbolic procedure ratexp(p,ftem)$
if null ftem then t
             else if rationalp(p,car ftem) then ratexp(p,cdr ftem)
                                           else nil$

symbolic procedure polyp(p,f)$
%  tests if p is a polynomial in f and its derivatives
%    p: expression
%    f: function
if my_freeof(p,f) then t
else
begin scalar a$
if atom p then a:=t
else
if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,'QUOTIENT,'DF)) then
                                        %  erlaubte Funktionen
        <<if (car p='PLUS) or (car p='TIMES) then
                <<p:=cdr p$
                while p do
                    if a:=polyp(car p,f) then p:=cdr p
                                         else p:=nil>>
        else if (car p='MINUS) then
                a:=polyp(cadr p,f)
        else if (car p='QUOTIENT) then
                <<if freeof(caddr p,f) then a:=polyp(cadr p,f)>>
        else if car p='EXPT then        %  Exponent
                <<if (fixp caddr p) then
                   if caddr p>0 then a:=polyp(cadr p,f)>>
        else if car p='DF then          %  Ableitung
                if (cadr p=f) or freeof(cadr p,f) then a:=t>>
else a:=(p=f)$
return a
end$

symbolic procedure starp(ft,n)$
%  yields T if all functions from ft have less than n arguments 
begin scalar b$
  while not b and ft do                % searching a fct of all vars
  if fctlength car ft=n then b:=t
			else ft:=cdr ft$
  return not b
end$

symbolic procedure stardep(ftem,vl)$
%  yields: nil, if a function (from ftem) in p depends
%               on all variables (from vl)
%          cons(v,n) otherwise, with v being the list of variables
%               which occur in a minimal number of n functions
begin scalar b,v,n$
  if starp(ftem,length vl) then
  <<n:=sub1 length ftem$
    while vl do                % searching var.s on which depend a
			       % minimal number of functions
    <<if n> (b:=for each h in ftem sum
				   if member(car vl,fctargs h) then 1
							       else 0)
      then <<n:=b$v:=list car vl>>       % a new minimum
      else if b=n then v:=cons(car vl,v)$
      vl:=cdr vl>> >>$
  return if v then cons(v,n)             % on each varible from v depend n
		  		         % functions
	      else nil
end$

%symbolic procedure no_of_sep_var(ftem)$
%% assuming ftem are all functions from an ise
%% How many are there indirectly separable variables?
%% If just two then the new indirect separation is possible
%begin scalar v,vs$
%  vl:=argset(ftem);
%  for each f in ftem do 
%  vs:=union(setdiff(vl,fctargs f),vs)$
%  return vs
%end$

symbolic operator parti_fn$
symbolic procedure parti_fn(fl,el)$
% fl ... alg. list of functions, el ... alg. list of equations
% partitions fl such that all functions that are somehow dependent on
% each other through equations in el are grouped in lists, 
% returns alg. list of these lists
begin
 scalar f1,f2,f3,f4,f5,e1,e2;
 fl:=cdr fl;
 el:=cdr el;
 while fl do <<
  f1:=nil;         % f1 is the sublist of functions depending on each other
  f2:=list car fl; % f2 ... func.s to be added to f1, not yet checked
  fl:=cdr fl;
  while f2 and fl do <<
   f3:=car f2; f2:=cdr f2;
   f1:=cons(f3,f1);
   for each f4 in  
   % smemberl will be all functions not registered yet that occur in 
   % an equation in which the function f3 occurs
   smemberl(fl,    % fl ... the remaining functions not known yet to depend
            <<e1:=nil;  % equations in which f3 occurs
              for each e2 in el do 
              if smember(f3,e2) then e1:=cons(e2,e1);
              e1
            >>
           )        do <<
    f2:=cons(f4,f2);
    fl:=delete(f4,fl)
   >>
  >>;
  if f2 then f1:=append(f1,f2);
  f5:=cons(cons('LIST,f1),f5)
 >>;
 return cons('LIST,f5)
end$

symbolic procedure plot_dependencies(pdes)$
begin scalar ps,fl$
  ps:=promptstring!*$ promptstring!*:=""$ 
  fl:=ftem_;
  if flin_ and yesp 
  "Shall only functions from the linear list flin_ be considered? " 
  then fl:=setdiff(fl,setdiff(fl,flin_))$
  promptstring!*:=ps$
  plot_dep_matrix(pdes,fl)
end$

symbolic procedure plot_dep_matrix(pdes,allf)$
begin scalar f,ml,lf,fl,h,lh,lc,n,m,h;
  terpri()$
  write "Horizontally: function names (each vertical),  ",
        "Vertically: equation indices"$
  terpri()$

  ml:=0;                % the maximal length of all variable names
  lf:=length allf$
  for each f in reverse allf do <<
    h:=explode f; 
    lh:=length h;
    if lh>ml then ml:=lh;
    lc:=cons(h,lc);
  >>$

  % print the variable names
  for n:=1:ml do <<
    terpri()$ write"     "$
    for m:=1:lf do write <<
      h:=nth(lc,m);
      if n>length h then " "
                    else nth(nth(lc,m),n)
    >>
  >>$

  m:=add1 add1 ml;
  terpri()$terpri()$
  for each p in pdes do 
  if m>=0 then <<
   h:=explode p;
   for n:=3:length h do write nth(h,n);
   for n:=(sub1 length(h)):5 do write" "$
   fl:=get(p,'fcts);
   if (not get(p,'fcteval_lin)) and
      (not get(p,'fcteval_nca)) and
      (not get(p,'fcteval_nli)) then fcteval(p,nil)$ % for writing "s"
   for each f in allf do
   if freeof(fl,f) then write" " else 
   if solvable_case(p,f,'fcteval_lin) or 
      solvable_case(p,f,'fcteval_nca) then write"s"
                       	              else write"+"$
   terpri()$
   m:=add1 m$
   if m=23 then if not yesp "Continue ?" then m:=-1
                                         else m:=0
  >>$
end$

symbolic procedure solvable_case(p,f,case)$
begin scalar fe;
 fe:=get(p,case);
 while fe and (cdar fe neq f) do fe:=cdr fe$
 return fe
end$

%symbolic procedure lin_check(pde,fl)$
%<<for each f in fl do pde:=subst({'TIMES,lin_test_const,f},f,pde);
%  freeof(reval {'QUOTIENT,pde,lin_test_const},lin_test_const)
%>>$

symbolic procedure lin_check(pde,fl)$
% This needs pde to have prefix form.
begin scalar a,f;
 a:=pde;
 for each f in fl do a:=err_catch_sub(f,0,a);
 return
 if a then <<
  for each f in fl do pde:=subst({'TIMES,lin_test_const,f},f,pde);
  freeof(reval {'QUOTIENT,{'DIFFERENCE,pde,a},lin_test_const},lin_test_const)
 >>   else nil
end$

symbolic procedure plot_non0_coeff_ld(s)$
begin scalar dv,dl,dlc,dr,fdl,avf;
 write " Leading derivatives with non-zero symbol: "$terpri()$
 dv:=get(s,'derivs);
 avf:=get(s,'allvarfcts);
 while dv do <<
  dr:=caar dv; dv:=cdr dv;
  if member(car dr,avf) then <<
   dlc:=dl;
   while dlc and ((caar dlc neq car dr) or 
                  which_deriv(car dlc,dr)  ) do dlc:=cdr dlc;
   if null dlc then dl:=cons(dr,dl);
   % which_deriv(a,b) takes two lists of derivatives and returns how
   % often you need to diff. a in order to get at least the
   % derivatives in b. e.g. which_deriv((x 2 y), (x y 2)) returns y
  >>
 >>;
 for each dr in dl do << 
  dr:=if null cdr dr then car dr
                     else cons('DF,dr);
  if get(s,'linear_) or 
     freeofzero(reval {'DF,get(s,'val),dr},get(s,'fcts),
                get(s,'vars),get(s,'nonrational)) then 
  fdl:=cons(dr,fdl)
 >>;
 mathprint cons('! ,fdl)
end$

%%%%%%%%%%%%%%%%%%%%%%%%%
%  leading derivatives  %
%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure listrel(a,b,l)$
%   a>=b  w.r.t list l; e.g. l='(a b c) ->  a>=a, b>=c 
member(b,member(a,l))$

symbolic procedure abs_dfrel(p,q,vl)$
% returns t if derivative of p is lower than derivative of q
%         0         "             equal          "
%        nil        "             higher         "
% p,q  : derivatives or functions from ftem like ((f x 2 y z 3) . 2)
% ftem : list of fcts
% vl   : list of vars
begin scalar a$
return 
if lex_df then dfrel2(p,q,vl) else
if zerop (a:=absdeg(cdar p)-absdeg(cdar q)) then dfrel2(p,q,vl)   
else a<0$
end$

symbolic procedure mult_derivs(a,b)$
% multiplies deriv. of a and b 
% a,b list of derivs of the form ((fct var1 n1 ...).pow)
begin scalar l$
 return 
  if not b then a
  else if not a then b
  else
   <<
   for each s in a do 
      for each r in b do
        if car s=car r then l:=union(list cons(car r,plus(cdr r,cdr s)),l)
                       else l:=union(list(r,s),l)$
   l>>$
end$

symbolic procedure all_deriv_search(p,ftem)$
%  yields all derivatives occuring polynomially in a pde  p 
begin scalar a$
if not pairp p then 
 <<if member(p,ftem) then a:=list cons(list p,1)>>
else
<<if member(car p,'(PLUS QUOTIENT EQUAL)) then
    for each q in cdr p do
	  a:=union(all_deriv_search(q,ftem),a)
  else if car p='MINUS then a:=all_deriv_search(cadr p,ftem)
  else if car p='TIMES then
    for each q in cdr p do
	  a:=mult_derivs(all_deriv_search(q,ftem),a)
  else if (car p='EXPT) and numberp caddr p then 
      for each b in all_deriv_search(cadr p,ftem) do
          <<if numberp cdr b then 
               a:=cons(cons(car b,times(caddr p,cdr b)),a)>>
  else if (car p='DF) and member(cadr p,ftem) then a:=list cons(cdr p,1) 
>>$
return a
end$

symbolic procedure abs_ld_deriv(p)$
if get(p,'derivs) then reval cons('DF,caar get(p,'derivs))$

symbolic procedure abs_ld_deriv_pow(p)$
if get(p,'derivs) then cdar get(p,'derivs)
                  else 0$

symbolic procedure which_first(a,b,l)$
if null l then nil else
if a = car l then a else
if b = car l then b else which_first(a,b,cdr l)$


symbolic procedure total_less_dfrel(a,b,ftem,vl)$
% = 0 if a=b, =t if a<b, = nil if a>b
begin scalar fa,ad,al,bl$
  fa:=caar a$
  return
  if a=b then 0 else
  if lex_fc then % lex. order. of functions has highest priority
  if fa=caar b then
  if (ad:=abs_dfrel(a,b,vl))=0 then             % power counts
  if cdr a < cdr b then t
                   else nil
                               else
  if ad then t
        else nil
               else
  if fa=which_first(fa,caar b,ftem) then nil
                                    else t
            else % order. of deriv. has higher priority than fcts.
                 % number of variables of functions has still higher priority
  if (al:=fctlength fa) > (bl:=fctlength caar b) then nil
                                                 else
  if bl>al then t
           else
  if (ad:=abs_dfrel(a,b,vl))=0 then 
  if fa=caar b then
  if cdr a < cdr b then t
                   else nil
               else
  if fa=which_first(fa,caar b,ftem) then nil
                                    else t
                               else
  if ad then t
        else nil
end$

symbolic procedure sort_derivs(l,ftem,vl)$
% yields a sorted list of all derivatives in l
begin scalar l1,l2,a$
 return
 if null l then nil
           else <<
  a:=car l$
  l:=cdr l$
  while l do <<
     if total_less_dfrel(a,car l,ftem,vl) then l1:=cons(car l,l1)
                                          else l2:=cons(car l,l2)$
     l:=cdr l
  >>$
  append(sort_derivs(l1,ftem,vl),cons(a,sort_derivs(l2,ftem,vl)))>>
end$

symbolic procedure dfmax(p,q,vl)$
%   yields the higher derivative
%   vl list of variables e.g. p=((x 2 y 3 z).2), q=((x y 4 z).1)
%                             df(f,x,2,y,3,z)^2, df(f,x,y,4,z)
if dfrel(p,q,vl) then q
		 else p$

symbolic procedure dfrel(p,q,vl)$
%   the relation "p is lower than q"
%   vl list of vars e.g. p=((x 2 y 3 z).2), q=((x y 4 z).1)
if cdr p='infinity then nil
else if cdr q='infinity then t
else begin scalar a$
 return 
  if lex_df then dfrel1(p,q,vl) else
  if zerop(a:=absdeg(car p)-absdeg(car q)) then dfrel1(p,q,vl)
                                           else if a<0 then t
                                                       else nil
end$

symbolic procedure diffrelp(p,q,v)$
% gives  t  when p "<" q
%       nil when p ">" q
%        0  when p  =  q
%   p, q Paare (liste.power), v Liste der Variablen
%   liste Liste aus Var. und Ordn. der Ableit. in Diff.ausdr.,
%   power Potenz des Differentialausdrucks
if cdr p='infinity then nil
else if cdr q='infinity then t
     else dfrel1(p,q,v)$

symbolic procedure dfrel1(p,q,v)$
% p,q like ((f x 2 y z 3) . 2)
if null v then 
%   if cdr p = t then if cdr q = t then 0 else nil 
%                else if cdr q = t then t else
   if cdr p>cdr q then nil else         %  same derivatives,
   if cdr p<cdr q then t   else 0       %  considering powers
   % for termorderings of non-linear problems the last 2 lines 
   % have to be extended
else begin
	scalar a,b$
	a:=dfdeg(car p, car v)$
	b:=dfdeg(car q, car v)$
	return if a<b then t
	else   if b<a then nil
	else dfrel1(p,q,cdr v)  %  same derivative w.r.t car v
end$

symbolic procedure dfrel2(p,q,v)$
% p,q like ((f x 2 y z 3) . 2)
if null v then 0 
else begin
	scalar a,b$
	a:=dfdeg(car p, car v)$
	b:=dfdeg(car q,car v)$
	return if a<b then t
	else   if b<a then nil
	else dfrel2(p,q,cdr v)  %  same derivative w.r.t car v
end$

symbolic procedure absdeg(p)$
if null p then 0
else eval cons('PLUS,for each v in p collect if fixp(v) then sub1(v)
                                                        else 1)$

symbolic procedure maxderivs(numberlist,deriv,varlist)$
if null numberlist then 
 for each v in varlist collect dfdeg(deriv,v)
else begin scalar l$
 for each v in varlist do
  <<l:=cons(max(car numberlist,dfdeg(deriv,v)),l)$
  numberlist:=cdr numberlist>>$
 return reverse l
end$
   
symbolic procedure dfdeg(p,v)$
%   yields order of deriv. wrt. v$
%   e.g p='(x 2 y z 3), v='x --> 2
if null(p:=member(v,p)) then 0           
else if null(cdr p) or not fixp(cadr p) 
        then 1                          %  v without order
        else cadr p$                    %  v with order

symbolic procedure lower_deg(p,v)$
%  reduces the order of the derivative p wrt. v by one
%   e.g p='(x 2 y z 3), v='z --> p='(x 2 y z 2)
%   e.g p='(x 2 y z 3), v='y --> p='(x 2 z 3)
%  returns nil if no v-derivative
begin scalar newp$
 while p and (car p neq v) do <<newp:=cons(car p,newp);p:=cdr p>>$
 if p then 
 if null(cdr p) or not fixp(cadr p) then p:=cdr p
                                    else <<
  newp:=cons(sub1 cadr p,cons(car p,newp));
  p:=cddr p
 >> else newp:=nil;
 while p do <<newp:=cons(car p,newp);p:=cdr p>>$
 return reverse newp
end$

symbolic procedure df_int(d1,d2)$
begin scalar n,l$
return 
 if d1 then
  if d2 then
   <<n:=dfdeg(d1,car d1)-dfdeg(d2,car d1)$
   l:=df_int(if cdr d1 and numberp cadr d1 then cddr d1
                                           else cdr d1 ,d2)$
   if n<=0 then l
   else if n=1 then cons(car d1,l)
   else cons(car d1,cons(n,l))
   >>
  else d1$
end$

symbolic procedure linear_fct(p,f)$
begin scalar l$
 l:=ld_deriv(p,f)$
 return ((car l=f) and (cdr l=1))
end$

% not used anymore:
%
%symbolic procedure dec_ld_deriv(p,f,vl)$
%%  gets leading derivative of f in p wrt. vars order vl
%%  result: derivative , e.g. '(x 2 y 3 z)
%begin scalar l,d,ld$
% l:=get(p,'derivs)$
% vl:=intersection(vl,get(p,'vars))$
% while caaar l neq f do l:=cdr l$
% ld:=car l$l:=cdr l$
% % --> if lex_df then dfrel1() else
% d:=absdeg(cdar ld)$
% while l and (caaar l=f) and (d=absdeg cdaar l) do
%   <<if dfrel1(ld,car l,vl) then ld:=car l$
%   l:=cdr l>>$
% return cdar ld$
%end$

symbolic procedure ld_deriv(p,f)$
%  gets leading derivative of f in p
%  result: derivative + power , e.g. '((DF f x 2 y 3 z) . 3)
begin scalar l$
 return if l:=get(p,'derivs) then 
    <<while caaar l neq f do l:=cdr l$
      if l then cons(reval cons('DF,caar l),cdar l)>>
 else cons(nil,0)$
end$

symbolic procedure ldiffp(p,f)$
%  liefert Liste der Variablen + Ordnungen mit Potenz
%  p Ausdruck in LISP - Notation, f Funktion
ld_deriv_search(p,f,fctargs f)$

symbolic procedure ld_deriv_search(p,f,vl)$
%  gets leading derivative of function f in expr. p w.r.t
%  list of variables vl
begin scalar a$
if p=f then a:=cons(nil,1)
else
<<a:=cons(nil,0)$
if pairp p then
  if member(car p,'(PLUS TIMES QUOTIENT EQUAL)) then
     <<p:=cdr p$
       while p do
	 <<a:=dfmax(ld_deriv_search(car p,f,vl),a,vl)$    
	   if cdr a='infinity then p:=nil
			      else p:=cdr p
	 >>
     >>
  else if car p='MINUS then a:=ld_deriv_search(cadr p,f,vl)
  else if car p='EXPT then 
     <<a:=ld_deriv_search(cadr p,f,vl)$    
       if numberp cdr a then
	  if numberp caddr p
	  then a:=cons(car a,times(caddr p,cdr a))
	  else if not zerop cdr a
	    then a:=cons(nil,'infinity)
	    else if not my_freeof(caddr p,f)
	            then a:=cons(nil,'infinity)
     >>
  else if car p='DF then 
	   if cadr p=f then a:=cons(cddr p,1) 
	   else if my_freeof(cadr p,f) 
		then a:=cons(nil,0)                       %  a constant
		else a:=cons(nil,'infinity)
  else if my_freeof(p,f) then a:=cons(nil,0)
		       else a:=cons(nil,'infinity)
>>$
return a
end$

symbolic procedure lderiv(p,f,vl)$
%  fuehrende Ableitung in LISP-Notation mit Potenz (als dotted pair)
begin scalar l$
l:=ld_deriv_search(p,f,vl)$
return cons(if car l then cons('DF,cons(f,car l))
		  else if zerop cdr l then nil
			  else f
		,cdr l)
end$

symbolic procedure splitinhom(q,ftem,vl)$
% Splitting the equation q into the homogeneous and inhom. part
% returns dotted pair qhom . qinhom
begin scalar qhom,qinhom,denm;
  vl:=varslist(q,ftem,vl)$
  if pairp q and (car q = 'QUOTIENT) then 
  if starp(smemberl(ftem,caddr q),length vl) then 
  <<denm:=caddr q; q:=cadr q>>        else return (q . 0)
                                     else denm:=1;

  if pairp q and (car q = 'PLUS) then q:=cdr q
				 else q:=list q;
  while q do <<
   if starp(smemberl(ftem,car q),length vl) then qinhom:=cons(car q,qinhom)
      	       	     	      	     else qhom  :=cons(car q,qhom);
   q:=cdr q
  >>;
  if null qinhom then qinhom:=0
		 else
  if length qinhom > 1 then qinhom:=cons('PLUS,qinhom)
		       else qinhom:=car qinhom;
  if null qhom then qhom:=0
	       else
  if length qhom   > 1 then qhom:=cons('PLUS,qhom)
		       else qhom:=car qhom;
  if denm neq 1 then <<qhom  :=list('QUOTIENT,  qhom,denm);
                       qinhom:=list('QUOTIENT,qinhom,denm)>>;
  return qhom . qinhom
end$

symbolic procedure search_den(l)$
% get all denominators and arguments of LOG,... anywhere in a list l
begin scalar l1$
if pairp l then 
 if car l='quotient then 
  l1:=union(cddr l,union(search_den(cadr l),search_den(caddr l)))
 else if member(car l,'(log ln logb log10)) then 
  if pairp cadr l and (caadr l='QUOTIENT) then 
   l1:=union(list cadadr l,search_den(cadr l))
  else l1:=union(cdr l,search_den(cadr l))
 else for each s in l do l1:=union(search_den(s),l1)$
return l1$
end$

symbolic procedure zero_den(l,ftem,vl)$
begin scalar cases$
 l:=search_den(l)$
 while l do <<
  if not freeofzero(car l,ftem,vl,ftem) then cases:=cons(car l,cases);
  l:=cdr l
 >>$
 return cases
end$

symbolic procedure forg_int(forg,fges)$
for each ex in forg collect
 if pairp ex and pairp cadr ex then forg_int_f(ex,smemberl(fges,ex))
                             else ex$

symbolic procedure forg_int_f(ex,fges)$
% try to integrate expr. ex of the form df(f,...)=expr.
begin scalar p,h,f$
 p:=caddr ex$
 f:=cadadr ex$
 if pairp p and (car p='PLUS)
    then p:=reval cons('PLUS,cons(list('MINUS,cadr ex),cdr p))
    else p:=reval list('DIFFERENCE,p,cadr ex)$
 p:=integratepde(p,cons(f,fges),nil,nil,nil)$
 if p and (car p) and not cdr p then 
    <<h:=car lderiv(car p,f,fctargs f)$
    p:=reval list('PLUS,car p,h)$
    for each ff in fnew_ do
      if not member(ff,ftem_) then ftem_:=fctinsert(ff,ftem_)$
    ex:=list('EQUAL,h,p)>>$
 return ex
end$

symbolic operator total_alg_mode_deriv$
symbolic procedure total_alg_mode_deriv(f,x)$
begin scalar tdf$ %,u,uli,v,vli$
 tdf:={'DF,f,x}$
% explicit program for chain rule of differentiation which is not used
% as currently f=f(u), u=u(x) gives df(f**2,x)=2*f*df(f,x)
%
% for each u in depl!* do 
% if not freeof(cdr u,x) then uli:=cons(car u,uli)$
% for each u in uli do <<
%  vli:=nil$
%  for each v in depl!* do 
%  if not freeof(cdr v,u) then vli:=cons(car v,vli)$
%  algebraic ( tdf:=tdf+df(f,v)*df(v,u)*df(u,x) )$  
% >>$
 return reval tdf
end$

symbolic procedure no_of_v(v,l)$
% v is a variable name, l a list of derivatives like (x 2 y z 3)
% it returns the order of v-derivatives
<<while l and car l neq v do l:=cdr l;
  if null l then 0 else
  if null cdr l or not fixp cadr l or (cadr l = 1) then 1 else
  cadr l
>>$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  general purpose procedures  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure memberl(a,b)$
%  member for a list
if a and b then
if member(car a,b) then cons(car a,memberl(cdr a,b))
                   else memberl(cdr a,b)$
                       
symbolic procedure smemberl(fl,ex)$
%  smember for a list
if fl and ex then
if smember(car fl,ex) then cons(car fl,smemberl(cdr fl,ex))
                      else smemberl(cdr fl,ex)$
                       
symbolic operator my_freeof$
symbolic procedure my_freeof(u,v)$
%  a patch for FREEOF in REDUCE 3.5
not(smember(v,u)) and freeofdepl(depl!*,u,v)$

lisp flag('(my_freeof),'BOOLEAN)$

symbolic procedure freeoflist(l,m)$
%   liefert t, falls kein Element aus m in l auftritt
if null m then t
else if freeof(l,car m) then freeoflist(l,cdr m)
                        else nil$

symbolic procedure freeofdepl(de,u,v)$
if null de then t
else if smember(v,cdar de) and smember(caar de,u) then nil
else freeofdepl(cdr de,u,v)$

symbolic procedure fctins(f,flen,ftem)$
if null ftem then list f else 
if fctlength car ftem < flen then cons(f,ftem)
else cons(car ftem,fctinsert(f,cdr ftem))$

symbolic procedure fctinsert(f,ftem)$
% isert a function f in the function list ftem
if freeof(ftem,f) then fctins(f,fctlength f,ftem)
                  else ftem$

symbolic procedure newfct(id,l,nfct)$
begin scalar f$
 % Only in the top level function names may be recycled otherwise
 % name clashes occur when passing back solutions with new functions
 % of integration but old used names
 if (null level_) and (id=fname_) and recycle_fcts then <<
  f:=car recycle_fcts$
  recycle_fcts:=cdr recycle_fcts
 >>                              else
 f:=mkid(id,nfct)$
 depl!*:=delete(assoc(f,depl!*),depl!*)$
 %put(f,'simpfn,'simpiden)$
 %if pairp l then f:=cons(f,l)$
 if pairp l then depl!*:=cons(cons(f,l),depl!*)$
 if print_ then
   <<terpri()$
   if pairp l then
     <<write "new function: "$
     fctprint list f>>
   else
     write "new constant: ",f>>$
 return f$
end$

symbolic procedure drop_fct(f)$
% check before that f is not one of the forg functions!
% check dropping f also from ftem_
<<if do_recycle_fnc then recycle_fcts:=f . recycle_fcts$
  depl!*:=delete(assoc(reval f,depl!*),depl!*)$
>>$

symbolic procedure varslist(p,ftem,vl)$
begin scalar l$
ftem:=argset smemberl(ftem,p)$
for each v in vl do
  if not my_freeof(p,v) or member(v,ftem) then l:=cons(v,l)$
return reverse l$
end$

symbolic procedure var_list(pdes,forg,vl)$
begin scalar l,l1$
for each p in pdes do l:=union(get(p,'vars),l)$
for each v in vl do
  if member(v,l) or not my_freeof(forg,v) then l1:=cons(v,l1)$
return reverse l1$
end$

symbolic procedure fctlist(ftem,pdes,forg)$
begin scalar fges,l$
 for each p in pdes do l:=union(get(p,'fcts),l)$
 for each f in ftem do
     if not freeof(forg,f) or member(f,l) then fges:=fctinsert(f,fges)$ 
 for each f in forg do
     if not pairp f and not member(f,fges) then fges:=fctinsert(f,fges)$
 for each f in l do
     if not member(f,fges) then fges:=fctinsert(f,fges)$
 l:=setdiff(ftem,fges);
 for each f in l do drop_fct(f)$
 return fges$
end$   

symbolic operator fargs$
symbolic procedure fargs f$
cons('LIST,fctargs f)$

symbolic procedure fctargs f$
%  arguments of a function
if (f:=assoc(f,depl!*)) then cdr f$

symbolic procedure fctlength f$
%  number of arguments
length fctargs f$

symbolic procedure fctsort(l)$
% list sorting
begin scalar l1,l2,l3,m,n$
return
if null l then nil
else
<<n:=fctlength car l$
  l2:=list car l$
  l:=cdr l$
  while l do
   <<m:=fctlength car l$
   if m<n then l1:=cons(car l,l1)
   else if m>n then l3:=cons(car l,l3)
   else l2:=cons(car l,l2)$
   l:=cdr l>>$
  append(fctsort reverse l3,append(reverse l2,fctsort reverse l1))>>
end$

symbolic procedure listprint(l)$
%   print elements of a lisp list
if pairp l then <<
 write car l$
 for each v in cdr l do <<write ",",v>> 
>>$ 

symbolic procedure fctprint1(f)$
%   print a function
begin scalar vl;
if f then
 if pairp f then <<
  write car f$
  if pairp cdr f then <<
   for each a in vl_ do 
   if not freeof(cdr f,a) then vl:=cons(a,vl);
   write "("$
%   listprint cdr f$
   listprint append(setdiff(cdr f,vl),reverse vl)$
   write ")">>
  >>
 else write f$
end$

symbolic procedure fctprint(fl)$
%   Ausdrucken der Funktionen aus fl
begin scalar l,f,a,n,nn$
  n:=0$
  while fl do <<
    f:=car fl$
    fl:=cdr fl$
    if pairp f then
      if car f='EQUAL then
        <<n:=no_of_terms caddr f$
        if (null print_) or (n>print_) then
          <<terpri()$write cadr f,"= expr. with ",n," terms"$
          if (l:=get(cadr f,'fcts)) then
             <<write " in "$
             myprin2l(l,", ")>>$
          terpri()>>
        else mathprint f$
        n:=0>>
      else
         <<if n = 4 then
              <<terpri()$n:=0>>$
         fctprint1 f$
         if fl then write ", "$
         n:=add1 n>>
   else <<
     nn:=reval {'PLUS,4,length explode f,
                for each a in fctargs f sum add1 length explode a};
     if nn+n > 79 then <<terpri()$n:=0>>$
     l:=assoc(f,depl!*)$
     fctprint1 if l then l
                    else f$
     if fl then write ", "$
     n:=nn+n>>
  >>$
end$

symbolic procedure deprint(l)$
%   Ausdrucken der Gl. aus der Liste l
if l and print_ then for each x in l do eqprint list('EQUAL,0,x)$

symbolic procedure eqprint(e)$
%   Ausdrucken der Gl. e
if print_ then
 begin scalar n$
 n:=delength e$
 if n>print_ then
        <<write %"expr. with ",
                n," factors in ",
                no_of_terms(e)," terms"$
          terpri()
        >>
 else
        mathprint e$
end$

symbolic procedure print_level(neu)$
if print_ and level_ then <<
 terpri()$
 if neu then write "New level : "
        else write "Current level : "$
 for each m in reverse level_ do write m,"."$
>>$

symbolic procedure print_statistic(pdes,fcts)$
if print_ then begin
 integer j,k,le,r,s$
 scalar n,m,p,el,fl,vl,pl$

 %--- printing the stats of equations:
 if pdes then <<
  terpri()$write "number of equations : ",length pdes$
  terpri()$write "total no of terms   : ",
  j:=for each p in pdes sum get(p,'terms)$
  k:=for each p in pdes sum get(p,'length)$
  if k neq j then <<terpri()$write "total no of factors :  ",k>>$

  while pdes do <<
   j:=0;
   el:=nil;
   for each p in pdes do <<
    vl:=get(p,'vars);
    if vl then le:=length vl
          else le:=0;
    if ((j=0) and null vl) or
       (j=le) then el:=cons(p,el)
              else if j<le then <<
     j:=le;
     el:=list(p)
    >>
   >>;
   pdes:=setdiff(pdes,el);
   if el then <<
    n:=length el$
    terpri()$write n," equation"$
    if n>1 then write"s"$write" in ",j," variable"$
    if j neq 1 then write"s"$
    write": "$

    if struc_eqn then el:=sort_deriv_pdes(el)$
    repeat <<
     if struc_eqn then <<
      pl:=first el; el:=cdr el;
      terpri()$
      write length cdr pl," equations with ",car pl," derivative",
            if car pl = 1 then ":" else "s:"$
      pl:=cdr pl
     >>           else <<pl:=el;el:=nil>>;
     terpri()$
     k:=0;
     while pl do <<
      if (k geq 70) then <<k:=0;terpri()>>$
      k:=k+4+length explode car pl + length explode get(car pl,'terms)$
      write car pl,"(",get(car pl,'terms)$
      r:=get(car pl,'val)$
      if (s:=get(car pl,'starde)) then <<
       for r:=1:(1+cdr s) do write"*"$
       k:=k+1+cdr s;
      >>$
      if (pairp r) and (car r='TIMES) then write"#"$
      if flin_ and freeoflist(r,flin_) then write"a"$
      write")"$
      pl:=cdr pl$
      if pl then write","$
     >>;

    >> until null el;

   >>$
   j:=add1 j;
  >>
 >>
 else <<terpri()$write "no equations">>$
 %--- printing the stats of functions:
 for each f in fcts do if not pairp f then fl:=cons(f,fl)$
 if fl then <<
  fl:=fctsort reverse fl$
  m:=fctlength car fl$
  while m>=0 do <<
   n:=0$
   el:=nil;
   while fl and (fctlength car fl=m) do <<
    n:=add1 n$
    el:=cons(car fl,el)$
    fl:=cdr fl
   >>$
   if n>0 then
   if m>0 then <<
    terpri()$
    write n," function"$
    if n>1 then write"s"$
    write" with ",m," argument",if m>1 then "s : "
                                       else "  : "
   >>     else <<
    terpri()$
    write n," constant"$
    if n>1 then write"s"$
    write" : "
   >>$
   k:=0;
   while el do <<
    if k=10 then <<k:=0;terpri()>>
            else k:=add1 k$
    write car el$
    el:=cdr el$
    if el then write","$
   >>$
   m:=if fl then fctlength car fl
            else -1
  >> 
 >>    else <<terpri()$write "no functions or constants">>$
 terpri()$
end$

symbolic procedure get_statistic(pdes,fcts)$
 % returns:    {time(),
 %              stepcounter_,
 %              number of pdes,
 %              number of terms,
 %              length of pdes,
 %              {{no of eq, no of var in eq}, ...}
 %              {{no of fc, no of var in fc}, ...}
 %             }
if contradiction_ then "contradiction" else
begin
 integer j,le$
 scalar n,p,el,fl,vl,li,stats$

 stats:={for each p in pdes sum get(p,'length),
         for each p in pdes sum get(p,'terms),
         length pdes,
         stepcounter_,
         time()}$

 %--- the statistics of equations:
 while pdes do <<
  % j is number of variables and el the list of equations
  j:=0;
  el:=nil;
  for each p in pdes do <<
   vl:=get(p,'vars);
   if vl then le:=length vl
         else le:=0;
   if ((j=0) and null vl) or
      (j=le) then el:=cons(p,el)
             else if j<le then <<
    j:=le;
    el:=list(p)
   >>
  >>;
  pdes:=setdiff(pdes,el);
  li:=cons({length el,j},li)
  % length el equations in j variables
 >>;
 stats:=cons(li,stats)$
 li:=nil;

 %--- the statistics of functions:
 for each f in fcts do if not pairp f then fl:=cons(f,fl)$
 if fl then <<
  fl:=fctsort reverse fl$
  j:=fctlength car fl$
  while j>=0 do <<
   n:=0$
   while fl and (fctlength car fl=j) do <<n:=add1 n$ fl:=cdr fl>>$
   li:=cons({n,j},li)$
   % n functions of j variables
   j:=if fl then fctlength car fl
            else -1
  >> 
 >>$

 return reverse cons(li,stats)
end$

symbolic procedure sort_deriv_pdes(pdes)$
begin scalar max_no_deri,cp,pl,res$
 max_no_deri:=0;
 cp:=pdes;
 while cp do << 
  if get(car cp,'no_derivs)>max_no_deri then 
  max_no_deri:=get(car cp,'no_derivs);
  cp:=cdr cp
 >>;
 repeat <<
  pl:=nil;
  cp:=pdes;
  while cp do << 
   if get(car cp,'no_derivs)=max_no_deri then pl:=cons(car cp,pl);
   cp:=cdr cp
  >>$
  if pl then res:=cons(cons(max_no_deri,reverse pl),res)$
  pdes:=setdiff(pdes,pl);
  max_no_deri:=if zerop max_no_deri then nil
                                    else sub1(max_no_deri);
 >> until (null max_no_deri) or (null pdes);
 return res
end$

symbolic procedure print_pdes(pdes)$
% print all pdes up to some size
begin scalar pl,ps,n,pdecp$
  terpri()$
  if pdes then <<
   if (null !*batch_mode) and 
      (batchcount_<stepcounter_) and
      (cdr pdes) then << % if more than one pde
    write"What is the maximal number of terms of equations to be shown? "$
    ps:=promptstring!*$
    promptstring!*:=""$
    terpri()$n:=termread()$
    promptstring!*:=ps$ 
    for each pl in pdes do 
    if get(pl,'terms)<=n then pdecp:=cons(pl,pdecp);
    pdecp:=reverse pdecp;
   >>          else pdecp:=pdes$

   write "equations : "$ 
   if struc_eqn then <<
    pl:=sort_deriv_pdes(pdecp)$
    while pl do <<
     terpri()$
     write length cdar pl," equations with ",caar pl," derivatives:"$
     typeeqlist(cdar pl)$
     pl:=cdr pl     
    >>
   >>           else typeeqlist(pdecp)
  >>      else <<write "no equations"$ terpri()>>$
end$

symbolic procedure print_ineq(ineqs)$
% print all ineqs
begin scalar a,b,c$
 terpri()$
 if ineqs then <<
  terpri()$write "non-vanishing expressions: "$
  for each a in ineqs do if pairp a then c:=cons(a,c) 
                                    else b:=cons(a,b);
  listprint b;terpri()$
  for each a in c do eqprint a
 >>
end$

symbolic procedure print_fcts2(pdes,fcts)$
% print all fcts and vars
begin scalar dflist,dfs,f,p,cp,h,hh,ps,showcoef$

for each h in fcts do if not pairp h then hh:=cons(h,hh);
 ps:=promptstring!*$ promptstring!*:=""$ 
% write "Which function out of "$terpri()$
% listprint(hh)$
% write"? "$terpri()$
% write"Enter the function name only or ; for all functions."$terpri()$
%
% h:=termread();
% if h neq '!; and not_included(list h,hh) then <<
%  write"This is not a function in the list."$
%  terpri();
%  h:=nil
% >>;
% if null h then return nil;
% if h='!; then fcts:=hh
%          else fcts:=list h;

 fcts:=select_from_list(hh,nil)$
 pdes:=select_from_list(pdes,nil)$

 write"Do you want to see the coefficients of all derivatives in all equations"$
 terpri()$
 write"in factorized form which may take relatively much time? y/n"$
 terpri()$
 repeat
  h:=termread()
 until (h='y) or (h='n);
 if h='n then showcoef:=nil else showcoef:=t;

 promptstring!*:=ps$

 while fcts do 
 if pairp car fcts then fcts:=cdr fcts
                   else <<
  f:=car fcts;  fcts:=cdr fcts;
  dflist:=nil;
  for each p in pdes do if not freeof(get(p,'fcts),f) then <<
   dfs:=get(p,'derivs);
   while dfs do <<
    if caaar dfs=f then <<
     cp:=dflist;
     while cp and (caar cp neq caar dfs) do cp:=cdr cp;
     if cdaar dfs then h:=cons('DF,caar dfs)
                  else h:=caaar dfs;
     if showcoef then 
     if null cp then dflist:=cons({caar dfs,{'LIST,p,
                                             factorize coeffn(get(p,'val),h,1)}},
                                  dflist)  
                else rplaca(cp,cons(caar cp,
                                    cons({'LIST,p,
                                          factorize coeffn(get(p,'val),h,1)},
                                         cdar cp)))
               else
     if null cp then dflist:=cons({caar dfs,p},dflist)  
                else rplaca(cp,cons(caar cp,cons(p,cdar cp)))
    >>;
    dfs:=cdr dfs
   >>;
  >>;
  while dflist do <<
   dfs:=car dflist;dflist:=cdr dflist;
   if cdar dfs then h:=cons('DF,car dfs)
               else h:=caar dfs;
   if showcoef then algebraic <<write h,": ",lisp cons('LIST,cdr dfs)>>
               else <<write h,": "$ print cdr dfs$ terpri()>>
  >>;
 >>;
end$

symbolic procedure print_fcts(fcts,vl)$
% print all fcts and vars
<<if fcts then
         <<terpri()$write "functions : "$
           fctprint(fcts);
           terpri()$write "with ",
           for each p in fcts sum no_of_terms(p)," terms">>$
 if vl then
         <<terpri()$write "variables : "$
         fctprint(vl)>>$
>>$

symbolic procedure print_pde_fct_ineq(pdes,ineqs,fcts,vl)$
% print all pdes, ineqs and fcts
if print_ then begin$
 print_pdes(pdes)$
 print_ineq(ineqs)$
 print_fcts(fcts,vl)$
 print_statistic(pdes,fcts)
end$ 

symbolic procedure no_of_terms(d)$
if not pairp d then if (null d) or (zerop d) then 0 
                                             else 1 else 
if car d='PLUS then length d - 1            else 
if car d='EQUAL then no_of_terms(cadr  d) +
                     no_of_terms(caddr d)   else
if (car d='MINUS) or (car d='QUOTIENT) then 
   no_of_terms(cadr d)                 else 
if car d='EXPT then 
if (not fixp caddr d) or (caddr d < 2) then 1 else
% number of terms of (a1+a2+..+an)**r = n+r-1 over r
begin scalar h,m,q$
 m:=no_of_terms(cadr d)-1;
 h:=1;
 for q:=1:caddr d do h:=h*(m+q)/q;
 return h
end            else 
if car d='TIMES then begin scalar h,r;
 h:=1;
 for each r in cdr d do h:=h*no_of_terms(r);
 return h
end             else 1$

symbolic procedure delength(d)$
%   Laenge eines Polynoms in LISP - Notation
if not pairp d then 
 if d then 1
      else 0
else
if (car d='PLUS) or (car d='TIMES) or (car d='QUOTIENT) 
   or (car d='MINUS) or (car d='EQUAL)
   then for each a in cdr d sum delength(a) 
else 1$

symbolic procedure pdeweight(d,ftem)$
%   Laenge eines Polynoms in LISP - Notation
if not smemberl(ftem,d) then 0
else if not pairp d then 1
else if (car d='PLUS) or (car d='TIMES) or (car d='EQUAL)
        or (car d='QUOTIENT) then
         for each a in cdr d sum pdeweight(a,ftem)
else if (car d='EXPT) then 
        if numberp caddr d then
         caddr d*pdeweight(cadr d,ftem)
                           else       
        pdeweight(caddr d,ftem)+pdeweight(cadr d,ftem)
else if (car d='MINUS) then pdeweight(cadr d,ftem)
else 1$

symbolic procedure desort(l)$
% sort expressions hat are the elements of the list l by size
for each a in idx_sort for each b in l collect cons(delength b,b)
collect cdr a$

symbolic procedure idx_sort(l)$
% All elements of l have a numerical first element and are sorted
% by that number
begin scalar l1,l2,l3,m,n$
return
if null l then nil
else
<<n:=caar l$
  l2:=list car l$
  l:=cdr l$
  while l do
   <<m:=caar l$
   if m<n then l1:=cons(car l,l1)
   else if m>n then l3:=cons(car l,l3)
   else l2:=cons(car l,l2)$
   l:=cdr l>>$
  append(idx_sort(l1),append(l2,idx_sort(l3)))>>
end$

symbolic procedure rat_idx_sort(l)$
% All elements of l have a rational number first element
% and are sorted by that number
% The rational number has to be reval-ed !
begin scalar l1,l2,l3,m,n$
return
if null l then nil
else
<<n:=caar l$
  l2:=list car l$
  l:=cdr l$
  while l do
   <<m:=caar l$
   if rational_less(m,n) then l1:=cons(car l,l1)
   else if rational_less(n,m) then l3:=cons(car l,l3)
   else l2:=cons(car l,l2)$
   l:=cdr l>>$
  append(rat_idx_sort(l1),append(l2,rat_idx_sort(l3)))>>
end$

symbolic procedure argset(ftem)$
%  List of arguments of all functions in ftem
if ftem then union(reverse fctargs car ftem,argset(cdr ftem))
        else nil$

symbolic procedure no_fnc_of_v$
begin scalar vl,v,nofu,f,nv$
  % How many functions do depend on each variable?
  vl:=argset(ftem_)$
  for each v in vl do <<
    nofu:=0;  % the number of functions v occurs in
    for each f in ftem_ do
    if not freeof(fctargs f,v) then nofu:=add1 nofu$
    nv:=cons((v . nofu),nv)$
  >>$
  return nv
end$

procedure push_vars(liste)$
for each x in liste collect 
if not boundp x then x else eval x$ % valuecell x$

symbolic procedure backup_pdes(pdes,forg)$
%  make a backup of all pdes
begin scalar allfl$
 return
 list(push_vars glob_var,
      for each p in pdes collect 
      list(p,
           for each q in prop_list collect cons(q,get(p,q)),
           <<allfl:=nil;
             for each q in allflags_ do 
             if flagp(p,q) then allfl:=cons(q,allfl);
             allfl>>),
      for each f in forg collect
           if pairp f then cons(f,get(cadr f,'fcts))
                      else cons(f,get(     f,'fcts)),
      for each id in idnties_ collect
      list(id,get(id,'val),flagp(id,'to_int),flagp(id,'to_subst))
     )
end$

%symbolic procedure backup_pdes(pdes,forg)$
%%  make a backup of all pdes
%begin scalar cop$
% cop:=list(nequ_,
%           for each p in pdes collect 
%               list(p,
%                    for each q in prop_list collect cons(q,get(p,q)),
%                    for each q in allflags_ collect if flagp(p,q) then q),
%           for each f in forg collect
%               if pairp f then cons(cadr f,get(cadr f,'fcts))
%                          else cons(f,get(f,'fcts)),
%           ftem_,
%           ineq_,
%           recycle_ens,
%           recycle_fcts)$
% return cop
%end$ 

symbolic procedure pop_vars(liste,altewerte)$
foreach x in liste do <<set (x,car altewerte);
                        altewerte := cdr altewerte>>$

symbolic procedure restore_pdes(bak)$
%  restore all data: glob_var, pdes, forg
begin scalar pdes,forg$

  % recover values of global variables
  pop_vars(glob_var,car bak)$
  % recover the pdes
  for each c in cadr bak do <<
    pdes:=cons(car c,pdes)$
    for each s in cadr  c do put(car c,car s,cdr s)$
    for each s in caddr c do flag1(car c,s)
  >>$

  % recover the properties of forg
  for each c in caddr bak do <<
    forg:=cons(car c,forg)$
    if pairp car c then put(cadar c,'fcts,cdr c)
  >>$

  % recover the properties of idnties_
  if cdddr bak then
  for each c in cadddr bak do <<
    put(car c,'val,cadr c);
    if caddr c then flag1(car c,'to_int)
               else if flagp(car c,'to_int) then remflag(car c,'to_int);
    if caddr c then flag1(car c,'to_subst)
               else if flagp(car c,'to_subst) then remflag(car c,'to_subst);
  >>$

  return {reverse pdes,reverse forg}$
end$

%symbolic procedure restore_pdes(cop)$
%%  restore the pde list cop 
%%  first element must be the old value of nequ_ 
%%  the elements must have the form (p . property_list_of_p)
%begin scalar pdes$
%  % delete all new produced pdes
%  for i:=car cop:sub1 nequ_ do setprop(mkid(eqname_,i),nil)$
%  nequ_:=car cop$
%  for each c in cadr cop do
%    <<pdes:=cons(car c,pdes)$
%    for each s in cadr c do
%        put(car c,car s,cdr s)$
%    for each s in caddr c do 
%        if s then flag1(car c,s)$
%    >>$
%  for each c in caddr cop do
%    put(car c,'fcts,cdr c)$
%  ftem_:=cadddr cop$
%  ineq_:=car cddddr cop$
%  recycle_eqns:= cadr cddddr cop$
%  recycle_fcts:= caddr cddddr cop$
%  return reverse pdes$
%end$

%symbolic procedure copy_from_backup(copie)$
%%  restore the pde list copie
%%  first element must be the old value of nequ_ 
%%  the elements must have the form (p . property_list_of_p)
%%  at least recycle_eqns should not work with this procedure
%begin scalar l,pdes,cop$
%  cop:=cadr copie$         
%  l:=cop$
%  for i:=1:length cop do
%    <<pdes:=cons(mkid(eqname_,nequ_),pdes)$
%    setprop(car pdes,nil)$
%    nequ_:=add1 nequ_>>$
%  pdes:=reverse pdes$
%  for each p in pdes do  
%    <<cop:=subst(p,caar l,cop)$
%    l:=cdr l>>$
%  for each c in cop do
%    <<for each s in cadr c do
%        put(car c,car s,cdr s)$
%    for each s in caddr c do 
%        if s then flag1(car c,s)$
%    >>$
%  for each c in caddr copie do
%    put(car c,'fcts,cdr c)$
%  ftem_:=cadddr copie$
%  return pdes$
%end$

symbolic procedure deletepde(pdes)$
begin scalar s,l,ps$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 terpri()$
 write "Equations to be deleted: "$
 l:=select_from_list(pdes,nil)$
 promptstring!*:=ps$
 for each s in l do
 if member(s,pdes) then pdes:=drop_pde(s,pdes,nil)$
 return pdes
end$

symbolic procedure new_pde()$
begin scalar s$
 if null car recycle_eqns then <<
  s:=mkid(eqname_,nequ_)$
  nequ_:=add1 nequ_$
 >>                       else <<
  s:=caar recycle_eqns$
  recycle_eqns:=(cdar recycle_eqns) . (cdr recycle_eqns)
 >>$
 setprop(s,nil)$
 return s
end$

symbolic procedure drop_pde_from_properties(p,pdes)$
begin
 for each q in pdes do if q neq p then <<
  drop_dec_with(p,q,t)$
  drop_dec_with(p,q,nil)$
  drop_rl_with(p,q)
 >>
end$

symbolic procedure drop_pde_from_idties(p,pdes,pval)$
% to be used whenever the equation p is dropped or changed and
% afterwards newly characterized in update,
% pval is the new value of p in terms of the previous value,
% if this is unknown then pval=nil
% to be done before setprop(p,nil)
begin scalar q,newidval,idnt$
 for each q in pdes do if q neq p then
 if not freeof(get(q,'histry_),p) then 
 put(q,'histry_, if null pval then q
                              else subst(pval,p,get(q,'histry_)))$

 if record_hist and (getd 'show_id) then <<
  idnt:=idnties_$
  while idnt do <<
   if not freeof(get(car idnt,'val),p) then 
   if null pval then drop_idty(car idnt)
                else <<
    % Once pdes_ is available as global variable then simplify 'val
    % before put()
    newidval:=reval subst(pval,p,get(car idnt,'val))$
    if trivial_idty(pdes,newidval) then drop_idty(car idnt)
                                   else <<
     put(car idnt,'val,newidval);
     flag1(car idnt,'to_subst)$
     flag1(car idnt,'to_int)
    >>
   >>;
   idnt:=cdr idnt
  >>;
  if pval and not zerop pval and (p neq get(p,'histry_)) then <<
   idnt:=reval num reval {'PLUS,get(p,'histry_),{'MINUS,pval}}$
   if not zerop idnt then 
   new_idty(idnt,pdes,if pdes then t else nil)
  >>
 >>
end$

symbolic procedure drop_pde(p,pdes,pval)$
% pval is the value of p in terms of other equations, 
% if pval=nil then unknown
% pdes should be a list of all currently used pde-names
begin 
 if do_recycle_eqn then
 recycle_eqns:=(car recycle_eqns) . union({p},cdr recycle_eqns)$
 depl!*:=delete(assoc(reval p,depl!*),depl!*)$
 drop_pde_from_idties(p,pdes,pval)$
 setprop(p,nil)$
 return delete(p,pdes)
end$

symbolic procedure change_pde_flag(pdes)$
begin scalar p,ty,h$
 repeat <<
  terpri()$
  write "From which PDE do you want to change a ",
        "flag or property, e.g. e_23?"$
  terpri()$
  p:=termread()$
 >> until not freeof(pdes,p)$
 terpri()$
 write"Type in one of the following flags that is to be flipped "$
 terpri()$
 write"(e.g. to_int <ENTER>): "$
 terpri()$terpri()$
 write allflags_;
 terpri()$terpri()$
 write"or type in one of the following properties that is to be changed"$
 terpri()$
 write"(e.g. vars <ENTER>): "$
 terpri()$terpri()$
 write prop_list;
 terpri()$terpri()$
 ty:=termread()$
 if member(ty,allflags_) then <<
  if flagp(p,ty) then remflag1(p,ty)
                 else flag1(p,ty)$
  write"The new value of ",ty,": ",flagp(p,ty)$
 >>                      else 
 if member(ty,prop_list) then <<
  terpri()$
  write"current value: ",get(p,ty)$
  terpri()$
  write"new value (e.g. (x y z) ENTER): "$
  terpri()$
  h:=termread()$ 
  put(p,ty,h)$
  write"The new value of ",ty,": ",get(p,ty)$
 >>                      else write"Input not recognized."$
 terpri()$ 
end$

symbolic procedure restore_backup_from_file(pdes,forg,nme)$
begin scalar ps,s,p,echo_bak$ 
  if nme=t then <<
    ps:=promptstring!*$
    promptstring!*:=""$ 
    terpri()$
    write"Please give the name of the file in double quotes"$terpri()$
    write"without `;' : "$
    s:=termread()$  
  >>     else 
  if nme then s:=nme
         else s:=level_string(session_)$
  % infile s$
  echo_bak:=!*echo; semic!*:='!$; in s$ !*echo:=echo_bak;
  %-- cleaning up:
  for each p in pdes do setprop(p,nil)$
  for each p in forg do if pairp p then put(cadr p,'fcts,nil)$
  %-- assigning the new values:
  promptstring!*:=ps$
  size_hist   :=car backup_; backup_:=cdr backup_;
  stepcounter_:=car backup_; backup_:=cdr backup_;
  level_      :=car backup_; backup_:=cdr backup_;
  nfct_       :=car backup_; backup_:=cdr backup_;
  time_limit  :=car backup_; backup_:=cdr backup_;
  limit_time  :=car backup_; backup_:=cdr backup_;
  history_    :=car backup_; backup_:=cdr backup_;
  s:=restore_pdes(backup_)$
  backup_:=nil;
  orderings_:=car orderings_;
  return s
end$

symbolic procedure level_string(s)$
 << for each m in reverse level_ do s := if s then bldmsg("%w%d.",s,m)
                                               else bldmsg("%d.",m  );
 s>>;

symbolic procedure backup_to_file(pdes,forg,nme)$
begin scalar s,ps$ %,levelcp$
  if nme=t then <<
    ps:=promptstring!*$
    promptstring!*:=""$ 
    terpri()$
    write"Please give the name of the file in double quotes"$terpri()$
    write"without `;' : "$
    s:=termread()$
    promptstring!*:=ps$
  >>     else 
  if nme then s:=nme
         else s:=level_string(session_)$
  out s;
  off nat$
  orderings_:=list orderings_;
  write"off echo$"$
  write "backup_:='"$terpri()$
  print cons(size_hist,cons(stepcounter_,cons(level_,cons(nfct_,
        cons(time_limit,cons(limit_time,cons(history_,
        backup_pdes(pdes,forg))))))))$
  write"$"$terpri()$
  write "end$"$terpri()$
  shut s;
  on nat;
end$

symbolic procedure delete_backup$
begin scalar s$
 s:=level_string(session_);
 delete!-file s;
end$

symbolic procedure restore_and_merge(soln,pdes,forg)$
% pdes, forg are cleaned up
% one could just use restore_pdes without assigning bak
% but then it would not be stored in a file, such that 
% rb can reload the file
begin scalar bak,newfdep,sol,f,h$

 % store ongoing global values in bak
 newfdep:=nil$
 for each sol in soln do
 if sol then <<
   for each f in caddr sol do 
   if h:=assoc(f,depl!*) then newfdep:=cons(h,newfdep);
 >>;
 bak:=append(list(size_hist,stepcounter_,level_,nfct_,
                  time_limit,limit_time,history_),
             newfdep);

 h:=restore_backup_from_file(pdes,forg,nil)$

 size_hist   :=car bak; bak:=cdr bak;
 stepcounter_:=car bak; bak:=cdr bak;
 level_      :=car bak; bak:=cdr bak;
 nfct_       :=car bak; bak:=cdr bak;
 time_limit  :=car bak; bak:=cdr bak;
 limit_time  :=car bak; bak:=cdr bak;
 history_    :=car bak; bak:=cdr bak;
 depl!*      :=append(depl!*,bak);
 return h
end$

symbolic procedure write_in_file(pdes,forg)$
begin scalar p,pl,s,h,ps,wn,vl,ftem$
  ps:=promptstring!*$
  promptstring!*:=""$ 
  terpri()$
  write "Enter a list of equations, like e2,e5,e35; from: "$terpri()$
  listprint(pdes)$
  terpri()$write "To write all equations just enter ; "$terpri()$
  repeat <<
    s:=termlistread()$
    h:=s;
    if s=nil then pl:=pdes else <<
      pl:=nil;h:=nil$
      if (null s) or pairp s then <<
        for each p in s do
        if member(p,pdes) then pl:=cons(p,pl);
        h:=setdiff(pl,pdes);
      >> else h:=s;
    >>;
    if h then <<write "These are no equations: ",h,"   Try again."$terpri()>>$
  >> until null h$
  write"Shall the name of the equation be written? (y/n) "$
  repeat s:=termread()
  until (s='y) or (s='Y) or (s='n) or (s='N)$
  if (s='y) or (s='Y) then wn:=t$
  write"Please give the name of the file in double quotes"$terpri()$
  write"without `;' : "$
  s:=termread()$
  out s;
  off nat$

  write"load crack$"$terpri()$
  write"lisp(nfct_:=",nfct_,")$"$terpri()$
  if wn then write"lisp(nequ_:=",nequ_,")$"$terpri()$
  write"off batch_mode$"$terpri()$
  for each p in pl do <<h:=get(p,'vars);if h then vl:=union(h,vl)>>$
  write"list_of_variables:="$
  algebraic write lisp cons('LIST,vl)$

  for each p in pl do <<h:=get(p,'fcts);if h then ftem:=union(h,ftem)>>$
  write"list_of_functions:="$
  algebraic write lisp cons('LIST,ftem)$

  for each h in ftem do 
  if assoc(h,depl!*) then <<
    p:=pl;
    while p and freeof(get(car p,'val),h) do p:=cdr p;
    if p then <<
      write "depend ",h$
      for each v in cdr assoc(h,depl!*) do <<
        write ","$print v
      >>$
      write "$"$terpri()$

%      for each v in cdr assoc(h,depl!*) do 
%      algebraic write "depend ",lisp h,",",lisp v$
    >>
  >>$
  if wn then <<
    for each h in pl do algebraic (write h,":=",lisp (get(h,'val)))$
    write"list_of_equations:="$
    algebraic write lisp( cons('LIST,pl) )
  >>    else <<
    write"list_of_equations:="$
    algebraic write lisp( cons('LIST,
       for each h in pl collect get(h,'val)) )$
  >>$

  write"list_of_inequalities:="$
  algebraic write lisp( cons('LIST,ineq_))$

  terpri()$ write"solution_:=crack(list_of_equations,"$
  terpri()$ write"                 list_of_inequalities,"$
  terpri()$ write"                 list_of_functions,"$
  terpri()$ write"                 list_of_variables)$"$
  terpri()$

  for each h in forg do <<
   terpri()$
   if pairp h and (car h = 'EQUAL) then 
   algebraic
   write lisp(cadr  h)," := sub(second first solution_,",
         lisp(caddr h),")"  
  >>$
  terpri()$
  write"end$"$terpri()$terpri()$
  write"These data were produced with the following input:"$terpri()$terpri()$
  write"lisp( old_history := "$terpri()$
  write"'",reverse history_,")$"$terpri()$
  shut s;
  on nat;
  promptstring!*:=ps$
end$

symbolic procedure give_low_priority(pdes,f)$
% It assumes that f is element of ftem_.
% It assumes that if f is element of flin_ then flin_ functions
% come first in each group of functions with the same number
% of independent variables.
% If f is element of flin_ then f is put at the end of the flin_
% functions with equally many variables but before the first functions
% that occur in ineq_ in order to change ftem_ as little as possible
% not to invalidate previous decoupling.

begin scalar ftemcp,ano,h,s,fli$
 ftemcp:=ftem_$
 while ftemcp and (car ftemcp neq f) do <<
  h:=cons(car ftemcp,h)$
  ftemcp:=cdr ftemcp
 >>$
 % Is there an element of the remaining ftemcp with the same no of
 % variables and that is not in ineq_ ?

 if ftemcp then <<
  ftemcp:=cdr ftemcp;
  ano:=fctlength f$
  if member(f,flin_) then fli:=t$
  while ftemcp do
  if (ano > (fctlength car ftemcp)) or
     (fli and (not member(car ftemcp,flin_))) then ftemcp:=nil else <<
   h:=cons(car ftemcp,h)$
   ftemcp:=cdr ftemcp$
   if not member(car h,ineq_) then <<
    while ftemcp and 
          (ano = (fctlength car ftemcp)) and 
          (not member(car ftemcp,ineq_)) and
          ((not fli) or member(car ftemcp,flin_)) do <<
     h:=cons(car ftemcp,h)$
     ftemcp:=cdr ftemcp
    >>$ 

    if print_ or tr_orderings then <<
     write"The lexicographical ordering of unknowns is changed"$
     terpri()$
     write"because ",f," has to be non-zero, giving ",f," a low priority."$
     terpri()$
     write "Old ordering: "$
     s:=ftem_;
     while s do <<write car s$ s:=cdr s$ if s then write",">>$
     terpri()$
     write "New ordering: "$
     s:=append(reverse h,cons(f,ftemcp));
     while s do <<write car s$ s:=cdr s$ if s then write",">>$
     terpri()$

    >>$
    change_fcts_ordering(append(reverse h,cons(f,ftemcp)),pdes,vl_)$
    ftemcp:=nil
   >>  % of not member(car h,ineq_)
  >>   % of ano > (fctlength car ftemcp)
 >>    % of ftemcp
end$

symbolic procedure addineq(pdes,newineq)$
% it assumes newineq involves functions of ftem_ 
begin scalar h1,h2,h3$
 newineq:=simp_ineq(newineq);
 h1:=cdr err_catch_fac(newineq)$          % h1 is a lisp list of factors
 if null cdr h1 then <<
  h2:=signchange(car h1);                 % only one factor
  if not member(h2,ineq_) then <<ineq_:=cons(h2,ineq_);h3:=cons(h2,h3)>>
 >>             else for each h2 in h1 do <<
  h2:=signchange(h2);
  if (not freeoflist(h2,ftem_)) and
     (not member(h2,ineq_)) then <<ineq_:=cons(h2,ineq_);h3:=cons(h2,h3)>>
 >>$

 if print_ and h3 then <<
  write"The list of inequalities got extended."$terpri()
 >>$
 % h3 is the list of all new non-zero expressions
 % if any one of these expressions is an element of ftem_ then it
 % should get a low priority in the lexicographical ordering for
 % non-linear problems
 if flin_ and null lin_problem then % maybe also for flin_=nil
 for each h2 in h3 do
 if atom h2 and member(h2,ftem_) then give_low_priority(pdes,h2);
 % h2 gets a low priority so that it is eliminated late in decoupling
 % to be available as non-zero coefficient as long as possible to
 % allow substitutions of other functions without case-distinctions

 if pdes then for each h2 in h3 do update_fcteval(pdes,h2);

 % If one term of the equation is non-zero then the sum of the 
 % remaining terms has to be non-zero too
 if h3 and pdes then for each h2 in pdes do 
 if get(h2,'terms)=2 then new_ineq_from_pde(h2,pdes)

end$

% symbolic procedure drop_factor(h,pro)$
% % This procedure drops a factor h or its negative from an expression pro
% begin scalar hs,newpro,mi;
%  hs:=signchange(h);
%  if pairp pro and (car pro='MINUS) then <<pro:=cadr pro; mi:=t>>;
%  if pro = h  then newpro:= 1 else
%  if pro = hs then newpro:=-1 else 
%  if pairp pro and (car pro = 'TIMES) then 
%  if member(h ,pro) then newpro:=reval delete(h ,pro) else
%  if member(hs,pro) then newpro:=reval list('MINUS,delete(hs,pro));
%  if mi and newpro then newpro:=reval list('MINUS,newpro)
%  return newpro
% end$

symbolic procedure update_fcteval(pdes,newineq)$
begin scalar p,pv,ps,hist,h1,h2$
 for each p in pdes do <<
  pv:=get(p,'val)$
  if pairp pv and (car pv='TIMES) and member(newineq,pv) then <<
   pv:=reval {'QUOTIENT,pv,newineq};
   if record_hist then hist:=reval {'QUOTIENT,get(p,'histry_),newineq}$
   update(p,pv,get(p,'fcts),get(p,'vars),t,list(0),pdes)$
   drop_pde_from_idties(p,pdes,hist)$
   drop_pde_from_properties(p,pdes)
  >> else <<
   ps:=get(p,'fcteval_nli)$
   while ps and
     <<h1:=caar ps;
       h2:=signchange(h1);
       (not ((newineq=h1              ) or
             (pairp h1            and
              (car h1 = 'TIMES)   and
               member(newineq,h1)     )    )) and
       (not ((newineq=h2              ) or
             (pairp h2            and
              (car h2 = 'TIMES)   and
               member(newineq,h2)     )    ))
     >> do ps:=cdr ps;
   if ps then << % We would have to check whether apart from the
    % new non-zero factor, the other factors can vanish for
    % specific values of ftem_ or not. Instead of programming
    % this again here we simply change flags to re-compute all
    % fct... properties later when a substitution is to be done.
    flag1(p,'to_eval)$
    put(p,'fcteval_lin,nil)$
    put(p,'fcteval_nca,nil)$
    put(p,'fcteval_nli,nil)$
    put(p,'fct_nli_lin,nil)$
    put(p,'fct_nli_nca,nil)$
    put(p,'fct_nli_nli,nil)$
    put(p,'fct_nli_nus,nil)$
   >>
  >>
 >>$
end$

symbolic procedure addfunction(ft)$
begin scalar f,ff,l,ps,ok$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 ff:=mkid(fname_,nfct_)$
 repeat <<
  ok:=t;
  terpri()$
  write "What is the name of the new function?"$
  terpri()$
  write "If the name is ",fname_,"+digits then use ",ff,". Terminate with <ENTER>: "$
  f:=termread()$
  if f=ff then nfct_:=add1 nfct_
          else if member(f,ft) then <<
   terpri()$
   write"Choose another name. ",f," is already in use."$
   ok:=nil
  >>$
 >> until ok;
 depl!*:=delete(assoc(f,depl!*),depl!*)$
 terpri()$
 write "Give a list of variables ",f," depends on, for example x,y,z;  "$ 
 terpri()$
 write "For constant ",f," input a `;'  "$
 l:=termxread()$
 if (pairp l) and (car l='!*comma!*) then l:=cdr l;
 if pairp l then depl!*:=cons(cons(f,l),depl!*) else
 if l then depl!*:=cons(list(f,l),depl!*)$
 ft:=fctinsert(f,ft)$
 ftem_:=fctinsert(f,ftem_)$
 promptstring!*:=ps$
 return (ft . f)
end$

symbolic procedure newinequ(pdes)$
begin scalar ps,ex$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 write "Input of a value for the new non-vanishing expression."$
 terpri()$
 write "You can use names of pds, e.g. 3*e_12 - df(e_13,x) + 8; "$
 terpri()$
 write "Terminate the expression with ; or $ : "$
 terpri()$
 ex:=termxread()$
 for each a in pdes do ex:=subst(get(a,'val),a,ex)$
 terpri()$
 promptstring!*:=ps$
 addineq(pdes,ex)$
end$

symbolic procedure replacepde(pdes,ftem,vl)$
begin scalar p,q,ex,ps,h,newfct,again$
 ps:=promptstring!*$ promptstring!*:=""$ 
 repeat <<
  terpri()$
  write "Is there a"$
  if again then write" further"$
  write" new function in the changed/new PDE that"$
  terpri()$
  write "is to be calculated (y/n)? "$
  p:=termread()$
  if (p='y) or (p='Y) then <<
   h:=addfunction(ftem)$
   ftem:=car h$ 
   if cdr h then newfct:=cons(cdr h,newfct)
  >>;
  again:=t
 >> until (p='n) or (p='N)$
 terpri()$
 write "If you want to replace a pde then type its name, e.g. e_23 <ENTER>."$
 terpri()$
 write "If you want to add a pde then type `new_pde' <ENTER>. "$
 p:=termread()$
 if (p='NEW_PDE) or member(p,pdes) then
  <<terpri()$write "Input of a value for "$
  if p='new_pde then write "the new pde."
                else write p,"."$
  terpri()$
  write "You can use names of other pds, e.g. 3*e_12 - df(e_13,x); "$
  terpri()$
  write "Terminate the expression with ; or $ : "$
  terpri()$
  ex:=termxread()$
  for each a in pdes do ex:=subst(get(a,'val),a,ex)$
  terpri()$
  write "Do you want the equation to be"$terpri()$
%  write "- left completely unchanged"$
%  terpri()$
%  write "  (e.g. to keep the structure of a product to "$
%  terpri()$
%  write "   investigate subcases)                        (1)"$
%  terpri()$
  write "- simplified (e.g. e**log(x) -> x) without"$
  terpri()$
  write "  dropping non-zero factors and denominators"$
  terpri()$
  write "  (e.g. to introduce integrating factors)       (1)"$
  terpri()$
  write "- simplified completely                         (2) "$
  h:=termread()$
%  if h=2 then ex:=reval ex$
%  if h<3 then h:=nil
%         else h:=t$
  if h=1 then <<ex:=reval ex$h:=nil>>
         else h:=t$
  if p neq 'NEW_PDE then 
  pdes:=drop_pde(p,pdes,{'QUOTIENT,{'TIMES,p,ex},get(p,'val)})$
  q:=mkeq(ex,ftem,vl,allflags_,h,list(0),nil,pdes)$
  % A new equation with a new function appearing linear and only
  % algebraically can only have the purpose of a transformation
  % in which case the new equation should not be solved for the
  % new function as this would just mean dropping the new equation:
  if (p='NEW_PDE) and newfct then
  put(q,'not_to_eval,newfct)$
  terpri()$write q$
  if p='NEW_PDE then write " is added"
                else write " replaces ",p$
  pdes:=eqinsert(q,pdes)>>
 else <<terpri()$
        write "A pde ",p," does not exist! (Back to previous menu)">>$
 promptstring!*:=ps$
 return list(pdes,ftem)
end$

symbolic procedure select_from_list(liste,n)$
begin scalar ps,s$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 terpri()$
 if n then write"Pick ",n," from this list:"
      else write"Pick from this list"$
 terpri()$
 listprint(liste)$write";"$terpri()$
 if null n then <<
  write"a sublist and input it in the same form. Enter ; to choose all."$
  terpri()$
 >>$
 s:=termlistread()$
 if n and n neq length s then <<
  write "Wrong number picked."$terpri()$
  s:=nil;
 >>                      else
 if null s then s:=liste else
 if not_included(s,liste) then <<
  write setdiff(s,liste)," is not allowed."$terpri()$
  s:=nil;
 >>;
 promptstring!*:=ps$ 
 return s
end$

symbolic procedure selectpdes(pdes,n)$
% interactive selection of n pdes  
% n may be an integer or nil. If nil then the
% number of pdes is free.
if pdes then
begin scalar l,s,ps,m$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 terpri()$
 if null n then <<
  write "How many equations do you want to select? "$terpri()$
  write "(number <ENTER>) : "$terpri()$
  n:=termread()$
 >>$ 
 write "Please select ",n," equation"$
 if n>1 then write "s"$write " from: "$
 write pdes$ % fctprint pdes$
 terpri()$
 m:=0$
 s:=t$
 while (m<n) and s do
  <<m:=add1 m$
  if n>1 then write m,". "$
  write "pde: "$
  s:=termread()$
  while not member(s,pdes) do
   <<write "Error!!! Please select a pde from: "$
   write pdes$ % fctprint pdes$
   terpri()$if n>1 then write m,". "$
   write "pde: "$
   s:=termread()>>$
  if s then
   <<pdes:=delete(s,pdes)$
   l:=cons(s,l)>> >>$
 promptstring!*:=ps$
 return reverse l$
end$

symbolic operator nodepnd$
symbolic procedure nodepnd(fl)$
for each f in cdr fl do
depl!*:=delete(assoc(reval f,depl!*),depl!*)$

symbolic procedure err_catch_readin$
begin scalar h,mode_bak,echo_bak$
 mode_bak:=!*mode; % as the _stop_ file has to start with 'lisp;'
 echo_bak:=!*echo; semic!*:='!$; 
 h:= errorset({'in,mkquote {"_stop_"}},nil,nil)
     where !*protfg=t;
 !*echo:=echo_bak; semic!*:='!; ; 
 erfg!*:=nil; !*mode:=mode_bak$
 return not errorp h 
end$

symbolic procedure err_catch_solve(eqs,fl)$
% fl='(LIST x y z);    eqs='(LIST expr1 expr2 .. )
begin scalar h$
 h:=errorset({'solveeval,mkquote{eqs, fl}},nil,nil)
 where !*protfg=t;
 erfg!*:=nil; 
 return if errorp h then nil
                    else cdar h    % cdr for deleting 'LIST
end$

symbolic operator err_catch_sub$
symbolic procedure err_catch_sub(h2,h6,h3)$
% do sub(h2=h6,h3) with error catching
begin scalar h4,h5;
 h4 := list('EQUAL,h2,h6);
 h5:=errorset({'subeval,mkquote{reval h4, 
                                reval h3 }},nil,nil)
     where !*protfg=t;
 erfg!*:=nil; 
 return if errorp h5 then nil
                     else car h5   
end$

symbolic operator err_catch_int$
symbolic procedure err_catch_int(h2,h3)$
% do int(h2,h3) with error catching
begin scalar h5;
 h5:=errorset({'simpint,mkquote{reval h2, 
                                reval h3 }},nil,nil)
     where !*protfg=t;
 erfg!*:=nil; 
 return if errorp h5 then nil
                     else 
 if not freeof(car h5,'INT) then nil
                            else prepsq car h5
end$

symbolic procedure beforegcsystemhook()$ 
my_gc_counter:=add1 my_gc_counter$

symbolic procedure aftergcsystemhook()$
if my_gc_counter > max_gc_counter then <<
 if print_ % and print_more (User must know that not all is computed.)
  then <<
  write "Stop of a subroutine."$terpri()$
  write "Number of garbage collections exceeds ",backup_,".";
  terpri()$
 >>$
 rederr "Heidadeife "
>>$

symbolic operator err_catch_fac$
symbolic procedure err_catch_fac(a)$
begin scalar h,bak,kernlist!*bak,kord!*bak,bakup_bak;
 bak:=max_gc_counter;
 max_gc_counter:=my_gc_counter+max_gc_fac;
 kernlist!*bak:=kernlist!*$
 kord!*bak:=kord!*$
 bakup_bak:=backup_;backup_:='max_gc_fac$
 h:=errorset({'reval,list('FACTORIZE,mkquote a)},nil,nil)
    where !*protfg=t;
 kernlist!*:=kernlist!*bak$
 kord!*:=kord!*bak;
 erfg!*:=nil; 
 max_gc_counter:=bak;
 backup_:=bakup_bak;
 return if errorp h then {'LIST,a}
                    else car h
end$

symbolic operator err_catch_gcd$
symbolic procedure err_catch_gcd(a,b)$
begin scalar h,bak,kernlist!*bak,kord!*bak,bakup_bak;
 bak:=max_gc_counter;
 max_gc_counter:=my_gc_counter+max_gc_fac;
 kernlist!*bak:=kernlist!*$
 kord!*bak:=kord!*$
 bakup_bak:=backup_;backup_:='max_gc_fac$
 h:=errorset({'reval,list('GCD,mkquote a,mkquote b)},nil,nil)
    where !*protfg=t;
 kernlist!*:=kernlist!*bak$
 kord!*:=kord!*bak;
 erfg!*:=nil; 
 max_gc_counter:=bak;
 backup_:=bakup_bak;
 return if errorp h then 1
                    else car h
end$

%symbolic operator err_catch_fac$
%symbolic procedure err_catch_fac(a)$
%begin scalar h,bak,bakup_bak;
% bak:=max_gc_counter;
% max_gc_counter:=my_gc_counter+max_gc_fac;
% bakup_bak:=backup_;backup_:='max_gc_fac$
% h:=errorset({'reval,list('FACTORIZE,mkquote a)},nil,nil)
%    where !*protfg=t;
% erfg!*:=nil;
% max_gc_counter:=bak;
% backup_:=bakup_bak;
% return if errorp h then {'LIST,a}
%                    else car h
%end$

symbolic procedure factored_form(a)$
% a is expected to be in prefix form
begin scalar b;
 if (pairp a) and (car a = 'PLUS) then <<
  b:=err_catch_fac a$
  if b and (length b > 2) then a:=cons('TIMES,cdr b)
 >>;
 return a
end$

symbolic lispeval '(putd 'countids 'expr
          '(lambda nil (prog (nn) (setq nn 0)
                (mapobl (function (lambda (x) (setq nn (plus2 nn 1)))))
                                  (return nn))))$

symbolic operator low_mem$
% if garbage collection recovers only 500000 cells then backtrace
% to be used only on workstations, not PCs i.e. under LINUX, Windows

symbolic procedure newreclaim()$
   <<oldreclaim();
     if (known!-free!-space() < 500000 ) then backtrace()
   >>$

symbolic procedure low_mem()$
if not( getd 'oldreclaim) then <<
    copyd('oldreclaim,'!%reclaim);
    copyd('!%reclaim,'newreclaim);
>>$

symbolic operator polyansatz$
symbolic procedure polyansatz(ev,iv,fn,degr,homo)$
% ev, iv are algebraic mode lists
% generates a polynomial in the variables ev of degree degr
% with functions of the variables iv
% if homo then a homogeneous polynomial
begin scalar a,fi,el1,el2,f,fl,p,pr;
 a:=reval list('EXPT,cons('PLUS,if homo then cdr ev
                                        else cons(1,cdr ev)),degr)$
 a:=reverse cdr a$
 fi:=0$
 iv:=cdr iv$
 for each el1 in a collect <<
  if (not pairp el1) or
     (car el1 neq 'TIMES) then el1:=list el1
                          else el1:=cdr el1;
  f:=newfct(fn,iv,fi);
  fi:=add1 fi;
  fl:=cons(f,fl)$
  pr:=list f$
  for each el2 in el1 do
  if not fixp el2 then pr:=cons(el2,pr);
  if length pr>1 then pr:=cons('TIMES,pr)
                 else pr:=car pr;
  p:=cons(pr,p)
 >>$
 p:=reval cons('PLUS,p)$
 return list('LIST,p,cons('LIST,fl))
end$

symbolic operator polyans$
symbolic procedure polyans(ordr,dgr,x,y,d_y,fn)$
% generates a polynom
% for i:=0:dgr sum fn"i"(x,y,d_y(1),..,d_y(ordr-1))*d_y(ordr)**i
% with fn as the function names and d_y as names or derivatives of y
% w.r.t. x
begin scalar ll,fl,a,i,f$
    i:=sub1 ordr$
    while i>0 do
	  <<ll:=cons(list(d_y,i),ll)$
	  i:=sub1 i>>$
    ll:=cons(y,ll)$
    ll:=reverse cons(x,ll)$
    fl:=nil$
    i:=0$
    while i<=dgr do
    <<f:=newfct(fn,ll,i)$
      fl:=(f . fl)$
      a:=list('PLUS,list('TIMES,f,list('EXPT,list(d_y,ordr),i)),a)$
      i:=add1 i>>$
    return list('list,reval a,cons('list,fl))
end$ % of polyans

symbolic operator sepans$
symbolic procedure sepans(kind,v1,v2,fn)$
% Generates a separation ansatz
% v1,v2 = lists of variables, fn = new function name + index added
% The first variable of v1 occurs only in one sort of the two sorts of
% functions and the remaining variables of v1 in the other sort of
% functios.
% The variables of v2 occur in all functions.
% Returned is a sum of products of each one function of both sorts.
% form: fn1(v11;v21,v22,v23,..)*fn2(v12,..,v1n;v21,v22,v23,..)+...
% the higher "kind", the more general and difficult the ansatz is
% kind = 0 is the full case
begin scalar n,vl1,vl2,h1,h2,h3,h4,fl$
  if cdr v1 = nil then <<vl1:=cdr v2$vl2:=cdr v2>>
		  else <<vl1:=cons(cadr v1,cdr v2)$
			 vl2:=append(cddr v1,cdr v2)>>$
  return
  if kind = 0 then <<vl1:=append(cdr v1,cdr v2)$
		     h1:=newfct(fn,vl1,'_)$
		     list('LIST,h1,list('LIST,h1))>>
  else
  if kind = 1 then <<h1:=newfct(fn,vl1,1)$
		     list('LIST,h1,list('LIST,h1))>>
  else
  if kind = 2 then <<h1:=newfct(fn,vl2,1)$
		     list('LIST,h1,list('LIST,h1))>>
  else
  if kind = 3 then <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     list('LIST,reval list('PLUS,h1,h2),
			  list('LIST,h1,h2))>>
  else
  if kind = 4 then <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     list('LIST,reval list('TIMES,h1,h2),
			  list('LIST,h1,h2))>>
  else
  if kind = 5 then <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     h3:=newfct(fn,vl1,3)$
		     list('LIST,reval list('PLUS,list('TIMES,h1,h2),h3),
			  list('LIST,h1,h2,h3))>>
  else
  if kind = 6 then <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     h3:=newfct(fn,vl2,3)$
		     list('LIST,reval list('PLUS,list('TIMES,h1,h2),h3),
			  list('LIST,h1,h2,h3))>>
  else
  if kind = 7 then <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     h3:=newfct(fn,vl1,3)$
		     h4:=newfct(fn,vl2,4)$
		     list('LIST,reval list('PLUS,
			  list('TIMES,h1,h2),h3,h4),
			  list('LIST,h1,h2,h3,h4))>>
  else
% ansatz of the form FN = FN1(v11,v2) + FN2(v12,v2) + ... + FNi(v1i,v2)
  if kind = 8 then <<n:=1$ vl1:=cdr v1$ vl2:=cdr v2$
		    fl:=()$
		     while vl1 neq () do <<
		       h1:=newfct(fn,cons(car vl1,vl2),n)$
		       vl1:=cdr vl1$
		       fl:=cons(h1, fl)$
		       n:=n+1 
		     >>$
		     list('LIST, cons('PLUS,fl), cons('LIST,fl))>>
		       

  else
		   <<h1:=newfct(fn,vl1,1)$
		     h2:=newfct(fn,vl2,2)$
		     h3:=newfct(fn,vl1,3)$
		     h4:=newfct(fn,vl2,4)$
		     list('LIST,reval list('PLUS,list('TIMES,h1,h2),
						 list('TIMES,h3,h4)),
			  list('LIST,h1,h2,h3,h4))>>
end$ % of sepans

%
% Orderings support!
%
% change_derivs_ordering(pdes,fl,vl) changes the ordering of the
% list of derivatives depending on the current ordering (this
% is detected "automatically" by sort_derivs using the lex_df flag to
% toggle between total-degree and lexicographic.
%
symbolic procedure change_derivs_ordering(pdes,fl,vl)$
begin scalar p, dl;
   for each p in pdes do <<
      if tr_orderings then <<
	 terpri()$
	 write "Old: ", get(p,'derivs)$
      >>$
      dl := sort_derivs(get(p,'derivs),fl,vl)$
      if tr_orderings then <<
	 terpri()$
	 write "New: ", dl$
      >>$
      put(p,'derivs,dl)$
      put(p,'dec_with,nil)$    % only if orderings are not
      % investigated in parallel (-->ord) 
      put(p,'dec_with_rl,nil)  % only if orderings are not ..
   >>$
   return pdes
end$

symbolic procedure sort_according_to(r,s)$
% all elements in r that are in s are sorted according to their order in s
begin scalar ss,h;
 for each ss in s do
 if member(ss,r) then h:=cons(ss,h);
 return reverse h
end$

symbolic procedure change_fcts_ordering(newli,pdes,vl)$
begin scalar s$
 ftem_ := newli$
 for each s in pdes do put(s,'fcts,sort_according_to(get(s,'fcts),ftem_))$
 pdes := change_derivs_ordering(pdes,ftem_,vl)$
 if tr_orderings then <<
  terpri()$
  write "New functions list: ", ftem_$
 >>
end$

symbolic procedure check_history(pdes)$
begin scalar p,q,h$
 for each p in pdes do <<
  h:=get(p,'histry_);
  for each q in pdes do
  h:=subst(get(q,'val),q,h)$
  if algebraic((lisp(get(p,'val)) - h) neq 0) then <<
   write"The history value of ",p," is not correct!"$
   terpri()
  >>
 >>
end$

symbolic procedure check_globals$
% to check validity of global variables at start of CRACK
begin scalar flag, var$

 % The integer variables
 foreach var in global_list_integer do
 if not fixp eval(var) then <<
  terpri()$
  write var, " needs to be an integer: ", eval(var)," is invalid"$
  flag := var
 >>$

 % Now for integer variables allowed to be nil
 foreach var in global_list_ninteger do 
 if not fixp eval(var) and eval(var) neq nil then <<
  terpri()$
  write var, " needs to be an integer or nil: ",
  eval(var)," is invalid"$
  flag := var
 >>$

 % Finally variables containing any number
 foreach var in global_list_number do
 if not numberp eval(var) then <<
  terpri()$
  write var, " needs to be a number: ", eval(var)," is invalid"$
  flag := var
 >>$

 return flag
end$

symbolic procedure search_li(l,care)$
% Find the cadr of all sublists which have 'care' as car (no nesting)
if pairp l then
if car l = care then {cadr l} 
                else begin
 scalar a,b,resul;
 for each a in l do
 if b:=search_li(a,care) then resul:=union(b,resul);
 return resul
end$

symbolic procedure search_li2(l,care)$
% Find all sublists which have 'care' as car (no nesting)
if pairp l then
if car l = care then list l
                else begin
 scalar a,b,resul;
 for each a in l do
 if b:=search_li2(a,care) then resul:=union(b,resul);
 return resul
end$

symbolic operator backup_reduce_flags$
symbolic procedure backup_reduce_flags$
% !*nopowers   = t  to have output of FACTORIZE like in Reduce 3.6
% !*allowdfint = t  moved here from crintfix, to enable simplification
%                   of derivatives of integrals
begin
 !*dfprint_bak   := cons(!*dfprint,!*dfprint_bak)$     
 !*exp_bak       := cons(!*exp,!*exp_bak)$             
 !*ezgcd_bak     := cons(!*ezgcd,!*ezgcd_bak)$         
 !*fullroots_bak := cons(!*fullroots,!*fullroots_bak)$ 
 !*gcd_bak       := cons(!*gcd,!*gcd_bak)$             
 !*mcd_bak       := cons(!*mcd,!*mcd_bak)$             
 !*ratarg_bak    := cons(!*ratarg,!*ratarg_bak)$       
 !*rational_bak  := cons(!*rational,!*rational_bak)$   

 if null !*dfprint   then algebraic(on  dfprint)$
 if null !*exp       then algebraic(on  exp)$    
 if null !*ezgcd     then algebraic(on  ezgcd)$  
 if null !*fullroots then algebraic(on  fullroots)$
 if      !*gcd       then algebraic(off gcd)$
 if null !*mcd       then algebraic(on  mcd)$
 if null !*ratarg    then algebraic(on  ratarg)$
 if null !*rational  then algebraic(on  rational)$

 !#if (neq version!* "REDUCE 3.6")
  !*nopowers_bak   := cons(!*nopowers,!*nopowers_bak)$     
  !*allowdfint_bak := cons(!*allowdfint,!*allowdfint_bak)$ 
  if null !*nopowers   then algebraic(on nopowers)$
  if null !*allowdfint then algebraic(on allowdfint)$
 !#endif

end$

symbolic operator recover_reduce_flags$
symbolic procedure recover_reduce_flags$
begin
 
 if !*dfprint neq car !*dfprint_bak then
 if !*dfprint then algebraic(off dfprint) else algebraic(on dfprint)$
 !*dfprint_bak:= cdr !*dfprint_bak$ 

 if !*exp neq car !*exp_bak then
 if !*exp then algebraic(off exp) else algebraic(on exp)$
 !*exp_bak:= cdr !*exp_bak$ 

 if !*ezgcd neq car !*ezgcd_bak then
 if !*ezgcd then algebraic(off ezgcd) else algebraic(on ezgcd)$
 !*ezgcd_bak:= cdr !*ezgcd_bak$ 

 if !*fullroots neq car !*fullroots_bak then
 if !*fullroots then algebraic(off fullroots) else algebraic(on fullroots)$
 !*fullroots_bak:= cdr !*fullroots_bak$ 

 if !*gcd neq car !*gcd_bak then
 if !*gcd then algebraic(off gcd) else algebraic(on gcd)$
 !*gcd_bak:= cdr !*gcd_bak$ 

 if !*mcd neq car !*mcd_bak then
 if !*mcd then algebraic(off mcd) else algebraic(on mcd)$
 !*mcd_bak:= cdr !*mcd_bak$ 

 if !*ratarg neq car !*ratarg_bak then
 if !*ratarg then algebraic(off ratarg) else algebraic(on ratarg)$
 !*ratarg_bak:= cdr !*ratarg_bak$ 

 if !*rational neq car !*rational_bak then
 if !*rational then algebraic(off rational) else algebraic(on rational)$
 !*rational_bak:= cdr !*rational_bak$ 

 !#if (neq version!* "REDUCE 3.6")
  if !*nopowers neq car !*nopowers_bak then
  if !*nopowers then algebraic(off nopowers) else algebraic(on nopowers)$
  !*nopowers_bak:= cdr !*nopowers_bak$ 
  if !*allowdfint neq car !*allowdfint_bak then
  if !*allowdfint then algebraic(off allowdfint) else algebraic(on allowdfint)$
  !*allowdfint_bak:= cdr !*allowdfint_bak$ 
 !#endif
end$

algebraic procedure maklist(ex)$
% making a list out of an expression if not already
if lisp(atom algebraic ex) then {ex} else
if lisp(car algebraic ex neq 'LIST) then ex:={ex}
                                    else ex$

symbolic procedure add_to_last_steps(h)$
if length last_steps < 20 then last_steps:=cons(h,last_steps) else
last_steps:=cons(h,reverse cdr reverse last_steps)$

symbolic procedure same_steps(a,b)$
if (car a = car b            ) and
   ((a = b) or
    ((car a neq 'subst) and
     (car a neq   27  ) and
     (car a neq   11  )     )) then t
                               else nil$

symbolic procedure in_cycle(h)$
begin scalar cpls1,cpls2,n,cycle;
 cpls1:=last_steps$
 if car h='subst then <<
  n:=0$
  while cpls1 do <<
   if same_steps(h,car cpls1) then n:=add1 n;
   cpls1:=cdr cpls1
  >>$
  cycle:=
  if n>2 then << % the subst. had been done already >=3 times
   write"A partial substitution has been repeated too often."$ terpri()$
   write"It will now be made rigorously."$ terpri()$
   t 
  >>     else nil
  % add_to_last_steps(h) is done outside for substitutions as it is not
  % clear at this stage whether the substitution will be performed
 >>                         else <<
  n:=1$
  while cpls1 and (not same_steps(h,car cpls1)) do 
  <<n:=add1 n;cpls1:=cdr cpls1>>$
  if null cpls1 or 
    ((reval {'PLUS,n,n})>length last_steps) then cycle:=nil 
                                            else <<
   cpls1:=cdr cpls1;  
   cpls2:=last_steps$
   while (n>0) and same_steps(car cpls2,car cpls1) do
   <<cpls1:=cdr cpls1;cpls2:=cdr cpls2;n:=sub1 n>>$
   if (n=0) and print_ then <<
     write if car h = 11 then "An algebraic length reduction" else
           if car h = 27 then "A length reducing simplification" else
                              "A step",
           " was prevented"$    terpri()$
     write"to avoid a cycle."$  terpri()$
   >>$
   cycle:=if n>0 then nil else t
  >>;
  if null cycle then add_to_last_steps(h)$
 >>;
 return cycle
end$

endmodule$

%********************************************************************
module solution_handling$
%********************************************************************
%  Routines for storing, retrieving, merging and displaying solutions
%  Author: Thomas Wolf Dec 2001

symbolic procedure save_solution(eqns,assigns,freef,ineq,file_name)$
% input data are algebraic mode
% eqns    .. list of remaining unsolved equations
% assigns .. list of computed assignments of the form `function = expression'
% freef   .. list of list of functiones either free or in eqns
% ineq    .. list of inequalities
begin scalar s,levelcp,h,p,conti$
  if file_name then s:=file_name
               else <<
   s:=session_;
   levelcp:=reverse level_;
   while levelcp do <<setq(s,bldmsg("%w%d.",s,car levelcp));
                      levelcp:=cdr levelcp>>;
   s:=explode s;
   s:=compress cons(car s,cons('s,cons('o,cdddr s)))$
  >>$

  sol_list:=union(list s,sol_list)$

  out s;
  write"off echo$ "$
  write"backup_:='("$terpri()$

  for each h in freef do 
  if p:=assoc(h,depl!*) then conti:=cons(p,conti);

  % The first sub-list is a  list of dependencies, like ((f x y) (g x))
  write"% A list of dependencies, like ((f x y) (g x))"$terpri()$
  print conti$write" "$terpri()$

  % The next sublist is a list of unsolved equations
  write"% A list of unsolved equations"$terpri()$
  print eqns$write" "$terpri()$
  
  % The next sublist is a list of assignments
  write"% A list of assignments"$terpri()$
  print assigns$write" "$terpri()$
  
  % The next sublist is a list of free or unresolved functions
  write"% A list of free or unresolved functions"$terpri()$
  print freef$write" "$terpri()$
  
  % The next sublist is a list of non-vanishing expressions
  write"% A list of non-vanishing expressions"$terpri()$
  print ineq$terpri()$
  
  write")$"$

  % shorter but less clear:  print list(conti,eqns,freef,ineq)$

  write "end$"$terpri()$
  shut s;

  return s
end$

symbolic procedure print_indexed_list(li)$
begin scalar a,h$
 terpri()$
 h:=0$
 for each a in li do <<
  h:=add1 h;
  write"[",h,"]";terpri()$
  mathprint a
 >>
end$

symbolic procedure merge_two(s1,sol1,s2,sol2,absorb)$
% Is sol1 a special case of sol2 ?
% If yes, return the new generalized solution sol2 with one less inequality.
% If absorb then modify s2 and sol2 if s1 can be absorbed

begin scalar eli_2,singular_eli,regular_eli,a,b,cond2,sb,remain_sb,
             singular_sb,regular_sb,c2,remain_c2,remain_num_c2,h,
             try_to_sub,try_to_sub_cp,num_sb,num_sb_quo,singular_ex,
             singular_ex_cp,ineq2,ine,ineqnew,ineqdrop,tr_merge,
             extra_par_in_s1,gauge_of_s2,gauge_of_s2_cp,did_trafo,n,
             remain_c2_cp,dropped_assign_in_s2,new_assign_in_s2$

 % tr_merge:=t$
 if tr_merge then <<write"s1=",s1$terpri()$
                    write"s2=",s2$terpri()$
                    write"*** sol1 ***:"$print_indexed_list(caddr sol1)$
                    write"*** sol2 ***:"$print_indexed_list(caddr sol2)$
                    write"free param in sol1: ",cadddr sol1$terpri()$
                    write"free param in sol2: ",cadddr sol2$terpri()>>$
 % At first we list all lhs y_i in assignments y_i=... in sol2
 eli_2:=setdiff(for each a in caddr sol2 collect cadr a,cadddr sol2); 

 % We use setdiff because of assignments, like a6=a6 in sol2 
 % where a6 is a free parameter.

 % writing assignments of solution 2 as expressions to vanish
 cond2:=for each a in caddr sol2 
        collect {'PLUS,cadr a,{'MINUS,caddr a}};

 % Do all substitutions a=... from sol1 for which 
 % there is an assignment a=... in sol2 and 
 % collecting the others as remain_sb

 cond2:=cons('LIST,cond2);
 sb:=caddr sol1; % all assignments of solution 1
 while sb do <<
  a:=car sb; sb:=cdr sb;
  if member(cadr a,eli_2) then <<
   eli_2:=delete(cadr a,eli_2)$
   cond2:=err_catch_sub(cadr a,caddr a,cond2)
  >>                      else remain_sb:=cons(a,remain_sb)
 >>$
 eli_2:=append(eli_2,cadddr sol2)$ 
 % eli_2 is now the list of new sol2 parameters

 % At this stage any sol2 conditions either become singular or zero.
 % The singular ones are collected in remain_c2.

 % The same again, now taking only numerators
 remain_c2:=cond2;  

 cond2:=cdr cond2;
 c2:=nil$
 h:=0$
 while cond2 and (null c2 or zerop c2) do <<
  c2:=car cond2;
  h:=add1 h;
  if tr_merge then <<write"[",h,"]"$terpri()$mathprint c2>>$

  % Is the numerator of c2 fulfilled by assignments of solution 1?
  sb:=remain_sb;               % all remaining assignments of solution 1
  while sb and c2 and not zerop c2 do <<
   a:=car sb; sb:=cdr sb;
   c2:=algebraic(num(lisp(c2)));
   if tr_merge then b:=c2;
   c2:=err_catch_sub(cadr a,caddr a,c2);
   if tr_merge and (b neq c2) then <<
    write"Sub: ";mathprint a$
    if c2 then <<write"c2="$mathprint c2>>
          else <<write"singular result"$terpri()>>
   >>
  >>$
  if null c2 then remain_num_c2:=cons(car cond2,remain_num_c2);
  cond2:=cdr cond2
 >>$

 if c2 and not zerop c2 then return nil; % sol1 is not special case of sol2
 if remain_num_c2 then << % can only occur if there were singular subst.
  write"Even substitutions in the numerator is giving "$terpri()$
  write"singularities like for log(0)."$                terpri()$
  return nil
 >>$

 write"Substitutions in numerators give all zero"$terpri()$

 % At this stage in the program either all is satisfied (remain_c2 = nil)
 % or only singular solution 2 assignments remain (remain_c2 <> nil)
 % but which vanish if numerators are taken.

 % We now want to find a different order of substitutions, especially
 % substituting for the free parameter functions of solution 2
 % based on remain_sb to be done in remain_c2.

 % At first we sort all solution 1 assignments into regular_sb and singular_sb.
 % remain_c2 is not changed in this
 cond2:=remain_c2;
 sb:=remain_sb;               % all remaining assignments of solution 1
 while sb do <<
  a:=car sb; sb:=cdr sb;
  h:=err_catch_sub(cadr a,caddr a,cond2);
  if null h then singular_sb:=cons(a,singular_sb)
            else regular_sb:=cons(a,regular_sb)
 >>$
 if tr_merge then <<terpri()$
                    write"regular_sb: "$mathprint cons('LIST,regular_sb)>>$
 if tr_merge then <<write"singular_sb: "$mathprint cons('LIST,singular_sb)>>$

 if singular_sb then <<
  write"Substitutions lead to singularities."$terpri()$
  write"Solution ",s2," has to be transformed."$terpri()
 >>$
   
 % We now make a list of vanishing expressions based on singular_sb
 % which when replaced by 0 in remain_c2 give singularities
 singular_ex:=for each a in singular_sb
              collect reval {'PLUS,cadr a,{'MINUS,caddr a}};
 if tr_merge then <<
  write"The following are expressions which vanish due to sol1 and"$
  terpri()$
  write"which lead to singularities when used for substitutions in sol2"$
  terpri()$
  mathprint cons('LIST,singular_ex)
 >>$

 if tr_merge then <<
  write"The following are all free parameters in sol2 for which there are"$
  terpri()$
  write"substitutions in sol1"$ terpri()$
 >>$
 singular_eli:=for each a in singular_sb collect cadr a; 
 regular_eli:=for each a in regular_sb collect cadr a; 
 if tr_merge then <<terpri()$
                    write"singular_eli: "$mathprint cons('LIST,singular_eli)>>;
 if tr_merge then <<write"regular_eli: "$mathprint cons('LIST,regular_eli)>>;

 % Before continuing we want to check whether the supposed to be more special
 % solution sol1 has free parameters which are not free parameters in the more
 % general solution sol2. That can cause problems, i.e. division through 0
 % and non-includedness when in fact sol1 is included in sol2.

 extra_par_in_s1:=setdiff(cadddr sol1,cadddr sol2);
 if tr_merge then <<write"Param in sol1 and not in sol2: ",extra_par_in_s1;
                    terpri()>>$

 for each a in extra_par_in_s1 do <<
  h:=caddr sol2$
  while h and cadar h neq a do h:=cdr h;
  if null h then write"ERROR, there must be an assignment of a in sol2!"
            else <<
   if tr_merge then <<
    write"Assignment in ",s2," of a variable that is a free parameter in ",
         s1," :"$
    terpri()$
    mathprint car h$
   >>$
   dropped_assign_in_s2:=cons(car h,dropped_assign_in_s2);
   gauge_of_s2:=cons(algebraic(num(lisp({'PLUS,cadr car h,
                                         {'MINUS,caddr car h}}))),
                     gauge_of_s2)
  >>
 >>$

 gauge_of_s2:=cons('LIST,gauge_of_s2);

 if tr_merge then <<write"gauge_of_s2="$mathprint gauge_of_s2>>$

 % We should not do all regular substitutions in gauge_of_s2 (tried that)
 % because some of them may set variables to zero which limits the
 % possibilities of doing transformations of remain_c2 

 % We now search for a substitution based on one of the equations
 % gauge_of_s2. The substitution is to be performed on remain_c2.

 % One sometimes has to solve for regular_eli as singular_eli
 % might appear only non-linearly.
 % try_to_sub:=append(regular_eli,singular_eli);
 try_to_sub:=append(singular_eli,regular_eli);
 n:=1;
 repeat <<
  did_trafo:=nil;
  gauge_of_s2_cp:=cdr gauge_of_s2;
  while gauge_of_s2_cp do <<
   sb:=reval car gauge_of_s2_cp$ 
   gauge_of_s2_cp:=cdr gauge_of_s2_cp$
   if not zerop sb then <<
    try_to_sub_cp:=try_to_sub;
    if tr_merge then <<write"next relation to be used: 0="$mathprint sb$
                       write"try_to_sub=",try_to_sub$terpri()>>$
    h:=err_catch_fac(sb);
    if h then <<
     sb:=nil;
     h:=cdr h;
     while h do <<
      if pairp cadar h then sb:=cons(cadar h,sb);
      h:=cdr h;
     >>
    >>$

    % From the next condition 0=sb we drop all factors which are 
    % single variables which set to zero would be a limitation
    if tr_merge then <<write"After dropping single variable factors ",
                             length sb," factor(s) remain"$terpri()>>$
    sb:=reval cons('TIMES,cons(1,sb));
    if tr_merge then <<write"New relation used for substitution: sb="$
                       mathprint sb$terpri()>>$
    
    while try_to_sub_cp do <<
     a:=car try_to_sub_cp; try_to_sub_cp:=cdr try_to_sub_cp;
     if tr_merge then <<write"try to sub next: ",a$terpri()>>$
     if not freeof(sb,a) and lin_check(sb,{a}) then <<
      num_sb:=reval {'DIFFERENCE, sb,{'TIMES,a,coeffn(sb,a,1)}};
      if tr_merge then <<write"num_sb="$mathprint num_sb>>$
      singular_ex_cp:=singular_ex;
      while singular_ex_cp do << 
       if tr_merge then <<write"car singular_ex_cp=",car singular_ex_cp$
                          terpri()>>$
       % Search for an expression causing a singular substitution
       % which is a factor of the substituted expression for a
       num_sb_quo:=reval {'QUOTIENT,num_sb,car singular_ex_cp};
       if tr_merge then <<write"num_sb_quo="$mathprint num_sb_quo>>$
       % if (not pairp num_sb_quo) or 
       %    (car num_sb_quo neq 'QUOTIENT) then <<
       if t then <<
        eli_2:=delete(a,eli_2);
        % i.e. num_sb is a multiple of one of members of singular_ex, HURRAY!
        % Do the substitution in remain_c2
        b:=cadr solveeval list(sb,a)$
        h:=err_catch_sub(cadr b,caddr b,remain_c2);
        if tr_merge and null h then <<
         write"Trafo "$mathprint b$write" was singular."$ terpri()
        >>$
        if h then <<
         gauge_of_s2:=err_catch_sub(cadr b,caddr b,gauge_of_s2);
         gauge_of_s2:=cons('LIST, for each gauge_of_s2_cp in cdr gauge_of_s2
                           collect algebraic(num(lisp(gauge_of_s2_cp))));
         gauge_of_s2_cp:=nil$
         new_assign_in_s2:=cons(b,new_assign_in_s2);
         did_trafo:=t$
         write"In order to avoid a singularity when doing substitutions"$
         terpri()$
         write"the supposed to be more general solution was transformed using:"$
         terpri()$
         mathprint b$
         if tr_merge then <<write"The new gauge_of_s2: "$
                            mathprint gauge_of_s2>>$

         remain_c2:=h;
         h:=append(regular_sb,singular_sb);
         while h and a neq cadar h do h:=cdr h;
         if h then remain_c2:=append(remain_c2,list {'DIFFERENCE,caddar h,caddr b});
	 if tr_merge then <<write"remain_c2="$print_indexed_list(cdr remain_c2)>>$
         singular_ex_cp:=nil;
         try_to_sub:=delete(a,try_to_sub);
         try_to_sub_cp:=nil;
         n:=n+1
        >>   else singular_ex_cp:=cdr singular_ex_cp
       >>                                else singular_ex_cp:=cdr singular_ex_cp
      >>    % while singular_ex_cp
     >>    % if car try_to_sub_cp passes first test
    >>$   % while try_to_sub_cp
   >>    % if not zerop sb
  >>$   % while gauge_of_s2_cp
 >> until (did_trafo=nil)$  

 if tr_merge then <<
  write"After completing the trafo the new list of parameters of"$
  terpri()$      
  write"sol2 is: ",eli_2$terpri()$      
  write"sol1 has free parameters: ",cadddr sol1$terpri()
 >>$

 if not_included(cadddr sol1,eli_2) then return <<
  write"Something seems wrong in merge_sol(): after the transformation of"$
  terpri()$
  write"sol2, all free parameters of sol1 should be free parameters of sol2."$
  terpri();
  nil
 >>                                else <<
  if tr_merge then <<
   write"All free parameters of sol1 are free parameters of sol2"$
   terpri()
  >>
 >>$


 % Now all in remain_c2 has to become zero by using first substitutions
 % from regular_sb and substituting parameters from sol2 such that
 % the substituted expression has one of the singular_ex as factor.

 % We seek global substitutions, i.e. substitutions based on sol1
 % which satisfy all sol2 conditions and not for each sol2 condition a
 % different set of sol1 based substitutions. Therefore substitutions
 % are done in the whole remain_c2.

 % try_to_sub are free parameters in sol2 that are contained in
 % regular_eli and which are therefore not in singular_eli and not free
 % parameters in sol1. They are to be substituted next because sol1 is
 % obviously singularity free, so we have to express sol2 in the same
 % free parameters, so we have to substitute for the free parameters fo
 % sol2 which are not free parameters of sol1. But we must not use the
 % same substitutions regular_sb which substitute for them as they lead
 % to singular substitutions afterwards.

% try_to_sub:=memberl(cadddr sol2,regular_eli);
%
% write"try_to_sub=",try_to_sub$terpri()$
%
% % We now search for a substitution in regular_sb which leads to a
% % substitution of a member of try_to_sub, say p, ...
% b:=regular_sb;
% for each sb in b do <<
%  sb_cp:=algebraic(num(lisp({'PLUS,cadr sb,{'MINUS,caddr sb}})));
%  try_to_sub_cp:=delete(cadr sb,try_to_sub); % ... but the substitution
%                                             % does not originally 
%                                             % have the form p=...  .
%  while try_to_sub_cp do <<
%   a:=car try_to_sub_cp; try_to_sub_cp:=cdr try_to_sub_cp;
%   if not freeof(sb_cp,a) and lin_check(sb_cp,{a}) then <<
%    num_sb:={'DIFFERENCE, sb_cp,{'TIMES,a,coeffn(sb_cp,a,1)}};
%
%    singular_ex_cp:=singular_ex;
%    while singular_ex_cp do << 
%     % Search for an expression causing a singular substitution
%     % which is a factor of the substituted expression for a
%     num_sb_quo:=reval {'QUOTIENT,num_sb,car singular_ex_cp};
%     if (not pairp num_sb_quo) or 
%        (car num_sb_quo neq 'QUOTIENT) then <<
%      % i.e. num_sb is a multiple of one of members of singular_ex, HURRAY!
%      % Do the substitution in remain_c2
%      h:=err_catch_sub(cadr sb,caddr sb,remain_c2);
%      if h then <<
%       write"In order to avoid a singularity when doing substitutions"$
%       terpri()$
%       write"the supposed to be more general solution was transformed:"$
%       terpri()$
%       mathprint sb$
%       remain_c2:=h;
%       singular_ex_cp:=nil;
%       regular_sb:=delete(sb,regular_sb);
%       try_to_sub:=delete(a,try_to_sub);
%       try_to_sub_cp:=nil;
%      >>   else singular_ex_cp:=cdr singular_ex_cp
%     >>                                else singular_ex_cp:=cdr singular_ex_cp
%    >>    % while singular_ex_cp
%   >>    % if car try_to_sub_cp passes first test
%  >>$   % while try_to_sub_cp
% >>$   % for each sb

 % Do the remaining regular_sb
 sb:=append(regular_sb,singular_sb); % all remaining assignments of solution 1
 while sb and remain_c2 do <<
  a:=car sb; sb:=cdr sb;
  remain_c2_cp:=remain_c2$
  remain_c2:=err_catch_sub(cadr a,caddr a,remain_c2);
  if tr_merge then 
  if null remain_c2 then 
  <<write"The following subst. was singular: "$mathprint a>>
                    else <<
   write"Remaining substitution: ";mathprint a$
   %write"remain_c2="$mathprint remain_c2
  >>
 >>$

 if null remain_c2 then remain_c2:=remain_c2_cp
                   else remain_c2_cp:=remain_c2;

 % Drop all zeros.
 remain_c2_cp:=cdr remain_c2_cp$
 while remain_c2_cp and zerop car remain_c2_cp do 
 remain_c2_cp:=cdr remain_c2_cp;

 if remain_c2_cp then <<     % s1 is NOT a special case of s2

  remain_c2_cp:=remain_c2$
  if tr_merge then <<write"remain_c2="$
                     print_indexed_list(cdr remain_c2_cp)>>$

  % Is there a contradiction of the type that the equivalence of two
  % assignments, a8=A (from sol1), a8=B (from sol2) requires 0=A-B
  %  which got transformed into an expression C which is built only
  % from free parameters of sol1 and therefore should not vanish?

  h:=cadddr sol1; % all free parameters in sol1
  while h and <<
   if tr_merge then write"Substitution of ",car h," by: "$
   repeat <<     % find a random integer for the free parameter
    a:=1+random(10000);   % that gives a regular substitution
    if tr_merge then <<write a$terpri()>>$
    a:=err_catch_sub(car h,a,remain_c2_cp)
   >> until a;
   remain_c2_cp:=a;
   while a and ((not numberp car a) or (zerop car a)) do a:=cdr a;
   not a
  >> do h:=cdr h;

  if h then return <<
   write"In the following S1 stands for ",s1,"and S2 stands for ",s2," . ",
        "Solution S1 fulfills all conditions of solution S2 when conditions",
        "are made denominator free. But, after rewriting solution S2 so that",
        "all free parameters of solution S1 are also free parameters of S2",
        "then the new solution S2 now requires the vanishing of an expression",
        "in these free parameters which is not allowed by S1. Therefore S1",
        "is not a special case of S2."$
   nil
  >>$

  if tr_merge and remain_c2_cp then 
  <<write"remain_c2_cp after subst = "$mathprint cons('LIST,remain_c2)>>$
  write"Solution ",s1," is not less restrictive than solution"$terpri()$
  write s2," and fulfills all conditions of solution ",s2," ."$terpri()$
  write"But it was not possible for the program to re-formulate solution "$
  terpri()$ write s2," to include both solutions in a single set of"$terpri()$
  write"assignments without vanishing denominators. :-( "$
  terpri()$
  return nil

 >>              else return <<  % return the new s2 as s1 IS a special case of s2

  % Which inequality is to be dropped?
  ineq2:=car cddddr sol2$

  while ineq2 do <<
   ine:=car ineq2;
   % ine should not have denominators, so no extra precautions for substitution:
   for each a in caddr sol1 do ine:=reval(subst(caddr a,cadr a,ine));
   if not zerop reval ine then ineqnew:=cons(car ineq2,ineqnew)
                          else ineqdrop:=cons(car ineq2,ineqdrop)$
   ineq2:=cdr ineq2
  >>$

  if absorb then <<
   % delete the redundant solution
   sol_list:=delete(s1,sol_list); % system bldmsg ("rm %s",s1);

   % transform the general solution if that was necessary and
   % updating the list of free parameters
   h:=cons('LIST,caddr sol2);
   b:=cadddr sol2;
   if tr_merge then <<
    write"h0="$print_indexed_list(h)$
    write"dropped_assign_in_s2="$print_indexed_list(dropped_assign_in_s2)$
    write"new_assign_in_s2="$print_indexed_list(new_assign_in_s2)$
   >>$
   for each a in dropped_assign_in_s2 do 
   <<h:=delete(a,h);b:=cons(reval cadr a,b)>>$
   if tr_merge then <<write"h1="$print_indexed_list(h)>>$
   for each a in reverse new_assign_in_s2 do if h then <<
    b:=delete(reval cadr a,b)$
    if tr_merge then <<write"a=",a$terpri()$write"h2="$print_indexed_list(h)>>$
    h:=err_catch_sub(cadr a,caddr a,h);
    if h then h:=reval append(h,list a)
   >>$
   if null h then 
   write"A seemingly successful transformation of ",s2,
        "went singular when performing the transformation ",
        "finally on the whole solution."
             else     % save the generalized solution
   save_solution(cadr sol2,cdr h,b,ineqnew,s2)$
  >>;

  if absorb and null h then nil 
                       else <<
   % report the merging
   if null ineqdrop then <<   
    write"Strange: merging ",s1," and ",s2," without dropping inequalities!"$ 
    terpri()$
    write"Probably ",s2," had already been merged with ",s1,
         " or similar before."$ terpri()
   >>               else
   if print_ then <<
    write"Solution ",s2," includes ",s1," by dropping "$
    if length ineqdrop = 1 then write"inequality"
                           else write"inequalities"$terpri()$
    for each ine in ineqdrop do mathprint ine
   >>;
   s2 % the more general solution
  >>
 >>
end$

symbolic operator merge_sol$
symbolic procedure merge_sol$
begin scalar s,sol_cp,sl1,sl2,s1,s2,s3,sol1,sol2,echo_bak$

 if null session_ then ask_for_session() else <<
  write "Do you want to merge solutions computed in this session,"$
  terpri()$
  if not yesp "i.e. since loading CRACK the last time? " then 
  ask_for_session()
 >>$

 % reading in sol_list
 setq(s,bldmsg("%w%w",session_,"sol_list"));
 in s; 

 % % At fist sort sol_list by the number of free unknowns
 % for each s1 in sol_list do <<
 %  in s1; 
 %  s2:=if null cadddr backup_ then 0 else length cadddr backup_;
 %  if cadr backup_ then s2:=s2 - length cadr backup_;
 %  sol_cp:=cons((s2 . s1),sol_cp)
 % >>$
 % sol_cp:=idx_sort(sol_cp)$
 % while sol_cp do <<sl1:=cons(cdar sol_cp,sl1);sol_cp:=cdr sol_cp>>$

 sol_cp:=sol_list$
 sl1:=sol_cp$

 if sl1 then
 while sl1 and cdr sl1 do <<
  s1:=car sl1; sl1:=cdr sl1;
  %infile s1; 
  echo_bak:=!*echo; semic!*:='!$; in s1$ !*echo:=echo_bak;
  sol1:=backup_;     
  if print_ then <<write"Comparing ",s1," with:"$terpri()>>$

  sl2:=sl1;
  while sl2 do <<
   s2:=car sl2; sl2:=cdr sl2;
   %infile s2$ 
   echo_bak:=!*echo; semic!*:='!$; in s2$ !*echo:=echo_bak;
   sol2:=backup_;    
   if print_ then <<write"  ",s2$terpri()>>$

   if (null car sol1) and (null car sol2) then % algebraic problem
   if length cadddr sol1 < 
      length cadddr sol2 then s3:=merge_two(s1,sol1,s2,sol2,t) else
   if length cadddr sol1 > 
      length cadddr sol2 then s3:=merge_two(s2,sol2,s1,sol1,t) else
   <<s3:=merge_two(s1,sol1,s2,sol2,t)$                             
    if s3 then <<                                                  
     write"Strange: ",s1," is contained in ",s2$terpri()$          
     write"but both have same number of free unknowns!"$terpri()$  
     write"One of them has probably undergone earlier merging"$
     terpri()$
    >>                                                             
   >>                                     else 
   if null (s3:=merge_two(s1,sol1,s2,sol2,t)) then 
   s3:=merge_two(s2,sol2,s1,sol1,t);
   if s3=s1 then sl1:=delete(s2,sl1) else % not to pair s2 later
   if s3=s2 then sl2:=nil                 % to continue with next element in sl1
  >>
 >>;

 save_sol_list()
end$

symbolic procedure save_sol_list$
begin scalar s$
 setq(s,bldmsg("%w%w",session_,"sol_list"));
 out s;
 write"off echo$ "$ terpri()$
 write"sol_list:='"$
 print sol_list$write"$"$terpri()$
 write "end$"$terpri()$
 shut s
end$

symbolic procedure ask_for_session$
begin scalar ps$
 ps:=promptstring!*$
 promptstring!*:="Name of the session in double quotes: "$ 
 terpri()$session_:=termread()$
 promptstring!*:=ps
end$

symbolic operator pri_sol$
symbolic procedure pri_sol(sin,assgn,crout,html,solcount,fname,prind)$
% print the single solution sin
begin scalar a,b,sout$
 in sin$

 if html then <<
  setq(sout,bldmsg("%w%w%d%w",fname,"-s",solcount,".html"));
  out sout;

  write"<html>"$terpri()$
  terpri()$
  write"<head>"$terpri()$
  write"<meta http-equiv=""Content-Type"" content=""text/html;"$terpri()$
  write"charset=iso-8859-1"">"$terpri()$
  write"<title>Solution ",solcount," to problem ",prind,"</title>"$terpri()$
  write"</head>"$terpri()$
  terpri()$
  write"<BODY TEXT=""#000000"" BGCOLOR=""#FFFFFF"">"$terpri()$
  terpri()$
  write"<CENTER><H2>Solution ",solcount," to problem ",prind,"</H2>"$terpri()$
  write"<HR>"$terpri()$
  if cadr backup_ then <<write"<A HREF=""#1"">Remaining equations</A> | "$
                         terpri()>>$
  write"<A HREF=""#2"">Expressions</A> | "$terpri()$
  write"<A HREF=""#3"">Parameters</A> | "$terpri()$
  write"<A HREF=""#4"">Relevance</A> | "$terpri()$
  write"<A HREF=",prind,".html>Back to problem ",prind,"</A> "$
  write"</CENTER>"$terpri()$
  terpri()
 >>$
 for each a in car backup_ do
 for each b in cdr a do
 algebraic(depend(lisp(car a),lisp b));
 backup_:=cdr backup_;
 terpri()$
 if html then write"<!-- "$
 write">>>=======>>> SOLUTION ",sin," <<<=======<<<"$
 if html then write" --> "$ 
 terpri()$terpri()$

 if assgn or html then <<

  if car backup_ then <<
   if html then <<
    write"<HR><A NAME=""1""></A><H3>Equations</H3>"$terpri()$
    write"The following unsolved equations remain:"$terpri()$
    write"<pre>"$
   >>      else write"Equations:"$
   for each a in car backup_ do mathprint {'EQUAL,0,a}$
   if html then <<write"</pre>"$terpri()>>
  >>$
  if html then <<
   write"<HR><A NAME=""2""></A><H3>Expressions</H3>"$terpri()$
   write"The solution is given through the following expressions:"$terpri()$
   write"<pre>"$terpri()$
   for each a in cadr backup_ do mathprint a$
   write"</pre>"$terpri()
  >>      else <<
   b:=nil;
   for each a in cadr backup_ do
   if not zerop caddr a then b:=cons(a,b);
   print_fcts(b,nil)$
  >>$
  terpri()$

  if html then <<
   write"<HR><A NAME=""3""></A><H3>Parameters</H3>"$terpri()$
   write"Apart from the condition that they must not vanish to give"$terpri()$
   write"a non-trivial solution and a non-singular solution with"$terpri()$
   write"non-vanishing denominators, the following parameters are free:"$terpri()$
   write"<pre> "$
   fctprint caddr backup_;
   write"</pre>"$terpri()
  >>      else <<
   write length caddr backup_," free unknowns: "$ listprint caddr backup_;
   print_ineq(cadddr backup_)$
  >>
 >>$

 if html then <<
  write"<HR><A NAME=""4""></A><H3>Relevance for the application:</H3>"$
  terpri()$
  % A text for the relevance should be generated in crack_out()
%  write"The solution given above tells us that the system {u_s, v_s}"$
%  terpri()$
%  write"is a higher order symmetry for the lower order system {u_t,v_t}"$
%  terpri()$
%  write"where u=u(t,x) is a scalar function, v=v(t,x) is a vector"$
%  terpri()$
%  write"function of arbitrary dimension and f(..,..) is the scalar"$
%  terpri()$
%  write"product between two such vectors:"$
%  terpri()$
  write"<pre>"
 >>$
 if crout or html then <<
  algebraic (
  crack_out(lisp cons('LIST,car backup_),
           lisp cons('LIST,cadr backup_),
           lisp cons('LIST,caddr backup_),
           lisp cons('LIST,cadddr backup_) ))$
 >>$
 if html then <<
  write"</pre>"$terpri()$
  write"<HR>"$terpri()$
  write"</body>"$terpri()$
  write"</html>"$terpri()$
  shut sout
 >>
end$

symbolic operator print_all_sol$
symbolic procedure print_all_sol$
begin scalar s,a,assgn,crout,natbak,print_more_bak,fname,solcount,
             html,prind,ps$

 write"This is a reminder for you to read in any file CRACK_OUT.RED"$
 terpri()$
 write"with a procedure CRACK_OUT() in case that is necessary to display"$
 terpri()$
 write"results following from solutions to be printed."$
 terpri()$ terpri()$

 if null session_ then ask_for_session() else <<
  write "Do you want to print solutions computed in this session,"$
  terpri()$
  if not yesp "i.e. since loading CRACK the last time? " then 
  ask_for_session()$
  terpri()
 >>$

 % reading in sol_list
 setq(s,bldmsg("%w%w",session_,"sol_list"));
 in s; 

 natbak:=!*nat$ print_more_bak:=print_more$ print_more:=t$
 if yesp "Do you want to generate an html file for each solution? "
 then <<  
  html:=t$
  solcount:=0$
  terpri()$
  write "What is the file name (including the path)"$
  terpri()$
  write "that shall be used (in double quotes) ? "$
  terpri()$
  write "(A suffix '-si'  will be added for each solution 'i'.) "$
  ps:=promptstring!*$
  promptstring!*:=""$
  fname:=termread()$terpri()$
  write "What is a short name for the problem? "$
  prind:=termread()$
  promptstring!*:=ps$
  terpri()$
 >> else <<
  if yesp "Do you want to see the computed value of each function? " 
  then assgn:=t$
  if yesp "Do you want procedure `crack_out' to be called? " then <<
   crout:=t;
   if flin_ and homogen_ then
   if yesp "Do you want to print less (e.g. no symmetries)? " 
   then print_more:=nil$
   if not yesp 
   "Do you want natural output (no if you want to paste and copy)? " 
   then !*nat:=nil$
  >>$
 >>$
 for each a in sol_list do <<
  if html then solcount:=add1 solcount$
  pri_sol(a,assgn,crout,html,solcount,fname,prind)$
 >>$
 !*nat:=natbak;
 print_more:=print_more_bak
end$

symbolic procedure sol_in_list(set1,set2,sol_list2)$
begin scalar set2cp,s1,s2,found,sol1,sol2,same_sets,echo_bak$
 while set1 do <<
  s1:=car set1; set1:=cdr set1;
  %infile s1; 
  echo_bak:=!*echo; semic!*:='!$; in s1$ !*echo:=echo_bak;
  sol1:=backup_;     
  set2cp:=set2$
  found:=nil$
  while set2cp and not found do <<
   s2:=car set2cp; set2cp:=cdr set2cp;
   %infile s2; 
   echo_bak:=!*echo; semic!*:='!$; in s2$ !*echo:=echo_bak;
   sol2:=backup_;     
   found:=merge_two(s1,sol1,s2,sol2,nil)$
  >>;
  if not found then <<
   same_sets:=nil;
   if print_ then <<
    write"Solution ",s1," is not included in ",sol_list2$
    terpri()
   >>
  >>
 >>$
 return same_sets
end$

symbolic operator same_sol_sets$
symbolic procedure same_sol_sets$
begin scalar session_bak,set1,set2,sol_list1,sol_list2,echo_bak$
 session_bak:=session_;
 write"Two sets of solutions are compared whether they are identical."$

 write"What is the name of the session that produced the first set of solutions?"$
 terpri()$
 write"(CRACK will look for the file `sessionname'+`sol_list'.)"$terpri()$ 
 ask_for_session()$

 % reading in sol_list
 setq(sol_list1,bldmsg("%w%w",session_,"sol_list"));
 %infile sol_list1; 
 echo_bak:=!*echo; semic!*:='!$; in sol_list1$ !*echo:=echo_bak;
 set1:=sol_list$

 write"What is the name of the session that produced the second set of solutions?"$
 terpri()$
 ask_for_session()$

 % reading in sol_list
 setq(sol_list2,bldmsg("%w%w",session_,"sol_list"));
 %infile sol_list2; 
 echo_bak:=!*echo; semic!*:='!$; in sol_list2$ !*echo:=echo_bak;
 set2:=sol_list$

 session_:=session_bak$
 
 % 1. Check that all solutions in set1 are included in set2.
 
 sol_in_list(set1,set2,sol_list2)$
 sol_in_list(set2,set1,sol_list1)$

end$

% some PSL specific hacks

!#if (memq 'psl lispsystem!*)
  symbolic procedure delete!-file(fi);
   if memq('unix,lispsystem!*) then
	system bldmsg("rm '%s'",fi) else
	system bldmsg("del ""%s""",fi);
 !#endif

endmodule$

end$