Artifact e73a5ed57025d3b0b534395012d09b18d1678d46f344f0adf3bf70fc9708374c:
- Executable file
r38/packages/crack/crshort.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: 36479) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** module shortening$ %******************************************************************** % Routines for algebraically combining de's to reduce their length % Author: Thomas Wolf % Jan 1998 % % $Id: crshort.red,v 1.4 1998/04/28 21:36:27 arrigo Exp $ % symbolic procedure alg_length_reduction(arglist)$ % Do one length-reducing combination of two equations begin scalar pdes,l,l1$ %,cpu,gc$ % cpu:=time()$ gc:=gctime()$ pdes:=car arglist$ if expert_mode then l:=selectpdes(pdes,2) else l:=pdes$ !*rational_bak := cons(!*rational,!*rational_bak)$ if !*rational then algebraic(off rational)$ if struc_eqn then << while l do << if is_algebraic(car l) then l1:=cons(car l,l1); l:=cdr l >>$ l:=reverse l1 >>$ if l and cdr l and (l1:=err_catch_short(l,caddr arglist,pdes)) then <<for each a in cdr l1 do pdes:=drop_pde(a,pdes,nil)$ for each a in car l1 do if a then pdes:=eqinsert(a,pdes)$ for each a in car l1 do if a then dec_fct_check(a,pdes)$ l:=nil; for each a in car l1 do if a then l:=cons(a,l); l:=list(pdes,cadr arglist,l)>> else l:=nil$ %if print_ and !*time then << % write " time : ", time() - cpu, % " ms GC time : ", gctime() - gc," ms " %>>$ if !*rational neq car !*rational_bak then if !*rational then algebraic(off rational) else algebraic(on rational)$ !*rational_bak:= cdr !*rational_bak$ return l$ end$ %------------------- symbolic procedure err_catch_short(a1,vl,pdes)$ begin scalar h,bak,kernlist!*bak,kord!*bak,mi,newp,p1,bakup_bak; bak:=max_gc_counter$ max_gc_counter:=my_gc_counter+max_gc_short; kernlist!*bak:=kernlist!*$ kord!*bak:=kord!*$ bakup_bak:=backup_; backup_:='max_gc_short$ h:=errorset({'shorten_pdes,mkquote a1},nil,nil) where !*protfg=t; kernlist!*:=kernlist!*bak$ kord!*:=kord!*bak; erfg!*:=nil; max_gc_counter:=bak; backup_:=bakup_bak; return if (errorp h) or (caar h=nil) then nil else << mi:=caar h; newp:=cdar h; h:=nil; p1:=0; for each pc in cdr newp do p1:=p1+get(pc,'terms); mi:=(<<h:=for each pc in car newp collect if zerop pc then <<nequ_:=add1 nequ_;nil>> else mkeq(pc,fctsort union(get(caddr mi,'fcts), get(cadddr mi,'fcts)), vl,allflags_,t,list(0),nil,pdes); for each pc in h do if pc then p1:=p1-get(pc,'terms); h >> . cdr newp); if print_ then << if tr_short then << for each h in cdr newp do <<write h,": "$typeeq h>>$ for each h in car mi do if null h then <<write "This gives identity 0=0."$terpri()>> else <<write h,": "$typeeq h>>$ >>$ write "shortening by ",p1," term"$ if p1 neq 1 then write"s"$ terpri()$ >>; for each pc in cdr newp do drop_pde(pc,nil,nil); mi >> end$ %------------------- symbolic procedure is_algebraic(p)$ % checks whether the leading derivative is algebraic % if true and if lex_fc:=nil then all allvar functions turn up % only algebraically begin scalar h; h:=get(p,'derivs)$ return if null h then t else <<h:=caar h; if (pairp h) and (cdr h) then nil else t >> end$ %------------------- symbolic procedure shorten_pdes(des)$ begin scalar mi,p1,p1rl,p1le,pc,pcc,newp, l0,l1,l2,l3,l4,version,p1_is_alg$ %,valcp$ if pairp des and pairp cdr des then << version:=1; repeat << % find the pair of pdes not yet reduced with each other % with the lowest product of their number of terms % printlength's mi:=nil; pc:=des; while cdr pc do << p1:=car pc;pc:=cdr pc; if flagp(p1,'to_eval) %and % ((get(p1,'terms)>1) or << % valcp:=get(p1,'val); % if car valcp='!*sq then valcp:=reval valcp; % freeof(valcp,'PLUS) % >>) then << p1rl:=get(p1,'rl_with); p1le:=get(p1,'terms); l1:=get(p1,'derivs); l0:=length l1; if struc_eqn then p1_is_alg:=is_algebraic(p1)$ pcc:=pc; while pcc do if flagp(car pcc,'to_eval ) and % ((get(car pcc,'terms)>1) or << % valcp:=get(car pcc,'val); % if car valcp='!*sq then valcp:=reval valcp; % freeof(valcp,'PLUS) % >>) and ((not member(car pcc, p1rl )) or (not member(p1 ,get(car pcc,'rl_with))) ) and ((null struc_eqn) or p1_is_alg or (is_algebraic(car pcc) and (p1le=get(car pcc,'terms)) ) ) and <<l2:=get(car pcc,'derivs)$ if version=1 then << l3:=length(setdiff(l1,l2))$ if (l0>l3 ) and % necessary requirement (( null mi ) or (( car mi) > l3) ) then t else nil >> else << l3:=length(setdiff(l1,setdiff(l1,l2)))$ l4:=length(union(l1,l2))$ if (l3>0) and ((null mi) or ((((car mi)*l4) > ((cadr mi)*l3)))) then t else nil >> >> then <<mi:=list(l3,l4,p1,car pcc); if ((version = 1) and (l3= 0)) or ((version neq 1) and (l3=l4)) then <<pcc:=nil;pc:={nil}>> else pcc:=cdr pcc >> else pcc:=cdr pcc; >> >>$ if mi then << newp:=shorten(caddr mi,cadddr mi); if null newp then add_rl_with(caddr mi,cadddr mi) >> >> until (null mi) or newp; % if not possible then already returned with nil >>; return (mi . newp) end$ %------------------- symbolic procedure partition_1(l,la)$ % l is an equation, % returning (l1 . l2) where % l1=partitioning of equation l into ((lpow1.lc1),(lpow2.lc2),...) % l2=(lpow1,lpow2,...) % This works currently only for l that are linear in elem. of la begin scalar l1,l3; l:=reorder !*a2f l; while pairp l and member(l3:=car lpow l,la) do << l1:=cons((l3 . !*f2a lc l), l1)$ l:= red l; >>; return if l then (append(l1,list(1 . !*f2a l)) . append(la,list(1))) % inhomogeneous case else (l1 . la) % homogeneous case end$ %------------------- symbolic procedure partition_2(de,l)$ % dropping from de all parts that can not be matched by the other % equation, a list of ftem-functions and their derivatives from % the other equation is l begin scalar newde,dropped,n; % dropped is the number of terms that can not be matched and % which are therefore dropped dropped:=0$ while de do << n:=no_of_terms cdar de$ if member(caar de,l) then newde:=cons(cons(n,car de),newde) else dropped:=dropped+n; de:=cdr de >>; return (dropped . newde) end$ %------------------- symbolic procedure strip(d)$ begin scalar h; d:= if not pairp d then list d else if car d='QUOTIENT then cadr d else if car d = 'PLUS then cdr d else list(d)$ return for each h in d collect !*a2f h end$ %------------------- symbolic procedure shorten(de1,de2)$ % shorten the two pdes with each other % returns a dotted pair, where car is a list of the values of new pdes % and cdr is a list of names of pdes to be dropped begin scalar a,b,h,r,s,cp,l1,l2,l1ul2,l1ml2,l2ml1,l1il2,oldorder, de1p,de2p,termsof1,termsof2,flip,n1,n2,ql,maxcancel, take_first,non_linear,homo,de2pnew,tr_short_local; non_linear:=t; % take_first:=t; % "=t is not so radical, --> eqn.s are longer and in total it is slower % tr_short_local:=t; if tr_short_local then deprint list(get(de1,'val),get(de2,'val))$ if homogen_ and (1=car get(de1,'hom_deg) ) and (1=car get(de2,'hom_deg) ) and ((cadr get(de1,'hom_deg)) neq (cadr get(de2,'hom_deg)) ) then homo:=t; if non_linear and null homo then << a:=sort_partition(de1,nil,get(de1,'fcts),nil)$ b:=sort_partition(de2,nil,get(de2,'fcts),nil)$ if tr_short_local then << write"a=",a$ terpri()$ write"b=",b$ terpri()$ >>; de1p:=nil; de2p:=nil; for each h in a do << s:=car h; cp:=b; % Does s occur in b? while cp and (s neq caar cp) do cp:=cdr cp; if cp then << r:=if (pairp s) or (numberp s) then gensym() else s; %--- dropping the ftem-depending factors once at the beginning de1p:=cons(cons(cadr h, cons(r, reval list('QUOTIENT, if cadr h>1 then cons('PLUS,caddr h) else caaddr h, s) )), de1p); de2p:=cons(cons(cadar cp, cons(r, reval list('QUOTIENT, if cadar cp>1 then cons('PLUS,caddar cp) else car caddar cp, s) )), de2p); % %--- not dropping the ftem-depending factors % de1p:=cons(cons(cadr h,cons(r,if cadr h>1 then cons('PLUS,caddr h) % else caaddr h )),de1p); % de2p:=cons(cons(cadar cp,cons(r,if cadar cp>1 then cons('PLUS,caddar cp) % else car caddar cp )),de2p); if tr_short_local then << write"de1p=",de1p$terpri()$ write"de2p=",de2p$terpri()$ >> >> >> >> else << de1p:=get(de1,'val)$ de2p:=get(de2,'val)$ if homo then << % multiplication with flin_ functions is forbidden a:=get(de1,'derivs)$ h:=nil$ while a do << if not freeoflist(car a,flin_) then h:=cons(car a,h); a:=cdr a >> >> else h:=get(de1,'derivs)$ l1:=for each a in h collect if length car a = 1 then caar a else cons('DF,car a)$ % all derivs of de1 if homo then << % multiplication with flin_ functions is forbidden a:=get(de2,'derivs)$ h:=nil$ while a do << if not freeoflist(car a,flin_) then h:=cons(car a,h); a:=cdr a >> >> else h:=get(de2,'derivs)$ l2:=for each a in h collect if length car a = 1 then caar a else cons('DF,car a)$ % all derivs of de2 l1ml2:=setdiff(l1,l2); % l1 - l2 l2ml1:=setdiff(l2,l1); % l2 - l1 l1il2:=setdiff(l1,l1ml2); % intersection l1ul2:=union(l1,l2); % union if tr_short_local then << write"before substitution:"$terpri()$ write"l1=",l1$ terpri()$ write"l2=",l2$ terpri()$ write"de1p=",de1p$terpri()$ write"de2p=",de2p$terpri()$ write"l1ml2=",l1ml2$terpri()$ write"l2ml1=",l2ml1$terpri()$ write"l1il2=",l1il2$terpri()$ write"l1ul2=",l1ul2$terpri()$ >>; % substituting derivatives by a new variable to become kernels for each a in l1ml2 do if pairp a then << b:=gensym()$ l1:=subst(b,a,l1)$ l1ul2:=subst(b,a,l1ul2)$ de1p:=subst(b,a,de1p) >>$ for each a in l2ml1 do if pairp a then << b:=gensym()$ l2:=subst(b,a,l2)$ l1ul2:=subst(b,a,l1ul2)$ de2p:=subst(b,a,de2p) >>$ for each a in l1il2 do if pairp a then << b:=gensym()$ l1:=subst(b,a,l1)$ l2:=subst(b,a,l2)$ l1ul2:=subst(b,a,l1ul2)$ de1p:=subst(b,a,de1p)$ de2p:=subst(b,a,de2p) >>$ if tr_short_local then << write"after substitution:"$terpri()$ write"l1=",l1$ terpri()$ write"l2=",l2$ terpri()$ write"de1p=",de1p$terpri()$ write"de2p=",de2p$terpri()$ write"l1ml2=",l1ml2$terpri()$ write"l2ml1=",l2ml1$terpri()$ write"l1il2=",l1il2$terpri()$ write"l1ul2=",l1ul2$terpri()$ >>; %--- writing both equations as polynomials in elements of l1ul2 oldorder:=setkorder l1ul2; de1p:=partition_1(de1p,l1); l1:=cdr de1p; de1p:=car de1p; de2p:=partition_1(de2p,l2); l2:=cdr de2p; de2p:=car de2p; setkorder oldorder; %--- l1,l2 can now have the element 1 in case of inhomogeneous de's l1ul2:=nil; l1il2:=nil; %--- Partitioning each equation into 2 parts, one part that can %--- be matched by the other equation and one that can not. % de1p:=partition_2(de1p,l2)$ dropped1:=car de1p; de1p:=cdr de1p; % de2p:=partition_2(de2p,l1)$ dropped2:=car de2p; de2p:=cdr de2p; de1p:=cdr partition_2(de1p,l2)$ de2p:=cdr partition_2(de2p,l1)$ >>$ if (null de1p) or (null de2p) then return nil; termsof1:=no_of_terms get(de1,'val)$ termsof2:=no_of_terms get(de2,'val)$ if tr_short_local then << write"---------"$terpri()$ write"de1:",de1," with ",termsof1," terms"$terpri()$ a:=de1p; while a do << write "caar =",caar a;terpri()$ write "cadar=",cadar a;terpri()$ write "cddar=", algebraic write lisp cddar a;terpri()$ a:=cdr a; >>;terpri()$ write"de2:",de2," with ",termsof2," terms"$terpri()$ a:=de2p; while a do << write "caar =",caar a;terpri()$ write "cadar=",cadar a;terpri()$ write "cddar=",algebraic write lisp cddar a;terpri()$ a:=cdr a; >>;terpri()$ >>; % One can do a stronger restriction: The maximum that can be % canceled is sum of min of terms of the de1p,de2p sublists % corresponding to the coefficients of different ftem functions/deriv. a:=de1p; b:=de2p; n2:=nil; while a do << n1:=if (caar a)<(caar b) then caar a else caar b; % n1 is min of terms of the coefficients of the same ftem function/der. n2:=cons(2*n1,n2); a:=cdr a; b:=cdr b; >>$ % maxcancel is the maximal number of cancellations in all the % remaining runs of short depending on the current run. maxcancel:=list(0); n1:=0; while n2 do << n1:=n1+car n2; n2:=cdr n2; maxcancel:=cons(n1,maxcancel); >>; if ((car maxcancel)<termsof1) and ((car maxcancel)<termsof2) then return nil; if homo and (cadr get(de1,'hom_deg)<cadr get(de2,'hom_deg)) then flip:=nil else if homo and (cadr get(de1,'hom_deg)>cadr get(de2,'hom_deg)) then flip:= t else if (termsof1<termsof2) or (struc_eqn and (n1=n2) and (null is_algebraic(de2)) ) then flip:=nil else flip:=t; if flip then << a:=de1p; de1p:=de2p; de2p:=a; n1:=termsof2; n2:=termsof1 >> else << n1:=termsof1; n2:=termsof2 >>; if (n1=1) and (length de1p = 1) and ((atom cddar de1p) or (caddar de1p neq 'PLUS)) then << % one equation has only a single term which is not a product of sums a:=cadar de1p; % e.g. g0030 b:=de2p; while b and (cadar b neq a) do b:=cdr b; if tr_short_local then << write"one is a 1-term equation"$terpri()$ write"a=",a$terpri()$ write"b=",b$terpri()$ write"de1p.1=",de1p$terpri()$ write"de2p.1=",de2p$terpri()$ >>$ a:=if null b then nil % that term does not turn up in other equation else << % it does turn up --> success de1p:=cddar de1p; de2p:=cddar b; if tr_short_local then << write"de1p.2=",de1p$terpri()$ write"de2p.2=",de2p$terpri()$ >>$ if homo then << if pairp de2p and car de2p='PLUS then de2p:= cdr de2p else de2p:=list de2p; for each a in de2p do << r:=algebraic(a/de1p); % otherwise already successful if freeoflist(algebraic den r,ftem_) then de2pnew:=cons(r,de2pnew) >>; de2p:=if null de2pnew then <<b:=nil;nil>> else if cdr de2pnew then cons('PLUS,de2pnew) else car de2pnew; de1p:=1 >>; de2p % does only matter whether nil or not >> >> else << repeat << % one shortening if tr_short_local then <<write"cadar de1p=",cadar de1p$terpri()>>$ b:=short(ql,strip cddar de1p,strip cddar de2p,n1, 2*(caar de1p),car maxcancel-cadr maxcancel, cadr maxcancel,take_first,homo)$ % take_first:=car b; b:=cdr b; % to activate see end of short() if b then << ql:=car b; a:=cdr b; if a and take_first then << % the result de1p:=!*f2a car a; de2p:=!*f2a cdr a; >> else << de1p:=cdr de1p; de2p:=cdr de2p; >>; maxcancel:=cdr maxcancel; >> else a:=nil; >> until (null b) or % failure a and take_first or % success null de1p$ % the case of exact balance if b and (null take_first) then << % search of the best shortening r:=0; % highest number of saved terms so far de1p:=nil; % numerator of the best quotient so far de2p:=nil; % denominator of the best quotient so far while ql do << s:=cdar ql; while s do << cp:=car s; h:=cdar cp; % nall in short() while cdr cp do << if (cdadr cp)+h>r then << r:=(cdadr cp)+h; de1p:=!*f2a caar cp; if caaadr cp neq 1 then de1p:=reval {'TIMES,de1p,caaadr cp}; de2p:=!*f2a caar ql; if cdaadr cp neq 1 then de2p:=reval {'TIMES,de2p,cdaadr cp}; >>; rplacd(cp,cddr cp) >>; s:=cdr s; >>; ql:=cdr ql; >> >> >>; return if null b or (take_first and null a) then nil else << % numerator and denominator are de1p, de2p %--- computing the shorter new equation if flip then <<a:=get(de2,'val); b:=get(de1,'val)>> else <<a:=get(de1,'val); b:=get(de2,'val)>>$ ql:=if termsof1>termsof2 then de1 else de2; if print_ then << if null car recycle_eqns then n1:=mkid(eqname_,nequ_) else n1:=caar recycle_eqns$ if tr_short then algebraic write"The new equation ",n1," = ", de2p*(if flip then de2 else de1) - de1p*(if flip then de1 else de2)," replaces " >>$ a:=reval list('PLUS, list('MINUS, if de1p=1 then b else list('TIMES,de1p,b)), if de2p=1 then a else list('TIMES,de2p,a) )$ if in_cycle(cons(11,if flip then { get(de2,'printlength),length get(de2,'fcts),de2p, get(de1,'printlength),length get(de1,'fcts),de1p} else { get(de1,'printlength),length get(de1,'fcts),de1p, get(de2,'printlength),length get(de2,'fcts),de2p})) then nil else (list a . list(ql)) >> end$ %------------------- symbolic procedure clean_num(qc,j)$ begin scalar qc1,nall$ return if 2*(cdaar qc)<=j then t else << qc1:=car qc; % the representative and list to proportional factors nall:=cdar qc1; while cdr qc1 do if (cdadr qc1)+nall<=j then rplacd(qc1,cddr qc1) else qc1:=cdr qc1; if qc1=car qc then t else nil % whether empty or not after cleaning >> end$ %-------------------- symbolic procedure clean_den(qc,j)$ begin scalar qcc$ qcc:=qc$ while cdr qc do if clean_num(cdr qc,j) then rplacd(qc,cddr qc) else qc:=cdr qc$ return null cdr qcc % Are there any numerators left? end$ %-------------------- symbolic procedure short(ql,d1,d2,n1,n1_now,max_save_now, max_save_later,take_first,homo)$ begin % d1,d2 are two subexpressions of two expressions with n1,n2 terms % ql is the list of quotients % drp is the number of terms dropped as they can not cancel anything % dne is the number terms of d1 already done, including those dropped % mi is the minimum of n1,n2 % homo=t then non-linear equations --> must check that d2 is not % multiplied with ftem_ dependent factor scalar nall,d1cop,d2cop,m,j,e1,q,qq,qc,dcl,nu,preqc,ldcl,lnu,mi,tr_short_local; %tr_short_local:=t; mi:=n1; m:=0; nall:=0; d1cop:=d1; % n1_now is the maximum number of terms cancelling each other % in this run of short based on 2*(number of remaining terms of d1 % still to check). % max_save_now is the maximum number of cancellations based % on 2*min(terms of d1, min terms of d2) j:=if n1_now<max_save_now then n1_now else max_save_now$ % The following j-value is the minimal number of cancellations % of a quotient by now in order to lead to a reduction. % mi is the minimal umber of cancelled terms at the end = number % of terms of the shorter equation. % max_save_later is the maximal number of cancelling terms in all % later runs of short. j:=mi-j-max_save_later$ repeat << % for each term of d1 n1_now:=n1_now-2; e1:=car d1cop; d1cop:=cdr d1cop; d2cop:=d2; while d2cop and (nall+m<=n1) do << % for each term of d2 q:=cancel(e1 ./ car d2cop); % otherwise already successful d2cop:=cdr d2cop; %--- dropping a numerical factors dcl:=cdr q; % dcl is the denominator of the current quotient if numberp dcl then <<ldcl:=dcl;dcl:=1>> else << ldcl:=dcl; repeat ldcl:=lc ldcl until numberp ldcl$% or car ldcl = '!:RN!:$ dcl:=car cancel(dcl ./ ldcl) >>; nu:=car q; % nu is the numerator of the current quotient if numberp nu then <<lnu:=nu;nu:=1>> else if homo and not freeoflist(nu,ftem_) then nu:=nil else << lnu:=nu; repeat lnu:=lc lnu until numberp lnu$% or car ldcl = '!:RN!:$ nu:=car cancel(nu ./ lnu) >>; if (lnu>1000000000) or (ldcl>1000000000) then if tr_short then << write" Num. factors grew too large in shortening."$ terpri() >> else else if nu then << % - ql is a list of denominator classes: (dcl1 dcl2 dcl3 ...) % - each denominator class dcli is a dotted pair (di . nclist) where % - di is the denominator and % - nclist is a list of numerator classes. % Each numerator class is a list with % - first element: (ncl . n) where ncl is the numerator % up to a rational numerical factor and n is the number of % occurences of ncl (up to a rational numerical factor) % - further elements: (nfi . ni) where nfi is the numerical % proportionality factor and ni the number of occurences % of this factor %---- search for the denominator class qc:=ql; while qc and (dcl neq caar qc) do qc:=cdr qc; if null qc then % denominator class not found if j <= 0 then % add denominator class, here nall,m are not % assigned as it would only play a role if % one equation had only one term but that % is covered as special case ql:=cons((dcl . list(list((nu . 1),((lnu . ldcl) . 1)))), ql) else % too late to add this denominator else << % denominator class has been found %---- now search of the numerator class qc:=cdar qc; % qc is the list of numerator classes nclist while qc and (nu neq caaar qc) do <<preqc:=qc; qc:=cdr qc>>; if null qc then % numerator class not found if j leq 0 then % add numerator class rplacd(preqc,list(list((nu . 1),((lnu . ldcl) . 1))) ) else % too late to add this numerator else <<% numerator class found nall:=cdaar qc + 1; % increasing the total number of occur. rplacd(caar qc,nall); %---- now search for the numerical factor qq:=(lnu . ldcl); qc:=cdar qc; while qc and (qq neq caar qc) do <<preqc:=qc;qc:=cdr qc>>; if null qc then << % numerical factor not found m:=1; rplacd(preqc,list((qq . 1))) >> else << m:=add1 cdar qc$ rplacd(car qc,m) >> >> % numerator class found >> % denominator class found >> % not (homo and ftem_ - dep. factor for d2) >>$ % all terms of d2 j:=if n1_now<max_save_now then n1_now else max_save_now$ j:=mi-j-max_save_later$ if j>0 then << while ql and clean_den(car ql,j) do ql:=cdr ql; if ql then << qc:=ql; while cdr qc do if clean_den(cadr qc,j) then rplacd(qc,cddr qc) else qc:=cdr qc >> >>; if tr_short_local then << terpri();write length ql," denominators"; >>; % If there is only one quotient left and no new one can be added % (because of j>0) then take_first:=t % The following lines need only be un-commented but a test % showed no speed up, only slight slowing down % % if (null take_first) and % (j > 0) and % no new quotients will be added % ql and (null cdr ql) and % only one denominator class % (null cddar ql) and % only one numerator class in cdar ql % % the numerator class is cadar ql % (1=cdar cadar ql) then take_first:=t >> % all terms of d1 until (null d1cop) or % everything divided (take_first and (nall+m>n1)) or % successful: saving > cost ((j > 0) and (null ql))$ % all quotients are too rare --> end return % cons(take_first, if ((j > 0) and (null ql)) then nil else if m+nall<=mi then (ql . nil) else (ql . q) % ) end$ % of short symbolic procedure drop_lin_dep(arglist)$ % drops linear dependent equations begin scalar pdes,tr_drop,p,cp,incre,newpdes,m,h,s,r,a,v, vli,indp,indli,conli,mli,success$ % the pdes are assumed to be sorted by the number of terms, % shortest come first % vli is the list of all `independent variables' v in this lin. algebra % computation, i.e. a list of all different products of powers of % derivatives of ftem functions and constants % format: ((product1, v1, sum1),(product2, v2, sum2),...) % where sumi is the sum of all terms of all equations multiplied % with the multiplier of that equation % indli is a list marking whether equations are necessarily lin % indep. because they involve a `variable' v not encountered yet % mli is the list of multipliers of the equations pdes:=car arglist$ % tr_drop:=t$ if pdes and cdr pdes then << while pdes do << p:=car pdes; pdes:=cdr pdes; newpdes:=cons(p,newpdes); m:=gensym()$ a:=sort_partition(p,nil,get(p,'fcts),nil); if tr_drop then <<write "new eqn:"$prettyprint a; write"multiplier for this equation: m=",m$terpri() >>$ indp:=nil; for each h in a do << s:=car h; % Does s occur in vli? % Adding the terms multiplied with the multiplier to the corresp. sum cp:=vli; while cp and (s neq caar cp) do cp:=cdr cp; if tr_drop then << write"searched for: s=",s$terpri(); if cp then write"found: car cp=",car cp$terpri()$ >>$ if cp then <<r:=cadar cp; incre:=reval {'QUOTIENT, {'TIMES,m,r, if cadr h>1 then cons('PLUS,caddr h) else caaddr h}, s}; rplaca(cddar cp,cons(incre,caddar cp)) >> else <<r:=if (pairp s) or (numberp s) then gensym() else s; indp:=s; incre:=reval {'QUOTIENT, {'TIMES,m,r, if cadr h>1 then cons('PLUS,caddr h) else caaddr h}, s}; vli:=cons({s,r,{incre}},vli) >>; if tr_drop then << write"corresponding symbol: r=",r$terpri()$ write"upd: incre=",incre$terpri()$ write"vli="$prettyprint vli >>$ >>; mli:=cons(m,mli); indli:=cons(indp,indli) >>$ % Formulating a list of conditions while vli do << v:=caddar vli; vli:=cdr vli; conli:=cons(if cdr v then cons('PLUS,v) else car v, conli) >>; % Now the investigation of linear independence pdes:=nil; % the new list of lin. indep. PDEs while cdr newpdes do << if tr_drop then << terpri()$ if car indli then write"lin. indep. without search of ",car newpdes, " due to the occurence of ",car indli else write"lin. indep. investigation for ",car newpdes$ >>; if car indli then pdes:=cons(car newpdes,pdes) else << s:=cdr solveeval {cons('LIST,subst(1,car mli,conli)),cons('LIST,cdr mli)}; if s then <<drop_pde(car newpdes,nil,nil)$ % lin. dep. success:=t$ if print_ then << terpri()$ write"Eqn. ",car newpdes, " has been dropped due to linear dependence."$ >> >> else <<pdes:=cons(car newpdes,pdes); % lin. indep. if tr_drop then << terpri()$ write"Eqn. ",car newpdes," is lin. indep."$ >> >>; >>; newpdes:=cdr newpdes; indli:=cdr indli; conli:=subst(0,car mli,conli); mli:=cdr mli >>; pdes:=cons(car newpdes,pdes) >>; return if success then list(pdes,cadr arglist) else nil end$ symbolic procedure find_1_term_eqn(arglist)$ % checks whether a linear combination of the equations can produce % an equation with only one term if not lin_problem then nil else begin scalar pdes,tr_drop,p,cp,incre,m,h,s,r,a,v, vli,indp,indli,conli,mli,mpli,success, sli,slilen,maxlen,newconli,newpdes,newp,fl,vl$ %tr_drop:=t$ if tr_drop then terpri()$ pdes:=car arglist$ newpdes:=pdes$ %--------------------------------- % if struc_eqn then << % cp:=pdes; % while cp do << % if is_algebraic(car cp) then r:=cons(car cp,r) % else s:=cons(car cp,s); % cp:=cdr cp % >>; % r:=nil; % s:=nil; % >>$ % Drop all PDEs which have at least two derivs which no other have %--------------------------------- if pdes and cdr pdes then << while pdes do << p:=car pdes; pdes:=cdr pdes; m:=gensym()$ if tr_drop then <<terpri()$write"multiplier m=",m$terpri()>>$ a:=sort_partition(p,nil,get(p,'fcts),nil); for each h in a do << s:=car h; % Does s occur in vli? % Adding the terms multiplied with the multiplier to the corresp. sum cp:=vli; while cp and (s neq caar cp) do cp:=cdr cp; if tr_drop then << write"searched for: s=",s$terpri(); if cp then <<write"found: car cp=",car cp$terpri()$>> >>$ if cp then <<r:=cadar cp; incre:=reval {'QUOTIENT, {'TIMES,m,%r, if cadr h>1 then cons('PLUS,caddr h) else caaddr h}, s}; rplaca(cddar cp,cons(incre,caddar cp)) >> else <<r:=if (pairp s) or (numberp s) then gensym() else s; indp:=s; incre:=reval {'QUOTIENT, {'TIMES,m,%r, if cadr h>1 then cons('PLUS,caddr h) else caaddr h}, s}; vli:=cons({s,r,{incre}},vli) >>; if tr_drop then << write"corresponding symbol: r=",r$terpri()$ write"upd: incre=",incre$terpri()$ write"vli="$prettyprint vli >>$ >>; mli:=cons(m,mli); mpli:=cons((m . p),mpli); indli:=cons(indp,indli) >>$ % Formulating a list of conditions while vli do << sli:=cons(caar vli,sli); v:=caddar vli; vli:=cdr vli; conli:=cons(if cdr v then cons('PLUS,v) else car v, conli) >>; % Now the investigation of the existence of equations with only one % term slilen:=length sli; mli :=cons('LIST, mli); conli:=cons('LIST,conli); if tr_drop then << write"sli=",sli$terpri()$ algebraic(write"mli=",mli)$ algebraic(write"conli=",conli)$ write"mpli=",mpli$terpri()$ >>; for h:=1:slilen do << % for each possible single term newp:=car sli;sli:=cdr sli; pdes:=newpdes; while pdes and ((get(car pdes,'terms)>1) or (not zerop reval {'DIFFERENCE,get(car pdes,'val),newp})) do pdes:=cdr pdes; if null pdes then << cp:=conli; for s:=1:h do cp:=cdr cp; rplaca(cp,reval {'PLUS,1,car cp}); if tr_drop then << write"h=",h$terpri()$ algebraic(write"new conli=",conli)$ >>; s:=cdr solveeval {conli,mli}; if (null s) and tr_drop then <<write"no success"$terpri()>>$ if s then << % found 1-term equation if null success then for each p in newpdes do << fl:=union(get(p,'fcts),fl); vl:=union(get(p,'vars),vl) >>$ success:=t$ maxlen:=0$ s:=cdar s; % first solution (lin. system), dropping 'LIST % now find the equation to be replaced by the 1-term-equation % find the non-vanishing m in s, such that the corresponding p in % mpli has a maximum number of terms while s do << if caddar s neq 0 then << r:=cadar s; cp:=mpli; while caar cp neq r do cp:=cdr cp; if get(cdar cp,'terms)>maxlen then << p:=cdar cp; % p will be the equation to be replaced m:=r; maxlen:=get(p,'terms); >> >>$ s:=cdr s >>$ % Now replacing the old equation p by the new 1-term equation in conli: r:=0; newconli:=nil$ while conli do << v:=subst(0,m,car conli)$ conli:=cdr conli$ if r=h then << v:=reval {'PLUS,{'TIMES,m,newp},v}$ >>$ newconli:=cons(v,newconli); r:=add1(r) >>$ conli:=reverse newconli$ % the new equation: newp:=mkeq(newp,fl,vl,allflags_,t,list(0),nil,nil); % last argument is nil as no new inequalities can result % if new equation has only one term newpdes:=cons(newp,newpdes); if print_ then << terpri()$ write"The new equation ",newp$ typeeq newp$ write" replaces ",p$ typeeq p$ >>; drop_pde(p,nil,nil)$ newpdes:=delete(p,newpdes); % update of mpli: mpli:=subst(newp,p,mpli)$ if tr_drop then << write"mpli=",mpli$terpri()$ >>; >>; % end of successful find cp:=conli; for s:=1:h do cp:=cdr cp; rplaca(cp,reval {'PLUS,-1,car cp}); >> % if the 1-term PDE is not already known >>$ % for each possible single term >>; return if success then list(newpdes,cadr arglist) else nil end$ endmodule$ end$ % moegliche Verbesserungen: % - auch subtrahieren, wenn 0 Gewinn (Zyklus!) % - kann Zyklus mit decoupling geben % - evtl erst alle Quotienten bestimmen, dann die Heuristik: % . erst wo nur die kleinere Gleichung mit ftem-Funktionen multipliziert % wird % . wo nur die kleinere Gleichung multipliziert wird % . alle Quotienten werden auf Hauptnenner gebracht und der mit der % groessten Potenz mit der die kuerzere Gleichung multipliziert wird, % ist der erste % - Erweiterung auf mehrere Gleichungen % - counter example: % 0 = +a+b+c % 0 = -b +d+e % 0 = -c-d +f % 0 = -a -e-f % combining any 2 gives a longer one % the sum of all 4 is even zero. % - In order not to run into a cycle with decouple one could use % dec_hist_list but that costs memory.