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

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