２．Galois 群
現在のオリジナルプログラムでは，Galois 群の元と分解体の絶対定義式の根との対応の把握に若干の検索を行っています．
これに対して，下記の nGG では，１．で得た RA における A が絶対定義式 DA の任意の根であるという点に着目して，それらの根と「入力 p の根の置換」との対応を直接求めています．すなわち，１．で得た DA の数値根 （絶対誤差は p の根差の半分より小とする）を
の A に代入したものは
の順列
であり，この
が絶対定義式の根
に対応する Galois 群の元，そして，
自体も１．で得た KK により，
と表せます．
nGG の引数は DA のみですが，mp が生成した RA，KK も上記のように内部で参照し，出力は Galois 群 GG2，また，その元と DA の根との対応のリスト GGRS も生成します．
nGG(DA):=block([h],
n:50,fpprintprec:n,ratepsilon:fpprec:2*n,er:0.1^n,
nRS:map('rhs,bfallroots(DA)),
nRAS:expand(map(lambda([s],subst(A=s,RA)),nRS)),
nRAS1:nRAS[1],
GG2:fullmapl(lambda([s],sublist_indices(nRAS1,lambda([t],abs(s-t)<er))[1]),nRAS),
RS:rat(map(lambda([s],KK.s),fullmapl(lambda([s],RA[s]),map(lambda([s],makelist(s[i],i,length(KK))),GG2)))),
GGRS:map("[",GG2,RS),
GG2)$
実行例．
(%i5) nGG(mp(x^5-2));
Evaluation took 0.2930 seconds (0.3980 elapsed) using 148.313 MB.
(%o5) [[1,2,3,4,5],[4,2,5,1,3],[2,4,1,5,3],[2,1,4,3,5],[3,5,1,4,2],
[5,3,4,1,2],[5,4,3,2,1],[3,1,5,2,4],[1,3,2,5,4],[4,5,2,3,1],
[2,5,3,1,4],[2,3,5,4,1],[1,4,5,3,2],[4,1,3,5,2],[3,4,2,1,5],
[5,1,2,4,3],[1,5,4,2,3],[4,3,1,2,5],[3,2,4,5,1],[5,2,1,3,4]]
(%i6) GGRS;
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o6) [[[1,2,3,4,5],A],[[4,2,5,1,3],(A^16-5*A^11+2350*A^6-3500*A)/11000],
[[2,4,1,5,3],-(A^16-5*A^11+2350*A^6-3500*A)/11000],[[2,1,4,3,5],-A],
[[3,5,1,4,2],-(A^16-10*A^11+2900*A^6-18000*A)/22000],
[[5,3,4,1,2],(A^16-10*A^11+2900*A^6-18000*A)/22000],
[[5,4,3,2,1],(A^16+10*A^11+2900*A^6+18000*A)/22000],
[[3,1,5,2,4],(A^16+5*A^11+2350*A^6+3500*A)/11000],
[[1,3,2,5,4],-(A^16+5*A^11+2350*A^6+3500*A)/11000],
[[4,5,2,3,1],-(A^16+10*A^11+2900*A^6+18000*A)/22000],
[[2,5,3,1,4],-(3*A^16+7600*A^6+11000*A)/22000],
[[2,3,5,4,1],-(A^16+5*A^11+2350*A^6+14500*A)/11000],
[[1,4,5,3,2],-(A^16-5*A^11+2350*A^6-14500*A)/11000],
[[4,1,3,5,2],(A^16-5*A^11+2350*A^6-14500*A)/11000],
[[3,4,2,1,5],(A^11+1800*A)/1100],
[[5,1,2,4,3],(3*A^16+7600*A^6-11000*A)/22000],
[[1,5,4,2,3],-(3*A^16+7600*A^6-11000*A)/22000],
[[4,3,1,2,5],-(A^11+1800*A)/1100],
[[3,2,4,5,1],(A^16+5*A^11+2350*A^6+14500*A)/11000],
[[5,2,1,3,4],(3*A^16+7600*A^6+11000*A)/22000]]
１．分解体の絶対定義多項式
今回のオリジナルプログラム https://github.com/YasuakiHonda/GaloisGroupSolver が行う主な処理は
１．分解体の絶対定義式
２．Galois 群
３．組成列
４．添加する元の定義式
５．分解体の相対定義式
の算出です．種々の方法の比較の意味も含め，私家版を作成しましたので，順次公開させて頂きます．
まず，オリジナルプログラムの１．では，入力の多項式 p の次数 が低い場合を想定して私が気軽に提案した（数値根と
次式の因数分解を用いる）方法を採用して頂いたのですが，これは明らかにハイコストであり，例えば，
の処理に数分を要します．他にも，primitive element と p の根との関係の調整はユーザーが行わねばなりません．
これに対し，下記の mp は有理数体 に p の根を１つずつ添加し，その都度，拡大体の
上の定義式（絶対定義式）の次数を上げてゆく（代数体上の
次以上の既約因子から 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基底の利用（その１）
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.
-----------------------------------------------------------------------------
