Artifact b7e9b30a39fafae40f7ae099d7edb416d7b09ede68625bc42d5aaacca1f4016a:
- Executable file
r37/packages/crack/crutil.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: 60399) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** module utilities$ %******************************************************************** % Routines for finding leading derivatives and others % Author: Andreas Brand % 1990 1994 % Updates by Thomas Wolf % % $Id: crutil.red,v 1.23 1998/06/25 18:20:43 tw Exp tw $ % %%%%%%%%%%%%%%%%%%%%%%%%% % 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 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 << 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(f,de1,de2,rl)$ % add (f 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(f,a)$ a:=delete(b,a)$ if b then b:=cons(f,cons(de1,cdr b)) else b:=list(f,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(f,de1,de2,rl)$ % add (f de1) to 'dec_with or 'dec_with_rl of de2 and % add (f de2) to 'dec_with or 'dec_with_rl of de1 begin add_dec_with(f,de1,de2,rl)$ add_dec_with(f,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)$ 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)$ 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 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 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_)$ return val end$ symbolic procedure mkeqlist(vallist,ftem,vl,flaglist,simp_flag,orderl)$ % 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 begin scalar l1$ for each a in vallist do l1:=eqinsert(mkeq(a,ftem,vl,flaglist,simp_flag,orderl),l1)$ return l1 end$ symbolic procedure mkeq(val,ftem,vl,flaglist,simp_flag,orderl)$ % 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 begin scalar s$ s:=mkid(eqname_,nequ_)$ nequ_:=add1 nequ_$ setprop(s,nil)$ for each a in flaglist do flag(list s,a)$ if not update(s,val,ftem,vl,simp_flag,orderl) then <<nequ_:=sub1 nequ_$ setprop(s,nil)$ s:=nil>>$ return s end$ symbolic procedure update(equ,val,ftem,vl,simp_flag,orderl)$ % 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 begin scalar l$ if val and not zerop val then << ftem:=reverse union(smemberl(ftem,val),nil)$ 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 % get(equ,'starde) then simplifypde(val,ftem,nil) else simplifypde(val,ftem,t)$ if val and not zerop val then << ftem:=reverse union(smemberl(ftem,val),nil)$ put(equ,'terms,no_of_terms(val))$ >> >>$ >>$ 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:=reverse l$ put(equ,'vars,vl)$ put(equ,'nvars,length vl)$ put(equ,'level,level_)$ put(equ,'derivs,sort_derivs(all_deriv_search(val,ftem),ftem,vl))$ put(equ,'fcteval_lin,nil)$ put(equ,'fcteval_nca,nil)$ put(equ,'fcteval_nli,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,'starde,stardep(ftem,vl))$ if (l:=get(equ,'starde)) then <<%remflag1(equ,'to_eval)$ remflag1(equ,'to_int)$ 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 % remflag1(equ,'to_decoup) 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 <<for each f in get(equ,'fcts) do val:=subst(0,f,val)$ if not zerop reval reval val then <<contradiction_:=t$ terpri()$ write "Contradiction in ",equ, "!!! PDE has to be homogeneous!!!">> >>$ remprop(equ,'full_int)$ return equ >> 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$ % if print_ and tr_main then % <<terpri()$write "inserting ",s," in the pde list: ">>$ l1:=proddel(s,l)$ if car l1 then <<n:=get(s,'length)$ l:=cadr l1$ l1:=nil$ % if print_ and tr_main then write s," is inserted."$ 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)$ % 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,go_on$ if pairp(lnew:=get(s,'val)) and (car lnew='TIMES) then lnew:=cdr lnew else lnew:=list lnew$ n:=length lnew$ go_on:=t$ while l do << % while l and go_on do << % in case one wants to stop if s % is found to be a consequence of l 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 setprop(car l,nil) >> else << if null not_included(l1,lnew) then % s is a consequence of car l <<if print_ and tr_main then <<terpri()$write s," is a consequence of ",car l,".">>$ % l2:=append(reverse l,l2); % one could stop here but continuation can still be useful go_on:=nil$ >>$ % else l2:=cons(car l,l2)$ >>; l:=cdr l >>$ if print_ and tr_main and l3 then << terpri()$listprint l3$ if cdr l3 then write " are consequences of ",s else write " is a consequence of ",s >>$ if null go_on then <<setprop(s,nil);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 typeeq(s)$ % print equation if get(s,'printlength)>print_ then begin scalar a,b$ write "expr. with ",(a:=get(s,'terms))," terms"$ if a neq (b:=get(s,'length)) then write", ",b," factors"$ write" in "$ myprin2l(get(s,'fcts),", ")$ terpri()$ end else mathprint list('EQUAL,0,get(s,'val))$ symbolic procedure typeeqlist(l)$ % print equations and its property lists <<terpri()$ for each s in l do <<%terpri()$ write s," : "$ typeeq(s)$%terpri()$ % write s," : weight = ",get(s,'length),", "$ % if print_ and (get(s,'printlength)>print_) then % <<write "expr. with ",get(s,'terms)," terms in "$ % myprin2l(get(s,'fcts),", ")$ % terpri()>> % else % mathprint list('EQUAL,0,get(s,'val))$ if print_all then << write " fcts : ",get(s,'fcts)$ terpri()$write " vars : ",get(s,'vars)$ terpri()$write " nvars : ",get(s,'nvars)$ terpri()$write " derivs : ",get(s,'derivs)$ terpri()$write " terms : ",get(s,'terms)$ 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 " 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_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 " not_to_eval: ",get(s,'not_to_eval)$ terpri()$write " used_ : ",flagp(s,'used_)$ terpri()$write " orderings : ",get(s,'orderings)$ 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) otherweise, 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$ % a new minimum v:=list car vl>> 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 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$ %%%%%%%%%%%%%%%%%%%%%%%%% % 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,ftem,vl)$ symbolic procedure abs_dfrel(p,q,vl)$ % the relation "p is a lower derivative than q" (total deg. ord.) % p,q : derivatives or functions from ftem % ftem : list of fcts % vl : list of vars begin scalar a$ return if lex_ then dfrel1(p,q,vl) else if zerop (a:=absdeg(cdar p)-absdeg(cdar q)) then dfrel1(p,q,vl) else a<0$ end$ %%symbolic procedure lower_lderiv(p,q,ftem,vl)$ %symbolic procedure lower_lderiv(p,q,vl)$ %% the relation "p has a lower (absolute) derivative than q" %if abs_ld_deriv(p) and abs_ld_deriv(q) then %% abs_dfrel(car get(p,'derivs),car get(q,'derivs),ftem,vl) % abs_dfrel(car get(p,'derivs),car get(q,'derivs),vl) %else if abs_ld_deriv(q) then t$ 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 sort_derivs(l,ftem,vl)$ % yields a sorted list of all derivatives in p begin scalar l1,l2,a,fa$ return if null l then nil else <<a:=car l$ fa:=caar a$ l:=cdr l$ while l do <<if fa=caaar l then % if abs_dfrel(a,car l,ftem,vl) then l1:=cons(car l,l1) % else l2:=cons(car l,l2) if abs_dfrel(a,car l,vl) then l1:=cons(car l,l1) else l2:=cons(car l,l2) else if fa=which_first(fa,caaar l,ftem) then l2:=cons(car l,l2) else l1:=cons(car l,l1)$ 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_ 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)$ % liefert t, falls p einfacherer Differentialausdruck, sonst nil % 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)$ if null v then % same derivatives, if cdr p>cdr q then nil % considering powers else t 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 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 degree 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 degree else cadr p$ % v with degree 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_ 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 funktion 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) 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)$ ftem_:=reverse ftem_$ for each ff in reverse fnew_ do if not member(ff,ftem_) then ftem_:=fctinsert(ff,ftem_)$ ftem_:=reverse ftem_$ ex:=list('EQUAL,h,p)>>$ return ex end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % general purpose procedures % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 fctinsert(f,ftem)$ % isert a function f in the function list ftem if null ftem then list f else if fctlength car ftem >= fctlength f then cons(f,ftem) else cons(car ftem,fctinsert(f,cdr ftem))$ symbolic procedure newfct(id,l,nfct)$ begin scalar f$ 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 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 my_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)$ return reverse 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 a funktion if pairp l then <<write car l$ for each v in cdr l do <<write ","$write v>> >>$ symbolic procedure fctprint1(f)$ % print a funktion if f then if pairp f then <<write car f$ if pairp cdr f then <<write "("$ listprint cdr f$ write ")">> >> else write f$ symbolic procedure fctprint(fl)$ % Ausdrucken der Funktionen aus fl begin scalar n,l,f$ n:=0$ while fl do <<f:=car fl$ fl:=cdr fl$ if pairp f then if car f='EQUAL then <<n:=delength 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 eq 4 then <<terpri()$n:=0>>$ fctprint1 f$ if fl then write ", "$ n:=add1 n>> else <<if n eq 4 then <<terpri()$n:=0>>$ l:=assoc(f,depl!*)$ fctprint1 if l then l else f$ if fl then write ", "$ n:=add1 n>> >>$ %if n neq 0 then terpri() 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 "$write n$write " terms"$ terpri() >> else mathprint e$ end$ symbolic procedure print_level$ if print_ and level_ then << terpri()$ write "New 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$ %--- 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": "$ k:=0; while el do << if k=4 then <<k:=0;terpri()>> else k:=add1 k$ write car el,"(",get(car el,'length)$ if (s:=get(car el,'starde)) then for r:=1:(1+cdr s) do write"*"$ write")"$ el:=cdr el$ if el then write","$ >> >>$ 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 print_pde_ineq(pdes,ineqs)$ % print all pdes and ineqs <<if pdes then <<terpri()$terpri()$write "equations : "$ typeeqlist(pdes)>> else <<terpri()$terpri()$write "no equations">>$ if ineqs then <<terpri()$write "non-vanishing expressions: "$ for each aa in ineqs do eqprint aa>>$ >>$ symbolic procedure print_fcts(fcts,vl)$ % print all fcts and vars <<if fcts then <<terpri()$write "functions : "$ fctprint(fcts)>>$ 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_pde_ineq(pdes,ineqs)$ print_fcts(fcts,vl)$ print_statistic(pdes,fcts) end$ symbolic procedure no_of_terms(d)$ if not pairp d then if d then 1 else 0 else if car d='PLUS then length d - 1 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 1 else 1$ symbolic procedure desort(l)$ % list sorting begin scalar l1,l2,l3,m,n$ return if null l then nil else <<n:=delength car l$ l2:=list car l$ l:=cdr l$ while l do <<m:=delength 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(desort(l1),append(l2,desort(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 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_)$ return cop 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$ 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) 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,ps$ ps:=promptstring!*$ promptstring!*:=""$ terpri()$ write "pde to be deleted: "$ s:=termread()$ promptstring!*:=ps$ if member(s,pdes) then setprop(s,nil)$ return delete(s,pdes)$ end$ symbolic procedure addfunction(ft)$ begin scalar f,l,ps$ ps:=promptstring!*$ promptstring!*:=""$ terpri()$ write "What is the name of the new function? "$ terpri()$ write "(terminate with Enter, no ; or $) "$ f:=termread()$ depl!*:=delete(assoc(f,depl!*),depl!*)$ terpri()$ write "Give a list of variables ",f," depends on, for example x,y,z; "$ terpri()$ write "If ",f," shall be constant then just input ; "$ 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)$ promptstring!*:=ps$ return ft end$ symbolic procedure change_pde_flag(pdes)$ begin scalar p,fla$ repeat << terpri()$ write "From which PDE do you want to change a flag, e.g. e23?"$ terpri()$ p:=termread()$ >> until not freeof(pdes,p)$ repeat << terpri()$ write "Type in one of the following flags which should be flipped:"$ terpri()$ write allflags_; terpri()$ fla:=termread()$ >> until not freeof(allflags_,fla)$ if flagp(p,fla) then remflag1(p,fla) else flag(list(p),fla) end$ symbolic procedure write_in_file(pdes,ftem)$ begin scalar p,pl,s,h,ps,wn$ ps:=promptstring!*$ promptstring!*:=""$ terpri()$ write "Enter a list of equations, like e2,e5,e35; from: "$terpri()$ fctprint pdes$terpri()$ write "To write all equations just enter ; "$terpri()$ repeat << s:=termxread()$ h:=s; if s=nil then pl:=pdes else if s and (car s='!*comma!*) then << s:=cdr s; 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 (without ;) : "$ s:=termread()$ off nat; out s; 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 ",",v$ write"$"$terpri() >> >>$ if wn then for each h in pl do algebraic write h,":=",lisp get(h,'val) else for each h in pl do algebraic write lisp get(h,'val)$ write"end$"$terpri()$ shut s; on nat; promptstring!*:=ps$ end$ symbolic procedure replacepde(pdes,ftem,vl)$ begin scalar p,q,ex,ps,h$ ps:=promptstring!*$ promptstring!*:=""$ repeat << terpri()$ write "Is there a (further) 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 ftem:=addfunction(ftem) >> until (p='n) or (p='N)$ terpri()$ write "If you want to replace a pde then type its name, e.g. e23."$ terpri()$ write "If you want to add a pde then type `new_pde'. "$ terpri()$ write "In both cases terminate with Enter (no ; or $) : "$ 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*e12 - df(e13,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) (2)"$ terpri()$ write "- simplified completely (3) "$ h:=termread()$ if h=2 then ex:=reval ex$ if h<3 then h:=nil else h:=t$ if p neq 'NEW_PDE then <<setprop(p,nil)$pdes:=delete(p,pdes)>>$ q:=mkeq(ex,ftem,vl,allflags_,h,list(0))$ 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 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: "$ 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: "$ 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 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 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,ordr)$ % ev, iv are algebraic mode lists % generates a polynomial in the variables ev of order ordr % with functions of the variables iv begin scalar a,fi,el1,el2,f,fl,p,pr; a:=reval list('EXPT,cons('PLUS,cons(1,cdr ev)),ordr)$ 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,vl,fl) changes the ordering of the % list of derivatives depending on the current ordering (this % is detected "automatically" by sort_derivs using the lex_ flag to % toggle between total-degree and lexicogrpahic. % 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)$ >>$ return pdes end$ % % check_globals() is used to check that the values with which CRACK is % being called are coherent and/or valid. % % !FIXME! would be nice if we had a list which we could use to store % the global variables and the associated valid values. % symbolic procedure check_globals()$ begin scalar flag$ flag := nil$ if !*batch_mode neq nil and !*batch_mode neq t then << terpri()$ write "!*batch_mode needs to be a boolean: ", !*batch_mode," is invalid"$ flag := t$ >>$ if printmenu_ neq nil and printmenu_ neq t then << terpri()$ write "printmenu_ needs to be a boolean: ", printmenu_," is invalid"$ flag := t$ >>$ if expert_mode neq nil and expert_mode neq t then << terpri()$ write "expert_mode needs to be a boolean: ", expert_mode," is invalid"$ flag := t$ >>$ if repeat_mode neq nil and repeat_mode neq t then << terpri()$ write "repeat_mode needs to be a boolean: ", repeat_mode," is invalid"$ flag := t$ >>$ if not fixp genint_ and genint_ neq nil then << terpri()$ write "genint_ needs to be an integer or nil: ", genint_," is invalid"$ flag := t$ >>$ if not fixp facint_ and facint_ neq nil then << terpri()$ write "facint_ needs to be an integer or nil: ", facint_," is invalid"$ flag := t$ >>$ if potint_ neq nil and potint_ neq t then << terpri()$ write "potint_ needs to be a boolean: ", potint_," is invalid"$ flag := t$ >>$ if safeint_ neq nil and safeint_ neq t then << terpri()$ write "safeint_ needs to be a boolean: ", safeint_," is invalid"$ flag := t$ >>$ if freeint_ neq nil and freeint_ neq t then << terpri()$ write "potint_ needs to be a boolean: ", potint_," is invalid"$ flag := t$ >>$ if not fixp odesolve_ then << terpri()$ write "odesolve_ needs to be an integer: ", odesolve_," is invalid"$ flag := t$ >>$ if not fixp max_factor then << terpri()$ write "max_factor needs to be an integer: ", max_factor," is invalid"$ flag := t$ >>$ if not fixp gensep_ then << terpri()$ write "gensep_ needs to be an integer: ", gensep_," is invalid"$ flag := t$ >>$ if not fixp new_gensep and new_gensep neq nil then << terpri()$ write "new_gensep needs to be an integer or nil: ", new_gensep," is invalid"$ flag := t$ >>$ if not fixp subst_0 then << terpri()$ write "subst_0 needs to be an integer: ", subst_0," is invalid"$ flag := t$ >>$ if not fixp subst_1 then << terpri()$ write "subst_1 needs to be an integer: ", subst_1," is invalid"$ flag := t$ >>$ if not fixp subst_2 then << terpri()$ write "subst_2 needs to be an integer: ", subst_2," is invalid"$ flag := t$ >>$ if not fixp subst_3 then << terpri()$ write "subst_3 needs to be an integer: ", subst_3," is invalid"$ flag := t$ >>$ if not fixp subst_4 then << terpri()$ write "subst_4 needs to be an integer: ", subst_4," is invalid"$ flag := t$ >>$ if not fixp cost_limit5 then << terpri()$ write "cost_limit5 needs to be an integer: ", cost_limit5," is invalid"$ flag := t$ >>$ if not fixp max_red_len then << terpri()$ write "max_red_len needs to be an integer: ", max_red_len," is invalid"$ flag := t$ >>$ if not fixp pdelimit_0 and pdelimit_0 neq nil then << terpri()$ write "pdelimit_0 needs to be an integer or nil: ", pdelimit_0," is invalid"$ flag := t$ >>$ if not fixp pdelimit_1 and pdelimit_1 neq nil then << terpri()$ write "pdelimit_1 needs to be an integer or nil: ", pdelimit_1," is invalid"$ flag := t$ >>$ if not fixp pdelimit_2 and pdelimit_2 neq nil then << terpri()$ write "pdelimit_2 needs to be an integer or nil: ", pdelimit_2," is invalid"$ flag := t$ >>$ if not fixp pdelimit_3 and pdelimit_3 neq nil then << terpri()$ write "pdelimit_3 needs to be an integer or nil: ", pdelimit_3," is invalid"$ flag := t$ >>$ if not fixp pdelimit_4 and pdelimit_4 neq nil then << terpri()$ write "pdelimit_4 needs to be an integer or nil: ", pdelimit_4," is invalid"$ flag := t$ >>$ if not numberp length_inc then << terpri()$ write "length_inc needs to be an number: ", length_inc," is invalid"$ flag := t$ >>$ if tr_main neq t and tr_main neq nil then << terpri()$ write "tr_main needs to be a boolean: ", tr_main," is invalid"$ flag := t$ >>$ if tr_gensep neq t and tr_gensep neq nil then << terpri()$ write "tr_gensep needs to be a boolean: ", tr_gensep," is invalid"$ flag := t$ >>$ if tr_genint neq t and tr_genint neq nil then << terpri()$ write "tr_genint needs to be a boolean: ", tr_genint," is invalid"$ flag := t$ >>$ if tr_decouple neq t and tr_decouple neq nil then << terpri()$ write "tr_decouple needs to be a boolean: ", tr_decouple," is invalid"$ flag := t$ >>$ if tr_redlength neq t and tr_redlength neq nil then << terpri()$ write "tr_redlength needs to be a boolean: ", tr_redlength," is invalid"$ flag := t$ >>$ if tr_orderings neq t and tr_orderings neq nil then << terpri()$ write "tr_orderings needs to be a boolean: ", tr_orderings," is invalid"$ flag := t$ >>$ if homogen_ neq t and homogen_ neq nil then << terpri()$ write "homogen_ needs to be a boolean: ", homogen_," is invalid"$ flag := t$ >>$ if solvealg_ neq t and solvealg_ neq nil then << terpri()$ write "solvealg_ needs to be a boolean: ", solvealg_," is invalid"$ flag := t$ >>$ if print_more neq t and print_more neq nil then << terpri()$ write "print_more needs to be a boolean: ", print_more," is invalid"$ flag := t$ >>$ if print_all neq t and print_all neq nil then << terpri()$ write "print_all needs to be a boolean: ", print_all," is invalid"$ flag := t$ >>$ if logoprint_ neq t and logoprint_ neq nil then << terpri()$ write "logoprint_ needs to be a boolean: ", logoprint_," is invalid"$ flag := t$ >>$ if poly_only neq t and poly_only neq nil then << terpri()$ write "poly_only needs to be a boolean: ", poly_only," is invalid"$ flag := t$ >>$ if time_ neq t and time_ neq nil then << terpri()$ write "time_ needs to be a boolean: ", time_," is invalid"$ flag := t$ >>$ if (not null print_) and (not fixp print_) then << terpri()$ write "print_ needs to be an integer: ", print_," is invalid"$ flag := t$ >>$ if not fixp dec_hist then << terpri()$ write "dec_hist needs to be an integer: ", dec_hist," is invalid"$ flag := t$ >>$ if not fixp maxalgsys_ then << terpri()$ write "maxalgsys_ needs to be an integer: ", maxalgsys_," is invalid"$ flag := t$ >>$ if adjust_fnc neq t and adjust_fnc neq nil then << terpri()$ write "adjust_fnc needs to be a boolean: ", adjust_fnc," is invalid"$ flag := t$ >>$ if not vectorp orderings_ and orderings_ neq nil then << terpri()$ write "orderings_ needs to be a vector: ", orderings_," is invalid"$ flag := t$ >>$ if simple_orderings neq t and simple_orderings neq nil then << terpri()$ write "simple_orderings needs to be a boolean: ", simple_orderings," is invalid"$ flag := t$ >>$ if lex_ neq t and lex_ neq nil then << terpri()$ write "lex_ needs to be a boolean: ", lex_," is invalid"$ flag := t$ >>$ if collect_sol neq t and collect_sol neq nil then << terpri()$ write "lex_ needs to be a boolean: ", lex_," is invalid"$ flag := t$ >>$ if flag then return nil else return t$ end$ endmodule$ end$