ATPとCASのこと

2017-03-24

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  ) ))  );

f:id:ehito:20170325144940p:image:w640


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

2017-02-11

Maxima さん,勘弁してください.

1と2は特別なのでしょうか?

k@k ~ $ rmaxima --init=
Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 21b (21B Unicode)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) newcontext();facts();
(%o1)                             context2947
(%o2)                                 []
(%i3) apply(assume,[x > 0,equal(1,sqrt(x))]);
(%o3)                     [x > 0, equal(1, sqrt(x))]
(%i4) apply(assume, [equal(-1,-sqrt(x))]);
*A2 gc_alloc_new_region failed, nbytes=8.
 CMUCL has run out of dynamic heap space (512 MB).
  You can control heap size with the -dynamic-space-size commandline option.

Imminent dynamic space overflow has occurred:
Only a small amount of dynamic space is available now.
Please note that you will be returned to the Top-Level without
warning if you run out of space while debugging.

Maxima encountered a Lisp error:

 Heap (dynamic space) overflow

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i5) 
Maxima encountered a Lisp error:

 Interrupted at #xF7FDB430.

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i5) facts();
(%o5)                     [x > 0, equal(1, sqrt(x))]
(%i6) newcontext();facts();
(%o6)                            context178628
(%o7)                                 []
(%i8) apply(assume,[x > 0,equal(2,sqrt(x))]);
(%o8)                     [x > 0, equal(2, sqrt(x))]
(%i9) apply(assume, [equal(-2,-sqrt(x))]);
*A2 gc_alloc_new_region failed, nbytes=8.
 CMUCL has run out of dynamic heap space (512 MB).

 Returning to top-level.
*A2 gc_alloc_new_region failed, nbytes=8.
 CMUCL has run out of dynamic heap space (512 MB).

 Returning to top-level.
*A2 gc_alloc_new_region failed, nbytes=8.
 CMUCL has run out of dynamic heap space (512 MB).

 Returning to top-level.
rlwrap: warning: maxima killed by SIGTRAP.
rlwrap has not crashed, but for transparency,
it will now kill itself (without dumping core)with the same signal

Trace/breakpoint trap
k@k ~ $ rmaxima --init=
Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 21b (21B Unicode)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) newcontext();facts();
(%o1)                             context2947
(%o2)                                 []
(%i3) apply(assume,[x > 0,equal(2,sqrt(x))]);
(%o3)                     [x > 0, equal(2, sqrt(x))]
(%i4) apply(assume, [equal(-2,-sqrt(x))]);

Maxima encountered a Lisp error:

 Interrupted at #x1006880B.

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i5) facts();
(%o5)                     [x > 0, equal(2, sqrt(x))]
(%i6) newcontext();facts();
(%o6)                            context101642
(%o7)                                 []
(%i8) apply(assume,[x > 0,equal(3,sqrt(x))]);
(%o8)                     [x > 0, equal(3, sqrt(x))]
(%i9) apply(assume, [equal(-3,-sqrt(x))]);
(%o9)                             [redundant]
(%i10) facts();
(%o10)                    [x > 0, equal(3, sqrt(x))]
(%i11) newcontext();facts();
(%o11)                           context101663
(%o12)                                []
(%i13) apply(assume,[x > 0,equal(4,sqrt(x))]);
(%o13)                    [x > 0, equal(4, sqrt(x))]
(%i14) apply(assume, [equal(-4,-sqrt(x))]);
(%o14)                            [redundant]

2017-02-07

cineqs4 更新

変更点

・名前が cineqs4g になりました.

・処理の効率を改善.数パーセント高速化しましたが,下記の処理のため結果として遅くなりました(苦).

