reduce-historical/r37/packages/crack/crsimp.red

%********************************************************************
module simplifications$
%********************************************************************
%  Routines for simplifications, contradiction testing
%  and substitution of functions
%  Author: Andreas Brand
%  1991 1993 1994
%  Updates, changes, additions by Thomas Wolf
%
% $Id: crsimp.red,v 1.17 1998/06/25 18:16:41 tw Exp tw $
%

symbolic procedure signchange(g)$
%  ensure, that the first term is positive
if pairp g then
 if (car g='MINUS) then cadr g
 else if (car g='PLUS) and (pairp cadr g) and (caadr g='MINUS)
      then reval list('MINUS,g)
      else g
else g$

symbolic procedure simplifyterm(p,ftem)$
%  simplify a single factor p of g=p*q*r*...=0
if (ftem:=smemberl(ftem,p)) then
  if pairp p and member(car p,'(MINUS SQRT QUOTIENT))  
  then simplifyterm(cadr p,ftem)
  else if pairp p and (car p='EXPT) then
       if smemberl(ftem,cadr p) then simplifyterm(cadr p,ftem)
                                else 1
  else if member((p:=signchange p),ineq_) then 1
     	      	       	                  else p
else if not p or zerop p then 0
                         else 1$

symbolic procedure contradictioncheck(s,pdes)$
begin scalar v,p$
if s then 
 while pdes do
   <<p:=car pdes$pdes:=cdr pdes$
   v:=get(p,'val)$
   if pairp v and (car v='TIMES) then
     (if member(s,cdr v) then
       <<v:=delete(s,cdr v)$
       update(p,if length v=1 then car v else cons('TIMES,v),
              get(p,'fcts),get(p,'vars),nil,list(0))>>)
   else if s=v then
    <<raise_contradiction(v,nil)$
    pdes:=nil>>
   >>$
return contradiction_$
end$

symbolic procedure raise_contradiction(g,text)$
<<contradiction_:=t$
  if print_ then
     <<terpri()$if text then write text
                        else write "contradiction : "$
  deprint list g>> >>$ 

symbolic procedure simplifypde(g,ftem,tofactor)$
%  simplify g=0
begin scalar h,l,ruli$
 if rulelist_ then g:=reval evalwhereexp list(rulelist_,g)$
 ruli:=start_let_rules()$
 g:=reval aeval g$
 stop_let_rules(ruli)$
 if g and not zerop g and not (ftem:=smemberl(ftem,g)) then 
   <<raise_contradiction(g,nil)$g:=1>>
 else if pairp g then
  if member(car g,'(EXPT QUOTIENT MINUS SQRT)) then 
     g:=simplifypde(cadr g,ftem,tofactor)
  else if member(car g,'(LOG LN LOGB LOG10)) then
     g:=simplifypde(reval list('PLUS,cadr g,-1),ftem,tofactor)
  else if tofactor then
   <<if car g='TIMES 
     then l:=for each a in cdr g join 
                 cdr reval list('FACTORIZE,a)
     else l:=cdr reval list('FACTORIZE,g)$
     while l do 
       <<if not member(car l,cdr l) then
            h:=union(list simplifyterm(car l,ftem),h)$
       l:=cdr l>>$
     h:=delete(1,h)$
     if null h then 
        <<raise_contradiction(g,nil)$
        g:=1>>
     else
       if pairp cdr h then g:=cons('TIMES,reverse h)
                      else g:=car h>>$
return g$
end$ 
     
symbolic procedure fcteval(p,less_vars)$
%  looks for a function which can be eliminated
%  if one is found, it is stored in p as (coeff_of_f.f)
%  if less_vars neq nil then the expr. to be substituted 
%  must have only fcts. of less vars
%  'to_eval neq nil iff not checked yet for substitution
%                   of if subst. possible
%               i.e. 'to_eval=nil if checked and not possible
%  'fcteval_lin includes subst. with coefficients that do not
%               include ftem functions/constants
%  'fcteval_nca includes subst. with non-vanishing coefficients
%               and therefore no case distinctions (linearity)
%  'fcteval_nli includes subst. with possibly vanishing coefficients
%               and therefore case distinctions (non-linearity)
begin scalar ft,a,b,fl,li,nc,nl,f,cpf,fv,fc$
  if flagp(p,'to_eval) then <<  
    b:=get(p,'not_to_eval)$   % functions that replace a derivative
    if (not get(p,'fcteval_lin)) and
       (not get(p,'fcteval_nca)) and
       (not get(p,'fcteval_nli)) then << 
      ft:=get(p,'allvarfcts)$                     
      if null ft then <<
        ft:=get(p,'rational)$
        % drop all functions f from ft for which there is another
        % function which is a function of all variables of f + at
        % least one extra variable
        for each f in ft do <<
          cpf:=get(p,'fcts)$
          fv:=fctargs f$
          while cpf and 
                (not_included(fv,fc:=fctargs car cpf) or
                 (length fv >= length fc)                ) do 
          cpf:=cdr cpf; 
          if null cpf then fl:=cons(f,cpf)
        >>$
        ft:=fl$
      >>$
      if ft then <<
	if (not less_vars)       or 
	   (not cdr ft)          or
	   (zerop get(p,'nvars)) then <<  
	  % either all functions allowed or only one fnc of all vars
	  for each f in ft do        
	  if (f neq b) and linear_fct(p,f) then <<
            % only linear algebr. fcts
	    a:=reval coeffn(get(p,'val),f,1)$
	    if fl:=smemberl(delete(f,get(p,'fcts)),a) then 
	    if freeofzero(a,fl,get(p,'vars)) then nc:=cons(cons(a,f),nc) 
					     else nl:=cons(cons(a,f),nl) 
						      else 
						  li:=cons(cons(a,f),li)
	  >>$
	  if li then put(p,'fcteval_lin,reverse li); % else
	  if nc then put(p,'fcteval_nca,reverse nc); % else
	  if nl then put(p,'fcteval_nli,reverse nl);
	  if not (li or nc or nl) then remflag1(p,'to_eval)
	>> 
      >>$
    >>$
    return (get(p,'fcteval_lin) or 
            get(p,'fcteval_nca) or 
	    get(p,'fcteval_nli)    )
  >>
end$

symbolic procedure freeofzero(p,ft,vl)$
%   gets p (factorized), if p does not vanish identically
if null ft then p
else 
begin scalar a,b,fr,pri$
 pri:=print_$
 print_:=nil$
 for each s in cdr reval list('FACTORIZE,p) do 
    a:=union(list simplifyterm(s,ft),a)$
 if length a>1 then p:=cons('TIMES,a)$
 while a do
  if null smemberl(ft,car a) or member(car a,ineq_) then a:=cdr a
  else if pairp cdr
    (b:=union(for each s in separ(car a,ft,vl) collect cdr s,nil)) then
      <<fr:=nil$
      while b do if freeofzero(car b,ft,vl) then <<b:=nil$fr:=t>>
      	       	     	      	       	    else b:=cdr b$
      if fr then a:=cdr a
            else <<a:=nil$p:=nil>> >>
    else <<a:=nil$p:=nil>>$
 print_:=pri$
return p
end$

symbolic procedure get_subst(l,length_limit,less_vars,no_df)$
%
% get the most simple pde from l which leads to a function substitution
% if less_vars neq nil: the expr. to subst. has only fcts. of less vars
% if no_df neq nil:     the expr. to subst. has no derivatives
begin scalar p,l1,l2,m,ntms,mdu$
  % mdu=(1:lin, 2:nca, 3:nli)

  % drop all equations longer than length_limit
  if length_limit then <<
    while l do
    if get(car l,'length)>length_limit then l:=nil
                                       else <<
      l1:=cons(car l,l1)$
      l:=cdr l
    >>$
    l:=reverse l1
  >>$
  % l is now the list of equations <= length_limit

  % next: substitution only if no_df=nil or 
  %       no derivative of any function occurs
  if no_df then <<
    for each s in l do 
      <<l2:=get(s,'derivs)$
      while l2 do 
 	if pairp(cdaar l2) then
 	   <<l2:=nil;
           l1:=cons(s,l1)>>
        else l2:=cdr l2>>$
    l:=setdiff(l,l1)>>$

  % next: restrict to substitutions, if any, 
  %       that have a coefficient without ftem-dependence
  l1:=nil;
  for each s in l do 
  if fcteval(s,less_vars) and
     get(s,'fcteval_lin) then l1:=cons(s,l1)$
  if l1 then <<l:=l1;mdu:=1>>
        else <<
   % at least substitutions that do not lead to subcases?
   for each s in l do 
   if get(s,'fcteval_nca) then l1:=cons(s,l1)$
   if l1 then <<l:=l1;mdu:=2>>
         else mdu:=3 
  >>$

  % next: find an equation with as many as possible variables 
  %       and few as possible terms for substitution    
  m:=-1$
  for each s in l do <<
   l1:=get(s,'nvars);
   if get(s,'starde) then l1:=sub1 l1;
   if l1>m then m:=l1$
  >>$

  while m>=0 do <<
    l1:=l$
    ntms:=10000000$
    while l1 do
    if ((get(car l1,'nvars) -
         if get(car l1,'starde) then 1 
                                else 0) = m )   and 
       fcteval(car l1,less_vars)                and
       (get(car l1,'terms) < ntms)              then <<
      p:=car l1$
      l1:=cdr l1$
      ntms:=get(p,'terms)$
    >>                             else l1:=cdr l1$
    m:=if p then -1
            else sub1 m
  >>$

  if p then return if mdu=1 then {mdu,p,car get(p,'fcteval_lin)} else
                   if mdu=2 then {mdu,p,car get(p,'fcteval_nca)} else
                                 {mdu,p,car get(p,'fcteval_nli)}
end$

symbolic procedure ineqsplit(q,ftem)$
% q into factors and
% drop quotients
begin scalar l$
 if pairp q and (car q='QUOTIENT) then q:=cadr q$
 q:=cdr reval list('FACTORIZE,q)$
 for each s in q do
     if smemberl(ftem,s) then
        <<s:=signchange s$
        if not member(s,l) then l:=cons(s,l)>>$
return l$
end$

symbolic procedure ineqsubst(new,old,ftem)$
% tests all q's in ineq_ for subst(new, old,q)=0
% result: nil, if 0 occurs
%         otherwise list of the subst(car p,...)
begin scalar l,a$
while ineq_ do
 if not my_freeof(car ineq_,old) then
 <<a:=simplifyterm(reval reval subst(new, old,car ineq_),ftem)$
   if zerop a then
   <<if print_ then
     <<terpri()$write "contradiction from the substitution:"$
       eqprint list('EQUAL,old,new)$
       write "because of the non-vanishing expression:"$
       eqprint car ineq_>>$
     contradiction_:=t$
     l:=nil$
     ineq_:=nil>>
   else
   <<l:=union(ineqsplit(a,ftem),l)$
     ineq_:=cdr ineq_>> >>
 else
 <<l:=cons(car ineq_,l)$
   ineq_:=cdr ineq_>>$
ineq_:=reverse l$
end$

symbolic procedure do_one_subst(ex,f,a,ftem,vl,level,eqn)$
% substitute f by ex in a
begin scalar l,l1,p,oldstarde$
 l:=get(a,'val)$
 oldstarde:=get(a,'starde)$
 if pairp l and (car l='TIMES) then l:=cdr l
                               else l:=list l$
 while l do <<  % for each factor
  if smember(f,car l) then <<
   p:=reval reval subst(ex,f,car l)$
   if not p or zerop p then <<l:=list nil$l1:=list 0>>
                       else <<
    if pairp p and (car p='QUOTIENT) then p:=cadr p$
    l1:=if no_of_terms(p)>max_factor then cons(p,l1)
                                     else
    append(reverse cdr reval list('FACTORIZE,p),l1)
   >> 
  >>                  else l1:=cons(car l,l1)$
 l:=cdr l
 >>$
 l:=nil$
 while l1 do <<
  if not member(car l1,cdr l1) then
  l:=union(list simplifyterm(car l1,ftem),l)$
  l1:=cdr l1
 >>$
 l:=delete(1,l)$
 if null l then <<
  if print_ then <<
   terpri()$  % new
   write"Substitution of "$
   fctprint list f$
   if cdr get(eqn,'fcts) then <<
    write " by an expression in "$terpri()$
    fctprint delete(f,get(eqn,'fcts))
   >>$
   write " found in ",eqn," : "$
   eqprint(list('EQUAL,f,ex))
  >>$
  raise_contradiction(get(a,'val),
                      "leads to a contradiction in : ")$
  a:=nil
 >>        else <<
  if pairp cdr l then l:=cons('TIMES,l)
                 else l:=car l$
  if get(a,'level) neq level then 
     a:=mkeq(l,ftem,vl,allflags_,nil,list(0))
  else <<
   for each b in allflags_ do flag(list a,b)$
   if null update(a,l,ftem,vl,nil,list(0)) then <<
    setprop(a,nil)$
    a:=nil
   >> 
  >>$
  put(a,'level,level)
 >>$
 if oldstarde and not get(a,'starde) then put(a,'dec_with,nil);
 return a$
end$

symbolic procedure do_subst(md,p,l,pde,ftem,forg,vl,plim)$
% md is the mode of substitution, needed in case of an ISE
% Substitute a function in all pdes 
begin scalar f,fl,h,ex,res,slim,too_large,was_subst,ps,
             ruli,ise,cf,vl,nof,stde$ %,l$
% l:=get(p,'fcteval_lin)$
% if null l then l:=get(p,'fcteval_nca)$
% if null l then l:=get(p,'fcteval_nli)$
% if l then << % l:=car l$
  f:=cdr l$
  cf:=car l$
  if get(p,'starde) then ise:=t;
  slim:=get(p,'length)$
  ruli:=start_let_rules()$
  ex:=reval aeval list('QUOTIENT,
                       list('PLUS,list('MINUS,get(p,'val)),
                                  list('TIMES,cf,f)),
                       cf)$
  if not ise then ineqsubst(ex,f,ftem)$
  if not contradiction_ then <<
   if not ise then <<
    fl:=delete(f,smemberl(ftem_,ex))$   % functions occuring in ex
    forg:=for each h in forg collect
          if atom h then
          if f=h then <<put(h,'fcts,fl)$was_subst:=t$list('EQUAL,f,ex)>>
                 else h
                    else 
          if (car h='EQUAL) and member(f,get(cadr h,'fcts)) then <<
           put(cadr h,'fcts,union(fl,delete(f,get(cadr h,'fcts))))$
           was_subst:=t$
           reval subst(ex,f,h)
          >>                                                else h$
   >>$
   if expert_mode then <<
    ps:=promptstring!*$
    promptstring!*:="pdes selected for substitution: "$
    l:=xread nil$ 
    promptstring!*:=ps$
    if idp l then <<if not member(l,pde) then l:=nil>>
             else if pairp l then
    if car l='!*comma!* then l:=intersection(pde,cdr l)
                        else l:=nil
                             else l:=nil$
    if not_included(pde,l) then too_large:=t$
    if not l then <<
     terpri()$
     write "This is not a pde (e.g. e2) or a list of pdes (e.g. e1,e2,e3)"
    >>
   >>             else l:=delete(p,pde)$
   % Do the substitution in all suitable equations
   if ise then <<
    h:=nil;
    vl:=get(p,'vars)$
    fl:=get(p,'fcts)$
    nof:=cdr get(p,'starde)$
    while l do <<
     if (stde:=get(car l,'starde)) and
        (nof<=cdr stde) and
        (not not_included(vl,get(car l,'vars))) and
        (not not_included(fl,get(car l,'fcts))) then h:=cons(car l,h);
     l:=cdr l
    >>$
    l:=h;
   >>$
   while l and not contradiction_ do <<
    if member(f,get(car l,'fcts)) then
    if not expert_mode and plim and (slim*get(car l,'length)>plim)
    then too_large:=t
    else <<
     pde:=eqinsert(do_one_subst(ex,f,car l,ftem,vl,get(p,'level),p),
                   delete(car l,pde))$
     for each h in pde do drop_rl_with(car l,h);         
     put(car l,'rl_with,nil);
     for each h in pde do drop_dec_with(car l,h,'dec_with_rl);     
     put(car l,'dec_with_rl,nil);
     flag(list car l,'to_int);
     was_subst:=t
    >>$
    l:=cdr l
   >>$
   if print_ and (not contradiction_) and was_subst then <<
    terpri()$write "Substitution of "$
    fctprint list f$
    if cdr get(p,'fcts) then <<
     write " by an "$
     if ise then write"(separable) "$
     write "expression in "$terpri()$
     fctprint delete(f,get(p,'fcts))
    >>$
    write " found in ",p," : "$
    eqprint(list('EQUAL,f,ex))
   >>$
   % To avoid using p repeatedly for substitutions of different
   % functions in the same equations:
   if ise then <<
    put(p,'fcteval_lin,nil);
    put(p,'fcteval_nca,nil);
    put(p,'fcteval_nli,nil);
    remflag1(p,'to_eval)$ % otherwise 'fcteval_??? would be computed again
    md:=md;   % only in order to do something with md if the next
              % statement is commented out
    % if too_large then
    % if md=1 then put(p,'fcteval_lin,list((cf . f))) else
    % if md=2 then put(p,'fcteval_nca,list((cf . f))) else
    %              put(p,'fcteval_nli,list((cf . f)))$  
    % could probably unnecessarily be repeated
   >>;
   % delete f and p if not anymore needed
   if (not ise) and 
      (not too_large) and 
      (not contradiction_) then <<
    if not assoc(f,depl_copy_) then 
    depl!*:=delete(assoc(f,depl!*),depl!*)$
    was_subst:=t$            % in the sense that pdes have been updated
    ftem_:=delete(f,ftem_)$
    pde:=delete(p,pde)$
    setprop(p,nil)
   >>$
%   if was_subst then 
   res:=list(pde,forg)  
   % also if not used to delete the pde if the function to be
   % substituted does not appear anymore
  >>$
  stop_let_rules(ruli)$
% >>$
 if not contradiction_ then return cons(was_subst,res)$
end$

symbolic procedure make_subst(pdes,forg,vl,l1,length_limit,pdelimit,
                              less_vars,no_df,no_cases,lin_subst,
                              mingrowth,cost_limit)$
% make a subst. 
% l1 is the list of possible "candidates"
begin scalar p,q,r,l,h,cop,cases_,w,md$   % ,ineq
  if expert_mode then l1:=selectpdes(pdes,1)$
again:
  if (mingrowth and (w:=search_subs(l1,cost_limit,no_cases))) or
     ((null mingrowth) and
      (w:=get_subst(l1,length_limit,less_vars,no_df))) then 

  if (   car w = 1)                             or
     ((lin_subst=nil)                     and
      ( (car w = 2)                  or
       ((car w = 3)            and 
        member(caaddr w,ineq_)     )    )     ) then <<
   l:=do_subst(car w,cadr w,caddr w,pdes,ftem_,forg,vl,pdelimit)$
   if l and null car l then << % not contradiction but not used
    l1:=delete(cadr w,l1);
    if l1 then <<
     pdes:=cadr l;
     forg:=caddr l;     
     l:=nil;
     goto again
    >>    
   >>;
   if l then l:=cdr l
  >>                                            else 
  if (null lin_subst) and (null no_cases) then <<
    md:=car w;   % md = type of substitution, needed in case of ISE
    p:=cadr  w;  % p = the equation
    w:=caddr w;  % w = (coeff . function)
    % make an equation from the coefficient
    q:=mkeq(car w,get(p,'fcts),get(p,'vars),allflags_,t,list(0))$
    % and an equation from the remainder
    r:=mkeq(list('PLUS,get(p,'val),
                 list('TIMES,car w,
                             list('MINUS,cdr w))),
            get(p,'fcts),get(p,'vars),allflags_,t,list(0))$
    if contradiction_ then <<
      contradiction_:=nil$
      h:=get(q,'val)$
      if pairp h and (car h='TIMES) then ineq_:=union(cdr  h,ineq_)
                                    else ineq_:=union(list h,ineq_)$      
      setprop(q,nil)$
      setprop(r,nil)$
      l:=do_subst(md,p,w,pdes,ftem_,forg,vl,pdelimit)$
      if print_ then  write"therefore no special investigation"$

      if l and null car l then << % not contradiction but not used
       l1:=delete(p,l1);
       if l1 then <<
	pdes:=cadr l;
	forg:=caddr l;     
        l:=nil;
	goto again
       >>    
      >>;
      if l then l:=cdr l

    >>                else <<
      cop:=backup_pdes(pdes,forg)$
      remflag1(p,'to_eval)$
      if print_ then <<
        terpri()$
        write "for the substitution of ",cdr w," by ",p$
        write " we have to consider the case : "$
        eqprint list('EQUAL,0,car w)
      >>$
      pdes:=eqinsert(q,delete(p,pdes))$
      pdes:=eqinsert(r,pdes)$      
      level_:=cons(1,level_)$
      print_level()$
      h:=get(q,'val)$    % to add it to ineq_ afterwards
      l:=crackmain(pdes,ineq_,forg,vl)$
      level_:=cons(2,level_)$  
      if print_ then <<
	print_level()$
        terpri()$
	write "now back to the substitution of ",cdr w," by ",p$
      >>$
      pdes:=restore_pdes(cop)$
      if pairp h and (car h='TIMES) then ineq_:=union(cdr  h,ineq_)
                                    else ineq_:=union(list h,ineq_)$      
      setprop(q,nil);
      setprop(r,nil);

      if contradiction_ then <<
	 contradiction_:=nil$
	 l:=do_subst(md,p,w,pdes,ftem_,forg,vl,pdelimit)$

	 if l and null car l then << % not contradiction but not used
	  l1:=delete(p,l1);
	  if l1 then <<
	   pdes:=cadr l;
	   forg:=caddr l;     
           l:=nil;
	   goto again
	  >>    
	 >>;
	 if l then l:=cdr l

      >>                else <<
        cases_:=t$
        cop:=nil;     % to save memory
	h:=crackmain(pdes,ineq_,forg,vl)$
	if contradiction_ then contradiction_:=nil
			  else l:=union(h,l)
      >>
    >>$
  >>$

  return if contradiction_  then nil % list(nil,nil) 
                            else if cases_ then list l
                                           else l$
end$


symbolic procedure get_fact_pde(pdes)$
% look for pde in pdes which can be factorized
begin scalar m,p,pv,fdegl,pdeg,h,f,fcc,fcl,tf,pmdeg,fmdeg,tfm,fm,pm$
 m:=0;   % counts how often the factor occurs that occurs most
 while pdes do <<
  p:=car pdes$  pdes:=cdr pdes$
  pv:=get(p,'val)$
  if pairp pv and (car pv='TIMES) then <<
   pv:=cdr pv$  % the list of factors in p
   fdegl:=for each f in pv collect 
          pde_degree(f,smemberl(get(p,'rational),f))$
   pdeg:=for each h in fdegl sum h$

   for each f in pv do <<

    % counting how often f has occured as factor
    fcc:=fcl$
    while fcc and (caar fcc neq f) do fcc:=cdr fcc$
    if fcc then <<
     h:=add1 cadar fcc;
     rplaca(cdar fcc,h);
    >>     else <<
     fcl:=cons({f,1,p},fcl);
     h:=1
    >>;
    tf:=nil;
    % Is f the best factor so far and p the best equation?
    if ( h>m                                                       ) or 
       ((h=m                                                 ) and 
        ((  pdeg < pmdeg                                ) or
         (( pdeg = pmdeg                          ) and
          ((car fdegl < fmdeg                ) or
           ((car fdegl = fmdeg         ) and 
            ((tf:=no_of_terms(f)) < tfm)     )    )     )    )     ) then <<
     m:=h;
     fm:=f;
     pm:=p;
     pmdeg:=pdeg;
     fmdeg:=car fdegl; 
     tfm:=if tf then tf
                else no_of_terms(f)$
    >>$
    fdegl:=cdr fdegl
   >>
  >>$
 >>$
 return if m=0 then nil 
               else <<
  put(pm,'val,cons('TIMES,cons(fm,delete(fm,cdr get(pm,'val)))))$
  pm
 >>
end$

algebraic procedure start_let_rules$
begin scalar ruli;
  ruli:={};
  let explog_$
  if sin(!%x)**2+cos(!%x)**2 neq 1         then <<ruli:=cons(1,ruli);let trig1_>> else ruli:=cons(0,ruli)$
  if cosh(!%x)**2 neq (sinh(!%x)**2 + 1)   then <<ruli:=cons(1,ruli);let trig2_>> else ruli:=cons(0,ruli)$
  if sin(!%x)*tan(!%x/2)+cos(!%x) neq 1    then <<ruli:=cons(1,ruli);let trig3_>> else ruli:=cons(0,ruli)$
  if sin(!%x)*cot(!%x/2)-cos(!%x) neq 1    then <<ruli:=cons(1,ruli);let trig4_>> else ruli:=cons(0,ruli)$
  if cos(2*!%x) + 2*sin(!%x)**2   neq 1    then <<ruli:=cons(1,ruli);let trig5_>> else ruli:=cons(0,ruli)$
  if sin(2*!%x) neq 2*cos(!%x)*sin(!%x)    then <<ruli:=cons(1,ruli);let trig6_>> else ruli:=cons(0,ruli)$
  if sinh(2*!%x) neq 2*sinh(!%x)*cosh(!%x) then <<ruli:=cons(1,ruli);let trig7_>> else ruli:=cons(0,ruli)$
  if cosh(2*!%x) neq 2*cosh(!%x)**2-1      then <<ruli:=cons(1,ruli);let trig8_>> else ruli:=cons(0,ruli)$
  return ruli;
end$

algebraic procedure stop_let_rules(ruli)$
begin
  clearrules explog_$
  if first ruli = 1 then clearrules trig8_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig7_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig6_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig5_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig4_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig3_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig2_$ ruli:=rest ruli$
  if first ruli = 1 then clearrules trig1_$ 
end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  procedures  for finding an optimal substitution  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure fbts(a,b)$
% fbts ... first better than second
(cadr   a <= cadr   b) and 
(caddr  a <= caddr  b) and 
(cadddr a <= cadddr b);


symbolic procedure list_subs(p,fevl,fli,mdu)$
% p is an equation, fevl a substitution list of p,
% fli is a list of lists (f,p1,p2,..) where 
%   f is a function, 
%   pi are lists (eqn,nco,nte,mdu) where
%   eqn is an equation that can be used for substituting f
%   nco is the number of terms of the coefficient of f in the eqn
%   nte is the number of terms without f in the eqn
%   mdu is the kind of substitution (1:lin, 2:nca, 3:nli)
begin
 scalar a,f,nco,nte,cpy,cc,ntry;
 for each a in fevl do <<
  f:=cdr a;
  nco:=no_of_terms(car a);
  nte:=get(p,'terms);
  nte:=if nte=1 then 0 
                else nte-nco$

  % Is there already any substitution list for f?
  cpy:=fli;
  while cpy and (f neq caar cpy) do cpy:=cdr cpy$
  ntry:={p,nco,nte,mdu}$

  if null cpy then fli:=cons({f,ntry},fli) % no, there was not
              else <<                         % yes, there was one
   cc:=cdar cpy$
   while cc and (null fbts(car cc,ntry)) do cc:=cdr cc$
   if null cc then << % ntry is at least in one criterium better 
                      % than a known one
    rplaca(cpy,cons(f,cons(ntry,cdar cpy)));
    cc:=cdar cpy$ % just the list of derivatives with ntry as the first
    while cdr cc do 
    if fbts(ntry,cadr cc) then rplacd(cc,cddr cc)
                          else cc:=cdr cc$
   >>
  >>
 >>;
 return fli
end$

algebraic procedure cwrno(n,r)$
% number of terms of (a1+a2+..+an)**r if ai are pairwise prime
% number of combinations of r factors out of n possible factors
% with repititions and without order = (n+r-1 over r)
<<n:=n+r-1;
  % The rest of the procedure computes binomial(n,r).
  if 2*r>n then k:=n-r;
  for i:=1:r product (n+1-i)/i
>>$

symbolic procedure besu(ic1,mdu1,ic2,mdu2)$
% Is the first substitution better than the second?
((mdu1<mdu2) and (ic1<=ic2)) or
((mdu1=mdu2) and (ic1< ic2)) or
% ########## difficult + room for improvement as the decision is
% actually dependent on how precious memory is 
% (more memory --> less cases and less time):
((mdu1=2) and (ic1<(ic2+ 4))) or
((mdu1=3) and (ic1<(ic2+25)))$

symbolic procedure search_subs(pdes,cost_limit,no_cases)$
begin
 scalar fli,p,el,f,fpl,dv,drf,d,ffl,hp,ff,nco,be,s,nte,ic,fp,
        rm,mc,subli,mdu,tr_search$

 % at first find the list of all functions that could be substituted
 % together with a list of pdes, the number of terms in the coeff and 
 % the type of substitution

% tr_search:=t$

 for each p in pdes do fcteval(p,nil)$

 fp:=pdes;
 while fp                                    and 
       ((get(car fp,'terms)>2)          or 
        (null get(car fp,'fcteval_lin))    ) do fp:=cdr fp;
 if fp then return {1,car fp,car get(car fp,'fcteval_lin)}$

 for each p in pdes do <<
  fli:=list_subs(p,get(p,'fcteval_lin),fli,1)$
  fli:=list_subs(p,get(p,'fcteval_nca),fli,2)$
  if null no_cases then fli:=list_subs(p,get(p,'fcteval_nli),fli,3)$
 >>$

 if tr_search then <<
  write"equations substitution: (eqn, no of coeff. t., no of other t., mdu)"$
  terpri()$
  for each el in fli do <<write el;terpri()>>$
 >>$

 if fli then
 if (null cdr   fli) and  % one function
    (null cddar fli) then % one equation, i.e. no choice
 return <<
  fli:=cadar fli;  % fli is now (eqn,nco,nte,mdu)
  mdu:=cadddr fli;
  {mdu,car fli,car get(car fli,if mdu = 1 then 'fcteval_lin else
                               if mdu = 2 then 'fcteval_nca else
                                               'fcteval_nli)     }
 >>                  else
 % (more than 1 fct.) or (only 1 function and more than 1 eqn.)
 for each el in fli do << % for any function to be substituted
                          % (for the format of fli see proc list_subs)

  f:=car el$ el:=cdr el$  
  % el is now a list of possible eqn.s to use for subst. of f

  fpl:=nil$ % fpl will be a list  of lists (p,hp,a1,a2,..) where 
            % p is an equation that involves f,
            % hp the highest power of f in p
            % ai are lists {ff,cdr d,nco} where ff is a derivative of f, 
            % cdr d its power and nco the number of coefficients
  for each p in pdes do << % for each equation in which f could be subst.
   dv:=get(p,'derivs)$    %   ((fct var1 n1 ...).pow)
   drf:=nil$
   for each d in dv do 
   if caar d = f then drf:=cons(d,drf)$
   % drf is now the list of powers of derivatives of f in p

   ffl:=nil$      % ffl will be a list of derivatives of f in p
                  % together with the power of f and number of 
                  % terms in the coeff.
   if drf then << % f occurs in this equation and we estimate the increase
    hp:=0$
    for each d in drf do <<
     if cdar d then ff:=cons('DF,car d)
               else ff:=caar d; 
     nco:=no_of_terms(coeffn(get(p,'val),ff,cdr d));
     if cdr d > hp then hp:=cdr d$
     ffl:=cons({ff,cdr d,nco},ffl);
    >>
   >>;

   if drf then fpl:=cons(cons(p,cons(hp,ffl)),fpl);
  >>$

  % now all information about all occurences of f is collected and for
  % all possible substitutions of f the cost will be estimated and the 
  % cheapest substitution for f will be determined

  be:=nil; % be will be the best equation with an associated min. cost mc
  for each s in el do << 
   % for each possible equation that can be used to subst. for f

   % number of terms of (a1+a2+..+an)**r = n+r-1 over r
   % f = (a1+a2+..+a_nte) / (b1+b2+..+b_nco)
   nco:=cadr s; 
   nte:=caddr s;
   ic:= - get(car s,'terms);  % ic will be the cost associated with
                              % substituting f by car s and car s
                              % will be dropped after the substitution
   for each fp in fpl do 
   if (car s) neq (car fp) then <<
    rm:=get(car fp,'terms);   % to become the number of terms without f
    hp:=cadr fp;
    ic:=ic - rm;              % as the old eqn. car fp will be replaced
    
    for each ff in cddr fp do << % for each power of each deriv. of f
     ic:=ic + (caddr ff)*           % number of terms of coefficient of ff
              cwrno(nte,cadr ff)*      % (numerator of f)**(power of ff)
              cwrno(nco,hp - cadr ff); % (denom. of f)**(hp - power of ff)
     rm:=rm - caddr ff;       % caddr ff is the number of terms with ff 
    >>;
    % Now all terms containing f in car fp have been considered. The
    % remaining terms are multiplied with (denom. of f)**hp
    ic:=ic + rm*cwrno(nco,hp)
   >>;

   % Is this substitution better than the best previous one?
   if (null be) or besu(ic,cadddr s,mc,mdu) then 
   <<be:=car s; mc:=ic; mdu:=cadddr s>>;

  >>;

  % It has been estimated that the substitution of f using the 
  % best eqn be has an additional cost of ic terms 

  if tr_search and (length el > 1) then <<
   terpri()$
   write"Best substitution for ",f," : ",{ic,f,be,mdu}$
  >>$

  if (null cost_limit) or (ic<cost_limit) then 
  subli:=cons({ic,mdu,f,be},subli)$
 >>$
 
 % Now pick the best substitution
 if subli then <<
  s:=car subli;
  subli:=cdr subli;
  for each el in subli do
  if besu(car el,cadr el,car s,cadr s) then s:=el$

  if tr_search then <<
   terpri()$
   write"Optimal substitution:"$terpri()$
   write"  replace ",caddr s," with the help of ",cadddr s,","$terpri()$
   if car s < 0 then write"  saving ", - car s," terms, "
                else write"  with a cost of ",car s," additional terms, "$
   terpri()$
   write if cadr s = 1 then "  linear substitution" else
         if cadr s = 2 then "  nonlinearity inceasing substitution" else
                            "  with case distinction"  $
  >>$

  el:=get(cadddr s,if (cadr s) = 1 then 'fcteval_lin else
                   if (cadr s) = 2 then 'fcteval_nca else
                                        'fcteval_nli);
  while (caddr s) neq (cdar el) do el:=cdr el;

  return {cadr s,cadddr s,car el}
     % = {mdu   ,p       ,car get(p,'fcteval_???)}
 >>$

end$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  procedures  for substitution of a derivative by a new function  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

symbolic procedure check_subst_df(pdes,forg)$
%  yields a list of derivatives which occur in all
%  pdes and in forg
begin scalar l,l1,l2,n,cp,not_to_substdf$
 if pdes then <<
  for each s in pdes do l:=union(for each a in get(s,'derivs) 
                                   collect car a,l)$    % all derivs
  for each s in forg do 
   if pairp s then l:=union(for each a in all_deriv_search(s,ftem_)
                                    collect car a,l)$
  l1:=df_min_list(l)$
  l:=nil$
  for each s in l1 do 
   if pairp s and not member(car s,not_to_substdf) then <<
    l:=cons(cons('DF,s),l)$
    not_to_substdf:=cons(car s,not_to_substdf)
  >> $

  % Derivatives of functions should only be substituted if the
  % function occurs in at least 2 equations
  while l do <<
   n:=0; % counter
   cp:=pdes;
   while cp and (n<2) do <<
    if member(car l,get(car cp,'fcts)) then n:=add1 n;
    cp:=cdr cp
   >>;
   if n=2 then l2:=cons(car l,l2);
   l:=cdr l
  >>
 >>$
 return l2$
end$

symbolic procedure df_min_list(dflist)$
%  yields the lowest derivative for each function in the list of
%  deriv. dflist.
%   e.g. dflist='((f x z) (g x) (g) (f y) (h x y) (h x z)) 
%   ==> result='(f g (h x))
if dflist then
begin scalar l,d,m,lmax$
while dflist do 
 <<m:=car dflist$
   dflist:=cdr dflist$
   l:=nil$
   while dflist do
    <<if (d:=df_min(car dflist,m)) then m:=d
       	     	      	       	   else l:=cons(car dflist,l)$
    dflist:=cdr dflist$
    >>$
   if pairp m and null cdr m then lmax:=cons(car m,lmax)
                             else lmax:=cons(m,lmax)$
   dflist:=l$
 >>$
return lmax$
end$

symbolic procedure df_min(df1,df2)$
%  yields the minimal derivative of d1,d2
%  e.g. df_min('(f x y),'(f x z))='(f x), df_min('(f x z),'(g x))=nil
<<if not pairp df1 then df1:=list df1$
if not pairp df2 then df2:=list df2$
if car df1=car df2 then 
  if (df1:=df_min1(cdr df1,cdr df2)) then cons(car df2,df1)
                                     else car df2>>$

symbolic procedure df_min1(df1,df2)$
begin scalar l,a$
while df1 do
 <<a:=car df1$
 if not zerop (a:=min(dfdeg(df1,car df1),dfdeg(df2,car df1))) then 
  l:=cons(car df1,l)$
 if a>1 then l:=cons(a,l)$
 df1:=cdr df1$
 if df1 and numberp car df1 then df1:=cdr df1>>$
return reverse l$
end$

symbolic procedure dfsubst_forg(p,g,d,forg)$
% substitute the function d in forg by an integral g
% of the function p
for each h in forg collect
  if pairp h and member(d,get(cadr h,'fcts)) then 
          <<put(cadr h,'fcts,
                cons(p,delete(d,get(cadr h,'fcts))))$
          reval subst(g,d,h)>>
  else h$

symbolic procedure expand_INT(p,varlist)$
if null varlist then p
else begin scalar v,n$
  v:=car varlist$
  varlist:=cdr varlist$
  if pairp(varlist) and numberp(car varlist) then
     <<n:=car varlist$
       varlist:=cdr varlist>>
  else n:=1$
  for i:=1:n do p:=list('INT,p,v)$
  return expand_INT(p,varlist)
end$


endmodule$

end$