1.分解体の絶対定義多項式
今回のオリジナルプログラム https://github.com/YasuakiHonda/GaloisGroupSolver が行う主な処理は
1.分解体の絶対定義式
2.Galois 群
3.組成列
4.添加する元の定義式
5.分解体の相対定義式
の算出です.種々の方法の比較の意味も含め,私家版を作成しましたので,順次公開させて頂きます.
まず,オリジナルプログラムの1.では,入力の多項式 p の次数 が低い場合を想定して私が気軽に提案した(数値根と 次式の因数分解を用いる)方法を採用して頂いたのですが,これは明らかにハイコストであり,例えば, の処理に数分を要します.他にも,primitive element と p の根との関係の調整はユーザーが行わねばなりません.
これに対し,下記の mp は有理数体 に p の根を1つずつ添加し,その都度,拡大体の 上の定義式(絶対定義式)の次数を上げてゆく(代数体上の 次以上の既約因子から primitive element を生成する),分解体の定義式の計算としては一般的なものです.これによれば,primitive element と p の根との関係(*)は自動的に得られ,また, の場合には,絶対定義式の次数が を超えた段階で,非可解と判定できます.
mp の入力 p は, を変数とする 上の 次以上の monic な既約多項式であり,出力は p の分解体(単純拡大)の primitive element A の 上の最小多項式の変数に A を代入した形の式 DA です.また,内部では,p の全ての根を A で表した式のリスト RA,および,(*)つまり,A=KK.makelist(RA[i],i,1,length(KK)) を満たすリスト KK も生成し,これらはすべてグローバルに定義されているので,随時参照できます.
mp(p):=block([],DEG:hipow(p,x), [Dx,DA]:divide(p,x-A,x),t:2,K:1,KK:[1],F:[x-A],F2:[], for i:0 unless K=0 do ( Dx:num(factor(Dx,DA)), printn("001"), /*Dx:num(algfac(Dx,DA)), printn("001"),*/ F3:[],if op(Dx)="*" then (map(lambda([s],if hipow(rat(s),x)=1 then push(s,F2) else push(s,F3)),args(Dx)),Dx:apply("*",F3),if length(F2)+length(F)=DEG then (DX:subst(A=X,DA),return(DA))), [Dx,DB]:divide(Dx,x-B,x),push(x-B,F), t:t+1,for j:0 while hipow(rat(gcd(f:resultant(DA,DBX:subst(B=(X-A)/(K:(-1)^t*floor(t/2)),DB),A),diff(f,X),X)),X)>0 do t:t+1, printn("002"), push(K,KK), /*print(f),*/ DX:factor(content(f,X)[2]), printx(DX), printn("003"), DX:if op(num(DX))="*" then args(DX)[1] else DX, if DEG=5 and hipow(DX,X)>20 then (DA:0,print(concat(string(p)," is not solvable.")),return()), printx(DX), tellrat(DX),algebraic:on,gcd:subres, AX:gcd(DBX,DA,A), printn("004"), ABX:solve([AX,A+K*B-X],[A,B]), [Dx,F,F2]:rat(subst(ABX,[Dx,F,F2])), algebraic:off,untellrat(X), [DA,Dx,F,F2]:subst(X=A,[DX,Dx,F,F2]), printn(i),printx([DA,Dx,F,F2]), if hipow(rat(Dx),x)<=1 then K:0 ), RA:sublist(flatten([F,Dx,F2]),lambda([s],not(numberp(s)))), RA:map(lambda([s],rhs(solve(s,x)[1])),RA), DX:subst(A=X,DA), return(DA) )$ PL:[x^2-2,x^3-3*x-1,x^4-2,x^4+x^2-1,x^4-2*x^3+2*x^2+2,x^4+2*x^3+3*x^2+4*x+5,x^4+x+1,x^5-2,x^5-5*x+12,x^5+20*x+32,x^5+11*x+44,x^5+x^4-4*x^3-3*x^2+3*x+1,x^5+100*x^2+1000,x^6+x^3+1,x^6-2,x^5+10*x^3-10*x+4,x^5+20*x+16]$ /* verbose mode on, off */ pon():=(printn(s):=print(s),printx(s):=print(s))$ poff():=(printn(s):=0,printx(s):=0)$
上記を例えば,mp.mac として保存し,load("mp.mac")$ とした後,for i:1 thru 15 do print([p:PL[i],mp(p)])$ と入力すれば,次のように出力されます.
[x^2-2,A^2-8] [x^3-3*x-1,A^3-3*A-1] [x^4-2,A^8+28*A^4+2500] [x^4+x^2-1,A^8+10*A^6+47*A^4+110*A^2+841] [x^4-2*x^3+2*x^2+2,A^12+4*A^10+24*A^8+48*A^6-560*A^4+3136] [x^4+2*x^3+3*x^2+4*x+5, A^24+24*A^23+336*A^22+3344*A^21+25740*A^20+159984*A^19+820856*A^18 +3519504*A^17+12721926*A^16+39075680*A^15+104485896*A^14+257189424*A^13 +603068156*A^12+1264487184*A^11+1484791560*A^10-3707413456*A^9 -23515353279*A^8-53513746296*A^7-7075256024*A^6+299352120960*A^5 +770653544880*A^4+869309952000*A^3+1145273500800*A^2+1451723788800*A +1818528595200] [x^4+x+1, A^24-80*A^20+340*A^18+7520*A^16+23120*A^14-973378*A^12-462400*A^10 +50899280*A^8+74190340*A^6+67773664*A^4+2114616240*A^2+266962921] [x^5-2,A^20+2500*A^10+50000] [x^5-5*x+12,A^10-10*A^8-75*A^6+1500*A^4-5500*A^2+16000] [x^5+20*x+32,A^10-20*A^8+100*A^6+2000*A^4-32000*A^2+128000] [x^5+11*x+44,A^10-22*A^8+77*A^6+4356*A^4-53724*A^2+189728] [x^5+x^4-4*x^3-3*x^2+3*x+1,A^5+A^4-4*A^3-3*A^2+3*A+1] [x^5+100*x^2+1000, A^20+250000*A^14+20000000*A^12+625000000*A^10-5300000000*A^8+700000000000*A^6 +18750000000000*A^4-598000000000000*A^2+4205000000000000] [x^6+x^3+1,A^6+A^3+1] [x^6-2,A^12-1012*A^6+19307236] Evaluation took 3.3130 seconds (3.3360 elapsed) using 1579.350 MB.
非可解な例を for i:16 thru 17 do print([p:PL[i],mp(p)])$ と入力すると...
x^5+10*x^3-10*x+4 is not solvable. [x^5+10*x^3-10*x+4,0] x^5+20*x+16 is not solvable. [x^5+20*x+16,0] Evaluation took 5.8240 seconds (5.8760 elapsed) using 3379.830 MB.
maxima の組み込み関数 splitfield との比較.
(%i3) for i:1 thru 15 do splitfield(PL[i],x)$ Evaluation took 9.5910 seconds (9.6560 elapsed) using 3121.559 MB. (%i4) for i:1 thru 15 do (mp(PL[i]),RA)$ Evaluation took 3.2100 seconds (3.2320 elapsed) using 1579.139 MB.
正規部分群の構成
今回のオリジナルプログラムでは各元の一点集合の conjugate closure から極大正規部分群を選んでいますが,それが必ずしも正しい結果を与えないことも,例えば
G:[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5], [2,1,3,4,5,6],[2,1,3,4,6,5],[2,1,4,3,5,6],[2,1,4,3,6,5]];
[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5]]
のように指摘されています.
一方,共役類に分割された群 とその部分群 について,次の2つは等価でした.
(1) は の正規部分群である.
(2) となる空でない がある.
従って,有限群 を共役類に分割すれば, の位数 の正規部分群の全体は,元の個数の合計が ,かつ,和が積閉である共役類の和の全体に一致します.
以下,この性質による組成列計算のコードを記します.なお,mnsg では powerset を用いましたが,入力の群のサイズが程々(笑)なら共役類も多くはないので,爆発しないと思われます(母関数を用いると,元の個数の合計が になる組のみを抽出できますが,…).
/* 置換の積,逆元 */ mul(A,B):=map(lambda([s],A[s]),B)$ inv(A):=map(second,sort(map("[",A,sort(A)),lambda([s,t],s[1]<t[1])))$ /* 共役類分割 */ CC(GG):=block([g,C:[],G:setify(GG),H], unless G={} do (g:listify(G)[1], H:setify(map(lambda([s],mul(s,mul(g,inv(s)))),GG)), G:setdifference(G,H),push(H,C)),map(listify,C))$ /* 積閉判定 */ isG(GG):=is(setify(apply('append,makelist(makelist(mul(i,j),i,GG),j,GG)))=setify(GG))$ /* 極大正規部分群 */ mnsg(GG):=block([S:[],PS:map(lambda([s],apply('append,listify(s))),listify(powerset(setify(CC(GG)))))], for d in rest(reverse(listify(divisors(length(GG))))) while S=[] do S:sublist(PS,lambda([s],length(s)=d and isG(s))),S[1])$ /* Galois 群の組成列 */ cs(GG):=block([S:[GG],m:GG],unless length(m)=1 do push(m:mnsg(m),S), [S:reverse(S),S:map(length,S),makelist(S[i-1]/S[i],i,2,length(S))])$
例えば,冒頭のGに対して,cs(G); と入力すれば
[[[[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5], [2,1,3,4,5,6],[2,1,3,4,6,5],[2,1,4,3,5,6],[2,1,4,3,6,5]], [[1,2,3,4,5,6],[1,2,3,4,6,5],[1,2,4,3,5,6],[1,2,4,3,6,5]], [[1,2,3,4,5,6],[1,2,3,4,6,5]],[[1,2,3,4,5,6]]],[8,4,2,1],[2,2,2]]
と出力(最後の2つは各部分群の位数,隣接するそれらの比の値)されます.
次の位数9の部分群は,冒頭の例と同じく,singleton の conjugate closure にはならないものです.
[x^6+x^5+x^4+x^3-4*x^2+5, [[[[1,2,3,4,5,6],[1,3,2,5,4,6],[3,4,1,6,2,5],[3,1,4,2,6,5],[2,4,1,3,6,5], [2,1,4,6,3,5],[4,2,3,5,1,6],[4,3,2,1,5,6],[3,5,1,2,6,4],[3,1,5,6,2,4], [2,5,1,6,3,4],[2,1,5,3,6,4],[5,2,3,1,4,6],[5,3,2,4,1,6],[3,5,4,6,2,1], [3,4,5,2,6,1],[2,5,4,3,6,1],[2,4,5,6,3,1],[1,6,3,5,4,2],[1,3,6,4,5,2], [1,6,2,4,5,3],[1,2,6,5,4,3],[6,4,1,2,3,5],[6,1,4,3,2,5],[4,6,3,1,5,2], [4,3,6,5,1,2],[4,6,2,5,1,3],[4,2,6,1,5,3],[6,5,1,3,2,4],[6,1,5,2,3,4], [5,6,3,4,1,2],[5,3,6,1,4,2],[5,6,2,1,4,3],[5,2,6,4,1,3],[6,5,4,2,3,1], [6,4,5,3,2,1]], [[1,2,3,4,5,6],[1,2,6,5,4,3],[1,3,2,5,4,6],[1,6,3,5,4,2],[4,2,6,1,5,3], [4,3,2,1,5,6],[4,6,3,1,5,2],[5,2,6,4,1,3],[5,3,2,4,1,6],[5,6,3,4,1,2], [1,3,6,4,5,2],[1,6,2,4,5,3],[4,2,3,5,1,6],[5,2,3,1,4,6],[4,3,6,5,1,2], [4,6,2,5,1,3],[5,3,6,1,4,2],[5,6,2,1,4,3]], [[1,2,3,4,5,6],[1,3,6,4,5,2],[1,6,2,4,5,3],[4,2,3,5,1,6],[5,2,3,1,4,6], [4,3,6,5,1,2],[5,6,2,1,4,3],[4,6,2,5,1,3],[5,3,6,1,4,2]], [[1,2,3,4,5,6],[1,3,6,4,5,2],[1,6,2,4,5,3]],[[1,2,3,4,5,6]]],[36,18,9,3,1], [2,2,3,3]]]
冪乗に分解できるタイプ
今回のオリジナルプログラムでは,例えば,SolveSolvable(x^2-2)$ の根の表示には [[alpha[1],alpha[1]^2-8]] が現れるので,SolveSolvable(x^2-8)$ と問うと,[[alpha[1],alpha[1]^2-32]] が現れるので,...となってしまいます.
プログラムの趣意に反する可能性もありますが,例えば,x^2-2,x^6+x^3+1,x^20+x^15+3*x^5+4 といった冪乗に分解できるタイプについては,それぞれ x-2,x^2+x+1,x^4+x^3+3*x+4 の根の表示を得た上で,その各根の2,3,5乗根のすべてを出力するコードを書いてみました.
まず,1の原始N乗根の最小多項式を出力する関数を用意します.
ROU(x,N):=block([D:reverse(listify(divisors(N)))], num(factor((x^pop(D)-1)/apply("*",map(lambda([s],x^s-1),D)))))$ /* テスト */ makelist(ROU(x,i),i,1,10);
次に SolveSolvable の定義の局所変数の宣言をコメントアウトした SolveSolvable2.mac をリロードして,以下を実行すると,上記の分解を挟んだ結果が得られます.
/* サンプル選択 */ p:[x^2-2,x^6+x^3+1,x^20+x^15+3*x^5+4][2]; /* 分解(polydecompの出力は気まぐれなので...) */ for i:(N:hipow(p,x)) step -1 unless polynomialp(q:rat(subst(x=x^(1/i),p)),[x]) do N:i-1$[q,x^N]; /* ベースの方程式の求根 */ if hipow(q,x)=1 then [C,SolN,x[1]]:[[],1,rhs(solve(q,x)[1])] else (SolveSolvable(q),C:SI[StageN]@ExtensionList)$ /* 各根x[i]のN乗根y[i]に1の原始N乗根の冪乗を掛ける(得られるもの全体はy[i]の選び方に拠らない) */ for i:1 thru SolN do C:push([y[i],y[i]^N-x[i]],C)$push([w,ROU(w,N)],C); RROOTS:flatten(makelist(makelist(y[i]*w^j,j,0,N-1),i,1,SolN)); /* 検算 */ ef_polynomial_reduction(P:apply("*",map(lambda([s],x-s),RROOTS)),C);
Groebner基底の利用(その1)
jurupapaさん( http://maxima.hatenablog.jp/ )が公開しておられる「可解な方程式を冪根で解く」プログラムの処理の一部に Groebner 基底を用いてみました.例は以下のものです.
http://maxima.hatenablog.jp/entry/2018/10/09/005917
http://maxima.hatenablog.jp/entry/2018/10/11/012019
http://maxima.hatenablog.jp/entry/2018/10/13/134605
load("SolveSolvable2.mac")$ load("grobner")$ gal_phi(xlist):=ev(sum(xlist[i]*cl[i],i,1,length(xlist)),cl:[1,-1,2,0,3,2,5])$ p1:x^4+2*x^3+3*x^2+4*x+5$ ldc(f,s):=coeff(f,s,hipow(f,s))$ tst():=( PS:gal_init_polynomial_info(p1), gal_minimal_polynomial_V2(PS), gal_sol_V4(PS), gal_galois_group2(PS), get_matrix_permutation_group(PS), FG:new(FiniteGroup), gr_gen_tables(FG,PS@solperm,PS@galois_group_index), NSGS:gr_subnormal_series(FG), g0v:PS@gal_minpoly, C:[], gr1:listify(NSGS[2]@index_elements), gr12:makelist(gr_mult(2,elem,NSGS[2]),elem,gr1), h0:product(x-remainder(ev(V[gr1[i]],PS@vncond),g0v),i,1,length(gr1)), h1:product(x-ev(V[gr12[i]],PS@vncond),i,1,length(gr12)), R:remainder(rat(h0-h1)/2,g0v,V), R:ldc(R,x), [g1v,Da]:poly_reduced_grobner([g0v,R-a],[V,a]), push([a,Da],C), gr2:listify(NSGS[3]@index_elements), gr22:makelist(gr_mult(4,elem,NSGS[3]),elem,gr2), gr23:makelist(gr_mult(4,elem,NSGS[3]),elem,gr22), h0:product(x-remainder(ev(V[gr2[i]],PS@vncond),g1v),i,1,length(gr2)), h1:product(x-remainder(ev(V[gr22[i]],PS@vncond),g1v),i,1,length(gr22)), h2:product(x-remainder(ev(V[gr23[i]],PS@vncond),g1v),i,1,length(gr23)), R:remainder(rat(h0+w*h1+w^2*h2)/3,g1v,V), push([w,Dw:w^2+w+1],C), R:ef_polynomial_reduction(R,C), R:ldc(R,x), [Da,Dw,Db,g2v]:poly_reduced_grobner([Da,Dw,g1v,R-b],[V,b,w,a]), push([b,Db],C), gr3:listify(NSGS[4]@index_elements), gr32:makelist(gr_mult(17,elem,NSGS[4]),elem,gr3), h0:product(x-ev(V[gr3[i]],PS@vncond),i,1,length(gr3)), h1:product(x-ev(V[gr32[i]],PS@vncond),i,1,length(gr32)), R:remainder(rat(h0-h1)/2,g2v,V), R:ef_polynomial_reduction(R,C), R:ldc(R,x), [Da,Dw,Db,Dc,g3v]:poly_reduced_grobner([Da,Dw,Db,g2v,R-c],[V,c,b,w,a]), push([c,Dc],C), gr4:listify(NSGS[5]@index_elements), gr42:makelist(gr_mult(8,elem,NSGS[5]),elem,gr4), h0:product(x-ev(V[gr4[i]],PS@vncond),i,1,length(gr4)), h1:product(x-ev(V[gr42[i]],PS@vncond),i,1,length(gr42)), R:remainder(rat(h0-h1)/2,g3v,V), R:ef_polynomial_reduction(R,C), R:ldc(R,x), [Da,Dw,Db,Dc,g4v,Dd]:poly_reduced_grobner([Da,Dw,Db,Dc,g3v,R-d],[V,d,c,b,w,a]), push([d,Dd],C), RROOTS:[], for i:1 thru 4 do ( x[i]:remainder(ev(solV[i](V),PS@solcond),g1v,V), for gv in [g2v,g3v,g4v] do ( x[i]:remainder(x[i],gv,V)), push((x[i]:ratexpand(ef_polynomial_reduction(x[i],C))),RROOTS)))$ showtime:all$ SolveSolvable(p1)$ tst()$[RROOTS,C];
実行結果です.
Minimal polynomial of V V^24+24*V^23+336*V^22+3344*V^21+25740*V^20+159984*V^19+820856*V^18 +3519504*V^17+12721926*V^16+39075680*V^15+104485896*V^14+257189424*V^13 +603068156*V^12+1264487184*V^11+1484791560*V^10-3707413456*V^9 -23515353279*V^8-53513746296*V^7-7075256024*V^6+299352120960*V^5 +770653544880*V^4+869309952000*V^3+1145273500800*V^2+1451723788800*V +1818528595200 Galois Group of x^4+2*x^3+3*x^2+4*x+5 matrix([a,b,c,d],[a,b,d,c],[a,c,b,d],[a,c,d,b],[a,d,b,c],[a,d,c,b],[b,a,c,d], [b,a,d,c],[b,c,a,d],[b,c,d,a],[b,d,a,c],[b,d,c,a],[c,a,b,d],[c,a,d,b], [c,b,a,d],[c,b,d,a],[c,d,a,b],[c,d,b,a],[d,a,b,c],[d,a,c,b],[d,b,a,c], [d,b,c,a],[d,c,a,b],[d,c,b,a]) Subnormal series and quotients of orders FiniteGroup[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24] FiniteGroup[1,4,5,8,9,12,13,16,17,20,21,24] FiniteGroup[1,8,17,24] FiniteGroup[1,8] FiniteGroup[1] x^4+2*x^3+3*x^2+4*x+5 is solvable. Solutions 'x[1] = (alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000 -(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816 +(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416 +(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000 -(7*alpha[2]^2*alpha[3]*alpha[4])/54080+(alpha[2]*alpha[3]*alpha[4])/104 -(alpha[3]*alpha[4])/16-alpha[4]/4+alpha[3]/4-1/2 'x[2] = (-(alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000) +(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816 -(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416 -(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000 +(7*alpha[2]^2*alpha[3]*alpha[4])/54080-(alpha[2]*alpha[3]*alpha[4])/104 +(alpha[3]*alpha[4])/16+alpha[4]/4+alpha[3]/4-1/2 'x[3] = (-(alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000) +(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816 -(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416 -(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000 +(7*alpha[2]^2*alpha[3]*alpha[4])/54080-(alpha[2]*alpha[3]*alpha[4])/104 +(alpha[3]*alpha[4])/16-alpha[4]/4-alpha[3]/4-1/2 'x[4] = (alpha[1]*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/292032000 -(3*alpha[2]^2*Zeta[3]*alpha[3]*alpha[4])/10816 +(3*alpha[2]*Zeta[3]*alpha[3]*alpha[4])/416 +(23*alpha[1]*alpha[2]^2*alpha[3]*alpha[4])/292032000 -(7*alpha[2]^2*alpha[3]*alpha[4])/54080+(alpha[2]*alpha[3]*alpha[4])/104 -(alpha[3]*alpha[4])/16+alpha[4]/4-alpha[3]/4-1/2 with [[alpha[4], alpha[4]^2+(23*alpha[1]*alpha[2]^2*Zeta[3])/4563000 -(7*alpha[2]^2*Zeta[3])/845-(2*alpha[2]*Zeta[3])/13 +(11*alpha[1]*alpha[2]^2)/2281500+(8*alpha[2]^2)/845 +(6*alpha[2])/13+4], [alpha[3], alpha[3]^2-(11*alpha[1]*alpha[2]^2*Zeta[3])/2281500 -(8*alpha[2]^2*Zeta[3])/845+(8*alpha[2]*Zeta[3])/13 +(alpha[1]*alpha[2]^2)/4563000-(3*alpha[2]^2)/169+(2*alpha[2])/13 +4], [alpha[2], (14*alpha[1]*Zeta[3])/27-1440*Zeta[3]+alpha[2]^3-(19*alpha[1])/135-2120], [Zeta[3],Zeta[3]^2+Zeta[3]+1],[alpha[1],alpha[1]^2-38880000]] Verification Numeric calcuration of the above solutions [(-0.8578967583284902*%i)-1.287815479557649, 1.416093080171908*%i+0.2878154795576489, 0.2878154795576489-1.416093080171908*%i, 0.85789675832849*%i-1.287815479557649] Numeric solutions with allroots( x^4+2*x^3+3*x^2+4*x+5 ) [x = 1.416093080171908*%i+0.287815479557648, x = 0.287815479557648-1.416093080171908*%i, x = 0.8578967583284904*%i-1.287815479557648, x = (-0.8578967583284904*%i)-1.287815479557648] Evaluation took 16.7930 seconds (16.8840 elapsed) using 6984.509 MB. (%i8) Evaluation took 6.1570 seconds (6.1670 elapsed) using 2897.350 MB. Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes. (%o9) [[(a*b^2*c*d*w)/292032000-(3*b^2*c*d*w)/10816+(3*b*c*d*w)/416 +(23*a*b^2*c*d)/292032000-(7*b^2*c*d)/54080 +(b*c*d)/104-(c*d)/16+d/4-c/4-1/2, (-(a*b^2*c*d*w)/292032000)+(3*b^2*c*d*w)/10816-(3*b*c*d*w)/416 -(23*a*b^2*c*d)/292032000+(7*b^2*c*d)/54080 -(b*c*d)/104+(c*d)/16-d/4-c/4-1/2, (-(a*b^2*c*d*w)/292032000)+(3*b^2*c*d*w)/10816-(3*b*c*d*w)/416 -(23*a*b^2*c*d)/292032000+(7*b^2*c*d)/54080 -(b*c*d)/104+(c*d)/16+d/4+c/4-1/2, (a*b^2*c*d*w)/292032000-(3*b^2*c*d*w)/10816+(3*b*c*d*w)/416 +(23*a*b^2*c*d)/292032000-(7*b^2*c*d)/54080 +(b*c*d)/104-(c*d)/16-d/4+c/4-1/2], [[d, (-23*a*b^2*w)+37800*b^2*w+702000*b*w-4563000*d^2-22*a*b^2-43200*b^2 -2106000*b-18252000], [c, 22*a*b^2*w+43200*b^2*w-2808000*b*w-4563000*c^2-a*b^2+81000*b^2 -702000*b-18252000], [b,(-70*a*w)+194400*w-135*b^3+19*a+286200],[w,w^2+w+1], [a,a^2-38880000]]]
CAD-QE on PARI+Singular(実行例)
同じ例を開発中の PARI+Singular 上での CAD-QE で処理すると...
? tst11Sv2s01([ex,ex],[a,b,c,x,y],and,"x^2+y^2=1,a*x+b*y+c=0",2*7);Ans() [y,2] [x,4] [c,6] [b,4] [a,1] var.:a b c x y *** connecting adjacent 40/1608 cells. 1[a < [a,1],true,[a^2+(b^2-c^2),1] <= c <= [a^2+(b^2-c^2),2],true,true] 2[a = [a,1],b < [a^2+b^2,1],[a^2+(b^2-c^2),1] <= c <= [a^2+(b^2-c^2),2],true,true] 3[a = [a,1],b = [a^2+b^2,1],c = [c,1],true,true] 4[a = [a,1],[a^2+b^2,1] < b,[a^2+(b^2-c^2),1] <= c <= [a^2+(b^2-c^2),2],true,true] 5[[a,1] < a,true,[a^2+(b^2-c^2),1] <= c <= [a^2+(b^2-c^2),2],true,true] time = 208 ms.
? tst11Sv2s01([ex,ex,ex],[a,b,c,d,x,y,z],and,"x^2+y^2+z^2=1,a*x+b*y+c*z+d=0",2*7);Ans() [z,2] [y,3] [x,6] [d,11] [c,11] [b,8] [a,1] var.:a b c d x y z *** connecting adjacent 896/107615 cells. 1[a < [a,1],true,true,[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] 2[a = [a,1],b < [2*a^2-b^2,1],true,[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] 3[a = [a,1],b = [2*a^2-b^2,1],c < [a^2+(b^2-c^2),1],[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] 4[a = [a,1],b = [2*a^2-b^2,1],c = [a^2+(b^2-c^2),1],d = [d,1],true,true,true] 5[a = [a,1],b = [2*a^2-b^2,1],[a^2+(b^2-c^2),1] < c,[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] 6[a = [a,1],[2*a^2-b^2,1] < b,true,[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] 7[[a,1] < a,true,true,[a^2+(b^2+(c^2-d^2)),1] <= d <= [a^2+(b^2+(c^2-d^2)),2],true,true,true] time = 26,807 ms.
QEPCAD Version B 1.72
https://www.usna.edu/Users/cs/wcbrown/qepcad/B/QEPCAD.html
上記の通り,更新されていましたので,インストールスクリプトをUPします.
export CAS=~/CAS mkdir $CAS export qe=$CAS/qesource export saclib=$CAS/saclib2.2.7 cd $CAS wget https://www.usna.edu/Users/cs/wcbrown/qepcad/INSTALL/saclib2.2.7.tar.gz tar zxvf ./saclib2.2.7.tar.gz -C ./ cd $saclib/bin ./sconf ./mkproto ./mkmake ./mklib all cd $CAS wget https://www.usna.edu/Users/cs/wcbrown/qepcad/INSTALL/qepcad-B.1.72.tar.gz tar zxvf ./qepcad-B.1.72.tar.gz -C ./ cd $qe sed -i "s/csh/sh/g" $qe/Makefile make sed -i "s/\#SINGULAR/SINGULAR/g" $qe/default.qepcadrc cd $CAS wget https://www.usna.edu/Users/cs/wcbrown/qepcad/SLFQ/simplify-1.20.tar.gz tar zxvf ./simplify-1.20.tar.gz -C ./ cd ./simplify make sudo ln -s $CAS/qesource/bin/qepcad /usr/local/bin/qepcad sudo ln -s $CAS/simplify/slfq /usr/local/bin/slfq
実行例(1)
echo -e "[] \n (a,b,x,y) \n 2 \n (Ex)(Ey)[x^2+y^2<=1 /\ a x+b y=1]. \n finish" | qepcad
======================================================= Quantifier Elimination in Elementary Algebra and Geometry by Partial Cylindrical Algebraic Decomposition Version B 1.72, Sun Mar 18 22:49:08 EDT 2018 by Hoon Hong (hhong@math.ncsu.edu) With contributions by: Christopher W. Brown, George E. Collins, Mark J. Encarnacion, Jeremy R. Johnson Werner Krandick, Richard Liska, Scott McCallum, Nicolas Robidoux, and Stanly Steinberg ======================================================= Enter an informal description between '[' and ']': []Enter a variable list: (a,b,x,y)Enter the number of free variables: 2 Enter a prenex formula: (Ex)(Ey)[x^2+y^2<=1 /\ a x+b y=1]. ======================================================= Before Normalization > finish An equivalent quantifier-free formula: b^2 + a^2 - 1 >= 0 ===================== The End ======================= ----------------------------------------------------------------------------- 0 Garbage collections, 0 Cells and 0 Arrays reclaimed, in 0 milliseconds. 334662 Cells in AVAIL, 500000 Cells in SPACE. System time: 31 milliseconds. System time after the initialization: 17 milliseconds. -----------------------------------------------------------------------------
実行例(2)
echo -e "[] \n (a,b,x,y) \n 2 \n (Ex)(Ey)[x^4+y^3<=1 /\ a x+b y=1]. \n finish" | qepcad +L100000 +N10000000
======================================================= Quantifier Elimination in Elementary Algebra and Geometry by Partial Cylindrical Algebraic Decomposition Version B 1.72, Sun Mar 18 22:49:08 EDT 2018 by Hoon Hong (hhong@math.ncsu.edu) With contributions by: Christopher W. Brown, George E. Collins, Mark J. Encarnacion, Jeremy R. Johnson Werner Krandick, Richard Liska, Scott McCallum, Nicolas Robidoux, and Stanly Steinberg ======================================================= Enter an informal description between '[' and ']': []Enter a variable list: (a,b,x,y)Enter the number of free variables: 2 Enter a prenex formula: (Ex)(Ey)[x^4+y^3<=1 /\ a x+b y=1]. ======================================================= Before Normalization > finish An equivalent quantifier-free formula: a^8 + 6 a^4 - 3 > 0 \/ 256 b^12 - 768 b^9 + 1728 a^4 b^6 + 768 b^6 - 432 a^8 b^3 + 432 a^4 b^3 - 256 b^3 + 27 a^12 - 54 a^8 + 27 a^4 > 0 \/ [ b > 0 /\ 256 b^12 - 768 b^9 + 1728 a^4 b^6 + 768 b^6 - 432 a^8 b^3 + 432 a^4 b^3 - 256 b^3 + 27 a^12 - 54 a^8 + 27 a^4 = 0 ] ===================== The End ======================= ----------------------------------------------------------------------------- 638 Garbage collections, -1442980436 Cells and 12 Arrays reclaimed, in 25548 milliseconds. 881113 Cells in AVAIL, 5000000 Cells in SPACE. System time: 205804 milliseconds. System time after the initialization: 205697 milliseconds. -----------------------------------------------------------------------------
CAD-QE on Risa/Asir(特徴)
CAD の要は,projection と代数体上での根の分離です.
昨年の maxima における実装は,Lazard method と複素数値根のクラスタリングのための pscs の計算(Mathematica 似)との併用というスタイルであり,それは primitive elements の計算(squarefree norm や代数体上の多項式の GCD の計算)が maxima では非常に重いことによる選択でしたが,肝心の多倍長浮動小数点数の計算も maxima では高速とは言えず,求根には外部の MPSolve を用いていました.
そうした経緯から,Lazard method の軽さを損なわない,primitive elements 生成の実装を探した結果,PARI+Singular での実装(後日公開予定,かなり高速)に至る訳ですが,当初,多変数多項式の因数分解(PARI での実装は実験的かつ遅い)などを含むルーチンには,Risa を利用しており,今回公開したコードは,それと Risa の特長の一つである逐次拡大体上の各種関数とを組み合わせた,Risa/Asir(PARI ライブラリも含む)プロパーな CAD-QE 関数です(ので,速度は関心の外ということに).