Artifact cef01032e441073a790d0d8f1392ef6e9e04bf1b8dcf608ca4d6dd3e9879d350:
- Executable file
r37/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: 17474) [annotate] [blame] [check-ins using] [more...]
%******************************************************************** module shortening$ %******************************************************************** % Routines for algebraically combining de's to reduce their length % Author: Thomas Wolf % Jan 1998 % symbolic procedure length_reduction_1(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$ if (l1:=shorten_pdes(l,caddr arglist)) then <<for each a in cdr l1 do pdes:=delete(a,pdes)$ 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:=list(pdes,cadr arglist)>> else l:=nil$ %if print_ and !*time then << % write " time : ", time() - cpu, % " ms GC time : ", gctime() - gc," ms " %>>$ return l$ end$ %------------------- symbolic procedure shorten_pdes(des,vl)$ begin scalar mi,h,p1,p1rl,p1le,p2,pc,newp$ if pairp des and pairp cdr des then << 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) then << p1rl:=get(p1,'rl_with); % p1le:=get(p1,'printlength); p1le:=get(p1,'terms); for each p2 in pc do if flagp(p2,'to_eval) and (not member(p2,p1rl)) and ((null mi) or % (car mi > (p1le*get(p2,'printlength)))) % then mi:=list(p1le*get(p2,'printlength),p1,p2) (car mi > (p1le*get(p2,'terms)))) then mi:=list(p1le*get(p2,'terms),p1,p2) >> >>$ if mi then << newp:=shorten(cadr mi,caddr mi); if null newp then add_rl_with(cadr mi,caddr mi); >> >> until (null mi) or newp; % if not possible then already returned with nil if mi then << p1:=0; % for each pc in cdr newp do p1:=p1+get(pc,'length); 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(cadr mi,'fcts), get(caddr mi,'fcts)), vl,allflags_,t); % for each pc in h do if pc then p1:=p1-get(pc,'length); for each pc in h do if pc then p1:=p1-get(pc,'terms); h >> . cdr newp); if print_ 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 "Length reduction is ",p1," term"$ if p1 neq 1 then write"s" >>; for each pc in cdr newp do setprop(pc,nil); >>; >>; return mi 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,l1,l2,l1ul2,l1ml2,l2ml1,l1il2,oldorder, de1p,de2p,dropped1,dropped2,termsof1,termsof2,tr_short, dne,flip,n1,n2,ql,maxcancel; %tr_short:=t; l1:=for each a in get(de1,'derivs) collect if length car a = 1 then caar a else cons('DF,car a)$ % all derivs of de1 l2:=for each a in get(de2,'derivs) collect if length car a = 1 then caar a else cons('DF,car a)$ % all derivs of de2 de1p:=get(de1,'val)$ de2p:=get(de2,'val)$ 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 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 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; termsof1:=no_of_terms get(de1,'val)$ termsof2:=no_of_terms get(de2,'val)$ if tr_short then << write"---------"$terpri()$ write"de1:",de1," with ",termsof1," terms"$terpri()$ write"dropped:",dropped1$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()$ write"dropped:",dropped2$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 ofthe 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 (null de1p) or (null de2p) or (((car maxcancel)<termsof1) and ((car maxcancel)<termsof2) ) then return nil; if termsof1<termsof2 then << dne:=dropped1; n1:=termsof1; n2:=termsof2 >> else << flip:=t; a:=de1p; de1p:=de2p; de2p:=a; dne:=dropped2; n1:=termsof2; n2:=termsof1 >>$ if n1=1 then << % one equation has only a single term a:=cadar de1p; % e.g. g0030 b:=de2p; while b and (cadar b neq a) do b:=cdr b; 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; t >> >> else repeat << % one shortening if tr_short 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)$ if b then << ql:=car b; a:=cdr b; if a 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 or % success null de1p$ % the case of exact balance return if 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 << n1:=mkid(reval eqname_,reval nequ_)$ 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) )$ (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)$ 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 scalar nall,d1cop,d2cop,m,j,e1,q,qq,qc,dcl,nu,preqc,ldcl,lnu,tr_short,mi$ %tr_short:=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 << lnu:=nu; repeat lnu:=lc lnu until numberp lnu$% or car ldcl = '!:RN!:$ nu:=car cancel(nu ./ lnu) >>; % - 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 >>$ % 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 then << terpri();write length ql," denominators"; >> >> % all terms of d1 until (null d1cop) or % everything divided (nall+m>n1) or % successful: saving > cost ((j > 0) and (null ql))$ % all denominators are too rare and are canceled return if ((j > 0) and (null ql)) then nil else if m+nall<=mi then (ql . nil) else (ql . q) end$ % of short endmodule$ %tr length_reduction_1$ %tr shorten_pdes$ %tr shorten$ %tr short$ %tr partition_1,partition_2$ 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 nichtlin. Probleme % - 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.