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

%********************************************************************
module identities$
%********************************************************************
%  Routines for dealing with differential identities
%  Author: Thomas Wolf
%  May 1999

symbolic procedure drop_idty(id)$
% recycles a name of an identity
<<setprop(id,nil)$
  recycle_ids:=id . recycle_ids$
  idnties_:=delete(id,idnties_)>>$

symbolic procedure new_id_name$
% provides a name for a new identity
begin scalar id$
 if null recycle_ids then <<
  id:=mkid(idname_,nid_)$
  nid_:=add1 nid_$
 >>                  else <<
  id:=car recycle_ids$
  recycle_ids:=cdr recycle_ids
 >>$
 setprop(id,nil)$
 return id
end$

symbolic procedure replace_idty$
begin scalar p,ps,ex$
 ps:=promptstring!*$
 promptstring!*:=""$ 
 terpri()$
 write "If you want to replace an identity then ",
       "type its name, e.g. id_2 <ENTER>."$
 terpri()$
 write "If you want to add an identity then type `new_idty' <ENTER>. "$
 p:=termread()$
 if (p='NEW_IDTY) or member(p,idnties_) then
  <<terpri()$write "Input of a value for "$
  if p='NEW_IDTY then write "the new identity."
                 else write p,"."$
  terpri()$
  write "You can use names of other identities, e.g. 3*id_2 - df(id_1,x); "$
  terpri()$
  write "Terminate the expression with ; or $ : "$
  terpri()$
  ex:=termxread()$
  for each a in idnties_ do ex:=subst(get(a,'val),a,ex)$
  if p neq 'NEW_IDTY then drop_idty(p)$
  new_idty(reval ex,nil,nil)$

  terpri()$write car idnties_$
  if p='NEW_IDTY then write " is added"
                 else write " replaces ",p$
 >>
 else <<terpri()$
        write "An identity ",p," does not exist! (Back to previous menu)">>$
 promptstring!*:=ps$
end$


symbolic procedure trivial_idty(pdes,idval)$
if zerop idval or
   (pdes and
    null smemberl(pdes,search_li(idval,'DF))) % identity is purely alg.
   then t else nil$

symbolic procedure new_idty(idty,pdes,simp)$
% idty is the value of a new differential identity between equations
begin scalar id,idcp$
 if simp then idty:=simplifypde(idty,pdes,t,nil)$
 if not trivial_idty(pdes,idty) then <<
  idcp:=idnties_$
  while idcp and (get(car idcp,'val) neq idty) do idcp:=cdr idcp;
  if null idcp then <<
   id:=new_id_name();
   put(id,'val,idty)$
   flag1(id,'to_subst)$
   flag1(id,'to_int)$
   idnties_:=cons(id,idnties_)
  >>
 >>
end$

symbolic procedure show_id$
begin scalar l,n$
 terpri()$
 l:=length idnties_$
 write if l=0 then "No" else l, 
 if l=1 then " identity." else " identities"$ 
 if l=0 then terpri()
        else <<
  n:=1;
  for each l in reverse idnties_ do <<
   terpri()$
   algebraic write n,")  ",l," :  0 = ",lisp(get(l,'val));
   n:=add1 n;
   if print_all then <<
    terpri()$write "   to_int     : ",flagp(l,'to_int)$ 
    terpri()$write "   to_subst   : ",flap(l,'to_subst)$
   >>
  >>
 >>
end$

symbolic procedure del_red_id(pdes)$
begin scalar oldli,pl,s,idty,news,succ,p,l$ % ,r,newr$
 if idnties_ then <<
  oldli:=idnties_$
  while oldli do 
  if not flagp(car oldli,'to_subst) then oldli:=cdr oldli 
                                    else <<
   idty:=get(car oldli,'val)$
   pl:=smemberl(pdes,idty)$
   
   for each p in pl do l:=union(get(p,'vars),l)$
   if l then l:=length l else l:=0$

   pl:=setdiff(pl,search_li(idty,'DF));
   % now all pdes in pl are redundand, drop the longest
   if null pl then remflag1(car oldli,'to_subst)
              else <<
    drop_idty(car oldli);
    % find the longest equation s of those with the most variables
    s:=nil;
    while pl do <<
     if (null get(car pl,'starde)                     ) and
        (get(car pl,'nvars)=l                         ) and
        (null(% flagp(s,'to_int) or
              % flagp(s,'to_fullint)  or                              
              % flagp(s,'to_sep) or
              % flagp(s,'to_gensep) or
              % flagp(s,'to_decoup) or
              flagp(s,'to_eval))                      ) and
        ((null s                                 ) or
         (get(car pl,'nvars)>get(s,'nvars)       ) or
         ((get(car pl,'nvars)=get(s,'nvars)) and 
          (get(car pl,'terms)>get(s,'terms))     )    ) then s:=car pl;
     pl:=cdr pl
    >>;
    if null s then remflag1(car oldli,'to_subst)
              else <<
     if print_ then <<
      write "Equation ",s," is dropped as it is a consequence of others: "$
      algebraic write "0 = ",lisp(idty)$
     >>$
     % assuming s occurs linearly:
     pl:=coeffn(idty,s,1)$
     news:=reval {'QUOTIENT,{'DIFFERENCE,{'TIMES,pl,s},idty},pl};
     %for each r in idnties_ do
     %if not freeof(get(r,'val),s) then <<
     % newr:=reval subst(news,s,get(r,'val));
     % newr:=simplifypde(newr,pdes,t,nil)$
     % put(r,'val,newr)$
     % flag1(r,'to_subst)$
     % flag1(r,'to_int)$
     %>>$
     succ:=t$
     pdes:=drop_pde(s,pdes,news)$
     oldli:=cdr oldli 
    >>
   >>
  >>
 >>;
 if succ then return pdes
end$

symbolic procedure del_redundant_de(argset)$
begin scalar pdes;
 if pdes:=del_red_id(car argset) then return {pdes,cadr argset}$
end$

symbolic procedure write_id_to_file(pdes)$
begin scalar s,p,h,pl,ps$
 if idnties_ then <<
  ps:=promptstring!*$
  promptstring!*:=""$ 
  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(nequ_:=",nequ_,")$"$terpri()$
  write"off batch_mode$"$terpri()$
  write"list_of_variables:="$
  algebraic write lisp cons('LIST,vl_)$

  write"list_of_functions:="$
  algebraic write lisp cons('LIST,pdes)$

  for each h in pdes do 
  if pl:=assoc(h,depl!*) then 
  for each p in cdr pl do 
  algebraic write "depend ",lisp h,",",lisp p$
  
  write"list_of_equations:="$
  algebraic write 
  lisp( cons('LIST,for each h in idnties_ collect get(h,'val)));

  terpri()$ write"solution_:=crack(list_of_equations,{},"$
  terpri()$ write"                 list_of_functions,"$
  terpri()$ write"                 list_of_variables)$"$
  terpri()$
  terpri()$
  write"end$"$terpri()$
  shut s;
  on nat;
  promptstring!*:=ps$ 
 >>
end$

symbolic procedure remove_idl$
<<for each h in idnties_ do setprop(h,nil);
  idnties_:=nil>>$

symbolic procedure start_history(pdes)$
begin scalar l,ps$
 ps:=promptstring!*$
 promptstring!*:=""$
 write"For recording the history of equations all currently"$ terpri()$
 write"recorded histories would be deleted as well as all"$ terpri()$
 write"present decoupling information, i.e. `dec_with'"$ terpri()$
 write"would be set to nil. Please confirm (y/n). "$
 l:=termread()$
 if (l='y) or (l='Y) then <<
  record_hist:=t;
  for each l in pdes do put(l,'histry_,l)$
  for each l in pdes do put(l,'dec_with,nil)$
 >>;
 promptstring!*:=ps$
end$

symbolic procedure stop_history(pdes)$
<<record_hist:=nil;
  for each l in pdes do put(l,'histry_,l)>>$

%  write"Do you want to delete all dec_with information? (y/n) "$
%  l:=termread()$
%  if (l='y) or (l='Y) then 
%  for each l in pdes do put(l,'dec_with,nil)$

symbolic procedure idty_integration(argset)$
begin scalar l,pdes,idcp;
 pdes:=car argset;
 idcp:=idnties_;
 while idcp do
 if not flagp(car idcp,'to_int) then idcp:=cdr idcp else
 if l:=integrate_idty(car idcp,pdes,%cadr argset,
                      ftem_,vl_) then <<
  pdes:=l;idcp:=nil>>                                      else <<
  remflag1(car idcp,'to_int);
  idcp:=cdr idcp;
 >>;
 if l then return {pdes,cadr argset}
end$

symbolic procedure integrate_idty(org_idty,allpdes,%forg,
                                  fl,vl)$
% idty is a differential identity between equations
% allpdes, fl, vl are lisp lists of equation names, functions and variables
% ways to optimize: use conlaw instead of the initial intcurrent2
%                   use more general methods to take advantage of
%                   non-conservation laws
if idnties_ then
begin scalar cl,ncl,vlcp,xlist,eql,a,f,newpdes,ftem_bak,
             nx,dl,l,k,ps,idty,pdes,extrapdes,newidtylist$ %nclu
 if null org_idty then 
 if null cdr idnties_ then org_idty:=car idnties_ 
                      else <<
  show_id()$
  ps:=promptstring!*$
  promptstring!*:=""$ 
  write"Which of the identities shall be integrated? (no) "$
  k:=length(idnties_);
  repeat 
   l:=termread()
  until (fixp l) and (0<l) and (l<=k);
  org_idty:=nth(idnties_,k+1-l)$ 
  promptstring!*:=ps
 >>$

 idty:=reval num reval get(org_idty,'val)$
 if trivial_idty(allpdes,idty) then return nil$

 pdes:=smemberl(allpdes,idty)$
 a:=all_deriv_search(idty,pdes)$
 xlist:=smemberl(vl,a)$
 cl:=intcurrent3(idty,cons('LIST,pdes),cons('LIST,xlist))$
 % intcurrent3 is only successful if only 2 derivatives found
 if (not zerop caddr cl) and inter_divint then 
 cl:=intcurrent2(idty,cons('LIST,pdes),cons('LIST,xlist))$
 if zerop caddr cl then <<
  cl:=cdadr cl;   
  vlcp:=xlist;
  xlist:=nil;
  while vlcp do <<
   if not zerop car cl then <<
    ncl:=cons(car cl,ncl);
    xlist:=cons(car vlcp,xlist)
   >>;
   cl:=cdr cl;
   vlcp:=cdr vlcp
  >>;
%  ncl:=reverse ncl;
%  xlist:=reverse xlist;

  cl:=ncl;

%  % Now try to get a divergence in less differentiation variables.
%  % Each component of the divergence is tried to be written as
%  % a divergence in the other (right) variables
%  while ncl do <<
%   a:=intcurrent2(car ncl,cons('LIST,pdes),cons('LIST,cdr xlist))$
%   if not zerop caddr a then <<
%    cl:=cons(car ncl,cl);       ncl:=cdr ncl;
%    vlcp:=cons(car xlist,vlcp); xlist:=cdr xlist
%   >>                   else <<
%    % It was possible to integrate car ncl to div(cdadr a,cdr xlist).
%    % distribute {'DF,car a,car xlist} to the divergence of cdr ncl
%    ncl:=cdr ncl;
%    a:=cdadr a;
%    nclu:=nil;
%    while ncl do <<
%     nclu:=cons(reval {'PLUS,car ncl,{'DF,car a,car xlist}}, nclu);
%     ncl:=cdr ncl;
%     a:=cdr a
%    >>;
%    ncl:=reverse nclu;
%    xlist:=cdr xlist
%   >>
%  >>$
%  ncl:=cl;
%  xlist:=vlcp;
  nx:=length xlist;
  while pdes do <<
   ncl:=subst(get(car pdes,'val),car pdes,ncl);
   pdes:=cdr pdes
  >>$
  ftem_bak:=ftem_;
  eql:=int_curl(reval cons('LIST,ncl),  cons('LIST,fl),
                      cons('LIST,xlist),cons('LIST,varslist(ncl,ftem_,vl)) )$
  % eql has the form {'LIST,reval cons('LIST,resu),cons('LIST,neweq)}
  if (null eql) or (null cdadr eql) or (zerop cadadr eql) then return nil;
  eql:=cdr eql;
  if print_ then <<
   ncl:=for i:=1:nx collect {'DF,nth(cl,i),nth(xlist,i)};
   ncl:=if cdr ncl then cons('PLUS,ncl)
                   else car ncl;
   terpri()$
   write"The identity "$
   % mathprint idty$
   mathprint reval ncl;
   write"can be integrated to "$terpri()$
   deprint(cdar eql)$
  >>$
  if nx < 3 then a:='y else 
  if (null inter_divint) or !*batch_mode then <<
   a:='n;
   if print_ then <<
    write"The integrated divergence is not used because it ",
          "has more than 2 terms and"$    terpri()$
    if !*batch_mode then write"`inter_divint' is nil."
                    else write"batch_mode is on."$
   >>$
   terpri()
  >>                   else <<
   ps:=promptstring!*$
   promptstring!*:=""$ 
   write"Shall this integration be used? (y/n) "$
   repeat a:=termread() until (a='y) or (a='n);
   promptstring!*:=ps 
  >>;
  if a='n then <<
   a:=setdiff(ftem_,ftem_bak);
   for each f in a do drop_fct(f)$
   ftem_:=ftem_bak
  >>      else <<
   % the extra conditions from the generalized integration:
   extrapdes:=cdadr eql$
   eql:=cdar eql; % eql are now the integrated curl conditions
   drop_idty(org_idty)$ 
   while eql do <<
    if not zerop car eql then <<
     a:=mkeq(car eql,ftem_,vl,allflags_,nil,list(0),nil,allpdes);
     newpdes:=cons(a,newpdes);
    >>;
    eql:=cdr eql;
   >>;
   newpdes:=reverse newpdes;
   % formulate the new identities
   for i:=1:nx do <<
    idty:=nth(cl,i);
    if nx=1 then a:=car newpdes
            else <<
     % at first sum over df(q^{ji},j),  j<i
     l:=i-1;
     dl:=nx-2;
     a:=for j:=1:(i-1) collect << 
      k:=l;
      l:=l+dl;
      dl:=sub1 dl;
      {'DF,nth(newpdes,k),nth(xlist,j)}
     >>;
     a:=if null a then 0 else
        if cdr a then cons('PLUS,a)
                 else car a;
     idty:={'PLUS,idty,a};
     % then sum over -df(q^{ij},j), j>i
     if i=1 then l:=1
            else l:=k+nx-i+1;
     a:=for j:=(i+1):nx collect << 
      k:=l;
      l:=l+1;
      {'DF,nth(newpdes,k),nth(xlist,j)}
     >>;
     a:=if null a then 0 else
        if cdr a then cons('PLUS,a)
                 else car a;
    >>$
    newidtylist:=cons({'DIFFERENCE,idty,a},newidtylist);
   >>;
   eql:=nil;
   for each a in extrapdes do <<
    a:=mkeq(a,ftem_,vl,allflags_,t,list(0),nil,allpdes);
    allpdes:=eqinsert(a,allpdes);
    to_do_list:=cons(list('subst_level_35,%allpdes,forg,vl_,
                          list a),
                     to_do_list);
    eql:=cons(a,eql)
   >>;
   if print_ then <<
    write"Integration gives: "$
    listprint(newpdes)$terpri()$
    if eql then <<
     write"with extra conditions: "$ 
     listprint(eql)
    >>$
   >>;
   for each a in newpdes do allpdes:=eqinsert(a,allpdes)$
   % now that allpdes is updated:
   for each a in newidtylist do new_idty(a,allpdes,t)$
   return allpdes
  >>
 >>
end$

symbolic procedure sortpermuli(a)$
% a is a list of numbers to be sorted and the exchanges of neighbours
% are to be counted
begin scalar flp,conti,newa;
 repeat <<
  newa:=nil;
  conti:=nil;
  while cdr a do 
  if car a < cadr a then <<newa:=cons(car a,newa); a:=cdr a>>
                    else <<
   conti:=t;
   flp:=not flp;
   newa:=cons(cadr a,newa); 
   a:=cons(car a,cddr a);
  >>$
  newa:=cons(car a,newa);
  a:=reverse newa
 >> until null conti;
 return flp . a
end$

symbolic procedure curlconst(xlist,vl)$
% generates a list q^ij=r^ijk,_k with r^ijk totally antisymmetric
% in the order q^(n-1)n,...
% xlist is the list of xi,xj,xk
% vl is the list of all variables new functions should depend on

begin scalar n,qli,i,j,k,qij,a,flp,f,resu,qlicp$
 n:=length xlist$
 for i:=1:(n-1) do 
 for j:=(i+1):n do << % generation of r^ijk,k
  qij:=nil;
  for k:=1:n do 
  if (k neq i) and (k neq j) then <<
   a:=sortpermuli({i,j,k});
   flp:=car a;
   a:=cdr a;
   qlicp:=qli;
   while qlicp and (caar qlicp neq a) do qlicp:=cdr qlicp;
   if qlicp then f:=cdar qlicp
            else <<
    f:=newfct(fname_,vl,nfct_);
    nfct_:=add1 nfct_;
    ftem_:=fctinsert(f,ftem_);
    qli:=cons(a . f,qli)
   >>;
   f:={'DF,f,nth(xlist,k)};
   if flp then f:={'MINUS,f};
   qij:=cons(f,qij)
  >>$
  if null qij then <<qij:=newfct(fname_,setdiff(vl,xlist),nfct_);
                     nfct_:=add1 nfct_;
                     ftem_:=fctinsert(qij,ftem_)>>
              else
  if cdr qij then qij:=reval cons('PLUS,qij)
             else qij:=car qij;
  resu:=cons(qij,resu)  
 >>$
 return resu
end$

symbolic procedure updt_curl(h2,rmdr,fl,done_xlist,x,cdrxlist,n,k)$
% a subroutine of int_curl
begin scalar i,h4,h5,h6,h7,rmdr,y,pint,succ$
 if (not zerop reval reval {'DF,rmdr,x}) then <<
  if print_ then <<terpri()$write"No success."$terpri()>>$
  succ:=nil
 >>                                      else <<
  succ:=t;
  if done_xlist then << % there is one computed curl component to be updated
   % integration wrt done_xlist
   h7:=intcurrent2(rmdr,fl,cons('LIST,done_xlist));
   rmdr:=caddr h7;
   h7:=cdadr h7;
   % update the already computed h2-components with the new h7-comp.
   h4:=nil;
   h5:=-1;
   for i:=1:(k-1) do <<
    h5:=add1 h5;
    for h6:=1:(n-k) do <<h4:=cons(car h2,h4);h2:=cdr h2>>;
    h4:=cons({'DIFFERENCE,car h2,car h7},h4);
    h2:=cdr h2;
    h7:=cdr h7;
    for h6:=1:h5 do <<h4:=cons(car h2,h4);h2:=cdr h2>>
   >>;
   h2:=reverse h4;
  >>$
  % now generalized integration of the remainder
  if zerop rmdr then pint:=cons(0,nil)
                else <<
   y:=if cdrxlist then car cdrxlist
                  else car done_xlist;
   fnew_:=nil$
   pint:=partint(rmdr,fl,vl_,y,genint_);
   % genint is max number of new terms
   if null pint then succ:=nil
                else for each h4 in fnew_ do ftem_:=fctinsert(h4,ftem_)
  >>
 >>;
 return if null succ then nil
                     else cons(h2,pint) 
 % pint=cons(generalized integral of rmdr,list of new eqn)
end$


symbolic procedure int_curl(pli,fl,xlist,vl)$
% given a vector p^i satisfying p^i,_i=0, find q^{ij}=-q^{ji}
% such that p^i=q^{ij},j
% car result: (q^{12}, q^{13},.., q^{1n}, q^{23},.., q^{2n},.., q^{(n-1)n})
%             each q^{ij} comes with r^{ijk},k
% cdr result: list of new conditions in fewer variables
% works only if identically satisfied, not modulo some equations
% vl is the list of all relevant variables
% during computation is h2 = 
% (q^{kn},.., q^{k(k+1)},q^{(k-1)n},.., q^{(k-1)k},.., 
%  q^{2n},.., q^{23},    q^{1n},.., q^{13},q^{12}}     )
begin scalar h1,h2,h3,resu,newpli,xcp,done_xlist,n,k,ok,neweq,ftem_bak$
 % conversion from algebraic mode lists to lisp lists:
 pli:=cdr pli$ xlist:=cdr xlist$ vl:=cdr vl; xcp:=xlist$
 n:=length(xlist);
 k:=0;
 ok:=t;
 ftem_bak:=ftem_;
 if n=1 then return {'LIST,reval cons('LIST,pli),{'LIST}}$

 while cdr pli and ok do <<
  k:=add1 k;
  % the integration has to be done first wrt cdr xlist. The resulting
  % curl will be used to change the remining pli to be integrated
  h3:=intcurrent2(reval car pli,fl,cons('LIST,cdr xlist));
  pli:=cdr pli;
  h1:=cdadr h3;   
  h3:=reval reval caddr h3;           
  % h3 now = the remainder of the integration wrt cdr xlist

  if not zerop h3 then <<
   % here the integration wrt the done_xlist. These curl updates will
   % not be used to update pli, because df(h3,car xlist)=0 is assumed
   h3:=updt_curl(h2,h3,fl,done_xlist,car xlist,cdr xlist,n,k)$
   if null h3 then ok:=nil
              else <<      % generalized integration of the remainder
    neweq:=append(cddr h3,neweq);
    h2:=car h3;
    h1:=cons({'PLUS,car h1,cadr h3},cdr h1);
    % because of cdr xlist neq nil here q^{k(k+1)} is updated
   >>
  >>$

  if ok then << % In the first round h2 is naturally nil --> use ok for test
   % append (q^{kn},.., q^{k(k+1)}) and h2
   h2:=append(reverse h1,h2);
   % update the remaining pli to be integrated
   newpli:=nil;
   while h1 do <<
    newpli:=cons({'PLUS,{'DF,car h1,car xlist},car pli},newpli);
    h1:=cdr h1;
    pli:=cdr pli
   >>;
   pli:=reverse newpli
  >>;
  done_xlist:=cons(car xlist,done_xlist);
  xlist:=cdr xlist
 >>$
 if ok then <<
  pli:=reval car pli;
  % to get the remainder of the last component of pli integrated
  if pli neq 0 then <<
   k:=k+1;
   h3:=updt_curl(h2,pli,fl,done_xlist,car xlist,nil,n,k)$
   if null h3 then ok:=nil
              else <<
    neweq:=append(cddr h3,neweq);
    h2:=car h3;
    h2:=cons({'DIFFERENCE,car h2,cadr h3},cdr h2)
    % because of null xlist here car h2=q^{n-1,n} is updated
   >>$
  >>
 >>;
 if null ok then << % drop all new functions
  h1:=setdiff(ftem_,ftem_bak);
  for each h2 in h1 do drop_fct(h2)$
  ftem_:=ftem_bak
 >>         else <<
  h1:=curlconst(xcp,vl)$
  while h1 do <<
   resu:=cons({'PLUS,car h2,car h1},resu);
   h1:=cdr h1; 
   h2:=cdr h2; 
  >>$
  return {'LIST,reval cons('LIST,resu),cons('LIST,neweq)}
 >>
end$

endmodule$

end$