ご冗談でしょう，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]
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)
)$
cineqs4 原理
cineqs4 は「真理値に関する中間値定理」により，原始式同士の結合の様子に立ち入ることなく，系を解いています．すなわち，．．．
一変数半代数系 F(x) と実数 a，b (a
例えば，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) が得られます．
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
nns シリーズ更新
・一変数半代数系の高速ソルバー cineqs2 の正式公開．
・solvex から呼ぶ不等式ソルバー sineq も一変数なら高速になりました（cineqs2 の一世代前のアルゴリズムなので少し遅いですが）．
load(qepmax)$
showtime:on$
display2d:false$
prtoff():=kill(prt)$
prton():=(prt([x]):=apply(print,x))$
/*********/
/* 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(ne2inq, (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"),
if is(ncondtype=x) then
(
exs:fullmap(lambda([x],[E,x]),
sort(rest(full_listify(powerset(setify(ncondvars)))),comparelength2)),
for exk in exs do
(q:qe(exk,(f)),prt([L,q]),
if is(not(qe([],(L)%implies(q))=true)) then L:(L)%and(q))
)
else
(
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)) or (listofvars(x)=[] and imagpart(x)#0)$
/* del formulas containnig imginary unit */
delimg(f):=block([sel3],
if atom(f) then return(f),
sel3(g):=block([],
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., separate ng type */
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 (prt("verifing..."),prt(g),
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:(ng)%and(qe([],"%not"(qe(map(lambda([x],[E,x]),solvar),p2t(g))))),
prt("sel2 returned..."),prt(gg),gg
)
),
qq:apply(nnsolve,flatten([eqs,v])),prt(qq),
if ins=true then prt("no ins"),
ng:f,
qq:mapply(sel2,qq),
if is((ng)#false) then
(prt("other type..."),prt(ng),
/* ncondtype:full,ncondsimp:off, */
inner:true,
qq:(qq)%or(ora(apply(nnsolvexx,flatten([ng,solvar])))),
/* ncondtype:eq,ncondsimp:nn, */
inner:false),
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):=block([g],g:expand( lhs(f)-rhs(f)),
if atom(g) then tdeg0(g) else apply(max,map(tdeg0,args(g)) )
)$
/* new ver */
nnsolvexx(f,[vv]):=block([t0,g,v,gg,q,dele1q],
if not(inner=true) then (thisin:f,t0:elapsed_real_time()),
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")
elseif is(ncondtype=x) then prt("ncond-type:x")
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(applyb1(lineq(f),ne2inq)))
else
(prt("pre-simplified dnf"),g:dnf(applyb1(qe([],p2t(f)),ne2inq))),prt(g),
if atom(g) then return(bra(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),ncondvars: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),prt(elapsed_real_time()-t0),
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))),
*/
deleq:dele(q),
/* if length(Lminus(listofvars(dele1q),listofvars(facts())))=1 then */
if length(listofvars(deleq))=1 then
(thisout:mapply(lambda([g],(cineqs(dele1(g))) %and (deli1(g))),dnf(q)),
if is(deleq=dnf(q)) then thisout:(cineqs(thisout)),
thisout:bra(l2r(thisout)),thisout)
elseif 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((qex(thisin))%eq(qex(ora(thisout))))$
/* for scanmap */
sineq(ineq,[vv]):=block([p,v,xxx,sineqx,q,r,s,lc,sol,sols,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,sols:[],
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),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)),
sols:append(sols,makelist(lhs(sk)=rhs(sk),j,part(multiplicities,k))))),
s:factor(quotient(p,q,sineqx)),
lc:ratcoef(expand(s),sineqx,hipow(expand(s),sineqx)),/* lc:lc/abs(lc), */
prt("factor type"),prt([p,sineqx,r,s,lc]),
if listofvars(q)=[sineqx] then qq:cineq(sols,lc,op(ineq))
else
(
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
(inner:true,
qq:(qq) %or (ora(nnsolvexx((ineq)%and(%not(excond)),sineqx))),
inner:false)
),
prt("sineq returned..."),prt(qq),
qq
)$
/* with cineq */
sineq(ineq,[vv]):=block([p,v,xxx,sineqx,q,r,s,lc,sol,sols,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 listofvars(p)=[] then (prt("const."),return("%and"(apply(op(ineq),[p,0])))),
if (tdeg(num(factor(p)))=1 and listofvars(denom(factor(p)))=[]
and length(listofvars(expand(p)))=1) then
(prt("linear"),/*********************/
v:first(listofvars(p)),lc:ratcoef(expand(p),v,hipow(expand(p),v)),
if lc>0 then return(apply(op(ineq),[v,expand((lc*v-p)/lc)]))
else return(apply(op(ineq),[expand((lc*v-p)/lc),v]))),
/* 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,sols:[],
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),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)),
sols:append(sols,makelist(lhs(sk)=rhs(sk),j,part(multiplicities,k))))),
s:factor(quotient(p,q,sineqx)),
lc:ratcoef(expand(s),sineqx,hipow(expand(s),sineqx)),/* lc:lc/abs(lc), */
prt("factor type"),prt([p,sineqx,r,s,lc]),
if listofvars(q)=[sineqx] then qq:cineq(sols,lc,op(ineq))
/*
elseif length(delete(sineqx,listofvars(q)))=1 then qq:tineq(sols,lc,op(ineq))
*/
else
(
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
(inner:true,
qq:(qq) %or (ora(nnsolvexx((ineq)%and(%not(excond)),sineqx))),
inner:false)
),
prt("sineq returned..."),prt(qq),
qq
)$
defrule(delsame, ((aa) < (bb)) %and ((bb) < (aa)), false)$
defrule(eqsame, ((aa) <= (bb)) %and ((bb) <= (aa)), (bb)=(aa))$
sineq2(ineq,sineqx):=block([p,sol,sols,conds,cfs,cqs,lc,crs,css,cts],
/* sineq2(ineq,sineqx):=block([], */
p:(lhs(ineq)-rhs(ineq)),prt(p),
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol),
sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol),
sols:[],for c in sol do if not(last(c)=false) then
sols:append(sols,[
[(lhs(first(c))-rhs(first(c)))^(second(c)),
(lhs(first(c))-concat(rt,length(sols)+1))^(second(c)),
concat(rt,length(sols)+1)=rhs(first(c)),
last(c)]]),prt(sols),
conds:listify(setify(makelist(last(s),s,sols))),prt(conds),
conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)),
cfs:makelist([c,append([[[1,1],t0=t0]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then [rest(s,-2),part(s,3)]),s,sols)))],c,conds),prt(cfs),
cqs:map(lambda([s],([first(s),substpart("*",map(first,last(s)),0),delete(t0=t0,map(last,last(s)))])),cfs),prt(cqs),
lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)),
crs:map(lambda([s],[first(s), lc*second(second(s)) ,factor(quotient(p,first(second(s)),sineqx))/lc ,last(s)]),cqs),prt(crs),
css:map(lambda([s],apply(if op(ineq)="<" or op(ineq)=">" then "%and" else "%or",[(first(s)) %and (subst(last(s), dele( subst(if op(ineq)="<" or op(ineq)=">" then [] else [">"=">=","<"="<="], qex(apply(op(ineq),[second(s),0])))))), if op(ineq)="<" or op(ineq)=">" then ora(nnsolvexx((first(s))%and(third(s)#0))) else ora(nnsolvexx((first(s))%and(third(s)=0))) ])),crs),prt("raw sols:"),prt(css),
cts:map(nns,l2rx(map(nnscand,css),sineqx)),prt("reduced sols (cts):"),prt(cts),
substpart("%or",cts,0)
)$
/* new ver */
/* sineq3(ineq,sineqx):=block([], */
sineq3(ineq,sineqx):=block([p,k,c,sol,sols,conds,cfs,crs,cqs,lc,css,cts],
p:(lhs(ineq)-rhs(ineq)),prt(p),
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol),
sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol),
sols:[],for c in sol do if not(last(c)=false) then
sols:append(sols,[
[[(lhs(first(c))-rhs(first(c)))^(second(c))],
[if second(c)=1 then lhs(first(c))=rhs(first(c)) else
makelist(lhs(first(c))=rhs(first(c)),second(c))],
last(c)]]),prt(sols),
conds:listify(setify(makelist(last(s),s,sols))),prt(conds),
conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)),
cfs:makelist([c,append([append,[[1],[]]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then rest(s,-1)),s,sols)))],c,conds),
cqs:map(lambda([s],([first(s), substpart("*",first(apply(map,last(s))),0), last(apply(map,last(s)))])),cfs),prt(cqs),
lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)),
crs:map(lambda([s],[first(s), factor(quotient(p,second(s),sineqx))/lc, last(s)]),cqs),prt("crs:"),prt(crs),
ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off,
css:(apply(append,l2r(map(lambda([s],
if op(ineq)="<" or op(ineq)=">" then
block([a1],a1:nnssdn(p2t((first(s))%and(second(s)>0))),
brand(
tineq2(third(s),lc,op(ineq)),
sepaxa(if member(sineqx,listofvars(a1)) then bra(a1) else [[a1]]
,sineqx),
sineqx
)
)
else
block([a1],a1:nnss(p2t((first(s))%and(second(s)=0))),
append(
brand(tineq2(third(s),lc,op(ineq)), [[ first(s) ]], sineqx),
sepaxa(if member(sineqx,listofvars(a1)) then bra(a1) else [[a1]]
,sineqx)
)
)
),crs)))),
ncondtype:ncondtype_out,postsimp:postsimp_out,
prt("raw sols:"),prt(css),css
/*cts:map(nns,l2r(map(nnscand,css))),prt("reduced sols (cts):"),prt(cts),
substpart("%or",cts,0)*/
)$
/*
sineq3(x^4-2*a*x^2<b,x);
*/
sepaxa(braform,xxx):=
map(lambda([l],
block([L,M,e],L:[],M:[],for e in l do
if member(xxx,listofvars(e)) then L:append(L,[e]) else M:append(M,[e]),
[L,M])),map(flatten,braform))$
brand(bra1,bra2,v):=block([s1,s2],
sepaxa(apply(append,makelist(makelist(append(s1,s2),s2,bra2),s1,bra1)),v)
)$
braor(bra1,bra2):=append(bra1,bra2)$
/* sineq5(ineq,[vv]):=block([], */
sineq5(ineq,[vv]):=block([sineqx,p,k,c,sol,sols,conds,cfs,crs,cqs,lc,css,cts],
if length(listofvars(ineq))<=1 then
return(map(lambda([s],[s,[]]),cineqs2("%and"(ineq)))),
if vv#[] then sineqx:first(vv),p:(lhs(ineq)-rhs(ineq)),prt(p),
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol),
sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol),
sols:[],for c in sol do if not(last(c)=false) then
sols:append(sols,[
[[(lhs(first(c))-rhs(first(c)))^(second(c))],
[if second(c)=1 then lhs(first(c))=rhs(first(c)) else
makelist(lhs(first(c))=rhs(first(c)),second(c))],
last(c)]]),prt(sols),
conds:listify(setify(makelist(last(s),s,sols))),prt(conds),
conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)),
crs:makelist([c,append([],delete(false,makelist((if (fullall((c)%implies(last(s)))) then second(s)),s,sols)))],c,conds),
prt("crs"),prt(crs),
lc:ratcoef(expand(p),sineqx,hipow(expand(p),sineqx)),
ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off,
css:(apply(append,l2r(
map(lambda([s],brand(tineq2(second(s),lc,op(ineq)),[[first(s)]],sineqx))
,crs)))),
ncondtype:ncondtype_out,postsimp:postsimp_out,
prt("raw sols:"),prt(css),css
)$
/* sineq4(ineq,[vv]):=block([], */
sineq4(ineq,[vv]):=block([p,v,sineqx,hd,k,c,sol,sols,conds,cfs,cqs,lc,crs,css,cts],
/* pre-checks */
if atom(ineq) then return([[ineq]]),
if not(op(ineq)=">" or op(ineq)="<" or op(ineq)=">=" or op(ineq)="<=")
then return([[ineq]]),
p:(lhs(ineq)-rhs(ineq)),prt(p),
if imgp(p)=true then return([[ineq]]),
if listofvars(p)=[]
then (prt("const."),return([["%and"(apply(op(ineq),[p,0]))]])),
/* def. of sineqx */
v:vv,for xxx in Lminus(vv,listofvars(p)) do v:delete(xxx,v),
if v=[] then v:listofvars(p),sineqx:first(v),
/* linear */
hd:hipow(expand(p),sineqx),lc:ratcoef(expand(p),sineqx,hd),
if hd=1 and listofvars(lc)=[] then
(prt("linear"),
if lc>0 then return([[apply(op(ineq),[sineqx,expand((lc*sineqx-p)/lc)])]])
else return([[apply(op(ineq),[expand((lc*sineqx-p)/lc),sineqx])]])),
/* main */
sol:delimg(ratsimp(solve(lhs(ineq)-rhs(ineq),sineqx))),prt(sol),
if listofvars(p)=[sineqx] then
return(bra(cineqs(cineq(
delete(false,flatten(
makelist(if not(part(sol,k)=false) then
makelist(lhs(part(sol,k))=rhs(part(sol,k)),part(multiplicities,k)),
k,length(sol)))),
lc,op(ineq))))),
sol:makelist([part(sol,k),part(multiplicities,k),qex([[E,sineqx]],part(sol,k))],k,length(sol)),prt(sol),
sols:[],for c in sol do if not(last(c)=false) then
sols:append(sols,[
[[(lhs(first(c))-rhs(first(c)))^(second(c))],
[if second(c)=1 then lhs(first(c))=rhs(first(c)) else
makelist(lhs(first(c))=rhs(first(c)),second(c))],
last(c)]]),prt(sols),
conds:listify(setify(makelist(last(s),s,sols))),prt(conds),
conds:delete(false,map(qex, delete(false,flatten(subst("Or"="[",applyb2(applyb2(xreduce("And",map(lambda([f],(f) Or (qex(%not(f)))),append(conds,[true]))),ru05a,ru08a),ru00a)))))),prt("case conds (conds):"),prt(l2r(conds)),
cfs:makelist([c,append([append,[[1],[]]],delete(false,makelist((if (fullall((c)%implies(last(s)))) then rest(s,-1)),s,sols)))],c,conds),
cqs:map(lambda([s],([first(s), substpart("*",first(apply(map,last(s))),0), last(apply(map,last(s)))])),cfs),prt(cqs),
crs:map(lambda([s],[first(s), factor(quotient(p,second(s),sineqx))/lc, last(s)]),cqs),prt(crs),
ncondtype:off,postsimp:off,
css:l2r(map(lambda([s],
if op(ineq)="<" or op(ineq)=">" then
(ora(tineq2(third(s),lc,op(ineq)))) %and ora(nnsolvexx((first(s))%and(second(s)#0)))
else ( (first(s)) %and (ora(tineq2(third(s),lc,op(ineq)))) )
%or (ora(nnsolvexx((first(s))%and(second(s)=0))))),crs)),
kill(ncondtype,postsimp),prt("raw sols:"),prt(css),
map("[",css)
/*cts:map(nns,l2r(map(nnscand,css))),prt("reduced sols (cts):"),prt(cts),
substpart("%or",cts,0)*/
)$
/* one type */
cineq(sols,lc,gl):=block([L,M],
L:sort(sols,lambda([s,t],rhs(s)>=rhs(t))),
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:append([first(L)],L),
L:makelist(
if k=1 then rhs(part(L,k))<lhs(part(L,k))
elseif evenp(k) and k<length(L)
then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k)))
elseif evenp(k) and k=length(L)
then lhs(part(L,k))<rhs(part(L,k))
else false,
k,length(L)),
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:rest(L),
if gl=">=" or gl="<=" then L:subst("<"="<=",L),
L:reverse(L),prt(L),
apply1(substpart("%or",L,0),delsame,eqsame)
)$
cineqsd(l):=block([v,ll,lr],
if listofvars(l)=[] then return(first(l)),
v:first(listofvars(l)),ll:{},lr:{},
for f in l2r(l) do
(if rhs(f)=v then ll:append(ll,{f}) else lr:append(lr,{f})),
ll:sort(listify(ll),lambda([s,t],lhs(s)<lhs(t) or (lhs(s)=lhs(t) and op(s)="<=" and op(t)="<"))),
lr:sort(listify(lr),lambda([s,t],rhs(s)<rhs(t) or (rhs(s)=rhs(t) and op(s)="<" and op(t)="<="))),
prt([ll,lr]),
if ll=[] then return(first(lr))
elseif lr=[] then return(last(ll))
elseif lhs(last(ll))<rhs(first(lr)) then (last(ll)) %and (first(lr))
elseif (op(last(ll))="<=" and op(first(lr))="<=" and lhs(last(ll))=rhs(first(lr))) then lhs(first(lr))=rhs(first(lr)) /* (last(ll)) %and (first(lr)) */
else false)$
cineqsc(l):=block([v,ll,lr],
if listofvars(l)=[] then return(first(l)),
v:first(listofvars(l)),ll:{},lr:{},
for f in l2r(l) do
(if rhs(f)=v then ll:append(ll,{f}) else lr:append(lr,{f})),
ll:sort(listify(ll),lambda([s,t],lhs(s)<lhs(t) or (lhs(s)=lhs(t) and op(s)="<=" and op(t)="<"))),
lr:sort(listify(lr),lambda([s,t],rhs(s)<rhs(t) or (rhs(s)=rhs(t) and op(s)="<" and op(t)="<="))),
prt([ll,lr]),
if ll=[] then return(last(lr))
elseif lr=[] then return(first(ll))
elseif rhs(last(lr))<lhs(first(ll)) then (last(lr)) %or (first(ll))
elseif (op(last(lr))="<" and op(first(ll))="<" and rhs(last(lr))=lhs(first(ll))) then (last(lr)) %or (first(ll)) /* lhs(last(lr))#rhs(last(lr)) */
else true)$
cineqsdm(L):=(prt("cineqs-dnf returned..."),
substpart("%or",map(cineqsd,bra(L)),0))$
cineqscm(L):=(prt("cineqs-cnf returned..."),
substpart("%and",map(cineqsc,brac(ora(L))),0))$
defrule(eq2inq, (aa) = (bb), (((aa)<=(bb)) %and ((bb)<=(aa))) )$
cineqs(L):=cineqsdm(cineqscm(apply1(L,ne2inq,eq2inq)))$
cineqs2(F):=block([lv,v,imagunit,L,M,N,k,K,J],
lv:listofvars(F),if length(lv)#1 then return([[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),["<",">","<=",">=","=","#"]) then map(rhs,delete(false,delimg(solve(substpart("=",e,0),v))))
else e),
F)))),"<"),
if L=[] then return([[is(subst(v=0,F))]]),
L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"),
/* */
M:map(lambda([e],is(ratsimp((subst(v=e,F))))),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 return([[is(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]))
)$
ratprint:false$
/*
asksign
*/
cineqs2f(F):=block([lv,v,imagunit,L,M,N,k,K],
lv:listofvars(F),if length(lv)#1 then return(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),["<",">","<=",">=","=","#"]) 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([[is(subst(v=0,F))]]),
L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"),
M:map(lambda([e],is(subst(v=float(e),F))),L),prt([L,M]),
if length(setify(M))=1 then return([[is(subst(v=0,F))]]),
L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]),
N:[],for k:3 thru length(M)-2 do (K:[],
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:delete([],append(N,[K]))),
N)$
/* 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),[lhs(ineq)-rhs(ineq),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(mapply(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))$
l2rxsort(f,g):=member(x,listofvars(f)) or not(member(x,listofvars(g)))$
/* 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)
)$
brac(f):=block([g],
if listp(f) then return(f),
if atom(f) then return([[f]]),
g:cnf(f),
if op(g)="%and" then g:args(g) else g:[g],
map(lambda([t],if op(t)="%or" then args(t) else [t]),g)
)$
ora(f):=block([],
if not(listp(f)) then return(f),
if atom(f) then return(f)
else substpart("%or",map(lambda([g],substpart("%and",g,0)),f),0)
)$
orac(f):=block([],
if not(listp(f)) then return(f),
if atom(f) then return(f)
else substpart("%and",map(lambda([g],substpart("%or",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]))$
/* tineq */
sorted_cineq(sols,lc,gl):=block([L,M],
L:sols,
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:append([first(L)],L),
L:makelist(
if k=1 then rhs(part(L,k))<lhs(part(L,k))
elseif evenp(k) and k<length(L)
then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k)))
elseif evenp(k) and k=length(L)
then lhs(part(L,k))<rhs(part(L,k))
else false,
k,length(L)),
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:rest(L),
if gl=">=" or gl="<=" then L:subst("<"="<=",L),
L:reverse(L),prt(L),
apply1(substpart("%or",L,0),delsame,eqsame)
)$
sorted_cineq5(sols,lc,gl):=block([L,M],
L:sols,
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:append([first(L)],L),
L:makelist(
if k=1 then rhs(part(L,k))<lhs(part(L,k))
elseif evenp(k) and k<length(L)
then (rhs(part(L,k+1))<lhs(part(L,k+1))) %and (lhs(part(L,k))<rhs(part(L,k)))
elseif evenp(k) and k=length(L)
then lhs(part(L,k))<rhs(part(L,k))
else false,
k,length(L)),
if (lc>0 and (gl="<" or gl="<=")) or (lc<0 and (gl=">" or gl=">="))
then L:rest(L),
if gl=">=" or gl="<=" then L:subst("<"="<=",L),
L:reverse(L),prt(L),
apply1(substpart("%or",L,0),delsame,eqsame)
)$
tineq2(sols,lc,gl):=block([x,PLx,delx,PL,IL,CL,M,N],
/*tineq2(sols,lc,gl):=block([],*/
if sols=[] then return([[["%and"(apply(gl,[lc,0]))],[ ]]]),
x:lhs(first(sols)),
PLx:map(flatten,listify(permutations(sols))),prt(PLx),
delx(l):=map(lambda([xa],rhs(xa)),l),PL:map(delx,PLx),prt(PL),
IL:flatten(map(lambda([s],substpart("%and",makelist(part(s,k)-part(s,k+1)>=0,k,length(s)-1),0)),PL)),prt(IL),
ncondtype_out:ncondtype,postsimp_out:postsimp,ncondtype:off,postsimp:off,
CL:map(lambda([s],ora(l2r(solvex(s)))),IL),
ncondtype:ncondtype_out,postsimp:postsimp_out,prt("case condition:"),prt(CL),
CL:map(lambda([s],if atom(s) then s elseif op(s)="=" then false else s),CL),
if not(casecheck=off) then (
M:[],CL:makelist(if member(e,M) then false elseif is(substpart("%or",makelist(qex((e)%implies(m)),m,M),0)=true) then false else (M:append(M,[e]),e),e,CL),
M:[],CL:reverse(makelist(if member(e,M) then false elseif is(substpart("%or",makelist(qex((e)%implies(m)),m,M),0)=true) then false else (M:append(M,[e]),e),e,reverse(CL)))
),prt(CL),
N:delete(false,makelist(if part(CL,j)#false then [[sorted_cineq(part(PLx,j),lc,gl)],[part(CL,j)]],j,length(CL))),
/* N:sepaxa(map(flatten,bra(ora(N))),x), */
thistineq2:[PLx,PL,IL,CL,N],
prt("tineq2 returned..."),prt(N),N
/* M:map(lambda([l],sorted_cineq(l,lc,gl)),PLx),
prt(M),map("[",map("%and",M,CL))*/
)$