(%i15) for c:1 thru 100 do (ep(0,0),cineqs4g('(   x^3+abs(x+1)^(-1/4)-2>0   )))$
Evaluation took 5.2900 seconds (5.2900 elapsed) using 216.100 MB.
(%i16) for c:1 thru 100 do (ep(0,0),cineqs4('(   x^3+abs(x+1)^(-1/4)-2>0   )))$
Evaluation took 1.7900 seconds (1.7900 elapsed) using 145.687 MB.

・解集合の境界点の近似値の表示を5桁までにしました.

・より正確な値を知る(というか厳密解を特定する)ための関数 NP() を新設.境界点の近似値,対応する厳密解を根とする整数係数多項式のリストのリストを出力します.

・その多項式を得るために load("grobner") しているので,初回コール時は少しモタツキます.

(%i17) ep(0,0)$cineqs4g('(   sqrt(x)+sqrt(3-x)+sqrt(3+x)<22/5  ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes.
Evaluation took 0.0500 seconds (0.0600 elapsed) using 1.867 MB.
(%o18) [[0 <= x,x < 0.96473],[2.9379 < x,x <= 3]]
(%i19) NP();
Evaluation took 0.1800 seconds (0.1800 elapsed) using 4.819 MB.
(%o19) [[0,x],
        [0.96473,
         9765625*x^4+285875000*x^3+1561490000*x^2-10002150400*x+7930971136],
        [2.9379,
         9765625*x^4+285875000*x^3+1561490000*x^2-10002150400*x+7930971136],
        [3,3-x]]
(%i20) bfloat(realroots( (NP())[2][2] ,1.0e-100)),fpprec:100;
Evaluation took 0.1900 seconds (0.1900 elapsed) using 6.588 MB.
(%o20) [x = 9.647334915176576305835677154346565758106422109360258068951819981502359285249692172963134559965398747b-1,
        x = 2.937906436935922003266179147367196899336830443976134830225254296066863763304387097967652650705038142b0]
/* 17.02.07.Tue. 22:19:01 */
showtime:on$
display2d:false$
prtoff():=kill(prt)$
prton():=block([],prt([x]):=apply(print,x),return(done))$
ep(e,p):=(
if e=1 then algexact:true else algexact:false,
if p=1 then prton() else (ratprint:off,prtoff()),
usersimp:simp,simp:off
)$
matchdeclare([aa,bb,cc],true)$
kill("implies")$
infix(implies,60,40)$
"implies"(x,y):=((not(x)) or (y))$
Powalgsyssingle2g(f,v):=block([F,eqs,inqs,db,pr0,prk,pr,el],
local(Powalgsyssinglelocal2g),
Powalgsyssinglelocal2g(t):=block([],
if atom(t) then return(t) 
elseif op(t)=Pow3 then
 (X:first(t),P:second(t),NU:num(P),DN:denom(P),prt([X,P,NU,DN]),
  if DN=1 and P>=0 then return(X^P)
  elseif evenp(DN) then (pslc:pslc+1,prk:concat(pr,pslc),
                     eqs:append(eqs,[prk^DN=X^NU]),
                     inqs:append(inqs,[concat(pr,pslc)>=0]),
                     return(prk))
  elseif oddp(DN) then
                    (pslc:pslc+1,prk:concat(pr,pslc),
                     eqs:append(eqs,[prk^DN=X^NU]),
                     return(prk))
  else return(X^P)
 )
else return(t)),
prt("Powalgsyssingle2g:"),pslc:0,eqs:[],inqs:[],
F:scanmap(Powalgsyssinglelocal2g,f,bottomup),prt(F),
if pslc>0 then
       (F:map(num,factor(map(lambda([p],lhs(p)-rhs(p)),
                             append([F],eqs)
             ))),prt([F,inqs]),
        el:[F, append(makelist(concat(pr,k),k,1,pslc),[v])],
        %PL:append(%PL,[el]),
        F:apply(algsys,el),prt(F),
        F:map(lambda([l],subst(l,[append([v],inqs)])),F),
        F:if is(%rnum_list=[]) then F
          else map( lambda([l],if listofvars(l)=[] then l else [[]]), F),
        F:ev(map(lambda([l],[substpart("and",first(l),0)]),F),eval)),
prt(F),return(F)
)$
kill("Lo","Lc2","Go","Gc","Eq")$
infix(Lo,60,40)$
infix(Lc2,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 Lc2 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)$
/*
float_approx_equal_tolerance
float_approx_equal(,)
bfloat_approx_equal(,)
*/
usersimp:simp$simp:off$
defrule(dev2pow, (aa) / (bb), (aa) * (bb) ^ (-1))$
defrule(sqrt2pow, sqrt((aa)), (aa) ^ (1/2))$
simp:usersimp$
defrule(pow2Pow, (aa) ^ (bb),
 if rat(aa)=0 then
  (if rat(bb)>0 then 0 elseif rat(bb)=0 then 1 else Pow3(false0,bb,Pow1))
 elseif listofvars(aa)=[] or (denom(rat(bb))=1 and rat(bb)>=0)
  then (aa) ^ (bb)
 else Pow3(aa,bb,Pow1))$
defrule(abs2Pow, abs((aa)), Pow3((aa)^2,1/2,Pow1))$
defrule(cineqs4n2N1, (aa) # (bb), Not( equal((aa),(bb)) ))$
defrule(cineqs4n2N2, notequal((aa), (bb)), Not( equal((aa),(bb)) ))$
defrule(defbound, Pow3((aa),(bb),(cc)),
 (prt("defbound:"),
  if is(bb<0) then Lnz:append(Lnz,{Not( equal((aa),0) )})
  elseif evenp(denom(rat(bb))) then Lnn:append(Lnn,{equal((aa),0)}),
  Pow3((aa),(bb),(cc))))$
realonly:true$
algexact:true$
/*
algepsilon:10^7$
rootsepsilon:1.0e-17$
*/
cineqs4g(f):=block([L0,F0,F,lv,v,L,sl,func,M,N,k,K,J,el],local(Pow3),
%PL:[],
lv:listofvars(f),if length(lv)#1 then 
 (simp:usersimp,return([[(true) and (f)]])) else v:first(lv),
F:apply1(f,dev2pow,sqrt2pow,pow2Pow),
simp:on,
F:subst(["<"="Lo",">"="Go","<="="Lc2",">="="Gc","="="Eq","not"=Not],F),
F:apply1(F,abs2Pow,cineqs4n2N1,cineqs4n2N2),
prt(F),
Lnn:{},
F:if member(Pow0,listofvars(F)) or member(Pow1,listofvars(F)) then
scanmap(lambda([e],
  if atom(e) then e
  elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal]) then
   (if member(Pow0,listofvars(e)) then false
    elseif member(Pow1,listofvars(e))
     then (Lnz:{},(applyb1(e,defbound)) and (substpart("and",Lnz,0)))
    else e)
  else e)
,F,bottomup)
else F,
if is(listofvars(F)=[]) then
 (F:ev(subst([Not="not"],F),eval),simp:usersimp,
  return([[(true) and (F)]])),
prt([F,Lnn]),
L:sort(delete(true,delete(false,listify(setify(flatten(
map(
 lambda([e],
if atom(e) then e
elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal,notequal])
 then (if member(Pow1,listofvars(e)) then (Powalgsyssingle2g(e,v))
       else
       (el:apply(num,[factor(lhs(e)-rhs(e))]),%PL:append(%PL,[[[el],[v]]]),
        sl:map(rhs,flatten(algsys([el], [v]))),
        delete(first(append(%rnum_list,[gensym()])),sl)))
else e),
flatten([subst(["%and"="[","%or"="[","and"="[","or"="[",mand="[",mor="[",Not="[","implies"="["],[F,listify(Lnn)])])
)
))))),"<"),
F:subst([Not="not"],F),
prt([L,F]),
Pow3(x,y,z):=if (is(x<0) and evenp(denom(y))) or (is(x Eq 0) and is(y<0)) then extd else x^y,
func(t):=subst(v=t,F),
/*prt(func(gensym("x"))),*/
if L=[] then (simp:usersimp,return([[apply(is,[func(0)])]])),
L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"),
M:map(is,(map(func,rat(L)))),prt([L,M]),
if length(setify(M))=1 then (simp:usersimp,return([[apply(is,[func(0)])]])),
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:[]),
simp:usersimp,
%BL:delete([],append(N,[J])),
eval_string(ev(string(%BL),fpprintprec:5))
)$
NP():=block([pl,v,bl,bp],
if not ?boundp ('poly_monomial_order) then load("grobner"),
v:first(listofvars(%BL)),
pl:factor(listify(setify(flatten(map( lambda([l], map(lambda([p],if polynomialp(p,[v]) then p else []),l) ),map(lambda([l],apply(poly_grobner,l)),%PL)))))),
pl:map(lambda([p],[p,map(rhs,flatten(algsys([p],[v])))]),pl),
prt(pl),
bl:delete(v,flatten(subst(["<"="[","<="="[",">"="[",">="="[","="="["],%BL))),
[pl,bl]:eval_string(ev(string([pl,bl]),fpprintprec:5)),
makelist(
 flatten(map( lambda([p],if member(bp,second(p)) then [bp,first(p)] else []), pl)),
bp,bl)
)$

2017-02-05

cineqs4 原理

cineqs4 は「真理値に関する中間値定理」により,原始式同士の結合の様子に立ち入ることなく,系を解いています.すなわち,...

変数代数系 F(x) と実数 a,b (a<b) についての

 {F(a),F(b)}={true,false}

F(x) のある原始式 rel(f(x),0) (rel∈{>,<,>=,<=,=,#}) について,{rel(f(a),0),rel(f(b),0)}={true,false}

例えば,f(a)>0はtrue,f(b)>0はfales(つまり,f(b)<=0はtrue)なので

F(x) のある原始式 rel(f(x),0) (rel∈{>,<,>=,<=,=,#}) について,f(c)=0,a<=c<=b となる実数 c がある

という性質(つまり,連続関数 f の零点に触れなければ F の真理値は変わらないこと)により,系に現れるすべての f のすべての相異なる実数の零点 n (>=0) 個とその間,および,両脇との合計 2n+1 個の代表点についての F の真理値は,各点,各区間における真理値に一致するので,true となった代表点に対する点,区間を表す式の選言が解集合の式となる訳です.

具体例です.

(%i1) ep(1,0)$cineqs4('(   ((x-1)*(x^2-3)<=0 and x>0) or x^2-x-1=0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 856 bytes.
Evaluation took 0.0800 seconds (0.0800 elapsed) using 2.356 MB.
(%o2) [[x = -(sqrt(5)-1)/2],[1 <= x,x <= sqrt(3)]]

処理の様子は,コントローラー ep (algExact and Print) の第二引数を 1 にすれば確認できます.

(%i3) ep(1,1)$cineqs4('(   ((x-1)*(x^2-3)<=0 and x>0) or x^2-x-1=0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 992 bytes.
((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) 
((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) 
[((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0),{}] 
[[-sqrt(3),-(sqrt(5)-1)/2,0,1,(sqrt(5)+1)/2,sqrt(3)],
 ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0)]
  
[[((-2*sqrt(3))-1)/2,-sqrt(3),((-(sqrt(5)-1)/2)-sqrt(3))/2,-(sqrt(5)-1)/2,
  -(sqrt(5)-1)/4,0,1/2,1,((sqrt(5)+1)/2+1)/2,(sqrt(5)+1)/2,
  ((sqrt(5)+1)/2+sqrt(3))/2,sqrt(3),(2*sqrt(3)+1)/2],
 [false,false,false,true,false,false,false,true,true,true,true,true,false]]
  
Evaluation took 0.0300 seconds (0.0400 elapsed) using 2.036 MB.
(%o4) [[x = -(sqrt(5)-1)/2],[1 <= x,x <= sqrt(3)]]

この例では多項式なので最初の3行はあまり変化しません(Lc2 は less close,Go は greater open,Eq は equal).

次の [-sqrt(3),-(sqrt(5)-1)/2,0,1,(sqrt(5)+1)/2,sqrt(3)] が (x-1)*(x^2-3),x,x^2-x-1 の零点 6 個であり,最後に表示されたリストが,上で述べた代表点 13 個と対応する ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) の真理値です.

その結果において,左から 4 番目が true なので [x = -(sqrt(5)-1)/2],8から12番目までがtrue なので [1 <= x,x <= sqrt(3)],他は false,ということで,(%o4) が得られます.

2017-02-03

cineqs4 使用例

作動確認バージョン.

Maxima 5.39.0
using Lisp CMU Common Lisp 21b (21B Unicode)

高次不等式.

(%i1) ep(1,0)$cineqs4('(   (x-1)*(x^2-2)*(x^3-3)*(x^2+4)<=0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 856 bytes.
Evaluation took 0.0600 seconds (0.0600 elapsed) using 1.337 MB.
(%o2) [[-sqrt(2) <= x,x <= 1],[sqrt(2) <= x,x <= 3^(1/3)]]

代数的な表示が出来ない場合は,近似値(内部では小数第6位以下を切り捨てています)で答えます.

(%i3) ep(1,0)$cineqs4('(   (x-1)*(x^2-2)*(x^3-3)*(x^2+4)<=1/10   ));
Evaluation took 0.0000 seconds (0.0100 elapsed) using 720 bytes.
Evaluation took 0.0100 seconds (0.0100 elapsed) using 1.091 MB.
(%o4) [[-1.414631932468004 <= x,x <= 1.010331702011963],
       [1.374072765807135 <= x,x <= 1.473559962228517]]

分数関数で表された不等式.

(%i5) ep(1,0)$cineqs4('(   (x-1)/x<2 ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0000 seconds (0.0100 elapsed) using 520.242 KB.
(%o6) [[x < -1],[0 < x]]

連言.

(%i7) ep(1,0)$cineqs4('(   (x-1)*(x^2-2)*(x^3-3)<=0 and (x-1)/x<2 ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0400 seconds (0.0500 elapsed) using 2.571 MB.
(%o8) [[-sqrt(2) <= x,x < -1],[0 < x,x <= 1],[sqrt(2) <= x,x <= 3^(1/3)]]

選言.

(%i9) ep(1,0)$cineqs4('(   (x-1)*(x^2-2)*(x^3-3)<=0 or (x-1)/x<2 ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0300 seconds (0.0300 elapsed) using 2.804 MB.
(%o10) [[true]]

こんな場合は...

(%i11) ep(1,0)$cineqs4('(   x^4-x-1<=0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.2500 seconds (0.2500 elapsed) using 24.248 MB.
(%o12) [[-(sqrt(sqrt(6)*sqrt((sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(1/3)
                                                *(2^(8/3)*18^(2/3)*sqrt(849)
                                                 -9*2^(8/3)*18^(2/3))
                              +(sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(2/3)
                                                  *(2^(1/3)*18^(2/3)*sqrt(849)
                                                   -9*2^(1/3)*18^(2/3))
                              +(sqrt(849)+9)^(1/3)
                               *(32*18^(2/3)*sqrt(849)-16*18^(5/3)))
                 +3*2^(11/3)*(3^(3/2)*sqrt(283)-27)^(1/3)
                 -8*18^(2/3)*(sqrt(3)*sqrt(283)+9)^(1/3))
          -sqrt(6)*sqrt(2^(1/3)*(3*sqrt(849)-27)^(1/3)*(3*sqrt(849)+155)^(1/3)
                         -2^(8/3)*(3*sqrt(849)-27)^(1/3)))
          /24
           <= x,
         x <= (sqrt(sqrt(6)*sqrt((sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(1/3)
                                                    *(2^(8/3)*18^(2/3)
                                                             *sqrt(849)
                                                     -9*2^(8/3)*18^(2/3))
                                  +(sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(2/3)
                                                      *(2^(1/3)*18^(2/3)
                                                               *sqrt(849)
                                                       -9*2^(1/3)*18^(2/3))
                                  +(sqrt(849)+9)^(1/3)
                                   *(32*18^(2/3)*sqrt(849)-16*18^(5/3)))
                     +3*2^(11/3)*(3^(3/2)*sqrt(283)-27)^(1/3)
                     -8*18^(2/3)*(sqrt(3)*sqrt(283)+9)^(1/3))
           +sqrt(6)*sqrt(2^(1/3)*(3*sqrt(849)-27)^(1/3)
                                *(3*sqrt(849)+155)^(1/3)
                          -2^(8/3)*(3*sqrt(849)-27)^(1/3)))
           /24]]

コントローラー ep の第一引数を 0 にすると,近似値で答えます.

(%i13) ep(0,0)$cineqs4('(   x^4-x-1<=0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes.
Evaluation took 0.0000 seconds (0.0100 elapsed) using 175.711 KB.
(%o14) [[-0.7244919786096257 <= x,x <= 1.220744081172491]]

無理関数で表された不等式.

(%i15) ep(1,0)$cineqs4('(   sqrt(x)<x-1  ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0100 seconds (0.0100 elapsed) using 674.156 KB.
(%o16) [[(sqrt(5)+3)/2 < x]]

含意.

(%i17) ep(1,0)$cineqs4('(   1<sqrt(x) implies 1<x  ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0100 seconds (0.0100 elapsed) using 428.797 KB.
(%o18) [[true]]

根号が複数あっても.

(%i19) ep(1,0)$cineqs4('(   sqrt(x)+sqrt(3-x)<2  ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0300 seconds (0.0200 elapsed) using 1.524 MB.
(%o20) [[0 <= x,x < -(2^(3/2)-3)/2],[(2^(3/2)+3)/2 < x,x <= 3]]
(%i21) ep(0,0)$cineqs4('(   sqrt(x)+sqrt(3-x)+sqrt(3+x)<22/5  ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes.
Evaluation took 0.0200 seconds (0.0200 elapsed) using 2.500 MB.
(%o22) [[0 <= x,x < 0.9647335423197492],[2.937906564163217 < x,x <= 3]]

ネストされても.

(%i23) ep(1,0)$cineqs4('(   1/sqrt(sqrt(1+x)-x)<2   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0700 seconds (0.0800 elapsed) using 4.835 MB.
(%o24) [[-1 <= x,x < 5/4]]

有理数冪でも.

(%i25) ep(1,0)$cineqs4('(   x^4-x^(1/3)-1>0   ));
Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes.
Evaluation took 0.0300 seconds (0.0400 elapsed) using 1.039 MB.
(%o26) [[x < -0.6196702671972711],[1.198342827550492 < x]]

一変数有理冪半代数系ソルバー cineqs4

cineqs2 を分数関数や有理冪乗(fractional power)関数で表された方程式,不等式まで扱えるように拡張しました.

コードの大半は,Maxima の非論理的挙動を封じるためのものです(笑).

作動確認は Maxima 5.39.0 on Ubuntu 及び windows で行っています.

/* 17.02.05.Sun. 12:17:08 */
showtime:on$
display2d:false$
prtoff():=kill(prt)$
prton():=block([],prt([x]):=apply(print,x),return(done))$
ep(e,p):=(
if e=1 then algexact:true else algexact:false,
if p=1 then prton() else (ratprint:off,prtoff()),
usersimp:simp,simp:off
)$
matchdeclare([aa,bb,cc],true)$
kill("implies")$
infix(implies,60,40)$
"implies"(x,y):=((not(x)) or (y))$
powalgsyssingle2(f,v):=block([powalgsyssinglelocal2,F,eqs,inqs,db,paral,pr0,prk,pr],
powalgsyssinglelocal2(t):=block([],
if atom(t) then return(t) 
elseif op(t)=Pow3 then
 (X:first(t),P:second(t),NU:num(P),DN:denom(P),prt([X,P,NU,DN]),
  if is(DN=1) and is(P>=0) then return(X^P)
  elseif evenp(DN) then (pslc:pslc+1,prk:concat(pr,pslc),
                     eqs:append(eqs,[prk^DN=X^NU]),
                     inqs:append(inqs,[concat(pr,pslc)>=0]),
                     db:append(db,map(rhs,flatten(algsys([X],[v])))),
                     paral:append(paral,%rnum_list),
                     return(prk))
  elseif oddp(DN) then
                    (pslc:pslc+1,prk:concat(pr,pslc),
                     eqs:append(eqs,[prk^DN=X^NU]),
                     db:append(db,map(rhs,flatten(algsys([X],[v])))),
                     return(prk))
  else return(X^P)
 )
else return(t)),
prt("powalgsyssingle2:"),pslc:0,eqs:[],inqs:[],db:[],paral:[],
F:scanmap(powalgsyssinglelocal2,apply1(f,cineqs4r1),bottomup),prt(F),
if is(pslc>0) then
       (F:map(num,factor(map(lambda([p],lhs(p)-rhs(p)),
             append([F],eqs)
             ))),prt([F,inqs,db]),
        F:algsys(F, append([v],makelist(concat(pr,k),k,1,pslc))  ),prt(F),
        F:map(lambda([l],subst(l,[append([v],inqs),db])),F),
        F:if is(%rnum_list=[]) then F
          else subst( map( lambda([e],e=false), %rnum_list) ,F),
        F:ev(map(lambda([l],[substpart("and",first(l),0),last(l)]),F),eval)),
prt(F),return(F)
)$
kill("Lo","Lc2","Go","Gc","Eq")$
infix(Lo,60,40)$
infix(Lc2,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 Lc2 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)$
/*
float_approx_equal_tolerance
float_approx_equal(,)
bfloat_approx_equal(,)
*/
defrule(cineqs4r1, (aa) ^ (bb),
 if is(listofvars(aa)=[]) or ( is(denom(bb)=1) and is(bb>=0) )
 then (aa) ^ (bb) else Pow3(aa,bb,Pow0))$
defrule(cineqs4r2, (aa) / (bb),
 if is(listofvars(bb)=[]) then (aa) / (bb) else (aa) * Pow3(bb,-1,Pow0))$
defrule(cineqs4r3, (aa) # (bb), Not( equal((aa),(bb)) ))$
defrule(cineqs4r4, notequal((aa), (bb)), Not( equal((aa),(bb)) ))$
defrule(cineqs4r5, Div((aa), 0), Div((aa),0,Div0))$
defrule(cineqs4r6, Div((aa), (bb)), Div((aa),(bb),Div1))$
defrule(cineqs4r7, abs((aa)), Pow3((aa)^2,1/2,Pow0))$
defrule(defbound, Pow3((aa),(bb),(cc)),
 (L0:append(L0,{equal((aa),0)}),Pow3((aa),(bb),(cc))))$
defrule(defbound2, Div((aa),(bb),(cc)),
 (prt("defbound2:"),L0:append(L0,{Not( equal((bb),0)) }),(aa)/(bb)))$
defrule(defbound3, Div((aa),(bb)), ((aa)/(bb)))$
realonly:true$
algexact:true$
algepsilon:10^17$
rootsepsilon:1.0e-17$
cineqs4(f):=block([L0,F0,F,lv,v,L,sl,func,M,N,k,K,J],local(Pow3),
lv:listofvars(f),
if length(lv)#1 then 
(simp:usersimp,return([[(true) and (f)]])) else v:first(lv),
F:subst(["<"="Lo",">"="Go","<="="Lc2",">="="Gc", "="="Eq","not"=Not,"/"=Div,0^0=1],f),
F:apply1(F,cineqs4r3,cineqs4r4,cineqs4r5,cineqs4r6,cineqs4r7),
prt(F),
F0:if member(Div0,listofvars(F)) then
scanmap(lambda([e],
  if atom(e) then e
  elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal]) and
         member(Div0,listofvars(e))
   then false
  else e)
,F,bottomup)
else F,
if is(listofvars(F0)=[]) then
 (F0:ev(subst([Not="not"],F0),eval),simp:usersimp,
  return([[(true) and (F0)]])),
F1:if member(Div1,listofvars(F)) then
scanmap(lambda([e],
  if atom(e) then e
  elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal])
   and  member(Div1,listofvars(e))
   then (L0:{},(applyb1(e,defbound2)) and (substpart("and",L0,0)))
  else e)
,F0,bottomup)
else F0,
F:F1,
prt(F),
simp:on,L0:{},
F:apply1(F,cineqs4r1,defbound),
prt([F,L0]),
L:sort(delete(false,listify(setify(flatten(scanmap(
 lambda([e],
if atom(e) then e
elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal,notequal])
 then (if member(Pow0,listofvars(e))
       then (powalgsyssingle2(e,v))
       else
       (sl:map(rhs,flatten(algsys(
        [apply(num,[factor(lhs(e)-rhs(e))])], [v]))),
        delete(first(append(%rnum_list,[gensym()])),sl))
      )
elseif member(op(e),["%and","%or","and","or",mand,mor,Not,"implies"])
 then substpart("[",e,0)
else e),
mand(F, substpart(mand,L0,0)) ))))),"<"),
F:subst([Not="not"],F),
F:apply1(F,defbound3),
prt([L,F]),
Pow3(x,y,z):=if (is(x<0) and evenp(denom(y))) or (is(x Eq 0) and is(y<0)) then extd else x^y,
func(t):=subst(v=t,F),
/*prt(func(gensym("x"))),*/
if L=[] then (simp:usersimp,return([[apply(is,[func(0)])]])),
L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"),
M:map(is,ratsimp(map(func,L))),prt([L,M]),
/*
M:map(is,(map(func,L))),prt([L,M]),
fpprec:64,
M:map(lambda([e],is(rationalize(bfloat(subst([v=e],F))))),L),prt([L,M]),
fpprec:16,
M:map(lambda([e],is(subst(v=float(e),F))),L),prt([L,M]), 
*/
if length(setify(M))=1 then (simp:usersimp,return([[apply(is,[func(0)])]])),
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:[]),
simp:usersimp,
delete([],append(N,[J]))
)$

Maxima さん,bug ってますよ.

Maxima 5.39.0 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 21b (21B Unicode)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) assume(a>0,b>0);
(%o1)                           [a > 0, b > 0]
(%i2) assume(a*b<0);
(%o2)                           [inconsistent]
(%i3) forget(a>0,b>0);
(%o3)                           [a > 0, b > 0]
(%i4) assume(a*b<0);
(%o4)                              [a b < 0]
(%i5) assume(a>0,b>0);
(%o5)                           [a > 0, b > 0]
(%i6) facts();
(%o6)                       [0 > a b, a > 0, b > 0]
(%i7) forget(a>0,b>0,a*b<0);
(%o7)                       [a > 0, b > 0, a b < 0]
(%i8) sqrt(x)>0;
(%o8)                             sqrt(x) > 0
(%i9) is(%);
(%o9)                               unknown
(%i10) sqrt(x)>=0;
(%o10)                           sqrt(x) >= 0
(%i11) is(%);
(%o11)                               true