solvex 概要
solvex は冪乗根関数で表された方程式，不等式の系を指定された変数についての１次式を原始式とする系に変換する関数です．
例えば，３数の相加平均と相乗平均とが一致する条件を得ようとしても，qex では．．．
(%i1) qex(a>=0 %and b>=0 %and c>=0 %and (a*b*c)^(1/3)=(a+b+c)/3);
Evaluation took 1.0500 seconds (1.4300 elapsed) using 117.327 MB.
(%o1) (a >= 0) %and (b >= 0) %and (c >= 0)
%and (c^3+3*b*c^2+3*a*c^2+3*b^2*c-21*a*b*c+3*a^2*c+b^3+3*a*b^2+3*a^2*b+a^3 = 0)
一方，solvex は
(%i2) solvex(a>=0 %and b>=0 %and c>=0 %and (a*b*c)^(1/3)=(a+b+c)/3); Evaluation took 8.3100 seconds (16.3700 elapsed) using 761.807 MB. (%o2) [[0 <= a,b = a,c = a]]
と答えます．
solvex の書式は solvex(formula,var1,...,varN)であり，formula を var1,...,varN のうち listofvars(formula) の要素について順に解こうとし，それら以外の変数は定数と見做します．
ただし，solvex(formula) は listofvars(formula) の全ての要素について順に解こうとし，それら以外の変数は定数と見做します．
出力は，bracket dnf，つまり，２重のリストであり，外側のリストのヘッドは %or，内側のリストのヘッドは %and を表わしています．
%or，%and の表示に変換するには，ora コマンドで
(%i3) ora(%); Evaluation took 0.0000 seconds (0.0100 elapsed) using 817.102 KB. (%o3) (0 <= a) %and (b = a) %and (c = a)
元に戻すには，bra コマンドで
(%i4) bra(%); Evaluation took 0.0200 seconds (0.0200 elapsed) using 1.058 MB. (%o4) [[0 <= a,b = a,c = a]]
といった具合です．
また，この書式による入力も可能で
(%i35) solvex([[x^2=a,x>-1]],x); Evaluation took 8.4900 seconds (18.3800 elapsed) using 469.256 MB. (%o35) [[a < 1,x = -sqrt(a)],[x = sqrt(a)]]
と出来，ご覧のように冪乗根を含んだ出力に対応しています．
qex 概要
以下は Maxima 5.39.0 using Lisp CMU Common Lisp 21b (21B Unicode) 上での実行です．
qex は qe コマンドを有理数冪入力対応に拡張したもの（qe extended version）であり，qex([],formula) は qex(formula) のように省略可能です．また，出力は qe と同じく整数係数の多変数多項式による表示です．
平方根．
(%i1) qex(y=sqrt(x)); Evaluation took 0.4700 seconds (0.8600 elapsed) using 28.694 MB. (%o1) (y >= 0) %and (x-y^2 = 0)
複雑な式．
(%i2) qex(y=2^(2/3)/sqrt(2-2^(1/3))); Evaluation took 0.6300 seconds (1.0100 elapsed) using 35.511 MB. (%o2) (y > 0) %and (3*y^6-6*y^4-12*y^2-8 = 0)
不等式．
(%i3) qex(y<sqrt(x)); Evaluation took 0.5900 seconds (0.9900 elapsed) using 40.485 MB. (%o3) (x >= 0) %and ((y < 0) %or (y^2-x < 0))
否定．
(%i4) qex(y#sqrt(x)); Evaluation took 0.4400 seconds (0.8400 elapsed) using 25.491 MB. (%o4) (y < 0) %or (x-y^2 # 0)
絶対値．
(%i5) qex(y=abs(x)); Evaluation took 0.7100 seconds (1.2000 elapsed) using 52.369 MB. (%o5) (y-x >= 0) %and (y+x >= 0) %and ((y-x = 0) %or (y+x = 0))
不等式．
(%i6) qex(x<=abs(x)); Evaluation took 0.4700 seconds (0.8500 elapsed) using 32.271 MB. (%o6) true
元が２つの集合の最大値．
(%i7) qex(x=max2(a,b)); Evaluation took 1.1800 seconds (1.7000 elapsed) using 88.835 MB. (%o7) ((b-x = 0) %or (x-a = 0)) %and (b-x <= 0) %and (x-a >= 0)
３つの場合．
(%i8) qex(x=min2(min2(a,b),c)); Evaluation took 2.4100 seconds (3.5600 elapsed) using 249.620 MB. (%o8) ((b-x = 0) %or (c-x = 0) %or (x-a = 0)) %and (b-x >= 0) %and (c-x >= 0) %and (x-a <= 0)
値域．
(%i9) qex([[E,x]],y=2/sqrt(x)+abs(x)); Evaluation took 1.0100 seconds (1.7000 elapsed) using 69.622 MB. (%o9) y-3 >= 0
値域．
(%i10) qex([[E,x]],y=abs(x)+abs(x-1)+abs(x-2)+abs(x-3)); Evaluation took 3.1500 seconds (25.8000 elapsed) using 237.202 MB. (%o10) y-4 >= 0
値域．
(%i11) qex([[E,x]],y=max2(x,min2(2-x,3-x))); Evaluation took 1.4100 seconds (2.0700 elapsed) using 89.527 MB. (%o11) y-1 >= 0
不等式．
(%i12) qex(abs(a+b)<=abs(a)+abs(b)); Evaluation took 0.9400 seconds (1.4900 elapsed) using 80.652 MB. (%o12) true
不等式．
(%i13) qex(min2(a,b)<=max2(x,a-x)+min2(x,b-x)); Evaluation took 1.3800 seconds (1.7700 elapsed) using 134.412 MB. (%o13) true
相加平均と相乗平均との大小関係．
(%i14) qex([[A,a],[A,b]],(a>=0 %and b>=0) %implies sqrt(a*b)<=(a+b)/n); Evaluation took 2.5200 seconds (3.1200 elapsed) using 252.836 MB. (%o14) (n-2 <= 0) %and (n > 0)
等式となる条件．
(%i15) qex("%and"(a>=0,b>=0,sqrt(a*b)=(a+b)/2));
Evaluation took 0.6700 seconds (1.1000 elapsed) using 57.813 MB.
(%o15) (a >= 0) %and (b-a = 0)
一様連続性．
(%i16) qex([[A,e]],(0<e) %implies (qex([[E,d]],(0<d) %and (qex([[A,x],[A,y]],(abs(x-y)<d) %implies (abs(x^(2/3)-y^(2/3))<e)))))); Evaluation took 6.2100 seconds (24.2100 elapsed) using 586.517 MB. (%o16) true
拡張された量化子．
(%i17) qex([[X2,x]],abs(abs(abs(x^2-1)-a)-1)=x+1);
Evaluation took 4.8000 seconds (10.2200 elapsed) using 333.388 MB.
(%o17) (((a-3 < 0) %and (4*a-5 > 0)) %or ((a-1 < 0) %and (4*a+3 > 0))
%or (a+1 < 0) %or (4*a-13 > 0))
%and (a-3 # 0) %and (a-1 # 0) %and (a+1 # 0) %and (a+3 > 0)
%and (4*a-13 # 0) %and (4*a-5 # 0) %and (4*a+3 # 0)
簡約．
(%i18) nnss(%);
Evaluation took 4.8300 seconds (8.3800 elapsed) using 186.479 MB.
(%o18) ((a-3 < 0) %and (4*a-5 > 0)) %or ((a-1 < 0) %and (4*a+3 > 0))
%or (4*a-13 > 0)
%or ((a+1 < 0) %and (a+3 > 0))
整理．(13/4 < a) が最初に来るのは %or が勝手に sort するためです．
(%i19) l2r(%);
Evaluation took 0.1400 seconds (0.1400 elapsed) using 4.450 MB.
(%o19) (13/4 < a) %or ((-3 < a) %and (a < -1)) %or ((-3/4 < a) %and (a < 1))
%or ((5/4 < a) %and (a < 3))
nns シリーズ更新
新規の関数
・qex
qeの拡張版，有理数冪の入力に対応．
・solvex
多変数の方程式，不等式を原始式とする quantifier free な論理式の solver，有理数冪の入力，出力に対応．
/* 17.01.09.Mon. 07:50:51 */
load(qepmax)$
showtime:on$
display2d:false$
prt(x):=print(x)$
prt(x):=null$
/*********/
/* tools */
/*********/
/* compare length function */
comparelength(f,g):=length(f)<=length(g)$
comparelength2(f,g):=length(g)<=length(f)$
/* cnf dnf */
infix(Or,60,40)$
infix(And,60,40)$
matchdeclare([aa,bb,cc,dd,ee,ff,gg,aa1,bb1,cc1,dd1,ee1,ff1,gg1],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)$
cnfb(f):=applyb2(applyb2(scanmap(orand2orand,f,bottomup),ru05r,ru08r),ru00r,ru00a)$
dnfb(f):=applyb2(applyb2(scanmap(orand2orand,f,bottomup),ru05a,ru08a),ru00r,ru00a)$
/* necessary condition */
defrule(ru00n, (aa) # (bb), "%or"((aa)<(bb),(bb)<(aa)))$
defrule(ru00n, (aa) # (bb), "%not"((aa)=(bb)))$
deli2(f):=(applyb2(f,ru001,ru002,ru003,ru004))$
ncondeq(f):=block([exs,L,exk],
if is(ncondtype=off) then return(f),
econd:[],L:true,
prt("necessary eq. condition"),
exs:fullmap(lambda([x],[E,x]),
sort(rest(full_listify(powerset(setify(listofvars(f))))),comparelength2)),
for exk in exs do
(if is(ncondtype=full) then q:qe(exk,(f)) else q:deli2(qe(exk,(f))),
prt([L,q]),
if is(not(qe([],(L)%implies(q))=true)) then L:(L)%and(q)),
prt(L),
if is(not(ncondsimp=off)) then L:nnscandc(L),
/* if is(qe([],(L)%implies(f))=true) and (tdeg(L)=1) then (prt("equivalent condition"),econd:[dnf(L)],return(dnf(L))) else */
return((L)%and(f))
)$
ncond(f):=block([exs,L,exk],L:true,
prt("necessary condition"),
exs:fullmap(lambda([x],[E,x]),
sort(rest(full_listify(powerset(setify(listofvars(f))))),comparelength2)),
for exk in exs do
(q:qe(exk,f),prt([L,q]),
if not(qe([],(L)%implies(q))=true) then L:(L)%and(nnscan(q))),
prt(L),
return((L)%and(f))
)$
/***************/
/* nns section */
/***************/
/* non-nonsense */
nns(f):=block([F,L,M],
if atom(f) then return(f)
elseif op(f)="%and" then
(F:sort(args(f),comparelength2),L:setify(F),
for g in F do
(M:disjoin(g,L),prt(M),
if p2emt((substpart("%and",M,0))%implies(g))=true then L:M),
return(substpart("%and",L,0)))
elseif op(f)="%or" then
(F:sort(args(f),comparelength2),L:setify(F),
for g in F do
(M:disjoin(g,L),prt(M),
if p2emt((g)%implies(substpart("%or",M,0)))=true then L:M),
return(substpart("%or",L,0)))
else f)$
/* non-nonsense scan */
nnscan(f):=(prt("nnscanning..."),scanmap(nns,f))$
nnscanb(f):=(prt("bottomup nnscanning..."),scanmap(nns,f,bottomup))$
/* non-nonsense scan with cnf */
nnscanc(f):=(prt("cnf nnscanning..."),scanmap(nns,cnf(f)))$
nnscancb(f):=(prt("cnf bottomup nnscanning..."),scanmap(nns,cnf(f),bottomup))$
/* non-nonsense scan with dnf */
nnscand(f):=(prt("dnf nnscanning..."),scanmap(nns,dnf(f)))$
nnscandb(f):=(prt("dnf bottomup nnscanning..."),scanmap(nns,dnf(f),bottomup))$
nnscand2(f):=(prt("dnf nnscanning..."),nns(map(nns,dnf(f))))$
/* non-nonsense scan with cnf,dnf */
nnscancd(f):=nnscanc(nnscand(f))$
nnscancdb(f):=nnscancb(nnscandb(f))$
nnscandc(f):=nnscand(nnscanc(f))$
nnscandcb(f):=nnscandb(nnscancb(f))$
/* non-nonsense set */
nnss(f):=block([f1,f2,f3,q,nns4],
if mode=d then nns4(g):=nnscand(g) else nns4(g):=nns(g),
prt("pre-simplified formula"),f1: dnf (qe ([],f) ) ,
if atom(f1) then return(f1),
prt(length(f1)),prt(f1),prt("necessary condition"),
if op(f1)="%or" then
(f2:[],
for k:1 thru length(f1) do
(q:ncond(part(f1,k)),prt(k),prt(q),f2:append(f2,[q])),
f2:substpart("%or",f2,0))
else
(f2:dnf(ncond(f1)),prt(f2)),
prt("simplifing..."),
if op(f2)="%or" then
(f3:[],
for k:1 thru length(f2) do
(q:qe([],part(f2,k)),prt(k),prt(q),f3:append(f3,[q])),
prt("nnscanning..."),f3:substpart("%or",f3,0),prt(f3),return(nns4(f3)))
else q:qe([],f2),prt(q),prt("nnscanning..."),return(nns4(q))
)$
/* nnss with double negation */
nnssdn(f):=qe([],%not(nnss(%not(f))))$
/**************/
/* p2t section */
/**************/
defrule(powertk1, (aa)^(bb), makeprtk(aa,num(bb),denom(bb)))$
defrule(powertk2, abs(aa), makeprtk(aa^2,1,2))$
domain:complex$
defrule(powertk2, abs(aa), (aa^2)^(1/2))$
defrule(powertk3, max2(aa,bb), (aa+bb+abs(aa-bb))/2)$
defrule(powertk4, min2(aa,bb), (aa+bb-abs(aa-bb))/2)$
p2e0(f):=block([L,makeprtk,g,q,rt,gg],L:[],
makeprtk(x,a,b):=block([rtk],
if b=1 then return(x^a)
else (rtk:concat(rt,length(L)+1),
if oddp(b) then L:append(L,[(rtk^b=x^a)])
else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]),
return(rtk))),
g:applyb1(f,powertk4,powertk3,powertk2,powertk1),
q:makelist([E,concat(rt,k)],k,length(L)),
if L=[] then return(f)
else (gg:(substpart("%and",L,0))%and(g),
ev(qe(q,gg),noeval))
)$
/* 1-1 type */
p2e1(f):=block([L,M,makeprtk,g,q,rt,gg],M:[],L:[],
makeprtk(x,a,b):=block([rtk],
if b=1 then return(x^a)
elseif is(assoc(x^(a/b),M)#false) then return(concat(rt,assoc(x^(a/b),M)))
else (M:append(M,[[x^(a/b),length(L)+1]]),rtk:concat(rt,length(L)+1),
if oddp(b) then L:append(L,[(rtk^b=x^a)])
else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]),
return(rtk))),
g:applyb1(f,powertk4,powertk3,powertk2,powertk1),
q:makelist([E,concat(rt,k)],k,length(L)),
if L=[] then return(f)
else (gg:(substpart("%and",L,0))%and(g),
ev(qe(q,gg),noeval))
)$
/* base64 ver */
p2e64(f):=block([sL,L,makeprtk,g,q,gg,rt],sL:[],L:[],
makeprtk(x,a,b):=block([strg,rtk],
if b=1 then return(x^a)
else (
strg:ssubst("equal","=",base64(string(x^(a/b)))),
strg:ssubst("plus","+",strg),
strg:ssubst("slash","/",strg),
rtk:eval_string(strg),
if not(member(rtk,sL)) then sL:append(sL,[rtk]),
if oddp(b) then L:append(L,[(rtk^b=x^a)])
else L:append(L,[(rtk^b=x^a)%and(rtk>=0)]),
return(rtk))),
g:applyb1(f,powertk4,powertk3,powertk2,powertk1),
q:makelist([E,s],s,sL),/* print([g,sL,q]),
if is(is(p2etime>0)=true) then p2etime:p2etime+1 else p2etime:1, */
if L=[] then return(f)
else (gg:(substpart("%and",L,0))%and(g),prt(sL),
subst(makelist(part(sL,k)=concat(rt,k),k,length(sL)),
ev(qe(q,gg),noeval))
)
)$
p2e(f):=if is(p2etype=64) then p2e64(f)
elseif is(p2etype=0) then p2e0(f) else p2e1(f)$
/* p2e map */
p2em(f):=block([p2etlocal],
if atom(f) then return(f)
else
p2etlocal(g):=block([],
if atom(g) then g
else
if atom(g) then g
elseif op(g)="#" then "%not"(p2e(substpart("=",g,0)))
elseif (op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=")
then p2e(g)
else g),
scanmap(p2etlocal,f)
/*qe([],scanmap(p2etlocal,f))*/
)$
/* p2e map tarski */
p2emt(f):=qe([],ev(p2em(f),eval))$
p2t(f):=p2emt(f)$
/*******************/
/* nnsolve section */
/*******************/
/* nnsolvexx */
algexact:true$
/* is imaginary */
imgp(x):=not(is(x>=0) or is(x<0))$
/* del formulas containnig imginary unit */
delimg(f):=block([sel3],
if atom(f) then return(f),
sel3(g):=if atom(g) then g
elseif (op(g)="#" or op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=") and (ratsimp(imagpart(lhs(g)))=0) and (ratsimp(imagpart(rhs(g)))=0) then ratsimp(realpart(g))
elseif (op(g)="#" or op(g)="=" or op(g)="<" or op(g)=">" or op(g)="<=" or op(g)=">=") and (imgp(lhs(g))=true or imgp(rhs(g))=true) then false else g,
scanmap(sel3,f))$
Lintersection(l,m):=listify(intersection(setify(l),setify(m)))$
Lminus(l,m):=listify(setdifference(setify(l),setify(m)))$
/* nnsolve: solver for cnj. of eq. */
/* new ver */
delp(L):=block([LL,ek],LL:L,
for ek in L do
(if member(rhs(ek),%rnum_list) then LL:subst(rhs(ek)=lhs(ek),LL)),
LL
)$
/* new ver */
nnsolve(f,[vv]):=block([fs,v,g,gg],
linsolvewarn:false,
if atom(f) then return(f)
elseif op(f)="%and" then fs:args(f)
elseif op(f)="=" then fs:[f] else return(f),
prt("nnsolving..."),prt(fs),%rnum_list:[],
v:Lintersection(vv,listofvars(fs)),
if v=[] then v:listofvars(fs),
g:delimg(algsys(fs,v)),prt([v,g]),
if g=[] then (v:listofvars(fs),g:delimg(algsys(fs,v)),prt([v,g])),
if g=[] or setify(g)={[false]} then return(false),
solvar:v,
return(substpart("%or",map(lambda([l],substpart("%and",l,0)),map(delp,g)),0))
)$
defrule(ru000, (aa) # (bb), true)$
defrule(ru001, (aa) > (bb), true)$
defrule(ru002, (aa) >= (bb), true)$
defrule(ru003, (aa) < (bb), true)$
defrule(ru004, (aa) <= (bb), true)$
defrule(ru006, (aa) = (bb), true)$
/* delete ineq. if f has no eq. then return true */
deli1(f):=(applyb2(f,ru000,ru001,ru002,ru003,ru004))$
deli(f):=block([g],
if atom(f) then return(f),
g:dnf(f),
if op(g)="%or" then substpart("%or",delete(true,map(deli1,args(g))),0) else deli1(g))$
/* delete eq. if f has no ineq. then return true */
dele1(f):=(applyb2(f,ru006))$
dele(f):=block([g],
if atom(f) then return(f),
g:dnf(f),
if op(g)="%or" then substpart("%or",delete(true,map(dele1,args(g))),0) else dele1(g))$
mapply(s,t):=if atom(t) then s(t) elseif op(t)="%or" then map(s,t) else s(t)$
/* solution selector for eq.s and ineq.s */
defrule(ru007, (aa) = (bb),
if Lintersection(v,listofvars(expand(aa-bb)))=[] then true else (aa)=(bb))$
defrule(ru008, (aa) = (bb),
if Lintersection(v,listofvars(expand(aa-bb)))=[] then (aa)=(bb) else true)$
delec1(f):=(applyb1(f,ru007))$
delec(f):=block([g],
if atom(f) then return(f),
g:dnf(f),
if op(g)="%or" then substpart("%or",delete(true,map(delec1,args(g))),0) else delec1(g))$
delev1(f):=(applyb1(f,ru008))$
delev(f):=block([g],
if atom(f) then return(f),
g:dnf(f),
if op(g)="%or" then substpart("%or",delete(true,map(delev1,args(g))),0) else delev1(g))$
/* new ver for ineq. and conj. of eq. */
solsel(f,[v]):=block([eqss,eqs,ins,sel2,ng],
eqss:deli(f),
if v=[] or length(listofvars(eqss))=1 then v:listofvars(eqss),
eqs:delec(eqss),
prt(eqs),
if atom(eqs) then return(f),
if op(eqs)="%or" then (prt("I can't solve disj. of eq."),return(f)),
ins:(dele(f))%and(delev(eqss)),
sel2(g):=block([gg],
/*(g)%and(qe(map(lambda([x],[E,x]),listofvars(g)),p2et((g)%and(ins))))*/
if /* ins=true then g
elseif */ solvar=[] or atom(g) then (g)%and(ins)
else (solvar:map(lhs,(if op(g)="%and" then args(g) else [g])),
gg:(g)%and(qe(map(lambda([x],[E,x]),solvar),p2t((g)%and(ins)))),
ng:qe([],"%not"(qe(map(lambda([x],[E,x]),solvar),p2t(g)))),
prt((ng)%and(f)),
insolsel:true,
gg:(gg)%or(ora(apply(nnsolvexx,flatten([(ng)%and(f),solvar])))),
insolsel:false,
prt("sel2 returned..."),prt(gg),gg
)
),
qq:apply(nnsolve,flatten([eqs,v])),prt(qq),
if ins=true then prt("no ins"),
qq:mapply(sel2,qq),
qq /* for nnsolvexx */
/* nnscan(qq) */ /* for stand alone */
)$
/* nnsolvexx: solver for semi-algebraic system */
nnsolvexx0(f,[v]):=block([g,gg],
prt("pre-simplified formula"),
g:dnf(qe([],p2t(f))),prt(g),
if atom(g) then return(g)
else (prt("necessary condition"),g:dnf(ncond(g)),prt(length(g)),prt(g)),
if op(g)="%or" then (prt("simplifing..."),gg:dnf(map(ncond,g)),prt(length(gg)),prt(gg),nnscan(map(lambda([t],apply(solsel,flatten([t,v]))),gg)))
else nnscan(apply(solsel,flatten([g,v])))
)$
/* total deg */
tdeg0(f):=block([tdegx],
hipow(expand(subst(map(lambda([s],s=tdegx),listofvars(f)),f)),tdegx)
)$
tdeg(f):=apply(max,map(tdeg0,args(expand( lhs(f)-rhs(f) ))))$
/* new ver */
nnsolvexx(f,[vv]):=block([g,v,gg,q],
if not(insolsel=true) then thisin:f,
if is(presimp=off) then prt("pre-simplify:off")
elseif is(presimp=l) then prt("pre-simplify:lineq with nnscan")
else prt("pre-simplify:p2t->qe"),
if is(ncondtype=off) then prt("ncond-type:off")
elseif is(ncondtype=full) then prt("ncond-type:full formula")
else prt("ncond-type:eq. only"),
if is(ncondsimp=off) then prt("ncond-post-simplify:off")
else prt("ncond-post-simplify:nnscandc"),
if is(sineqsimp=nn) then prt("sineq-simplify:nnscan")
else prt("sineq-simplify:lineq"),
if is(postsimp=off) then prt("post-simplify:off")
else prt("post-simplify:nnscanc->nnscand->nnscan"),
prt(""),prt("given formula"),prt(f),
if is(presimp=off) then g:dnf(f)
elseif is(presimp=l) then (prt("pre-simplified dnf"),g:dnf(lineq(f)))
else (prt("pre-simplified dnf"),g:dnf(qe([],p2t(f)))),prt(g),
prt("var for solve"),
v:vv,for xxx in Lminus(vv,listofvars(g)) do v:delete(xxx,v),
if v=[] then v:listofvars(f),prt(v),
if atom(g) then return(bra(g))
else
(if is(ncondsimp=off) then gg:dnf(mapply(ncondeq,g))
else gg:nnscan(dnf(mapply(ncondeq,g))),
prt(length(gg)),prt(gg),prt("solving..."),
q:scanmap(lambda([t],apply(fineq,flatten([t,last(v)]))),
mapply(lambda([t],apply(solsel,flatten([t,v]))),gg)
,bottomup),
prt("solved formula"),prt(q),
collect(ff):=expand(ff),
/*
collect(ff):=block([fff,vc],fff:expand(f),vc:listofvars(fff),
sum(factor(coeff(fff,vc,d))*x^d,d,0,hipow(fff,vc))),
*/
if is(postsimp=off) then
(thisout:scanmap(collect,l2rx(bra(l2r(q)),last(v))),thisout)
else (q:mapply(nnscanc,q),prt("reduced formula"),prt(q),
q:mapply(nnscand,q),prt("reduced formula"),prt(q),
q:nnscan(q),
thisout:scanmap(collect,l2rx(bra(l2r(q)),last(v))),thisout)
)
)$
thischeck():=fullall((thisin)%eq(ora(thisout)))$
/* for scanmap */
sineq(ineq,[vv]):=block([p,v,xxx,sineqx,q,r,s,lc,sol,sk,rt,qq,excond],
if atom(ineq) then return(ineq),
if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return(ineq),
p:ratsimp(lhs(ineq)-rhs(ineq)),
if imgp(p)=true then return(ineq),
if tdeg(num(factor(p)))<=1 and listofvars(denom(factor(p)))=[] then return(ineq),
if length(vv)=1 then (if hipow(p,last(vv))=1 then return(ineq)),
v:vv,for xxx in Lminus(vv,listofvars(p)) do v:delete(xxx,v),
if v=[] then v:listofvars(p),sineqx:last(v),
q:1,r:1,
sol:delimg(ratsimp(solve(p,sineqx))),
for k:1 thru length(sol) do
(sk:part(sol,k),if not(sk=false) then
(q:q*(lhs(sk)-rhs(sk))^(part(multiplicities,k)),
r:r*(lhs(sk)-concat(rt,k))^(part(multiplicities,k)))),
s:factor(quotient(p,q,sineqx)),lc:ratcoef(s,sineqx,hipow(s,sineqx)),lc:lc/abs(lc),
prt("factor type"),prt([p,sineqx,r,s,lc]),
if s=lc then
qq:qe([],qe([],apply(op(ineq),[r*lc,0])))
elseif ((op(ineq)=">") or (op(ineq)="<")) then
(prt("nnsolving..."),
qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %and nnsolvexx0(s#0)))
elseif ((op(ineq)=">=") or (op(ineq)="<=")) then
(prt("nnsolving..."),
qq:(qe([],qe([],apply(op(ineq),[r*lc,0]))) %or nnsolvexx0(s=0))),
prt(qq),
/* prt(qe([],cnf(qq))),*/
for k:1 thru length(sol) do
(sk:part(sol,k),if not(sk=false) then
qq:subst(concat(rt,k)=rhs(sk),qq)),
if is(sineqsimp=off) then qq:qq
elseif is(sineqsimp=ex) then qq:(map(lambda([t],(t)%and(fullex(t))),qq))
elseif is(sineqsimp=l) then qq:lineq(qq)
else qq:nnscan(qq),
excond:qe([[E,sineqx]],p2t(qq)),prt(excond),
if not(is(excond=true)) then
(insolsel:true,
qq:(qq) %or (ora(nnsolvexx((ineq)%and(%not(excond)),sineqx))),
insolsel:false),
prt("sineq returned..."),prt(qq),
qq
)$
/* factor for polynomial inequation */
/* new ver */
fineq(ineq,[vv]):=block([pp,v,xx,yy,dg,lc],
if atom(ineq) or listofvars(ineq)=[] then return("%and"(ineq)),
if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=") then return(ineq),
pp:ratsimp(lhs(ineq)-rhs(ineq)),
if atom(pp) then return(ineq),
prt(vv),
v:vv,for yy in Lminus(vv,listofvars(pp)) do v:delete(yy,v),
if v=[] then v:listofvars(pp),xx:last(v),dg:hipow(pp,xx),prt(xx),
if imgp(pp)=true or dg=0 then return(ineq),lc:ratcoef(pp,xx,dg),/*prt(lc),*/
if listofvars(lc)=[] then (prt("case1"),sineq(apply(op(ineq),[pp,0]),xx))
elseif is(fineqfull=on) then
((prt("case2-1"),(nnsolvexx(lc=0) %and fineq(apply(op(ineq),[expand(pp-lc*xx^dg),0]),xx)))
%or (prt("case2-2"),(fineq(lc>0) %and (sineq(apply(op(ineq),[pp/lc,0]),xx))))
%or (prt("case2-3"),(fineq(lc<0) %and (sineq(apply(op(ineq),[(-1)*pp/lc,0]),xx)))))
else ((prt("case2-1"),((lc=0) %and fineq(apply(op(ineq),[expand(pp-lc*xx^dg),0]),xx)))
%or (prt("case2-2"),((lc>0) %and (sineq(apply(op(ineq),[pp/lc,0]),xx))))
%or (prt("case2-3"),((lc<0) %and (sineq(apply(op(ineq),[(-1)*pp/lc,0]),xx)))))
)$
/* automatic quantify */
fullex(f):=qe(map(lambda([x],[E,x]),listofvars(f)),p2t(f))$
fullall(f):=qe(map(lambda([x],[A,x]),listofvars(f)),p2t(f))$
/* lineq */
defrule(tolineq1, (aa)^(bb), lineqdec(1,aa^bb,bb))$
defrule(tolineq2, (aa)/(bb), lineqdec(2,aa/bb,bb))$
lineq(f):=block([L,lineqdec,q],
if is(lsimp=off) then (prt("lineq off"),return(f)),
prt("linear simplifing..."),
L:[],
lineqdec(m,g,bb):=block([strg],
if (m=1 and denom(bb)=1) or (m=2 and listofvars(bb)=[]) then g
else (
strg:ssubst("equal","=",base64(string(g))),
strg:ssubst("plus","+",strg),
strg:ssubst("slash","/",strg),
L:append(L,[eval_string(strg)=g]),eval_string(strg))),
q:qe([],apply1(f,tolineq1,tolineq2)),L:listify(setify(L)),
prt([q,L]),
/* nnscan(nnscan(subst(L,q)))*/
/* nnscan(subst(L,q)) */
nnscan(map(lambda([t],(t)%and(fullex(t))),dnf(subst(L,q))))
/* subst(L,q) */
)$
/* l2r */
defrule(l2r1, (aa) > (bb),block([p,l,v],p:expand(bb-aa),l:listofvars(p),
if length(l)=1 then
(v:first(l),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v
else bb<aa)
else bb<aa
))$
defrule(l2r2, (aa) >= (bb),block([p,l,v],p:expand(bb-aa),l:listofvars(p),
if length(l)=1 then
(v:first(l),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v
else bb<=aa)
else bb<=aa
))$
defrule(l2r3, (bb) < (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p),
if length(l)=1 then
(v:first(l),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v
else bb<aa)
else bb<aa
))$
defrule(l2r4, (bb) <= (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p),
if length(l)=1 then
(v:first(l),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v
else bb<=aa)
else bb<=aa
))$
defrule(l2r5, (bb) = (aa),block([p,l,v],p:expand(bb-aa),l:listofvars(p),
if length(l)=1 then
(v:first(l),
if hipow(p,v)=1 then v=-coeff(p,v,0)/coeff(p,v,1)
else bb=aa)
else bb=aa
))$
defrule(l2r1x, (aa) > (bb),block([p,v],p:expand(bb-aa),
if length(vv)=1 then
(v:first(vv),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v
else bb<aa)
else bb<aa
))$
defrule(l2r2x, (aa) >= (bb),block([p,v],p:expand(bb-aa),
if length(vv)=1 then
(v:first(vv),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v
else bb<=aa)
else bb<=aa
))$
defrule(l2r3x, (bb) < (aa),block([p,v],p:expand(bb-aa),
if length(vv)=1 then
(v:first(vv),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<v
else bb<aa)
else bb<aa
))$
defrule(l2r4x, (bb) <= (aa),block([p,v],p:expand(bb-aa),
if length(vv)=1 then
(v:first(vv),
if hipow(p,v)=1 and coeff(p,v,1)>0 then v<=-coeff(p,v,0)/coeff(p,v,1)
elseif hipow(p,v)=1 and coeff(p,v,1)<0 then -coeff(p,v,0)/coeff(p,v,1)<=v
else bb<=aa)
else bb<=aa
))$
defrule(l2r5x, (bb) = (aa),block([p,v],p:expand(bb-aa),
if length(vv)=1 then
(v:first(vv),
if hipow(p,v)=1 then v=-coeff(p,v,0)/coeff(p,v,1)
else bb=aa)
else bb=aa
))$
defrule(l2r6, (bb) < (aa) %and (cc) < (bb), (cc) < (bb) %and (bb) < (aa))$
l2r(f):=block([],apply1(f,l2r1,l2r2,l2r3,l2r4,l2r5))$
l2rx(f,[vv]):=block([],apply1(f,l2r1x,l2r2x,l2r3x,l2r4x,l2r5x,l2r6))$
/* bra, ora */
bra(f):=block([g],
if listp(f) then return(f),
if atom(f) then return([[f]]),
g:dnf(f),
if op(g)="%or" then g:args(g) else g:[g],
map(lambda([t],if op(t)="%and" then args(t) else [t]),g)
)$
ora(f):=block([],
if atom(f) then return(f)
else substpart("%or",map(lambda([g],substpart("%and",g,0)),f),0)
)$
/* lfactor */
lfactor0(ff):=block([f,L,M,cf],L:[],M:[],
if not(listp(ff)) then f:bra(ff) else f:ff,
cf:delete({},map(lambda([t],apply(intersection,args(t))),sort(rest(listify( setdifference(powerset(setify(map(setify,f))),map(lambda([s],{s}),fullsetify(f))) )),comparelength2))),
if cf=[] then return(f),
cf:first(cf),
for fk in f do
if subsetp(cf,setify(fk)) then
L:append([substpart("%and",setdifference(setify(fk),cf),0)],L)
else M:append([fk],M),
append(M,[flatten([substpart("%or",L,0),args(cf)])])
)$
lfactor(ff):=block([f,afr,bfr],
if not(listp(ff)) then f:bra(ff) else f:ff,
afr:f,bfr:[],
for cc:1 unless length(bfr)=length(afr) do (bfr:afr,afr:lfactor0(bfr)),
factor(afr))$
/* plineq */
plineq(f,[p]):=sort(append(Lminus(f,part(f,p)),bra(lineq(ora(part(f,p))))))$
/* addex */
addex(L,[v]):=l2r(map(
lambda([Lk],append(Lk,first(bra(qe(makelist([E,vk],vk,v),p2t(substpart("%and",Lk,0))))))),
L))$
/* qex, solvex */
qex(L,[f]):=if listp(L) then qe(L,apply(p2t,f)) else p2t(L)$
solvex(f,[v]):=apply(nnsolvexx,flatten([ora(bra(f)),v]))$
ncond 動かない式を動かす
以前にも述べましたが，変数の個数を低減させた必要条件を付加する ncond は強力で，nnss，nnsolvexx の求解能力はこの関数によるヒントの賜物とも言えます．
(%i98) f:(a^2+b^2+c^2<=3)%and(a*b*c>=1)$ Evaluation took 0.0100 seconds (0.0200 elapsed) using 538.359 KB. (%i99) qe([],f); Evaluation took 0.3500 seconds (0.5500 elapsed) using 17.918 MB. (%o99) (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (a*b*c-1 = 0) (%i100) nnss(f); Evaluation took 17.1500 seconds (23.0500 elapsed) using 952.370 MB. (%o100) ((a-1 = 0) %and (b-1 = 0) %and (c-1 = 0)) %or ((a-1 = 0) %and (b+1 = 0) %and (c+1 = 0)) %or ((a+1 = 0) %and (b-1 = 0) %and (c+1 = 0)) %or ((a+1 = 0) %and (b+1 = 0) %and (c-1 = 0))
verbose モード．
(%i101) prt(x):=print(x)$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 200 bytes.
(%i102) nnss(f);
pre-simplified formula
2
(a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (a*b*c-1 = 0)
necessary condition
a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0
a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0
b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0
((a-1 = 0) %and (a-1 <= 0) %and (a+1 >= 0))
%or ((a-1 <= 0) %and (a+1 = 0) %and (a+1 >= 0))
((b-1 = 0) %and (b-1 <= 0) %and (b+1 >= 0))
%or ((b-1 <= 0) %and (b+1 = 0) %and (b+1 >= 0))
((c-1 = 0) %and (c-1 <= 0) %and (c+1 >= 0))
%or ((c-1 <= 0) %and (c+1 = 0) %and (c+1 >= 0))
true
8
((a-1 = 0) %and (a-1 <= 0) %and (a+1 >= 0) %and (b-1 = 0) %and (b-1 <= 0)
%and (b+1 >= 0) %and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0)
%and (c-1 = 0) %and (c-1 <= 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 = 0) %and (a-1 <= 0) %and (a+1 >= 0) %and (b-1 = 0)
%and (b-1 <= 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 <= 0)
%and (c+1 = 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 = 0) %and (a-1 <= 0) %and (a+1 >= 0) %and (b-1 <= 0)
%and (b+1 = 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 = 0)
%and (c-1 <= 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 = 0) %and (a-1 <= 0) %and (a+1 >= 0) %and (b-1 <= 0)
%and (b+1 = 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 <= 0)
%and (c+1 = 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 <= 0) %and (a+1 = 0) %and (a+1 >= 0) %and (b-1 = 0)
%and (b-1 <= 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 = 0)
%and (c-1 <= 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 <= 0) %and (a+1 = 0) %and (a+1 >= 0) %and (b-1 = 0)
%and (b-1 <= 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 <= 0)
%and (c+1 = 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 <= 0) %and (a+1 = 0) %and (a+1 >= 0) %and (b-1 <= 0)
%and (b+1 = 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 = 0)
%and (c-1 <= 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
%or ((a-1 <= 0) %and (a+1 = 0) %and (a+1 >= 0) %and (b-1 <= 0)
%and (b+1 = 0) %and (b+1 >= 0)
%and (a^2*b^4+a^4*b^2-3*a^2*b^2+1 = 0) %and (c-1 <= 0)
%and (c+1 = 0) %and (c+1 >= 0) %and (a*b*c-1 = 0)
%and (a^2*c^4+a^4*c^2-3*a^2*c^2+1 = 0)
%and (b^2*c^4+b^4*c^2-3*b^2*c^2+1 = 0))
simplifing...
{(a-1 = 0) %and (b-1 = 0) %and (c-1 = 0),
(a-1 = 0) %and (b+1 = 0) %and (c+1 = 0),
(a+1 = 0) %and (b-1 = 0) %and (c+1 = 0)}
{(a-1 = 0) %and (b-1 = 0) %and (c-1 = 0),
(a-1 = 0) %and (b+1 = 0) %and (c+1 = 0),
(a+1 = 0) %and (b+1 = 0) %and (c-1 = 0)}
{(a-1 = 0) %and (b-1 = 0) %and (c-1 = 0),
(a+1 = 0) %and (b-1 = 0) %and (c+1 = 0),
(a+1 = 0) %and (b+1 = 0) %and (c-1 = 0)}
{(a-1 = 0) %and (b+1 = 0) %and (c+1 = 0),
(a+1 = 0) %and (b-1 = 0) %and (c+1 = 0),
(a+1 = 0) %and (b+1 = 0) %and (c-1 = 0)}
{a-1 = 0,b-1 = 0}
{a-1 = 0,c-1 = 0}
{b-1 = 0,c-1 = 0}
{a-1 = 0,b+1 = 0}
{a-1 = 0,c+1 = 0}
{b+1 = 0,c+1 = 0}
{a+1 = 0,b-1 = 0}
{a+1 = 0,c+1 = 0}
{b-1 = 0,c+1 = 0}
{a+1 = 0,b+1 = 0}
{a+1 = 0,c-1 = 0}
{b+1 = 0,c-1 = 0}
Evaluation took 18.5100 seconds (24.4700 elapsed) using 954.651 MB.
(%o102) ((a-1 = 0) %and (b-1 = 0) %and (c-1 = 0))
%or ((a-1 = 0) %and (b+1 = 0) %and (c+1 = 0))
%or ((a+1 = 0) %and (b-1 = 0) %and (c+1 = 0))
%or ((a+1 = 0) %and (b+1 = 0) %and (c-1 = 0))
nnssdn 二重否定の利用
次の不等式を簡約します．
(%i90) f:x^3+y^3+z^3>=3*x*y*z; Evaluation took 0.0000 seconds (0.0000 elapsed) using 1.859 KB. (%o90) z^3+y^3+x^3 >= 3*x*y*z
まずは，qepmax から．．．．
(%i91) qe([],f); Evaluation took 0.0800 seconds (0.2700 elapsed) using 1.674 MB. (%o91) (z+y+x >= 0) %or (z^2-y*z-x*z+y^2-x*y+x^2 = 0)
一方，nnss は．．．．いい感じです．
(%i92) nnss(f); Evaluation took 1.4100 seconds (5.1900 elapsed) using 59.782 MB. (%o92) ((y-x = 0) %and (z-x = 0)) %or (z+y+x >= 0)
では，次はどうでしょう．
(%i93) g:x^3+y^3+z^3>3*x*y*z; Evaluation took 0.0000 seconds (0.0000 elapsed) using 1.859 KB. (%o93) z^3+y^3+x^3 > 3*x*y*z (%i94) qe([],g); Evaluation took 0.0700 seconds (0.2900 elapsed) using 1.510 MB. (%o94) (z+y+x > 0) %and (z^2-y*z-x*z+y^2-x*y+x^2 > 0) (%i95) nnss(g); Evaluation took 1.6900 seconds (3.7200 elapsed) using 92.056 MB. (%o95) (z+y+x > 0) %and (z^2-y*z-x*z+y^2-x*y+x^2 > 0)
共倒れです．．．．ということで，nnssdn の出番です．
(%i96) nnssdn(g); Evaluation took 1.8000 seconds (5.7800 elapsed) using 78.075 MB. (%o96) ((y-x # 0) %or (z-x # 0)) %and (z+y+x > 0)
勿論，万能ではありません．
(%i97) nnssdn(f); Evaluation took 2.5600 seconds (4.8100 elapsed) using 147.524 MB. (%o97) (z+y+x >= 0) %or (z^2-y*z-x*z+y^2-x*y+x^2 = 0)
nnss 見易い式を求めて
次の不等式を簡約します．
(%i86) f:(x^2-1)*(x^2-4)*(x^2-9)<0; Evaluation took 0.0000 seconds (0.0000 elapsed) using 1.344 KB. (%o86) (x^2-9)*(x^2-4)*(x^2-1) < 0
まずは，qepmax + nns です．悪くはありませんが，ちょっと解り難いですね．
(%i87) qe([],f);
Evaluation took 0.1000 seconds (0.3300 elapsed) using 2.794 MB.
(%o87) (x-3 < 0) %and (x-2 # 0) %and (x-1 # 0) %and (x+1 # 0) %and (x+2 # 0)
%and (x+3 > 0)
%and ((x-2 > 0) %or (x+2 < 0) %or ((x-1 < 0) %and (x+1 > 0)))
(%i88) nns(%);
Evaluation took 4.9000 seconds (6.1500 elapsed) using 232.032 MB.
(%o88) (x-3 < 0) %and (x+3 > 0)
%and ((x-2 > 0) %or (x+2 < 0) %or ((x-1 < 0) %and (x+1 > 0)))
ということで，nnss の出番です．
(%i89) nnss(f);
Evaluation took 4.4000 seconds (7.2900 elapsed) using 194.955 MB.
(%o89) ((x-3 < 0) %and (x-2 > 0)) %or ((x-1 < 0) %and (x+1 > 0))
%or ((x+2 < 0) %and (x+3 > 0))
dnf から cnf へ
式の簡約は，まず選言形に直し，得られた連言項毎に処理するのが基本です．よって結果は
(%i70) qe([],(a^2=1)%and(b^2=1));
Evaluation took 0.3400 seconds (0.5200 elapsed) using 12.210 MB.
(%o70) (a-1 <= 0) %and (a+1 >= 0)
%and (((a-1 = 0) %and (b-1 = 0))
%or ((a-1 = 0) %and (b+1 = 0))
%or ((a+1 = 0) %and (b-1 = 0))
%or ((a+1 = 0) %and (b+1 = 0))) %and (b-1 <= 0)
%and (b+1 >= 0)
(%i71) nns(%);
Evaluation took 4.5500 seconds (5.5200 elapsed) using 185.040 MB.
(%o71) ((a-1 = 0) %and (b-1 = 0)) %or ((a-1 = 0) %and (b+1 = 0))
%or ((a+1 = 0) %and (b-1 = 0))
%or ((a+1 = 0) %and (b+1 = 0))
のようになります．これは連立方程式の解の表示として標準的ですが，連言形のほうが見易そうです．実際
(%i72) cnf(%);
Evaluation took 0.2900 seconds (0.3300 elapsed) using 9.648 MB.
(%o72) ((a-1 = 0) %or (a+1 = 0)) %and ((a-1 = 0) %or (a+1 = 0) %or (b-1 = 0))
%and ((a-1 = 0) %or (a+1 = 0) %or (b-1 = 0)
%or (b+1 = 0))
%and ((a-1 = 0) %or (a+1 = 0) %or (b+1 = 0))
%and ((a-1 = 0) %or (b-1 = 0) %or (b+1 = 0))
%and ((a+1 = 0) %or (b-1 = 0) %or (b+1 = 0))
%and ((b-1 = 0) %or (b+1 = 0))
(%i73) nns(%);
Evaluation took 3.9800 seconds (5.3500 elapsed) using 177.244 MB.
(%o73) ((a-1 = 0) %or (a+1 = 0)) %and ((b-1 = 0) %or (b+1 = 0))
といった具合です．勿論，初めからこの結果を目指して
(%i74) map(nnsolve,(a^2=1)%and(b^2=1)); Evaluation took 0.0100 seconds (0.0100 elapsed) using 494.719 KB. (%o74) ((a = -1) %or (a = 1)) %and ((b = -1) %or (b = 1))
とすることも出来ます．
