CAD on maxima
SARAG https://github.com/andrejv/maxima/tree/master/share/contrib/sarag を待ちきれない方のために作ってみました.ただし,lifting のネックである根の分離は realroots に任せ,符号判定も近似的です.
使用例(wxmaxima 等で実行すると見易いと思います)
(vpeds ( [x,a,b],0,0,3,0 ),cnnm(cad ( '( a*x^2+b<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,11,0 ),cnnm(cad ( '( a*x^2+b<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,11,0 ),cnnm(cad ( '( a*x^5+b*x+1<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,3,0 ),cnnm(cad ( '( sqrt(x)<a*x+b ) )) );
(%i1) (vpeds ( [x,a],0,0,11,0 ),cnnm(cad ( '( x*(x-2)*(x-3)*(x-5)*(x-7)=a ) )) ); connecting cells. connecting cells. matrix([a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 2), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,2), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,2) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 3), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,5)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,3), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,3) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 4), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,4), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,4) < a, x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1)]) Evaluation took 2.3600 seconds (2.3600 elapsed) using 540.776 MB.
なお,例によって,コードの大半は結果の簡約に費やされています.
/* 20170324 22:04:46 */ showtime:on$ display2d:false$ prton2():=(prt:print,prt2:print)$ prton():=(prt:print)$ prtoff():=kill(prt,prt2)$ vpeds(v,p,e,d,s):=block([], %vv:v, algexact:true, bsolve:false, dtype:0, if p=1 then (prtoff(),prton()) elseif p=2 then prton2() else (ratprint:off,prtoff()), if e=0 then delv0 (L,vv):=delv0s (L,vv) else delv0 (L,vv):=delv0m (L,vv), if d=0 then getsample (a,b,c):=getsample0 (a,b,c) elseif d=1 then getsample (a,b,c):=getsample1 (a,b,c) elseif d=11 then (dtype:11,getsample (a,b,c):=getsample11 (a,b,c)) elseif d=2 then getsample (a,b,c):=getsample2 (a,b,c) else (getsample (a,b,c):=getsample1 (a,b,c), bsolve:true), if s=1 then dispsample:true else dispsample:false, usersimp:simp, simp:false )$ matchdeclare([aa,bb,cc],true)$ kill("implies")$ infix("implies",60,40)$ (x implies y):=(not(x) or y)$ %P:[]$ getlcP(f):=block([ff:expand(f),d,lc], d:hipow(ff,%v),lc:coeff(ff,%v,d),/*print(%v),*/ if member(%v,listofvars(ff)) then %P:append(%P,[ff]) elseif listofvars(ff)#[] then %C:append(%C,[ff]), if not(constantp(lc)) then (%C:append(%C,[lc]),getlcP(ff-lc*(%v)^d)) )$ gcd2(f,g,v):=block([p1:f,p2:g,p3:1,k], for k:0 unless p3=0 do (p3:remainder (p1,p2,v),p1:p2,p2:p3),factor(content(p1,v)[2]))$ sqfree(f,v):=block([], if member(v,listofvars (f)) then content( divide (f,gcd2 (f,diff(f,v),v),v) [1] ,v) [2] else f )$ sqfreem(L,v):=block([], if is(L=[]) then [] else listify(setify( map(lambda([f],sqfree (f,v)),L) )) )$ resultant1(f,v):=block ([g:sqfree(f,v)], resultant (g,diff (g,v),v ) )$ resultant2(f1,f2,v):=block ( [g1: quotient (f1,gcd2 (f1,f2,v),v), g2: quotient (f2,gcd2 (f1,f2,v),v)], resultant (g1,g2,v ) )$ delv0m(L,v):=block([P:[],PP,i,lp], delv0 (LL,vv):=delv0s (LL,vv), %C:[],%v:v, map(lambda([f],if member(%v,listofvars(f)) then P:append(P,[f]) else %C:append(%C,[f])),L), PP:map("[",P),prt(PP), /* PP:map(lambda([f],(%P:[],getlcP(f),%P)),P),prt(PP), */ map(lambda([f],%C:append(%C,[resultant1(f,%v)])),flatten(PP)), lp:length(PP),/*prt([%C,P]),*/ %C:append(%C,flatten(makelist(makelist(makelist(makelist(resultant2(PP[i][ii],PP[j][jj],%v),ii,1,length(PP[i])),jj,1,length(PP[j])),j,i+1,lp),i,1,lp-1))), listify(setify(delete(false,map(lambda([f],if constantp(f) then false else f),%C)))) )$ delv0s(L,v):=block([P:[],PP,i,lp], %C:[],%v:v, map(lambda([f],if member(%v,listofvars(f)) then P:append(P,[f]) else %C:append(%C,[f])),L), /* PP:map("[",P),print(PP), */ PP:map(lambda([f],(%P:[],getlcP(f),%P)),P),prt(PP), map(lambda([f],%C:append(%C,[resultant1(f,%v)])),flatten(PP)), lp:length(PP),/*prt([%C,P]),*/ %C:append(%C,flatten(makelist(makelist(makelist(makelist(resultant2(PP[i][ii],PP[j][jj],%v),ii,1,length(PP[i])),jj,1,length(PP[j])),j,i+1,lp),i,1,lp-1))), listify(setify(delete(false,map(lambda([f],if constantp(f) then false else f),%C)))) )$ load ("grobner")$ delv(L,vl):=block([p:[sqfreem(L,vl[1])],k], for k:1 thru length(vl)-1 do p:append(p,[sqfreem( delv0(last(p),vl[k]) ,vl[k+1])]), p:append(rest (p,-1), [listify(setify( poly_normalize_list (last (p), [last (vl) ] )))]), reverse(p) )$ factors(f):=block ([g:factor(f)], if atom(g) then g elseif op(g)="+" then g else (if op(g)="/" then g:args(g)[1], if atom(g) then g elseif op(g)="-" then g:-g, if atom(g) then g elseif op(g)="*" then delete(0, map(lambda([e], if constantp(e) then 0 elseif atom(e) then e elseif op(e)="^" then args(e)[1] else e), args(g))) elseif op(g)="^" then args(g)[1] else g ) )$ factorsm(L):=listify(setify(flatten(map(factors,L))))$ delv(L,vl):=block([p:[factorsm(L)],k], for k:1 thru length(vl)-1 do p:append(p,[factorsm( delv0(last(p),vl[k]) )]), p:append(rest (p,-1),[factorsm(last(p))]), reverse(p) )$ realonly:on$ getsample0(sub,p,v):=block([p1,p2,q,k], p1:listify(setify(flatten(map(lambda([f],algsys([f],listofvars(f))),subst(sub[2],p))))), if p1=[] then p2:[0] else (p1:sort(map(rhs,p1),"<"),q:map("[",p1), p2:sort(listify(setify(append(p1,(append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<")), makelist([[0.5*k],[v=p2[k]]],k,1,length(p2)) )$ getsample1(sub,p,v):=block([p1,sf,rt,T:[],S:[],p2,q,k], p1:apply(append,map(lambda([f], (sf: subst(sub[2],f), rt:if constantp (sf) then [] else flatten(solve( sf ,v)), rt:if rt=[] then [] else flatten(map(lambda([e],if is(imagpart(rhs(e))=0) then e else [] ), rt )), prt2([p,sf,v,rt]), rt:if rt=[] then [] else (rt:sort (map (rhs,rt),"<" ), makelist([[f,k ],rt[k]],k,1,length (rt) )) )),p)),prt2("p1:"),prt2(p1), if p1=[] then p2:[0] else (map (lambda ( [e], if not(member (e [2],T )) then (T:append(T,[e[2]]),S:append(S,[e])) ),p1 ), prt2 ("S:"), prt2(S), if S=[] then (p1:[],p2:[0]) else ( S:sort (S,lambda ( [u,t], u[2]<t[2] ) ), [q,p1]:[map (first,S),map (last,S)], prt2 ("q:"), prt2(q),prt2 ("p1:"), prt2(p1), p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ) ), prt2 ("p2:"), prt2(p2), makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[v=p2[k]]],k,1,length(p2)) )$ load ("grobner") $ getrt(sub,f,v):=block([eqv,subv:listofvars (sub),p1,sf,rt00:[],rt01:[],rt0,rt:[],sl,k],prt2(sub,f,v), p1:makelist ( if atom(sub [1] [k]) then sub[2][k][1] elseif is (op (sub [1] [k] )="=") then args(rhs(sub [1] [k] ))[1] else sub[2][k][1], k,1,length (sub [2] ) ), p1:append (p1, [f] ),prt2 ("p1:",p1), sf:poly_reduced_grobner (p1,append (subv,[v] ) ),prt2 ("sf:",sf), map (lambda ([e],if listofvars (e)=[v] then rt00:append (rt00, [e] ) else rt01:append (rt01, [e] ) ),sf ), if rt00=[] then rt:[] else (rt0:apply(realroots,rt00), if rt0=[] then rt:[] else (sub2:apply("+",map(lambda ( [e], abs (lhs (e[2])-rhs (e[2]) ) ),sub [2] )), /* cf. apply("+",[]) = 0 */ for eqv in rt0 do (sl:algsys(subst (eqv,rt01),subv), ans:apply(min, map (lambda ( [e], subst (e,sub2) ),sl )), if is(ans<10^(-5)) then rt:append (rt, [eqv] ), prt2(eqv,sl,"min:",float(ans)) ) ) ), rt:if rt=[] then [] else (rt:sort (map (rhs,rt),"<" ), makelist([root(f,k),[last(rt00),rt[k]]],k,1,length (rt) )), prt2("rt:",rt),rt )$ rootsepsilon:1.0e-10$ getsample11(sub,p,v):=block([p1,S:[],T:[],p2,p0,q,k], p1:apply(append,map(lambda([f],if member (v,listofvars (f) ) then getrt(sub,f,v) else []),p)), prt2("p1 in 11:",p1), if p1=[] then p2:[0] else (map (lambda ( [e], if not(member (e[2][2],T )) then (S:append(S,[e]),T:append(T,[e[2][2]])) ),p1 ), prt2 ("S:",S), if S=[] then (p1:[],p2:[0]) else ( S:sort (S,lambda ( [u,t], u[2][2]<t[2][2] ) ), [q,p0]:[map(first,S),map(last,S)], p0p:map(first,p0),p1:map(last,p0), prt2 ("q:"), prt2(q),prt2 ("p1:"), prt2(p1), p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ) ), prt2 ("p2:"), prt2(p2), makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[ [(if oddp(k) then v-p2[k] else p0p[k/2]), (if oddp(k) then v=p2[k] else v=p1[k/2])] ]],k,1,length(p2)) )$ getsample2(sub,p,v):=block([sub1,rt,q,p1,p2,k], sub1:delete([],map(lambda ( [e], if atom (e) then [] elseif op(e)="=" then e else []), sub[1] ) ), rt:subst(sub1,p), rt:flatten(map (lambda ( [e], if constantp (e) then [] else e ), rt )), if rt=[] then rt:[] else rt:flatten(map (lambda ( [e],solve (e,v) ),rt ) ), if rt=[] then rt:[] else (rt:map (lambda ( [e],[e,ratsimp(subst (sub [2],e )) ] ), map(rhs,rt) ), rt:delete([],map(lambda([e],if imagpart(e [2])=0 then e else [] ), rt )), rt:listify(setify(rt)) ), if rt=[] then (p1:[],p2:[0]) else ( rt:sort (rt,lambda ( [u,t], u[2]<t[2] )),/*print (rt),*/ [q,p1]:[map (first,rt),map (last,rt)], p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ), /*prt2 ("p2:"), prt2(p2),*/ makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[v=ratsimp( p2[k] )]],k,1,length(p2)) )$ is2(f):= if atom (f) then f elseif member (op (f), ["or","and"] ) then map (is,f) else is (f)$ kill("Lo","Lc","Go","Gc","Eq")$ infix("Lo",60,40)$ infix("Lc",60,40)$ infix("Go",60,40)$ infix("Gc",60,40)$ infix("Eq",60,40)$ (x Lo y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x<y)$ (x Lc y):=if member(extd,listofvars([x,y])) then false else float(x<y+10^(-5))$ (x Go y):=if member(extd,listofvars([x,y])) then false else float(x>y+10^(-5))$ (x Gc y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x>y)$ (x Eq y):=if member(extd,listofvars([x,y])) then false else abs(float(x-y))<10^(-5)$ kill("implies")$ infix(implies,60,40)$ "implies"(x,y):=((not(x)) or (y))$ defrule(imp2m, (aa) implies (bb), mor ( mnot((aa)),(bb) )); defrule(noteq2m, (aa) # (bb), mnot (equal((aa),(bb)))); defrule(symneq2m, notequal((aa),(bb)), mnot (equal((aa),(bb)))); cad_monic(F0):=block ( [F,mpp:[],p,sx,sxl,i,k], F:addD(subst(["="=equal,"and"=mand,"or"=mor,"not"=mnot],F0)), F:apply1(F,imp2m,noteq2m,symneq2m), simp:on, for vk in %vv do F:sepaF(F,vk), for vk in %vv do F:sepaS(sepaS(F,vk),vk), for vk in %vv do F:sepaF(F,vk), F:ev(subst([mand="and",mor="or",mnot="not"],F),eval),prt(F), F:scanmap(lambda([f],if atom(f) then f elseif member(op(f),["<","<=",">",">=",equal,notequal]) then (p:dnf(sepalcv(f,%vv[1])),prt(p), if op(p)="or" then (mpp:append (mpp, [args(p)] ), concat(mp,length (mpp) )) else f ) else f),F), print([F,mpp]), if mpp=[] then [F] else ( i:0, sxl:map(lambda([n],(i:i+1,product(concat(sx,i)-k,k,1,n))),map(length,mpp)), sxl:solve(sxl), map(lambda([sxl],subst(ev(subst(mppp=mpp,subst(sxl,makelist(concat(mp,i)=mppp[i][concat(sx,i)],i,1,length(mpp)))),eval),F)),sxl) ) ) $ /* epsilon sgn-det ver. */ cad (F0) :=block ([F,vk,P2:{},pp,vv:reverse(%vv),pn:[[[],[]]],pn0], F:addD(subst(["="=equal,"and"=mand,"or"=mor,"not"=mnot],F0)), F:apply1(F,imp2m,noteq2m,symneq2m), simp:on, for vk in %vv do F:sepaF(F,vk), for vk in %vv do F:sepaS(sepaS(F,vk),vk), for vk in %vv do F:sepaF(F,vk), F1:subst(["<"="Lo",">"="Go","<="="Lc",">="="Gc",equal="Eq"],F), F:ev(subst([mand="and",mor="or",mnot="not"],F),eval), F1:ev(subst([mand="and",mor="or",mnot="not"],F1),eval), prt2(F1), prt(F), /* F:subst(["<"="Lo",">"="Go","<="="Lc",">="="Gc",equal="Eq"],F), */ scanmap (lambda ( [f], if atom (f) then f elseif member (op (f), ["<","<=",">",">=",equal,notequal]) then P2:append (P2, {lhs (f)-rhs (f) } ) else f ),F ), simp:on,P2:listify (P2), /*P2:map(lambda ([f], sqfree (f,(%vv)[1])),P2), */ pp:map(factor,delv(P2,%vv)), prt(pp),%pp:pp, for n:1 thru length (vv) do pn:apply(append,map(lambda([e],map(lambda([t],map(append,e,t)),getsample(e,pp[n],vv[n]))),pn)),prt(length (pn)),pn00:pn, pn: if is(dtype=11) then delete(false,map(lambda([l],if is2( (subst( map(last,l[2]), F1)) ) then (if is(dispsample=true) or is(bsolve=true) then l else l[1])),pn)) else delete(false,map(lambda([l],if is2( (subst(l[2], F1)) ) then (if is(dispsample=true) or is(bsolve=true) then l else l[1])),pn)), /*print(pn),*/ simp:usersimp, if bsolve=true and dispsample=true then cad2s (pn) elseif bsolve=true and dispsample=false then cad2 (pn) else pn )$ cadm(f):=block([usr_display2d], thisout:cad(f), usr_display2d:display2d, set_display('xml), print(apply(matrix,thisout)), display2d:usr_display2d,return(done) )$ addD0(f):=block([DN], if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal]) then (DN:{}, scanmap(lambda([t],if atom(t) then t else (DN:ev(append(DN,{notequal(ratsimp(denom(t)),0)}),simp),t)), f),prt2([f,DN]), DN:ev(substpart("and",DN,0),eval),prt2([f,DN]), if is(DN=false) then return(false) else return( mand((f) , (apply1(DN,symneq2m))) )) else return(f) )$ addD(F):=block([addD0,S], prt2("addD:"), S:scanmap(addD0,F,bottomup),prt2(S), return(S) )$ sepaF0(f,v):=block([p,DN,NU], if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal]) then (p:factor(lhs(f)-rhs(f)), NU:num(p),DN:denom(p),prt2([p,NU,DN]), if member(v,listofvars(DN)) and member(op(f),["<",">"]) then return(apply(op(f),[NU*DN,0])) elseif member(v,listofvars(DN)) and member(op(f),["<=",">="]) then return( mand( (apply(op(f),[NU*DN,0])) , (mnot(equal(DN,0))) ) ) elseif member(v,listofvars(DN)) and op(f)=equal then return( mand( (equal(NU,0)) , (mnot(equal(DN,0))) ) ) else return(f)) else return(f))$ sepaF(F,v):=block([S], prt2("sepaF for",v,":"), S:scanmap(lambda([s],sepaF0(s,v)),F), prt2(S),return(S) )$ sepaSlocal(ff,v):=block([sepaSlocal0,getS,getSqrt,p,Sqrt,A,B,C], sepaSlocal0(f):=block([], if atom(f) or is(getS=true) then return(f) elseif op(f)=sqrt and member(v,listofvars(f)) then (getS:true,getSqrt:f,return(Sqrt)) else return(f)), if atom(ff) then return(ff) elseif member(op(ff),["<",">","<=",">=",equal]) then ( getS:false,p:scanmap(sepaSlocal0,lhs(ff)-rhs(ff)), (if is(getS=true) then (p:subst(getSqrt=Sqrt,p), [A,B]:[coeff(expand(p),Sqrt,1),getSqrt^2],C:expand(A*Sqrt-p), return(mand( [A,B,C,op(ff)] , (getSqrt^2 >= 0) ))) else return(ff)) ) else return(ff))$ defrule(sepaS1,[aa,bb,cc,equal], mand(mor(mand(aa >= 0,cc >= 0),mand(aa < 0,cc <= 0)),/*bb >= 0,*/ equal(cc^2-aa^2*bb,0)) )$ defrule(sepaS2,[aa,bb,cc,"<="], mand(mor(mand(aa < 0,aa^2*bb-cc^2 >= 0), mand(cc >= 0,aa^2*bb-cc^2 <= 0))/*,bb >= 0*/) )$ defrule(sepaS3,[aa,bb,cc,">="], mand(mor(mand(aa > 0,cc^2-aa^2*bb <= 0), mand(cc <= 0,cc^2-aa^2*bb >= 0))/*,bb >= 0*/) )$ defrule(sepaS4,[aa,bb,cc,"<"], mand(mor(mand(aa < 0,cc > 0),mand(aa < 0,aa^2*bb-cc^2 > 0), mand(cc > 0,aa^2*bb-cc^2 < 0))/*,bb >= 0*/) )$ defrule(sepaS5,[aa,bb,cc,">"], mand(mor(mand(aa > 0,cc < 0),mand(aa > 0,cc^2-aa^2*bb < 0), mand(cc < 0,cc^2-aa^2*bb > 0))/*,bb >= 0*/) )$ sepaS(F,v):=block([S], prt2("sepaS:"), S:scanmap(lambda([s],sepaSlocal(s,v)),expand(F)), prt2(S), S:expand(apply1(S,sepaS1,sepaS2,sepaS3,sepaS4,sepaS5)), prt2(S),return(S) )$ sepalcv(f,v):=block([pp,dg,lc], pp:expand(lhs(f)-rhs(f)),dg:hipow(pp,v),lc:ratcoef(pp,v,dg), if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal,notequal]) then (if dg=0 then return(f) else return( (notequal(lc,0) and f) or (equal(lc,0) and sepalcv(apply(op(f),[expand(pp-lc*v^dg),0]),v)) )) else return(f))$ /* cnf dnf */ kill("Or"); kill("And"); infix(Or,60,40)$ infix(And,60,40)$ matchdeclare([aa,bb,cc],true)$ defrule(ru05r, (aa) Or ((bb) And (cc)), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru08r, ((bb) And (cc)) Or (aa), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru05a, (aa) And ((bb) Or (cc)), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru08a, ((bb) Or (cc)) And (aa), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru00r, (aa) Or (bb), (aa) or (bb))$ defrule(ru00a, (aa) And (bb), (aa) and (bb))$ orand2orand(f):=block([g],g:f, if atom(g) then return(g), if op(g)="or" then g:xreduce("Or",args(g)), if op(g)="and" then g:xreduce("And",args(g)),g)$ cnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05r,ru08r),ru00r,ru00a)$ dnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05a,ru08a),ru00r,ru00a)$ cad2(F):=map(lambda([L], (U:[], map(lambda([f], (if atom (f) then f elseif op(f)="=" then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)=(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1], U:append (U, [sol] ),sol) elseif op(f)="<" and listp(rhs (f)) then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)<(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1]) elseif op(f)="<" and listp(lhs (f)) then (sol:solve(subst(U,[first(lhs(f))]),rhs(f)),prt2(sol), sol:(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs(f))][1]<rhs(f)) elseif op(f)="and" then (sol1:solve(subst(U,[first(lhs((args(f))[1]))]),rhs((args(f))[1])), sol1:map(rhs,sol1), sol2:solve(subst(U,[first(rhs((args(f))[2]))]),lhs((args(f))[2])), sol2:map(rhs,sol2), prt2([sol1,sol2]), (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol1 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs((args(f))[1]))][1] < rhs((args(f))[1]) and lhs((args(f))[2]) < (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol2 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs((args(f))[2]))][1] ) else f) ),L[1])) ),F)$ cad2s(F):=map(lambda([L], [(U:[], map(lambda([f], (if atom (f) then f elseif op(f)="=" then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)=(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1], U:append (U, [sol] ),sol) elseif op(f)="<" and listp(rhs (f)) then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)<(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1]) elseif op(f)="<" and listp(lhs (f)) then (sol:solve(subst(U,[first(lhs(f))]),rhs(f)),prt2(sol), sol:(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs(f))][1]<rhs(f)) elseif op(f)="and" then (sol1:solve(subst(U,[first(lhs((args(f))[1]))]),rhs((args(f))[1])), sol1:map(rhs,sol1), sol2:solve(subst(U,[first(rhs((args(f))[2]))]),lhs((args(f))[2])), sol2:map(rhs,sol2), prt2([sol1,sol2]), (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol1 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs((args(f))[1]))][1] < rhs((args(f))[1]) and lhs((args(f))[2]) < (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol2 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs((args(f))[2]))][1] ) else f) ),L[1])),L [2] ] ),F)$ cnn(LL0):=block([n,k,LL:subst("="=equal,LL0),LA,LB,LB0,s,t,L0,M], print("connecting cells."), if is(dispsample=true) then (n:length(LL[1][1]), for k:n step -1 thru 1 do ( LL:append(LL,[[append(makelist(0,n),[1]),makelist([0,0],n)]]), L0:LL[1][1],t:delete(L0[k],L0),s:false, LL:delete(0,map(lambda([L], (LA:L[1], if is(t=delete(LA[k],LA)) then (s:(s or LA[k]),LB:L[2],0) else (M:subst(L0[k]=s,L0),L0:LA,t:delete(LA[k],LA),s:LA[k],LB0:LB,LB:L[2],append([M],[LB0])) ) ),LL)) ), thisout0:LL, thisout:map(lambda([L],[map(lambda([e],if atom(e) then e elseif op(e)="or" then cnn0(e) else e),L[1]),L[2]] ),LL) ) else (n:length(LL[1]), for k:n step -1 thru 1 do ( LL:append(LL,[append(makelist(0,n),[1])]), L0:LL[1],t:delete(L0[k],L0),s:false, LL:delete(0,map(lambda([L], if is(t=delete(L[k],L)) then (s:(s or L[k]),0) else (M:subst(L0[k]=s,L0),L0:L,t:delete(L[k],L),s:L[k],M) ),LL)) ), thisout0:LL, thisout:map(lambda([L],map(lambda([e],if atom(e) then e elseif op(e)="or" then cnn0(e) else e),L) ),LL) ), thisout:subst([equal="="],thisout) )$ distb(f):=block([], if atom(f) then f elseif member(op(f),["or","and"]) then map(distb,f) elseif member(op(f),[equal,"<",">","<=",">="]) then apply(op(f),[{lhs(f)},{rhs(f)}]) else f )$ cnn0(f):=block( [f1,S:[],v,Sb,T,f01], if atom(f) then f else ( map(lambda([e], (if not(member(lhs(e),S)) then S:append(S,[lhs(e)]), if not(member(rhs(e),S)) then S:append(S,[rhs(e)]))), args(subst("and"="or",f))), v:if is("and"=op((args(f))[1])) then S[2] else S[1], Sb:map("{",delete(v,S)),T:makelist({10^5*k},k,length(Sb)), prt2(S,Sb,T,map("=",Sb,T),f,subst(map("=",Sb,T),distb(f))), f01:ev(subst("{"="(",subst(map("=",Sb,T),distb(f))),eval),prt2(f01), f01:ora(subst(["="=equal],cineqs2(f01))), if atom (f01) then f01 else subst(map("=",append([{v}],T),append([v],delete(v,S))),distb(f01)) ) )$ cnnm(LL):=block([k,L0:LL,L1:cnn(LL),L2,usr_display2d], for k:0 unless is(L0=L1) do (L2:cnn(L1),L0:L1,L1:L2), usr_display2d:display2d, set_display('xml), print(apply(matrix,L1)), display2d:usr_display2d,return(done) )$ infix ("==",50,50)$ nary ("&&",50,50)$ nary ("||",50,50)$ sqrtdispflag:off$ mma(F):=subst (["="="==","and"="&&","or"="||"],ora(subst (["="="==",%i=I],F)))$ sqrtdispflag:on$ ora(f):=block([], if not(listp(f)) then return(f), if atom(f) then return(f) else apply("or",map(lambda([e],apply("and",e)),subst("="=equal,f))) )$ cineqs2(F):=block([lv,v,imagunit,L,M,N,k,K,J], lv:listofvars(F),if length(lv)#1 then return([[is(F)]]) else v:first(lv), L:sort(listify(setify(flatten(scanmap( lambda([e], if atom(e) then e elseif member(op(e),["and","or"]) then substpart("[",e,0) elseif member(op(e),["<",">","<=",">=","=","#",equal]) then map(lambda([t],if member(imagunit,listofvars(t)) then [] else rhs(t)),subst(%i=imagunit,solve(substpart("=",e,0),v))) else e), F)))),"<"), if L=[] then return([[is2(subst(v=0,F))]]), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), M:map(lambda([e],if member(op(F),["and","or"]) then map(is2,rat(subst(v=e,F))) else is2(rat(subst(v=e,F)))),L), if length(setify(M))=1 then return([[is2(subst(v=0,F))]]), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), K:[],for k:3 thru length(M)-2 do ( if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K), N:[],J:[],for k in K do if rhs(k)=v then J:[k] else (N:append(N,[append(J,[k])]),J:[]), delete([],append(N,[J])) )$