Artifact aeb51c54dfc734314c0b407b623f84747e60e67a6aed0e02925935b5bed9e88a:
- Executable file
r36/src/ineq.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: 31135) [annotate] [blame] [check-ins using] [more...]
module ineq; % Inequalities and linear optimization. % Author: Herbert Melenk <melenk@zib-berlin.dbp.de> % Driver for solving inequalities and inequality systems. % Implemented methods: % % - linear multivariate system % - polynomial/rational univariate inequality and system % Common algebraic interface: % % ineq_solve(<ineq/ineqlist> [,<variable/variablelist>]) create!-package('(ineq linineq liqsimp1 liqsimp2 polineq),'(solve)); load!-package 'solve; % Some routines from solve are needed. fluid '(solvemethods!*); if not memq('ineqseval,solvemethods!*) then solvemethods!* := 'ineqseval!*!* . solvemethods!*; if not get('geq,'simpfn) then <<mkop 'leq; mkop 'geq; mkop 'lessp; mkop 'greaterp>>; if not get('!*interval!*,'simpfn) then <<mkop '!*interval!*; infix !*interval!*; put('!*interval!*,'prtch," .. ") >>; symbolic procedure ineqseval!*!* u; % Interface to solve. (if null w then nil else if w='(failed) then if smemql('(leq geq lessp greaterp),u) then w else nil else w) where w=ineqseval u; symbolic procedure ineqseval!* u; % Interface to ineq_solve. (if null w or w='(failed) then car u else w) where w=ineqseval u; put('ineq_solve,'psopfn,'ineqseval!*); symbolic procedure ineqseval(u); begin scalar s,s1,v,v1,l,w1,w2,err,ineqp,str; integer n; s:=reval car u; s:=if eqcar(s,'list) then cdr s else {s}; if cdr u then <<v:= reval cadr u; v:=if eqcar(v,'list) then cdr v else {v}; >>; % test for linearity, collect variables. l:=t; s1:=for each q in s join if not err then <<if atom q or not memq(car q,'(leq geq lessp greaterp equal)) then err:=t else << if not(car q eq 'equal) then ineqp := t; n:=n#+1; str := str or memq(car q,'(lessp greaterp)); w1:=simp cadr q; w2:=simp caddr q; v1:=union(v1,solvevars{w1,w2}); if not domainp denr w1 or not domainp denr w2 then l:=nil; {numr w1,denr w1,numr w2,denr w2} >> >>; if err or not ineqp then return nil; if null v then v:=v1; l := l and not nonlnrsys(s1,v); if length v1 > length v or not subsetp(v,v1) or not l and cdr v1 then return '(failed); % Too many indeterminates in inequality system; % if not l and n#>1 then % return '(failed); % Nonlinear system not implemented. if l and str then return '(failed); % No strict linear system. return if l then linineqseval u else polineqeval u; end; endmodule; module linineq; % Linear inequalities and linear optimization. % Author: Herbert Melenk <melenk@zib-berlin.dbp.de> % Version 1 January 1990 % Version 1.1 February 1990 % added parameter "record=t" % Version 2 May 1991 % added Branch-and-Bound for Integer Prgramming % Version 3 Dec 1994 % added formal simplifier for MAX/MIN expressions. % Changed "inf" to "infinity". % Operator linineq_solve new. % Version 4 Jan 95 % use polytope points for the simplification of MAX/MIN % expressions. % % Solution of linear inequalities & equations with numerical % coefficients. % % Fourier(1826) /Motzkin(1936): George. B. Dantzig, % Linear Programming and Extensions. put('linineq,'psopfn, function (lambda(u); rederr "USE simplex (package linalg) instead")); global '(!*trlinineq !*trlinineqint !*prlinineq); switch trlinineq,prlinineq,trlinineqint; fluid '(linineqinterval!* linineqrecord!*); fluid '(!*ineqerr); % error code symbolic procedure linineqeval(u); % Interface for algebraic mode. begin scalar prob,equa,requa,vars,oldorder,res,u1,x,y,p,e,msg; scalar direction,rec,linineqrecord!*,r,intvars,w1,w2,op; msg := !*prlinineq or !*trlinineq; !*ineqerr :=nil; u1 := getrlist reval car u; u := cdr u; if u and eqcar(x:= reval car u,'list) then <<vars := getrlist x; u:=cdr u>>; while u do <<x := reval car u; u := cdr u; if eqcar(x,'equal) and ((cadr x='record and (rec:=t)) or (cadr x='int and (intvars:=getrlist caddr x))) then t else <<!*ineqerr:=2; typerr(x,"illegal parameter")>>; >>; x := nil; for each u in vars do <<u := reval u; if eqcar(u,'equal) then if member(caddr u,'(min max)) then <<direction:=(cadr u. caddr u) . direction; u := cadr u>> else <<!*ineqerr:=2; rederr "illegal form in 2nd parameter">>; if smember(u,u1) and not member(u,x) then x:=u.x>>; x := vars := reversip x; while u1 do <<u := reval car u1; u1:=cdr u1; if not pairp u or not (car u memq '(geq leq equal)) then <<!*ineqerr:=2; typerr(u,"inequality")>>; op:=car u; w1:=reval cadr u; w2:=reval caddr u; if op='geq then if smemq('infinity,w2) then nil else if eqcar(w2,'max) then for each q in cdr w2 do u1:=append(u1,{{'geq,w1,q}}) else prob := (simp w1 . simp w2) . prob else if op='leq then if smemq('infinity,w2) then nil else if eqcar(w2,'min) then for each q in cdr w2 do u1:=append(u1,{{'leq,w1,q}}) else prob := (simp w2 . simp w1) . prob else if op='equal then if eqcar(w2,'!*interval!*) then u1:=append(u1,{{'geq,w1,cadr w2},{'leq,w1,caddr w2}}) else equa := (simp w1 . simp w2) . equa else <<!*ineqerr:=1; typerr(u,"inequality")>> >>; % control the linearity for each p in append(equa,prob) do <<if not domainp denr car p or not domainp denr cdr p then<<!*ineqerr:=1; rederr "unable to process nonlinear system">>; vars := linineqevaltest(numr car p, linineqevaltest(numr cdr p,vars)); >>; if msg then <<prin2 "variables:"; prin2t vars;>>; oldorder := setkorder vars; prob := for each p in prob collect (reorder numr car p./denr car p). (reorder numr cdr p./denr cdr p); equa:= for each p in equa collect (reorder numr car p./denr car p). (reorder numr cdr p./denr cdr p); % eliminate variables from equations while equa do <<e := car equa; equa := cdr equa; e := addsq(car e,negsq cdr e); if domainp numr e then <<if numr e then % nonzero constant equated to 0 <<!*ineqerr:=0; rederr "equation part inconsistent">> >> else <<u := list((x := mvar numr e). prepsq(y:=multsq(negf red numr e ./ 1, invsq(lc numr e ./ 1)))); if member(x,intvars) then % Dont eliminate integer variables - represent % equation by double inequality instead. <<x:=simp x; prob:=append({x.y,y.x},prob)>> else << prob := for each p in prob collect subsq(car p,u).subsq(cdr p,u); equa := for each p in equa collect subsq(car p,u).subsq(cdr p,u); requa := append(u,requa); if msg then <<prin2 " ";prin2 x; prin2 " eliminated by equation"; terpri()>>; vars := delete(x,vars); >> >> >>; res := if intvars then linineqint(prob,vars,msg,direction,rec,intvars) else linineq1(prob,vars,msg,direction,rec); % backsubstitution in equations; if null res then return '(list) else if res=t then res:=nil; for each e in requa do <<x := prepsq subsq(y:=simp cdr e,res); res := (car e . x) . res; if rec then <<x := prepsq y; linineqrecord!* := list(x,x) . linineqrecord!*>>; >>; setkorder oldorder; r := if rec then for each p in liqsimp!-maxmin pair(res,linineqrecord!*) collect list('list,list('equal,caar p,cdar p),cadr p,caddr p) else for each p in res collect list('equal,car p,cdr p); return 'list . r; end; % put('linineq_solve,'psopfn,'linineqseval); symbolic procedure linineqseval(u); 'list . reversip for each q in cdr linineqeval append(u,'((equal record t))) collect {'equal,cadr cadr q, if caddr q = cadddr q then caddr q else '!*interval!* . cddr q}; symbolic procedure linineqevaltest(f,v); % Collect the variables in standard form f and control linearity. if domainp f then v else if not(ldeg f=1) then <<!*ineqerr:=1; rederr "unable to process nonlinear system">> else if member(mvar f,v) then linineqevaltest(red f,v) else linineqevaltest(red f,mvar f.v); symbolic procedure linineq0(prob,vars,dir,rec); % Interface for symbolic mode. % Prob is a list (e1,e2,..) of algebraic expressions without % relational operators, which are interpreted as % set of inequalities ei >= 0. They are linear in the % variables vars. % Silent operation: result=nil if the system is inconsistent. begin scalar oldorder,res; linineqrecord!* := nil; oldorder := setkorder vars; prob := for each u in prob collect simp u . (nil ./ 1); res := linineq1(prob,vars,nil,dir,rec); setkorder oldorder; return res; end; symbolic procedure linineqint(prob,vars,msg,dir,rec,intvars); begin scalar x,x0,y,y0,y1,z,w,problems,best,z,z0,zbest,zf,bestr; % test integer variables and adjust order; for each x in vars do if member(x,intvars) then<<w:=x.w;intvars:=delete(x,intvars)>>; if intvars then <<!*ineqerr:=2;typerr('list.intvars,"int variables")>>; intvars:=reversip w; % select primary optimization principle. if dir then<<z:=caar dir;zf:=if cdar dir='max then 1 else -1>>; problems:=list (nil.prob); % macro loop. while problems do <<z0:=caar problems; prob:=cdar problems; problems:=cdr problems; if msg or !*trlinineqint then linineqprint2("=== next integer subproblem",prob); w:=if best and not evalgreaterp({'times,zf,z0},{'times,zf,zbest}) then nil % skip problem with suboptimal bound. else linineq1(prob,vars,msg,dir,rec); if !*trlinineqint then linineqprint3("=== subresult",w); if w and dir then <<% is better than best so far? z0:=cdr assoc(z,w); if best and evalgreaterp({'times,zf,zbest},{'times,zf,z0}) then w:=nil; >>; if w then <<% test feasability; y:=list prob; for each x in intvars do <<x0:=cdr assoc(x,w); if not fixp x0 then % branch and bound <<x:=simp x; y0:=simp{'ceiling,x0}; y1:=simp {'floor,x0}; y:= for each q in y join {(x.y0).q, (y1.x).q}; if msg or !*trlinineqint then <<writepri("branch and bound with",'first); writepri(mkquote{'list,{'geq,x:=prepsq x,prepsq y0}, {'leq,x,prepsq y1}},'last); >>; >>; >>; if cdr y then problems:=append(problems,for each q in y collect z0.q) else <<zbest:=z0; best:=w; bestr:=linineqrecord!*; if !*trlinineqint then prin2t "===> is feasable";>>; >>; % if w % without target dont need additional result. if best and null dir then problems:=nil; >>; % while problems linineqrecord!*:=bestr; return best; end; symbolic procedure linineq1(prob,vars,msg,dir,rec); % Algebraic evaluation of a set of inequalities: % prob is a list of pairs of standard quotients, % (( p1 . q1)(p2 . q2) .. (pn . qn)) % which are interpreted as inequalities: % pi >= qi ; % vars is the list of (linear) variables. % dir the direction of final optimization % rec switch; if t, the record of inequatlities is produced % Result is NIL if the system has no solution; otherwise % the solution has the form of an association list % ((v1 . val1)(v2 . val2) ... (vn . valn)), % where vi are the variables and vali are values in algebraic % form. NIL if the system has no solution. % begin scalar v,vq,lh,rh,x,y,z,prob1,prob2,prob3,prob4,nprob,sw,sol; if null vars then return linineq2(prob,msg); v := car vars; vars := cdr vars; vq := mksq(v,1); if !*trlinineq then linineqprint2(list("next variable:",v,"; initial system:"),prob); prob := linineqnormalize prob; for each p in prob do <<lh := car p; rh := cdr p; % if v appears on the lhs, isolate it if not domainp numr lh and mvar numr lh = v then <<x := invsq(lc numr lh ./ 1); sw := (numr x < 0); lh := multsq(lh,x); rh := multsq(rh,x); rh := addsq(rh,negf red numr lh ./ denr lh); if not sw then prob1 := (vq . rh) . prob1 else prob2 := (rh . vq) . prob2; >>else if domainp numr rh and domainp numr lh then prob4 := (lh . rh) . prob4 else prob3 := (lh . rh) . prob3; >>; if null prob1 and null prob2 and vars then << sol := linineq1(prob,vars,msg,dir,rec); if rec then linineqrecord!* := append(linineqrecord!*,'(((minus infinity) infinity))); return if sol then (v . 0) . sol else nil>>; if !*trlinineq then <<linineqprint2("class 1:",prob1); linineqprint2("class 2:",prob2); linineqprint2("class 3:",prob3); linineqprint2("class 4:",prob4); >>; if rec then << x := for each u in prob1 collect prepsq cdr u; y := for each u in prob2 collect prepsq car u; x := if null x then '(minus infinity) else if null cdr x then car x else 'max . x; y := if null y then ' infinity else if null cdr y then car y else 'min . y; linineqrecord!* := append(linineqrecord!*, list list(x,y)) >>; if not linineq2(prob4,msg) then return nil; nprob := append(prob3, for each x in prob1 join for each y in prob2 collect car y . cdr x); if vars then << if null (sol := linineq1(nprob,vars,msg,dir,rec)) then return nil>> else if not linineq2(nprob,msg) then return nil; % lower bound: x := if null prob1 then nil else linineqevalmax for each p in prob1 collect subsq(cdr p,sol); % upper bound: y := if null prob2 then nil else linineqevalmin for each p in prob2 collect subsq(car p,sol); if (z:=assoc(v,dir)) then z:= cdr z; if msg then <<writepri(" ",'first); writepri(mkquote if x then prepsq x else '(minus infinity),nil); writepri(" <= ",nil); writepri(mkquote v,nil); writepri(" <= ",nil); writepri(mkquote if y then prepsq y else 'infinity,nil); writepri("; ",nil)>>; linineqinterval!* := x . y; if z='min and null x or z='max and null y then <<if msg then writepri( " max/min cannot be resolved",'last); return nil>>; if not(x=y) then if z='min then y:=nil else if z='max then x:=nil; if msg then << writepri( if null x and null y then " completely free: " else if null y then " minimum: " else if null x then " maximum: " else if x=y then " zero length interval: " else " middle: ",nil); >>; if null x and null y then x := 0 else % completely free if null x then x := prepsq y else if null y then x := prepsq x else if sqlessp(y,x) then <<prin2 "system inconsistent:"; prin2 prepsq x; prin2 " not <= "; prin2t prepsq y; return nil>> else x := list('quotient,list('plus,prepsq x,prepsq y),2); x := aeval x; if msg then writepri(mkquote list('equal,v,x),'last); return (v . x) . sol; end; symbolic procedure linineq2(prob,msg); % All variables are elimitated. Control, if the % remaining numerical inequalities are consistent. begin scalar rh,lh; loop: if null prob then return t; lh := caar prob; rh := cdar prob; if not domainp numr rh or not domainp numr lh then <<!*ineqerr:=1; rederr{" non numeric:", rh, lh}; >>; if sqlessp(lh,rh) then <<if msg then <<writepri("system inconsistent: ",'first); writepri(mkquote prepsq lh,nil); writepri(" not >= ",nil); writepri(mkquote prepsq rh,'last);>>; return nil>>; prob := cdr prob; goto loop; end; symbolic procedure linineqnormalize prob; % Normalize system: reform all inequalities such that they have % the canonical form % polynomial >= constant % (canonical restriction: absolute term of lhs=0, % denominator of lhs = 1). % and remove those, which have same lhs, but smaller rhs % (the latter are superfluous). begin scalar r,lh,rh,d,ab,x; for each p in prob do <<lh := car p; rh := cdr p; % arithmetic normalizaion lh := addsq(lh,negsq rh); d := denr lh; lh := numr lh; ab := lh; x := if domainp lh then 1 else lc ab; while not domainp ab do <<x := gcdf(x,lc ab);ab := red ab>>; ab := negf ab; lh := multsq(addf(lh,ab)./1,1 ./ x); rh := multsq(ab ./ 1, 1 ./ x); % removal of redundant elements x := assoc(lh,r); if null x then r:=(lh.rh) . r else if sqlessp(cdr x,rh) then rplacd(x,rh); >>; if !*trlinineq then linineqprint2("normalized and reduced:",r); return r; end; symbolic procedure linineqevalmin u; % Compute the minimum among the list u with sq's. linineqevalmin1(car u,cdr u); symbolic procedure linineqevalmin1(q,u); if null u then q else (linineqevalmin1( if x and !:minusp x then q else car u, cdr u) ) where x=numr addsq(q,negsq car u); symbolic procedure linineqevalmax u; % compute the maximum among the list u with sq's linineqevalmax1(car u,cdr u); symbolic procedure linineqevalmax1(q,u); if null u then q else (linineqevalmax1( if x and !:minusp x then car u else q, cdr u) ) where x=numr addsq(q,negsq car u); symbolic procedure sqlessp(q1,q2); (x and !:minusp x) where x=numr addsq(q1,negsq q2); symbolic procedure liqsimp!-maxmin w; liqsimp2!-maxmin liqsimp1!-maxmin w; endmodule; module liqsimp1; % interval simplifcation level 1 by % inequality propagation. fluid'(liqsimp1bounds!*); symbolic procedure liqsimp1!-maxmin w; % W is a list of forms {x , l , u} where the interval [l,u] % has been assigned to the variable x. l and u may be formal % expressions dominated by an operator MAX or MIN and including % variables of the following intervals. I try to simplify the % bounds as far as possible by computing inequality chains. (liqsimp1!-maxmin1 w) where liqsimp1bounds!*=nil; symbolic procedure liqsimp1!-maxmin1 w; begin scalar x,l,u,r; x:=caaar w; l:=cadar w; u:=caddar w; if cdr w then % bottom reached? << r:=liqsimp1!-maxmin1 cdr w; l:=liqsimp1!-reducecases(l); u:=liqsimp1!-reducecases(u); >>; liqsimp1bounds!* := {x , liqsimp1!-maxmin1leaf(l), liqsimp1!-maxmin1leaf(u)} . liqsimp1bounds!*; return {caar w,l,u} . r; end; symbolic procedure liqsimp1!-reducecases(w); % M=t: upper bound, M=nil: lower bound. begin scalar op,l,tst,d,u; if atom w or not memq(car w,'(max min)) then return w; op:=car w; l:=for each u in cdr w collect u . simp u; % check whether an element is covered by another one. for each e in l do <<tst := nil; for each d in l do if d neq e and not tst then << % Can I prove that (with op=MAX) e <= d? % I compute u=d-e and check the thest u>= 0 down % the bounds for the other variables. u:=subtrsq(cdr d,cdr e); tst:=liqsimp1!-check(u,liqsimp1bounds!*,op); >>; % then delete it. if tst then l:=delete(e,l); >>; % collect the surviving elements. l:=for each u in l collect car u; return if cdr l then op. l else car l; end; symbolic procedure liqsimp1!-check(u,bounds,m); if m='min then liqsimp1!-check1(negsq u,bounds) else liqsimp1!-check1(u,bounds); symbolic procedure liqsimp1!-check1(u,bounds); % On this level I check whether u>=0 is true. if domainp numr u then not minusf numr u or null numr u else if null bounds then nil else if mvar numr u neq caar bounds then liqsimp1!-check1(u,cdr bounds) else begin scalar x,c,r,d,tst,bds; x:=caar bounds; % U = c*x + r, car bounds has lower and upper limits for x; % U >= 0 is then equivalent to % c >= 0 replace x by lower bounds % or % c <= 0 replace x by upper bounds. % d:=!*f2q denr u; c:=!*f2q lc numr u; r:=!*f2q red numr u; if liqsimp1!-check1(c,bounds) then % C>=0 bds:=cadar bounds else bds:=caddar bounds; for each b in bds do tst:=tst or liqsimp1!-check1(addsq(multsq(c,b),r),cdr bounds); return tst; end; symbolic procedure liqsimp1!-maxmin1leaf(q); if q='infinity or q='(minus infinity) then nil else for each w in (if pairp q and car q memq '(max min) then cdr q else {q}) collect simp w; endmodule; module liqsimp2; % interval simplification level2 by % removal of non-tight hyperplanes. fluid'(infinities!*); symbolic procedure liqsimp2!-maxmin w; % W is a list of forms {x.x0 , l , u} where the interval [l,u] % has been assigned to the variable x. l and u may be formal % expressions dominated by an operator MAX or MIN and including % variables of the following intervals. I try to simplify the % bounds as far as possible by computing inequality chains. begin scalar r; infinities!* := {simp 'infinity, simp '(minus infinity)}; w:= for each q in w collect {car q, minmax2ql cadr q, minmax2ql caddr q}; r:= liqsimp2!-maxmin1 w; return for each q in r collect {car q, ql2minmax('max,cadr q),ql2minmax('min,caddr q)}; end; symbolic procedure ql2minmax(m,l); <<l:=for each q in l collect prepsq q; if cdr l then m.l else car l>>; symbolic procedure minmax2ql(l); if pairp l and car l memq '(min max) then for each q in cdr l collect simp q else {simp l}; symbolic procedure liqsimp2!-maxmin1 w; if null w then nil else liqsimp2!-reducecases(car w,liqsimp2!-maxmin1 cdr w); symbolic procedure liqsimp2!-reducecases(w,ll); % ll is alreayd a simplified chain of intervals. begin scalar x,l,u,t1,e1,e2,pts,eqns,y; x:=caar w; l:=cadr w; u:=caddr w; if null cdr l and null cdr u then return w.ll; % If I have more than one inequality in the upper % or lower part, I compute all pairwise crossing point % because these form the new contribution to the edges. % An inequality which has no valid point can be excluded % from the set. I may ignore infinite points because each % line must have at least two points. eqns := for each q in delete(car infinities!*, delete(cadr infinities!*,append(l,u))) collect {q}; % Compute crossing points. t1:=eqns; while t1 and cdr t1 do <<e1 := car t1; t1:= cdr t1; for each e2 in t1 do <<pts :=liqsimp2_mk_edges(x,car e1,car e2,l,u,ll); cdr e1:=append(cdr e1,pts); cdr e2 := append(cdr e2,pts); >>>>; l:=for each q in l join if null (y:=assoc(q,eqns)) or cdr y then {q}; u:=for each q in u join if null (y:=assoc(q,eqns)) or cdr y then {q}; return{car w,l,u}.ll; end; symbolic procedure liqsimp2_mk_edges(x,e1,e2,l,u,ll); % x: current variable, % e1,e2: forms defining an edge contribution in x=e1,x=e2 at their % intersection points. e1,e2 free of x. % l: complete lower bounds for x, % u: complete upper bounds for x, % ll: simplified bounds for the lower variables. begin scalar form,pts,pl; form := subtrsq(e1,e2); pl := liqsimp2_mk_edges1(form,ll); pts := liqsimp2_mk_edges_check(pl,x,e1,l,u); return pts; end; symbolic procedure sfnegativep u; if domainp u then minusf u else if mvar u = 'infinity then sfnegativep lc u else typerr(prepf u,"numeric expression"); symbolic procedure liqsimp2_mk_edges1(f,ll); % check f=0 by substituting the hyperplanes in ll. if null ll and null numr f then '(nil) else if null ll then typerr (prepsq f,"soll nicht vorkommen") else begin scalar fx,fxx,t1,x,l,u,points,pl; x:=caaar ll; l:=cadar ll; u:=caddar ll; ll:=cdr ll; t1 := delete(car infinities!*, delete(cadr infinities!*,append(l,u))); if null t1 then t1:='((nil . 1)); fx:=liqsimp2_mk_edges2(f,x); fxx := '(nil . 1); points:= if null fx then % case 1: f does not depend of x. I must extend all % solutions of f wrt the remaining variables % by all possible edges from the interval bounds for x. <<pl :=liqsimp2_mk_edges1(fxx,ll); for each q in t1 join for each p in pl collect (x . prepsq subsq(q,p)) . p >> else if domainp numr fx then % case 2: f has the solution x=a where a does not depend % of any further variable. I must select those % extensions of x=a which are compatible under the local % inequalities. << pl:=liqsimp2_mk_edges1(fxx,ll); pl := liqsimp2_mk_edges_check(pl,x,fx,l,u); pl>> else % case 3: f depends of x and some more variables. % I compute all possible intrsection points with the % current interval bounds and extend the to solutions % with the remaining variables. for each p in t1 join <<fxx := subtrsq(fx,p); pl :=liqsimp2_mk_edges1(fxx,ll); % and remove points which don't satixfy the restrictions. pl := liqsimp2_mk_edges_check(pl,x,fx,l,u); pl>>; return points; end; symbolic procedure liqsimp2_mk_edges_check(pl,x,fx,l,u); % select those points of pl where sub(x=p,fx) is compatible % with the l and u bounds. Extend the compatible points by % a value for x. for each p in pl join begin scalar fine,x1; fine:=t; x1:=subsq(fx,p); for each l1 in l do if fine and sfnegativep numr subtrsq(x1,subsq(l1,p)) then fine:=nil; for each u1 in u do if fine and sfnegativep numr subtrsq(subsq(u1,p),x1) then fine:=nil; return if fine then {(x.prepsq x1).p}; end; symbolic procedure liqsimp2_mk_edges2(f,x); % solve f(x)=0 for linear standard quotient f. Return % nil if x does not occur in f. if not smemq(x,f) then nil else begin scalar w; w:=(reorder numr f) where kord!*={x}; return quotsq(negf red w./ 1,lc w ./ 1) ; end; % ============= printing ================================ symbolic procedure linineqprint1(text,lh,rh); <<writepri(text,'first); writepri(mkquote prepsq lh,nil); writepri(" >= ",nil); writepri(mkquote prepsq rh,'last)>>; symbolic procedure linineqprint2(text,prob); <<prin2t "--------------------------------"; if atom text then text:={text}; for each u in text do prin2 u; terpri(); writepri(mkquote('list . for each p in prob collect {'geq,prepsq car p,prepsq cdr p}),'last) >>; symbolic procedure linineqprint3(text,res); <<writepri(text,'first); writepri(mkquote('list . for each p in res collect {'equal,car p,cdr p}), 'last); >>; endmodule; module polineq; % Solve univariate polynomial inequality systems; % Author: Herbert Melenk, ZIB Berlin. % All rights reserved. % Method: compute the real roots of all numerators and denominators % and check the intervals between them. global '(!!arbint); if not get('arbreal,'simpfn) then mkop 'arbreal; symbolic procedure polineqeval u; begin scalar w,x; w:=reval car u; if eqcar(w,'list) then w:=for each q in cdr w collect reval q else w:={w}; if cdr u then x:=reval cadr u; if eqcar(x,'list) then if cddr x then typerr(x,"variable") else x:=cadr x; return polineq0(w,x); end; symbolic procedure polineq0(ul,x); begin scalar b,n,d,l,w,wl,op,u,r,s,x,y,z; loop: u:=car ul; ul:=cdr ul; if not pairp u or not((op:=car u) memq '(geq greaterp leq lessp)) then go to typerr; s:= s or op = 'greaterp or op = 'lessp; w:=simp if op='greaterp or op='geq then {'difference,cadr u,caddr u} else {'difference,caddr u,cadr u}; wl := w.wl; y:=(not domainp numr w and mvar numr w) or (not domainp denr w and mvar denr w); % check for a polynomial in a free variable. if null y or x and x neq y or pairp y and (get(car y,'!:rd!:) or get(car y,'opmtch)) then go to typerr; x:=y; n:= append(n,polineq!-realroots(numr w,x)); d:= append(d,polineq!-realroots(denr w,x)); if ul then go to loop; for each y in append(n,d) do if not(y member b) then b:=y.b; if null b then return if polineqcheck(wl,{x . 0}) then {'list,{'equal,x,{'arbreal,!!arbint := !!arbint+1}}} else '(list); b:=sort(b,'evallessp); % Create the intervals; while b do <<if null l then l:= {{{'difference,car b,1},nil,car b}}; z:= if cdr b then {'quotient,{'plus,car b,cadr b},2} else {'plus,car b,1}; l:={z,car b,if cdr b then cadr b}.l; b:=cdr b; >>; % check and collect the intervals; for each v in l do << if polineqcheck(wl,{x.car v}) then r:=(if null cadr v then {if s then 'lessp else 'leq, x, caddr v} else if null caddr v then {if s then 'greaterp else 'geq, x, cadr v} else {'equal,x, '!*interval!*.cdr v}) . r >>; return 'list.r; typerr: rederr("wrong arguments for polynomial inequality solver"); end; symbolic procedure polineqcheck(wl,p); null wl or not minusf numr subsq(car wl,p) and polineqcheck(cdr wl,p); symbolic procedure polineq!-realroots(u,x); % return real roots of u, if possible as rational numbers. if domainp u then nil else for each f in cdr fctrf u join <<f:= car f; if mvar f neq x then rederr "too many variables"; if ldeg f = 1 then {reval{'quotient,prepf negf red f,prepf lc f}} else for each y in cdr realroots list prepf f collect caddr y >>; endmodule; end;