一変数半代数系の高速ソルバー cineqs2

多変数半代数系のソルバーの基礎部分として作りました．かなり速いです．

以下のものは，パッケージのロードなしで動くように，真偽判定を簡素化，また，論理結合子も and, or になっています．なお，maxima の solve が，虚数なのに %i を使わなかったり，実数なのに %i を用いて根を表示する場合は，等価でない式を返すかも知れません．

cineqs2(F):=block([lv,v,imagunit,L,M,N,k,K,J],
lv:listofvars(F),if length(lv)#1 then return([[is(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],if member(op(F),["and","or"]) then map(is,rat(subst(v=e,F)))
 else is(rat(subst(v=e,F)))),L),
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]))
)$

実行例．

(%i21) showtime:on$display2d:false$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 8 bytes.
Evaluation took 0.0000 seconds (0.0000 elapsed) using 8 bytes.
(%i22) x*(x^2-2)*(x^2-3)<0;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 1.211 KB.
(%o22) x*(x^2-3)*(x^2-2) < 0
(%i23) cineqs2(%);
Evaluation took 0.0100 seconds (0.0100 elapsed) using 336.141 KB.
(%o23) [[x < -sqrt(3)],[-sqrt(2) < x,x < 0],[sqrt(2) < x,x < sqrt(3)]]

(%i24) (((x-1 > 0) and (x^2+x-3 < 0)) or ((x-1 < 0) and (x^2+x-3 > 0)) 
 or ((x < 0) and (x+2 > 0))) and (x-1 # 0) and (x # 0) and (x+2 # 0) and (x^2+x-3 # 0);
Evaluation took 0.0000 seconds (0.0000 elapsed) using 41.609 KB.
(%o24) x-1 > 0 and x^2+x-3 < 0 or x-1 < 0 and x^2+x-3 > 0 or x < 0 and x+2 > 0
(%i25) cineqs2(%);
Evaluation took 0.0100 seconds (0.0000 elapsed) using 388.031 KB.
(%o25) [[x <= -(sqrt(13)+1)/2],[-2 < x,x < 0],[1 < x,x < (sqrt(13)-1)/2]]

(%i26) product((x-sqrt(k))^(mod(k^3,5)),k,1,10)<=0;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 45.141 KB.
(%o26) (x-3)^4*(x-2)^4*(x-1)*(x-sqrt(2))^3*(x-2^(3/2))^2*(x-sqrt(3))^2
              *(x-sqrt(6))*(x-sqrt(7))^3
 <= 0
(%i27) cineqs2(%);
Evaluation took 0.0200 seconds (0.0200 elapsed) using 1.219 MB.
(%o27) [[1 <= x,x <= sqrt(2)],[x = sqrt(3)],[x = 2],
        [sqrt(6) <= x,x <= sqrt(7)],[x = 2^(3/2)],[x = 3]]