読者です 読者をやめる 読者になる 読者になる

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