solvex 原理
solvex の処理手順は,人間が計算機を用いて解く場合と同じです.すなわち,...
(ステップ1)presimp,つまり,入力を qex で半代数系に変換
(オプションスイッチ)presimp:l で qex ではなく,p2t を通さない簡約 lineq を利用.sqrt が係数のみに現れる場合など,それを残したまま処理できる.presimp:off で処理を無効化
(ステップ2)結果を dnf に変換
各連言項に対して
(ステップ3)ncondeq で必要条件のうち等式のみを添加
(オプションスイッチ)ncondtype:x でターゲット変数の存在条件のみを添加.ncondtype:full で等式以外も添加.ncondtype:off で処理を無効化.ncondsimp:off で ncondeq 内部での簡約を無効化
(ステップ4)solsel,nnsolve で等式系を処理
(ステップ5)fineq,sineq で不等式系を処理
(オプションスイッチ)sineqsimp:off で sineq 内部での簡約を無効化.sineqsimp:l で sineq 内部での簡約を lineq に変更
結果をまとめて
(ステップ6)postsimp,つまり,cnf,nnscan,dnf,nnscan,再度,nnscan で簡約
(オプションスイッチ)postsimp:off で処理を無効化
(ステップ7)bracket dnf 形式で出力
toy problems
最後の簡約は重いので無効化します(笑).
(%i1) postsimp:off$ Evaluation took 0.0000 seconds (0.0000 elapsed) using 32 bytes.
放物線と直線との交わり.
(%i2) solvex(y=x^2 %and x+y=1,x,y); Evaluation took 5.4600 seconds (11.3400 elapsed) using 315.180 MB. (%o2) [[x = (-sqrt(5)/2)-1/2,y = sqrt(5)/2+3/2], [x = sqrt(5)/2-1/2,y = 3/2-sqrt(5)/2]]
領域と直線との交わり.
(%i3) solvex(y>x^2 %and x+y=1,y); Evaluation took 2.1600 seconds (4.7800 elapsed) using 120.922 MB. (%o3) [[(-sqrt(5)/2)-1/2 < x,x < sqrt(5)/2-1/2,y = 1-x]]
x について解くと.
(%i4) solvex(y>x^2 %and x+y=1,x); Evaluation took 2.2100 seconds (4.9300 elapsed) using 122.540 MB. (%o4) [[3/2-sqrt(5)/2 < y,x = 1-y,y < sqrt(5)/2+3/2]]
放物線と領域との交わり.
(%i5) solvex(y=x^2 %and x+y<1,y); Evaluation took 2.4200 seconds (5.0200 elapsed) using 146.358 MB. (%o5) [[(-sqrt(5)/2)-1/2 < x,x < sqrt(5)/2-1/2,y = x^2]]
x について解くと.
(%i6) solvex(y=x^2 %and x+y<1,x); Evaluation took 6.7700 seconds (11.4300 elapsed) using 462.988 MB. (%o6) [[0 <= y,x = -sqrt(y),y < 3/2],[0 <= y,x = sqrt(y),y < 3/2-sqrt(5)/2], [3/2-sqrt(5)/2 < y,x = -sqrt(y),y < sqrt(5)/2+3/2]]
結果が見難いので,x を含む式とその他とに分離.
(%i7) sepaxa(%,x); Evaluation took 0.0000 seconds (0.0000 elapsed) using 8.836 KB. (%o7) [[[x = -sqrt(y)],[0 <= y,y < 3/2]], [[x = sqrt(y)],[0 <= y,y < 3/2-sqrt(5)/2]], [[x = -sqrt(y)],[3/2-sqrt(5)/2 < y,y < sqrt(5)/2+3/2]]]
解の選択.
(%i8) solvex(3*x^2=2 %and 1<x); Evaluation took 0.2100 seconds (0.5500 elapsed) using 9.020 MB. (%o8) [[false]]
空のようなので,改めて.
(%i9) solvex(2*x^2=3 %and 1<x); Evaluation took 2.9600 seconds (6.4500 elapsed) using 211.316 MB. (%o9) [[x = sqrt(3)/sqrt(2)]]
不等式の連言.
(%i10) solvex(2*x^2<=3 %and 1<x); Evaluation took 2.6700 seconds (3.7200 elapsed) using 174.558 MB. (%o10) [[1 < x,x <= sqrt(3)/sqrt(2)]]
いきなり 3 変数ですが.
(%i11) solvex(a*x^2=b %and a<x); Evaluation took 2.7200 seconds (7.5900 elapsed) using 215.724 MB. (%o11) [[0 < x-a,b = a*x^2]]
変数の指定が必要です.
(%i12) solvex(a*x^2=b %and a<x,x); Evaluation took 6.7700 seconds (19.4000 elapsed) using 487.579 MB. (%o12) [[0 < a,a^3-b < 0,x = sqrt(b/a)],[a < x,a = 0,b = 0], [a < 0,a^3-b < 0,b <= 0,x = -sqrt(b/a)],[a < 0,b <= 0,x = sqrt(b/a)]]
見易く整理.
(%i13) sepaxa(%,x); Evaluation took 0.0000 seconds (0.0000 elapsed) using 12.602 KB. (%o13) [[[x = sqrt(b/a)],[0 < a,a^3-b < 0]],[[a < x],[a = 0,b = 0]], [[x = -sqrt(b/a)],[a < 0,a^3-b < 0,b <= 0]], [[x = sqrt(b/a)],[a < 0,b <= 0]]]
http://www.cybernet.co.jp/maple/documents/pdf/product/maple/maple17/maple17_new_feature31.pdf
の 2 つ目の例.
(%i14) solvex(a*x^2<b %and a<x,x); Evaluation took 31.4200 seconds (47.6300 elapsed) using 2388.419 MB. (%o14) [[0 < a,a < x,a < 0,b < 0],[0 < a,a < x,-sqrt(b/a) < x,x < sqrt(b/a)], [0 < b,a < x,a = 0],[0 < b,a < x,a < 0],[x < -sqrt(b/a),a < x,a < 0], [a < x,sqrt(b/a) < x,a < 0]]
やはり簡約は必要ですね.
(%i15) plineq(%,1); Evaluation took 0.2200 seconds (0.7700 elapsed) using 10.936 MB. (%o15) [[0 < a,a < x,-sqrt(b/a) < x,x < sqrt(b/a)],[0 < b,a < x,a = 0], [0 < b,a < x,a < 0],[a < x,sqrt(b/a) < x,a < 0], [x < -sqrt(b/a),a < x,a < 0],[false]]
見易く整理.
(%i16) sepaxa(bra(ora(%)),x); Evaluation took 0.4100 seconds (0.4100 elapsed) using 30.975 MB. (%o16) [[[a < x,-sqrt(b/a) < x,x < sqrt(b/a)],[0 < a]], [[a < x],[0 < b,a = 0]],[[a < x],[0 < b,a < 0]], [[a < x,sqrt(b/a) < x],[a < 0]],[[a < x,x < -sqrt(b/a)],[a < 0]]]
一変数半代数系の高速ソルバー 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]]
qex 原理(その2)
p2t では abs,max2,min2 も
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)$
により,冪乗関数に変換して処理しています.流れを例示します.
次の式は恒真です.
(%i1) f:max2(min2(a,b),c)=min2(max2(a,c),max2(b,c)); Evaluation took 0.0000 seconds (0.0000 elapsed) using 560 bytes. (%o1) max2(min2(a,b),c) = min2(max2(a,c),max2(b,c))
実際
(%i2) qex(f); Evaluation took 1.9800 seconds (2.7400 elapsed) using 265.126 MB. (%o2) true
となります.前回同様,量化の具合を見ると
(%i3) p2em(f); Evaluation took 0.2200 seconds (0.2200 elapsed) using 22.855 MB. (%o3) qe([[E,rt1],[E,rt2],[E,rt3],[E,rt4],[E,rt5]], (rt1 >= 0) %and (rt1^2 = (b-a)^2) %and (rt2 >= 0) %and ((rt2+((-rt1)+b+a)/2+c)/2 = ((-rt5)+(rt4+c+b)/2+(rt3+c+a)/2)/2) %and (rt2^2 = (c-((-rt1)+b+a)/2)^2) %and (rt3 >= 0) %and (rt3^2 = (c-a)^2) %and (rt4 >= 0) %and (rt4^2 = (c-b)^2) %and (rt5 >= 0) %and (rt5^2 = ((rt4+c+b)/2-(rt3+c+a)/2)^2))
のように5変数です.
これに至る過程は,まず,powertk4 で,f における min2 を abs で
(%i4) applyb1(f,powertk4); Evaluation took 0.0000 seconds (0.0000 elapsed) using 88.875 KB. (%o4) max2(((-abs(b-a))+b+a)/2,c) = ((-abs(max2(b,c)-max2(a,c)))+max2(b,c)+max2(a,c))/2
powertk3 で max2 も abs で
(%i5) applyb1(f,powertk4,powertk3); Evaluation took 0.0000 seconds (0.0100 elapsed) using 809.602 KB. (%o5) (abs(c-((-abs(b-a))+b+a)/2)+c+((-abs(b-a))+b+a)/2)/2 = ((-abs((abs(c-b)+c+b)/2-(abs(c-a)+c+a)/2)) +(abs(c-b)+c+b)/2+(abs(c-a)+c+a)/2)/2
さらに,powertk2 で abs を sqrt で表し,前回述べた方法で束縛変数に置き換えています.
(%i6) applyb1(f,powertk4,powertk3,powertk2); Evaluation took 0.0200 seconds (0.0100 elapsed) using 1.960 MB. (%o6) (sqrt((c-((-sqrt((b-a)^2))+b+a)/2)^2)+c+((-sqrt((b-a)^2))+b+a)/2)/2 = ((-sqrt(((sqrt((c-b)^2)+c+b)/2-(sqrt((c-a)^2)+c+a)/2)^2)) +(sqrt((c-b)^2)+c+b)/2+(sqrt((c-a)^2)+c+a)/2)/2
なお,%o6 を見ると,sqrt の出現は8箇所ですが,sqrt((c-a)^2),sqrt((b-a)^2),sqrt((c-b)^2) は2箇所ずつあるので,5変数に納まっています.これが8変数だと
(%i7) p2em(f),p2etype:0; Evaluation took 0.2800 seconds (0.2900 elapsed) using 30.215 MB. (%o7) qe([[E,rt1],[E,rt2],[E,rt3],[E,rt4],[E,rt5],[E,rt6],[E,rt7],[E,rt8]], (rt1 >= 0) %and (rt1^2 = (b-a)^2) %and (rt2 >= 0) %and (rt2^2 = (b-a)^2) %and (rt3 >= 0) %and ((rt3+((-rt1)+b+a)/2+c)/2 = ((-rt8)+(rt5+c+b)/2+(rt4+c+a)/2)/2) %and (rt3^2 = (c-((-rt2)+b+a)/2)^2) %and (rt4 >= 0) %and (rt4^2 = (c-a)^2) %and (rt5 >= 0) %and (rt5^2 = (c-b)^2) %and (rt6 >= 0) %and (rt6^2 = (c-a)^2) %and (rt7 >= 0) %and (rt7^2 = (c-b)^2) %and (rt8 >= 0) %and (rt8^2 = ((rt7+c+b)/2-(rt6+c+a)/2)^2))
となり,これを評価するのは,...
qex 原理(その1)
qex は入力を関数 p2t を通して,qe に渡すだけの関数なので,以下は,p2t のお話です.
p2t (rational power to Tarski) とは,有理数乗を含む入力を特称量化することで多項式による系に変換する関数で,その原理は,原始式 A と x の既約分数 a/b 乗とに対して
A(x^(a/b)) ⇔ ∃y(y=x^(a/b) ∧ A(y)) ⇔ ∃y( (if b is even then y>=0 ∧) y^b=x^a ∧ A(y) )
となる性質であり,この変換後の式を qe に処理させています.
実際の量化を内部コマンド p2em (power to ex, scanmap ver.) で確認してみましょう.
(%i1) p2em(x<sqrt(2)); Evaluation took 0.0600 seconds (0.0600 elapsed) using 1.698 MB. (%o1) qe([[E,rt1]],(rt1 >= 0) %and (rt1^2 = 2) %and (x < rt1))
次の例では,sqrt(2) の出現は2箇所ですが,束縛変数は rt1 にまとめられています.出現が少ない場合は,出現毎に新たな変数を設ける方が処理が簡潔ですが,式が大きくなると後述のように変数の個数の節減は必須となるので,このようにしました.
(%i2) p2em(x<sqrt(2)/(3^(2/3)-sqrt(2))); Evaluation took 0.1000 seconds (0.1300 elapsed) using 7.201 MB. (%o2) qe([[E,rt1],[E,rt2]], (rt1 >= 0) %and (rt1^2 = 2) %and (rt2^3 = 9) %and (x < rt1/(rt2-rt1)))
一方,次の例では,変数の個数を減らすには量化を選言全体に施した方が良さそうです.
(%i3) p2em(x<-sqrt(3) %or sqrt(3)<x); Evaluation took 0.0700 seconds (0.0800 elapsed) using 2.432 MB. (%o3) (qe([[E,rt1]],(rt1 >= 0) %and (rt1 < x) %and (rt1^2 = 3))) %or (qe([[E,rt1]],(rt1 >= 0) %and (rt1^2 = 3) %and (x < -rt1))) (%i4) qe([],ev(%,eval)); Evaluation took 0.5200 seconds (1.1200 elapsed) using 28.026 MB. (%o4) x^2-3 > 0
実際,この場合には
(%i5) qe([[E,rt1]],(rt1 >= 0) %and (rt1^2 = 3) %and ((x < -rt1) %or (rt1 < x))); Evaluation took 0.2900 seconds (0.4700 elapsed) using 19.952 MB. (%o5) x^2-3 > 0
が得られます.しかし,...
(%i6) p2em(1+sqrt(3)<x %implies 4+2*sqrt(3)<x^2); Evaluation took 0.0400 seconds (0.0400 elapsed) using 1.239 MB. (%o6) qe([[E,rt1]],(rt1 >= 0) %and (rt1+1 < x) %and (rt1^2 = 3)) %implies qe([[E,rt1]],(rt1 >= 0) %and (2*rt1+4 < x^2) %and (rt1^2 = 3)) (%i7) qe([],ev(%,eval)); Evaluation took 0.6100 seconds (1.1800 elapsed) using 28.977 MB. (%o7) true
に対しては
(%i8) qe([[E,rt1]],(rt1 >= 0) %and (rt1^2 = 2) %and ((rt1+1 < x) %implies (2*rt1+4 < x^2))); Evaluation took 0.1700 seconds (0.3700 elapsed) using 3.884 MB. (%o8) (x-1 < 0) %or (x^2-2*x-1 <= 0) %or (x^4-8*x^2+8 > 0)
となり,これは
(%i9) fullall(%); Evaluation took 0.7300 seconds (1.1000 elapsed) using 31.596 MB. (%o9) false
のように恒真にはなりません.その理由は,特称量化が否定を潜ると全称量化となることであり,それが生じるか形か否かを判定させるより,各原始式に量化を map する(安易な?)設計になっています.
solvex 例
ソース http://d.hatena.ne.jp/ehito/20170106/1483705734 を見て戴くと判りますが,solvex には様々なスイッチがあり,説明が必要なのですが,それはまた後日とし(苦),ここではまずどんなものが解けるのか?をご紹介します.
係数にパラメータを含む不等式.
(%i1) solvex(x^2-2*x>a,x); Evaluation took 10.7700 seconds (19.8400 elapsed) using 595.634 MB. (%o1) [[x < 1-sqrt(a+1)],[sqrt(a+1)+1 < x],[a < -1]]
係数に冪乗根を含む不等式.
(%i2) solvex(x^2<=sqrt(3),x); Evaluation took 11.4000 seconds (19.3700 elapsed) using 809.069 MB. (%o2) [[-3^(1/4) <= x,x <= 3^(1/4)]]
次の例では計算時間短縮のため,幾つかの簡約を off にしています.
(%i3) solvex(x^2-sqrt(8)*x-1<0,x),presimp:off,ncondtype:off,postsimp:off; Evaluation took 6.2600 seconds (9.4200 elapsed) using 333.642 MB. (%o3) [[sqrt(2)-sqrt(3) < x,x < sqrt(3)+sqrt(2)]]
連立不等式.
(%i4) solvex([[x^2>a,x>1]],x); Evaluation took 16.9000 seconds (32.2600 elapsed) using 989.834 MB. (%o4) [[sqrt(a) < x,1 < x],[1 < x,a < 0]]
主係数にパラメータを含む方程式.
(%i5) solvex(a*x=b,x); Evaluation took 6.8400 seconds (25.5400 elapsed) using 424.958 MB. (%o5) [[a = 0,b = 0],[b = 0,x = b/a],[b # 0,x = b/a]]
簡約が不充分なので...
(%i6) nnscandc(ora(%)); Evaluation took 2.9800 seconds (5.6700 elapsed) using 252.291 MB. (%o6) ((a = 0) %and (b = 0)) %or (x = b/a)
(%i7) solvex([[a*x+b*y=c,x+y=1]],x,y); Evaluation took 32.3900 seconds (69.1900 elapsed) using 2624.115 MB. (%o7) [[b-a = 0,c-a = 0,y = 1-x],[x = b/(b-a)-c/(b-a),y = c/(b-a)-a/(b-a)]]
こちらの方が見易い?
(%i8) factor(%); Evaluation took 0.0000 seconds (0.0000 elapsed) using 101.938 KB. (%o8) [[b-a = 0,c-a = 0,y = -(x-1)],[x = -(c-b)/(b-a),y = (c-a)/(b-a)]]
絶対値を含む方程式.
(%i9) solvex(abs(x-a)+abs(x-b)=c,x); Evaluation took 111.1600 seconds (149.5500 elapsed) using 12527.014 MB. (%o9) [[0 <= c-b+a,0 <= c+b-a,x = (-c/2)+b/2+a/2], [0 <= c-b+a,0 <= c+b-a,x = c/2+b/2+a/2], [(-c/2)+b/2+a/2 <= x,c-b+a = 0,x <= c/2+b/2+a/2], [(-c/2)+b/2+a/2 <= x,c+b-a = 0,x <= c/2+b/2+a/2]]
このような式には logical factor 関数(開発中)が有効です.
(%i10) lfactor(%); Evaluation took 0.0200 seconds (0.0300 elapsed) using 3.229 MB. (%o10) [[(c-b+a = 0) %or (c+b-a = 0),(-c/2)+b/2+a/2 <= x,x <= c/2+b/2+a/2], [(x = (-c/2)+b/2+a/2) %or (x = c/2+b/2+a/2),0 <= c-b+a,0 <= c+b-a]]
こちらの方が見易い?
(%i11) factor(%); Evaluation took 0.0100 seconds (0.0100 elapsed) using 2.139 MB. (%o11) [[(c-b+a = 0) %or (c+b-a = 0),-(c-b-a)/2 <= x,x <= (c+b+a)/2], [(x = -(c-b-a)/2) %or (x = (c+b+a)/2),0 <= c-b+a,0 <= c+b-a]]
主係数にパラメータを含む不等式.
(%i12) solvex(a*x>b,x); Evaluation took 6.3000 seconds (14.1600 elapsed) using 505.691 MB. (%o12) [[0 < a,b/a < x],[x < b/a,a < 0],[a = 0,b < 0]]
無理関数で表された不等式.
(%i13) solvex(sqrt(x)>x-1,x); Evaluation took 14.7900 seconds (26.1800 elapsed) using 659.776 MB. (%o13) [[0 <= x,x < 1],[3/2-sqrt(5)/2 < x,x < sqrt(5)/2+3/2]]
簡約が不充分なので...
(%i14) nnscandc(ora(%)); Evaluation took 10.7700 seconds (16.4800 elapsed) using 428.886 MB. (%o14) (0 <= x) %and (x < sqrt(5)/2+3/2)
高次不等式.
(%i15) solvex(x*(x^2-2)*(x^2-3)<0,x); Evaluation took 45.8000 seconds (86.6800 elapsed) using 1963.461 MB. (%o15) [[-sqrt(2) < x,x < 0],[sqrt(2) < x,x < sqrt(3)],[x < -sqrt(3)]]
eqx の出力を入力にすると...
(%i16) qex(x=sqrt(2)-sqrt(3)); Evaluation took 0.6100 seconds (1.0100 elapsed) using 27.529 MB. (%o16) (x < 0) %and (x^2-3 < 0) %and (x^4-10*x^2+1 = 0) (%i17) solvex(%,x); Evaluation took 6.8200 seconds (15.1000 elapsed) using 284.881 MB. (%o17) [[x = sqrt(2)-sqrt(3)]]
方程式と不等式との系.
(%i18) solvex([[x+y=1,y>2*x,x^2+y^2=2]],y,x); Evaluation took 13.3600 seconds (24.1100 elapsed) using 652.512 MB. (%o18) [[x = 1/2-sqrt(3)/2,y = sqrt(3)/2+1/2]]
係数にパラメータを含む...
(%i19) solvex([[x+y>1,y>2*x,x^2+y^2=a]],y,x); Evaluation took 40.5400 seconds (67.2100 elapsed) using 2450.850 MB. (%o19) [[1/2-sqrt(2*a-1)/2 < x,x < sqrt(a/5),y = sqrt(a-x^2)]]
これまでの qe では
(%i20) qe([[E,x],[E,y]],(x^2+y^2<3*y) %and (y>1) %and (x+y=a)); Evaluation took 0.4100 seconds (0.6200 elapsed) using 29.980 MB. (%o20) (4*a^2-12*a-9 < 0) %and ((a > 0) %or (a^2-2*a-1 < 0))
ここまででしたが,solvex なら
(%i21) solvex(%); Evaluation took 24.3600 seconds (41.4200 elapsed) using 1047.936 MB. (%o21) [[0 < a,a < 3/sqrt(2)+3/2],[1-sqrt(2) < a,a < sqrt(2)+1]]
と解けます.でもまだ簡約が不充分なので
(%i22) nnscandc(ora(%)); Evaluation took 13.0000 seconds (20.4000 elapsed) using 545.835 MB. (%o22) (1-sqrt(2) < a) %and (a < 3/sqrt(2)+3/2)
合成して実行.
(%i23) solvex(qex([[X1,x]],sqrt(x)+a=abs(x-b)),a); Evaluation took 25.9000 seconds (43.3900 elapsed) using 1935.237 MB. (%o23) [[-b < a,b < a],[1/4 < b,a = -sqrt(b)],[a = (-b)-1/4,b <= 1/4]]
等価性
nns,nnss と同じく,qex,solvex の入力と出力とは等価,つまり,変数に任意の実数を代入したとき,それらの true,false は一致します.
例えば,sqrt を含む場合,定義:y=sqrt(x)⇔(y^2=x∧y>=0) により,x<0 のとき,y=sqrt(x) は false ですから,前回の最後の例の入力:[ [x^2=a,x>-1] ],出力:[ [a < 1,x = -sqrt(a)],[x = sqrt(a)] ] の等価性は,a<0 のとき入出力とも false,1<=a のとき出力は [ [false],[x = sqrt(a)] ],つまり,[ [x = sqrt(a)] ] といった具合に確かめられますし,xa 座標平面における入力,出力が表す図形が一致することからも判るでしょう.
システム自身に等価性をチェックさせるには,solvex の実行直後に
(%i36) thischeck(); Evaluation took 1.9700 seconds (3.1300 elapsed) using 117.474 MB. (%o36) true
とすればよく,入出力は
(%i37) thisin; Evaluation took 0.0000 seconds (0.0000 elapsed) using 48 bytes. (%o37) (x > -1) %and (x^2 = a) (%i38) thisout; Evaluation took 0.0000 seconds (0.0000 elapsed) using 48 bytes. (%o38) [[a < 1,x = -sqrt(a)],[x = sqrt(a)]]
のように参照出来ます.