Artifact b8641d751788fb87f69b4e61a7663f1eb374b36595704fda05c54e95e81ddddd:
- Executable file
r36/src/compact.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: 22088) [annotate] [blame] [check-ins using] [more...]
module compact; % Header module for compact code. % Author: Anthony C. Hearn. % Copyright (c) 1989 The RAND Corporation. All Rights Reserved. create!-package('(compact mv mvmatch reddom compactf comfac), '(contrib compact)); % These smacros are used in more than one module. symbolic smacro procedure mv!-!.!+(u,v); u . v; symbolic smacro procedure mv!-!.!*(u,v); u . v; symbolic smacro procedure mv!-lc u; cdar u; symbolic smacro procedure mv!-lpow u; caar u; symbolic smacro procedure mv!-lt u; car u; symbolic smacro procedure mv!-red u; cdr u; endmodule; module mv; % Operations on multivariate forms. % Author: Anthony C. Hearn. % Copyright (c) 1989 The RAND Corporation. All Rights Reserved. % These smacros are local to this module. symbolic smacro procedure mv!-term!-coeff u; cdr u; symbolic smacro procedure mv!-term!-pow u; car u; symbolic smacro procedure mv!-tpow u; car u; symbolic smacro procedure mv!-tc u; cdr u; symbolic procedure mv!-!+(u,v); if null u then v else if null v then u else if mv!-lpow u= mv!-lpow v then (lambda x; if x=0 then mv!-!+(mv!-red u,mv!-red v) else mv!-!.!+(mv!-!.!*(mv!-lpow u,x), mv!-!+(mv!-red u,mv!-red v))) (mv!-lc u + mv!-lc v) else if mv!-pow!-!>(mv!-lpow u,mv!-lpow v) then mv!-!.!+(mv!-lt u,mv!-!+(mv!-red u,v)) else mv!-!.!+(mv!-lt v,mv!-!+(u,mv!-red v)); symbolic smacro procedure domain!-!*(u,v); u*v; symbolic smacro procedure domain!-!/(u,v); u/v; symbolic procedure mv!-term!-!*(u,v); % U is a (non-zero) term and v a multivariate form. Result is % product of u and v. if null v then nil else mv!-!.!+(mv!-!.!*(mv!-pow!-!+(mv!-tpow u,mv!-lpow v), domain!-!*(mv!-tc u,mv!-lc v)), mv!-term!-!*(u,mv!-red v)); symbolic procedure mv!-term!-!/(u,v); % Returns the result of the (exact) division of u by term v. if null u then nil else mv!-!.!+(mv!-!.!*(mv!-pow!-!-(mv!-lpow u,mv!-tpow v), domain!-!/(mv!-lc u,mv!-tc v)), mv!-term!-!/(mv!-red u,v)); symbolic procedure mv!-domainlist u; if null u then nil else mv!-lc u . mv!-domainlist mv!-red u; symbolic procedure mv!-pow!-mv!-!+(u,v); if null v then nil else mv!-!.!+(mv!-pow!-mv!-term!-!+(u,mv!-lt v), mv!-pow!-mv!-!+(u,mv!-red v)); symbolic procedure mv!-pow!-mv!-term!-!+(u,v); mv!-!.!*(mv!-pow!-!+(u,mv!-term!-pow v), mv!-term!-coeff v); symbolic procedure mv!-pow!-!+(u,v); if null u then nil else (car u+car v) . mv!-pow!-!+(cdr u,cdr v); symbolic procedure mv!-pow!-!-(u,v); if null u then nil else (car u-car v) . mv!-pow!-!-(cdr u,cdr v); symbolic procedure mv!-pow!-!*(u,v); if null v then nil else (u*car v) . mv!-pow!-!*(u,cdr v); symbolic procedure mv!-pow!-minusp u; if null u then nil else car u<0 or mv!-pow!-minusp cdr u; symbolic procedure mv!-pow!-!>(u,v); if null u then nil else if car u=car v then mv!-pow!-!>(cdr u,cdr v) else car u>car v; symbolic procedure mv!-reduced!-coeffs u; % reduce coefficients of u to lowest terms. begin scalar x,y; x := mv!-lc u; y := mv!-red u; while y and x neq 1 do <<x := gcdn(x,mv!-lc y); y := mv!-red y>>; return if x=1 then u else mv!-!/(u,x) end; symbolic procedure mv!-!/(u,v); if null u then nil else mv!-!.!+(mv!-!.!*(mv!-lpow u,mv!-lc u/v),mv!-!/(mv!-red u,v)); % Functions that convert between standard forms and multivariate forms. symbolic procedure sf2mv(u,varlist); % Converts the standard form u to a multivariate form wrt varlist. sf2mv1(u,nil,varlist); symbolic procedure sf2mv1(u,powers,varlist); if null u then nil else if domainp u then list(append(powers,nzeros length varlist) . u) else if mvar u = car varlist % This should be eq, but seems to % need equal. then append(sf2mv1(lc u,append(powers,list ldeg u),cdr varlist), sf2mv1(red u,powers,varlist)) else sf2mv1(u,append(powers,list 0),cdr varlist); symbolic procedure nzeros n; if n=0 then nil else 0 . nzeros(n-1); symbolic procedure mv2sf(u,varlist); % converts the multivariate form u to a standard form wrt varlist. % This version uses addf to fold terms - there is probably a more % direct method. if null u then nil else addf(mv2sf1(mv!-lpow u,cdar u,varlist),mv2sf(cdr u,varlist)); symbolic procedure mv2sf1(powers,cf,varlist); if null powers then cf else if car powers=0 then mv2sf1(cdr powers,cf,cdr varlist) else !*t2f((car varlist .** car powers) .* mv2sf1(cdr powers,cf,cdr varlist)); endmodule; module mvmatch; % Side relation matching against expressions. % Author: Anthony C. Hearn. % Copyright (c) 1991 The RAND Corporation. All Rights Reserved. symbolic procedure mv!-compact(u,v,w); % Compares a multivariate form u with a multivariate form template v % and reduces u appropriately. % Previously, the content was removed from u. However, this does % not work well if the same content is in v. begin scalar x,y; % z; if null u then return mv!-reverse w; % if null w then <<z := mv!-content u; u := mv!-term!-!/(u,z)>> % else z := nzeros length mv!-lpow u . 1; % check first terms. if (x := mv!-pow!-chk(u,v)) and (y := mv!-compact2(u,mv!-pow!-mv!-!+(x,v))) % then return mv!-term!-!*(z,mv!-compact(y,v,w)) then return mv!-compact(y,v,w) % check second terms. else if (x := mv!-pow!-chk(u,mv!-red v)) and not mv!-pow!-assoc(y := mv!-pow!-!+(x,mv!-lpow v),w) and (y := mv!-compact2(mv!-!.!+(mv!-!.!*(y,0),u), mv!-pow!-mv!-!+(x,v))) % then return mv!-term!-!*(z,mv!-compact(y,v,w)) % else return mv!-term!-!*(z,mv!-compact(mv!-red u,v,mv!-lt u . w)) then return mv!-compact(y,v,w) else return mv!-compact(mv!-red u,v,mv!-lt u . w) end; symbolic procedure mv!-pow!-assoc(u,v); assoc(u,v); symbolic procedure mv!-reverse u; reversip u; symbolic procedure mv!-pow!-chk(u,v); % (u := mv!-pow!-!-(caar u,caar v)) and not mv!-pow!-minusp u and u; if v and (u := mv!-pow!-!-(caar u,caar v)) and not mv!-pow!-minusp u then u else nil; symbolic procedure mv!-compact2(u,v); % U and v are multivariate forms whose first powlists are equal. % Value is a suitable multiplier of v which when subtracted from u % results in a more compact expression. begin scalar x,y,z; x := equiv!-coeffs(u,v); z := mv!-domainlist v; y := reduce(x,z); return if y=x then nil else mv!-!+(mv!-coeff!-replace(v,mv!-domainlist!-!-(y,x)),u) end; symbolic procedure mv!-coeff!-replace(u,v); % Replaces coefficients of multivariate form u by those in domain % list v. if null u then nil else if car v=0 then mv!-coeff!-replace(mv!-red u,cdr v) else mv!-!.!+(mv!-!.!*(mv!-lpow u,car v), mv!-coeff!-replace(mv!-red u,cdr v)); symbolic procedure equiv!-coeffs(u,v); if null u then nzeros length v else if null v then nil else if mv!-lpow u = mv!-lpow v then cdar u . equiv!-coeffs(cdr u,cdr v) else if mv!-pow!-!>(mv!-lpow u,mv!-lpow v) then equiv!-coeffs(cdr u,v) else 0 . equiv!-coeffs(u,cdr v); endmodule; module reddom; % Reduction of domain elements. % Author: Anthony C. Hearn. % Copyright (c) 1989 The RAND Corporation. All Rights Reserved. fluid '(mv!-vars!*); global '(!*xxx !*yyy); % switch xxx,yyy; !*xxx := !*yyy := t; % Operations on domain elements. symbolic smacro procedure domain!-!+(u,v); u+v; symbolic smacro procedure domain!-!-(u,v); u-v; symbolic smacro procedure domain!-!*(u,v); u*v; symbolic smacro procedure domain!-divide(u,v); divide(u,v); % Operations on domain element lists. symbolic procedure mv!-domainlist!-!+(u,v); if null u then nil else domain!-!+(car u,car v) . mv!-domainlist!-!+(cdr u,cdr v); symbolic procedure mv!-domainlist!-!-(u,v); if null u then nil else domain!-!-(car u,car v) . mv!-domainlist!-!-(cdr u,cdr v); symbolic procedure mv!-domainlist!-!*(u,v); if null v then nil else domain!-!*(u,car v) . mv!-domainlist!-!*(u,cdr v); % Procedures for actually reducing domain elements. symbolic procedure reduce(u,v); % Reduce domain element list u with respect to an equal length domain % element list v. We assume that v has been reduced to lowest terms. begin scalar weightlist,x; % Look for equal ratios of elements. x := u; if !*yyy then x := reduce!-ratios(x,v); % Define weighting list. weightlist := set!-weights v; % Choose column elimination with lowest weight. if !*xxx then x := reduce!-columns(x,v,weightlist); % Look for a reduction in weight of the expression. if !*xxx then x := reduce!-weights(x,v,weightlist); return x end; symbolic procedure set!-weights v; % Define weights to be associated with the reduction test. % The current definition is pretty naive. begin integer n; % return reversip for each j in v collect (n := n+1) return reversip (0 . for each j in cdr v collect 1) end; symbolic procedure reduce!-ratios(u,v); begin scalar x; if null(x := red!-ratios1(u,v)) then return u; x := mv!-domainlist!-!-(mv!-domainlist!-!*(car x,u), mv!-domainlist!-!*(cdr x,v)); return if zeros u >= zeros x then u else reduce!-ratios(x,v) end; symbolic procedure zeros u; if null u then 0 else if car u = 0 then 1+zeros cdr u else zeros cdr u; symbolic procedure red!-ratios1(u,v); u and (red!-ratios2(cdr u,cdr v,car u,car v) or red!-ratios1(cdr u,cdr v)); symbolic procedure red!-ratios2(u,v,u1,v1); % The remainder check is needed for the example % reduce('(3 0 -3 0 0 0 0),(2 -1 -2 -1 3 -1 1)); begin integer n; return if null u then nil else if (n := u1*car v) = v1*car u and n neq 0 and remainder(gcdn(v1,u1),v1)=0 then red!-lowest!-terms(v1,u1) else red!-ratios2(cdr u,cdr v,u1,v1) end; symbolic procedure red!-lowest!-terms(u,v); begin scalar x; if u<0 then <<u := -u; v := -v>>; x := gcdn(u,v); % We must have x = u from call in red-ratios2. If % not, something is awfully wrong. if x neq u then errach list("red-lowest-terms",u,v); return 1 . (v/x) end; symbolic procedure reduce!-columns(u,v,weightlist); begin scalar w,x,y,z,z1; x := u; y := v; w := (u . red!-weight(u,weightlist)); a: if null x then return car w else if car x=0 or car y=0 then nil else if cdr(z := domain!-divide(car x,car y))=0 then <<z := mv!-domainlist!-!-(u,mv!-domainlist!-!*(car z,v)); z1 := red!-weight(z,weightlist); if red!-weight!-less!-p(z1,cdr w) and not more!-apartp(z . z1,w) then w := (z . z1)>>; x := cdr x; y := cdr y; go to a end; symbolic procedure more!-apartp(u,v); cadr u=2 and cadr u=cadr v and cadar u=0 and cadar v neq 0; symbolic procedure reduce!-weights(u,v,weightlist); begin scalar success,x,y,z; x := red!-weight(u,weightlist); a: y := mv!-domainlist!-!+(u,v); z := red!-weight(y,weightlist); if red!-weight!-less!-p(z,x) then <<success := t; u := y; x := z; go to a>>; if success then return u; b: y := mv!-domainlist!-!-(u,v); z := red!-weight(y,weightlist); if red!-weight!-less!-p(z,x) then <<u := y; x := z; go to b>>; return u end; symbolic procedure red!-weight(u,weightlist); nonzero!-length u . red!-weight1(u,weightlist); symbolic procedure red!-weight1(u,weightlist); if null u then 0 else abs car u*car weightlist + red!-weight1(cdr u,cdr weightlist); symbolic procedure nonzero!-length u; if null u then 0 else if car u=0 then nonzero!-length cdr u else add1 nonzero!-length cdr u; symbolic procedure red!-weight!-less!-p(u,v); if car u=car v then cdr u<cdr v else car u<car v; endmodule; module compactf; % Algorithms for compacting algebraic expressions. % Author: Anthony C. Hearn. % Copyright (c) 1991 The RAND Corporation. All Rights Reserved. fluid '(frlis!* mv!-vars!*); global '(!*trcompact); switch trcompact; % Interface to REDUCE simplifier. put('compact,'simpfn,'simpcompact); symbolic procedure simpcompact u; begin scalar bool; if null u or null cdr u then rerror(compact,1, list("Wrong number of arguments to compact")); if null !*exp then <<rmsubs(); bool := !*exp := t>>; u := errorset!*(list('simpcompact1,mkquote u),nil); if bool then !*exp := nil; if errorp u then rerror(compact,2,"Compact error"); return car u end; symbolic procedure simpcompact1 u; begin scalar v,x,y,w; v := simp!* car u; u := cadr u; if idp u then if eqcar(x := get(u,'avalue),'list) then u := cadr x else typerr(u,"list") else if getrtype u eq 'list then u := cdr u else typerr(u,"list"); u := for each j in u collect << w:=t; if eqcar(j,'equal) or eqcar(j,'replaceby) then << if eqcar(y:=caddr j,'when) then <<w:=compactbool formbool(caddr y,nil,'algebraic); y:=cadr y>>; j:= {'difference,cadr j,y}>>; % propagate free variables. if(y:=compactfmatch2 j) then <<j:=sublis(for each x in y collect x.cadr x,j); j:=sublis(for each x in y collect cadr x.x,j)>>; j.w>>; for each j in u do v := compactsq(v,simp!* car j,cdr j); return v end; symbolic procedure compactbool w; % Reform condtion w for later evaluation and substitution. % Without this reform (list (quote ~)(quote x)) would not % be substituted by subst('(((~ x).y)..)... . if atom w then w else if eqcar(w,'list) and cdr w and cadr w='(quote !~) then {'quote,{'!~,cadr caddr w}} else compactbool car w . compactbool cdr w; % True beginning of compacting routines. symbolic procedure compactsq(u,v,c); % U is a standard quotient, v a standard quotient for equation v=0. % Result is a standard quotient for u reduced wrt v=0. begin if denr v neq 1 then msgpri("Relation denominator",prepf denr v,"discarded", nil,nil); v := numr v; return multsq(compactf(numr u,v,c) ./ 1, 1 ./ compactf(denr u,v,c)) end; symbolic procedure compactf(u,v,c); % U is a standard form, v a standard form for an equation v=0. % C is a condition for applying v. % Result is a standard form for u reduced wrt v=0. begin scalar x; integer n; if !*trcompact then <<prin2t "*** Arguments on entering compactf:"; mathprint mk!*sq !*f2q u; mathprint mk!*sq !*f2q v>>; while x neq u do <<x := u; u := compactf0(u,v,c); n := n+1>>; if !*trcompact and n>2 then <<prin2 " *** Compactf looped ";prin2 n; prin2t " times">>; return u end; symbolic procedure compactf0(u,v,c); begin scalar x,y,w; x := kernels u; y := kernels v; if not smemq('!~,v) then return compactf1(u,v,x,y); for each p in compactfmatch(x,y) do if p and not smemq('!~,w:=sublis(p,c)) and eval w and not smemq('!~,w:=numr subf(v,p)) then u:=compactf1(u,w,x,kernels w); return u; end; symbolic procedure compactfmatch(x,y); % Finds all possible matches between free variables in % kernels of list x and pattern list y, including incomplete, % inconsistent and the empty match. if null x or null y then '(nil) else begin scalar y1,z,r; z:=compactfmatch(x,cdr y); if not smemq('!~,car y) then return z; y1:=car y; y:= cdr y; r:=for each x1 in x join for each w in compactfmatch1(x1,y1) join for each q in compactfmatch(delete(x1,x),sublis(w,y)) collect union(w,q); return union(r,z); end; symbolic procedure compactfmatch1(x,y); if car y = '!~ then {{y.x}} else if pairp x and car x=car y then mcharg(cdr x,cdr y,car y) where frlis!* =nconc(compactfmatch2 y,frlis!*); symbolic procedure compactfmatch2 y; if atom y then nil else if car y = '!~ then {y} else append(compactfmatch2(car y),compactfmatch2(cdr y)); symbolic procedure compactf1(u,v,x,y); begin scalar z; % x := kernels u; % y := kernels v; z := intersection(x,y); % find common vars. if null z then return u; % Unfortunately, it's too expensive in space to generate all perms. % as in this example: % l:={-c31*c21+c32*c22+c33*c23+c34*c24=t1}; % x:= -c31*c21+c32*c22+c33*c23+c34*c24; % compact(x,l); % out of heap space % for each j in permutations z do u := compactf11(u,v,x,y,j); return compactf11(u,v,x,y,z) % return u end; symbolic procedure compactf11(u,v,x,y,z); begin scalar w; if domainp u then return u; y := append(z,setdiff(y,z)); % vars in eqn. x := append(setdiff(x,z),y); % all vars. x := setkorder x; u := reorder u; % reorder expressions. v := reorder v; z := comfac!-to!-poly comfac u; u := quotf(u,z); u := remchkf(u,v,y); w := compactf2(u,mv!-reduced!-coeffs sf2mv(v,y),y); if termsf w < termsf u then u := w; % Now reduce z (required, e.g. for compact(u1*(h0+h1),{h0+h1=z1})) if not kernlp z then <<z := remchkf(z,v,y); w := compactf2(z,mv!-reduced!-coeffs sf2mv(v,y),y); if termsf w < termsf z then z := w>>; u := multf(z,u); setkorder x; u := reorder u; if !*trcompact then <<prin2t "*** Value on leaving compactf11:"; mathprint mk!*sq !*f2q u>>; return u end; symbolic procedure remchkf(u,v,vars); % This procedure returns u after checking if a smaller remainder % results after division by v. It is potentially inefficient, since % we check all the way down the list, term by term. However, the % process terminates when we no longer have any relevant kernels. (if domainp x or null intersection(kernels u,vars) then x else lt x .+ remchkf(red x,v,vars)) where x=remchkf1(u,v); symbolic procedure remchkf1(u,v); begin integer n; n := termsf u; v := xremf(u,v,n); if null v or termsf(v := car v)>=n then return u else if !*trcompact then prin2t "*** Remainder smaller"; return v end; symbolic procedure xremf(u,v,m); % Returns the quotient and remainder of U divided by V, or NIL if % the number of terms in the remainder exceeds M. % The goal is to keep terms u+terms z<=m. % There is some slop in the count, so one must check sizes on % leaving. begin integer m1,m2,n; scalar x,y,z; if domainp v then return list cdr qremd(u,v); m2 := termsf u; a: if m<= 0 then return nil else if domainp u then return list addf(z,u) else if mvar u eq mvar v then if (n := ldeg u-ldeg v)<0 then return list addf(z,u) else <<x := qremf(lc u,lc v); y := multpf(lpow u,cdr x); m := m+m1; z := addf(z,y); m1 := termsf z; m := m-m1+m2; u := if null car x then red u else addf(addf(u,multf(if n=0 then v else multpf(mvar u .** n,v), negf car x)), negf y); m2 := termsf u; m := m-m2; go to a>> else if not ordop(mvar u,mvar v) then return list addf(z,u); m := m+m1; x := xremf(lc u,v,m); if null x then return nil; z := addf(z,multpf(lpow u,car x)); m1 := termsf z; m := m-m1; u := red u; go to a end; symbolic procedure compactf2(u,v,vars); % U is standard form for expression, v for equation. W is ordered % list of variables in v. Result is a compacted form for u. if domainp u then u else if mvar u memq vars then compactf3(u,v,vars) else lpow u .* compactf2(lc u,v,vars) .+ compactf2(red u,v,vars); symbolic procedure compactf3(u,v,vars); begin scalar mv!-vars!*; mv!-vars!* := vars; return mv2sf(mv!-compact(sf2mv(u,vars),v,nil),vars) end; endmodule; module comfac; % Multivariate common factor/content routines. % Author: Anthony C. Hearn. % Copyright (c) 1989 The RAND Corporation. All Rights Reserved. symbolic smacro procedure domain!-gcd(u,v); gcdn(u,v); symbolic smacro procedure domain!-onep u; onep u; symbolic procedure mv!-pow!-zerop u; null u or zerop car u and mv!-pow!-zerop cdr u; symbolic procedure mv!-pow!-gcd(u,v); if null u then nil else min(car u,car v) . mv!-pow!-gcd(cdr u,cdr v); symbolic procedure mv!-content u; % Finds the term that is the content of u. if null u then nil else begin scalar x,y; x := mv!-lc u; y := mv!-lpow u; a: u := mv!-red u; if null u or domain!-onep x and mv!-pow!-zerop y then return mv!-!.!*(y,x); x := domain!-gcd(x,mv!-lc u); y := mv!-pow!-gcd(y,mv!-lpow u); go to a end; endmodule; end;