Artifact c1c98eb225effb84bc79f2fa0b87b1d3761ac0a5871fa6417e7e93d4b8c64a8d:
- Executable file
r37/packages/crack/crsimp.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 34797) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** 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$