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]))
)$