一変数有理冪半代数系ソルバー 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]))
)$
