4.添加する元の定義式,および,5.分解体の相対定義式(解説)

現在のオリジナルプログラムでは,すべての Lagrange Resolvent の冪乗,開冪,和を取って,表題の2つを算出しています.

一方,http://ehito.hatenablog.com/entry/20190123/1548255005 で提案した方法は,分解体の(基礎体 K 上の)定義式 P,新たに添加する冪根 ai の(K に分解体の primitive element A を添加した体 K(A) 上の)定義式 Q,そして,既に添加された元の定義式に対する Groebner 基底から,分解体の(拡大体 K(ai) 上の)定義式 R,冪根 ai の(K 上の)定義式 S を同時に得るもの(数学としては,P,Q を最初の2式,R,Sを終端の2式とする代数体上の互除法により,P=Q=0 ⇔ R=S=0 を保ちつつ,primitive element についての次数を下げる(特に S は 0 次)ことに当たります)ですが,maxima の Groebner 基底ルーチンによる実装 RR5 では

(%i2) for i:1 thru 15 do RR5(cs(nGG(mp(PL[i]))))$
PROU_5 is a member. 
PROU_5 is a member. 
PROU_3 is a member. 
PROU_3 is a member. 
Evaluation took 298.6010 seconds (299.8230 elapsed) using 29739.980 MB.

となり,Groebner 基底の計算に Risa/Asir の関数を用いる RR5a では

(%i3) for i:1 thru 15 do RR5a(cs(nGG(mp(PL[i]))))$
PROU_5 is a member. 
PROU_5 is a member. 
PROU_3 is a member. 
PROU_3 is a member. 
Evaluation took 7.3480 seconds (10.5720 elapsed) using 3442.276 MB.

となりました. ただ,今回のお話は maxima プロパーを旨とするもの(?)だと思いますので,前回掲載した RR6 での Groebner 基底の利用は無理のない(と思われる)範囲に留めています.

全体の処理の流れは,オリジナルプログラムと同様に,すべての Lagrange Resolvent の和として,5.を得るものですが,各 Resolvent LR[j] (j=1,...,pp) (pp は拡大次数)に冪根 ai を適当な回数(コードでは pp-1 回)乗じて,冪根の変数の冪乗 a[i]^(j-1) と基礎体上の多項式 ai^(pp-j)*LR[j] との積の形に整理しています.これは https://ikumi.que.jp/blog/wp-content/uploads/2018/09/galois-solution.pdf の 12.(4),https://ikumi.que.jp/blog/archives/717 の【2019, 1/21 追記】で間接,直接に指摘されているもので,コード内の

F:sum(a[i]^(j-1)*ai^(pp-j)*LR[j],j,1,pp)

が相対定義式の 0 でない定数倍 ai^(pp-1)*sum(LR[j],j,1,pp) となる訳ですが,この F の係数を冪根 ai の変数 a[i]と既に添加されているの元の変数で表すと,主係数に変数が残るので,有理化,つまり,代数体での逆元の計算が必要となり,Groebner 基底の出座に至ります(笑).例えば

(%i4) poly_reduced_grobner ([(a^2+b)*x^2+x+1,a^3+a+1,b^2+a*b+a^2],[x,b,a]);
Evaluation took 0.0020 seconds (0.0020 elapsed) using 159.969 KB.
(%o4) [a^3+a+1,b^2+a*b+a^2,(-x^2)+b*x-a^2*x+a*x+b-a^2+a]

といった感じです.

これに先立つもう1つのポイントが ai の検出です.ai は 0 でない Resolvent の主係数であり,いわゆる代数体での 0 判定が必要となり,その際,定義式の既約性が重要(「代入した値が0⇔定義式による整除」が成立)で,例えば,a^2-2 が定める Q(a) に,その上で既約ではない b^2-2 の根 b を添加するといった暴挙に出ると,a+b=0 か否かの判定が出来ません.実際,1の原始 pp 乗根 w[i] の(Q上の)定義式 Dw は,素数 pp が 2 以外ならば高次なので,基礎体上で既約でない可能性があり,因数分解が必要となり,自然なのは逐次拡大体としての各中間体上の分解ですが,処理時間を比較した結果,RR6 では分解体 Q(A) 上の分解 factor(Dw,DA) を利用しています.具体的には,低次の因数に分解されるならば,pp が素数であることから,それらの次数は全て等しく,1次ならば,w[i] は直接 A で表示し,次数が 2 以上ならば,0 でない ai の検出までは Q(A) 係数のまま用い(コード内の /* w[i] with A */ からの部分),ai^pp の簡約からは,分解体の(基礎体上の)定義式と既に添加された元の定義式により基礎体係数に直した形を用いています(コード内の /* w[i] with a[1],w[1],...,a[i-1],w[i-1] */ からの部分).参考として,2次式に分解される例を記しておきます.

(%i6) [p,Dw,DA];
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%o6) [x^5-5*x^3+5*x+5,w3^4+w3^3+w3^2+w3+1,
        A^20-50*A^18+1025*A^16-11250*A^14+72500*A^12-268125*A^10+456875*A^8
            +137500*A^6-1578125*A^4+1640625*A^2+1378125]
(%i7) Dwf:factor2list(factor(Dw,DA))[1];
Evaluation took 0.5380 seconds (0.5410 elapsed) using 91.638 MB.
(%o7) 240295101375*w3^2-54652*A^18*w3+2992335*A^16*w3-67936350*A^14*w3
                        +830924250*A^12*w3-5985718405*A^10*w3
                        +25902278250*A^8*w3-65067996000*A^6*w3
                        +83309465625*A^4*w3-24030550000*A^2*w3-202585732125*w3
                        +240295101375
(%i8) Dwf:poly_normalize(rrem(Dwf,[[(5*a2)/2-(A*(25*a1-75))/2+(A^3*(5*a1-25))/2+A^5,A],[a2^2-462*a1+1050,a2],[a1^2-5,a1]]),[w[3]]);
Evaluation took 0.0050 seconds (0.0050 elapsed) using 862.938 KB.
(%o8) w3^2-((a1-1)*w3)/2+1

他にも,RR シリーズでは
・remainder と algebraic モード との使い分け(後者の方が高速だが,tellrat は monic しか受け付けず,その為の変数変換はやや大掛かりになるので)
・元の変数に atom を使用(これも後者の方が高速なので)
・poly_reduced_grobner ではなく poly_grobner を利用(余分な式が含まれるが,やはり後者の方が高速なので)
といった工夫を施しています.

4.添加する元の定義式,および,5.分解体の相対定義式

取り敢えず,コードと例を記します.解説はまた後日.

load("grobner")$
factor2list(f0):=block([f:num(f0)],if op(f)="-" then args(-f) elseif op(f)="*" then args(f) else [f])$
rrem(S,T):=block([U:S],map(lambda([s],U:remainder(U,s[1],s[2])),T),U)$
RR6(csGG):=block([],[cs1,cs2,cs3]:csGG,kill(w),DA0:DA,Dws:[[DA,A]],DP:aiA:[],
for i:1 thru length(cs3) do (push([DA0,A],DP),
a[i]:concat(a,i),w[i]:concat(w,i),
/* w[i] with A */
Dw:rat((w[i]^(pp:cs3[i])-1)/(w[i]-1)),
Dwf:factor2list(factor(Dw,DA))[1],
if (degDwf:hipow(Dwf,w[i]))=1 then (w[i]:rhs(solve(Dwf,w[i])[1]),
  if pp>2 then printn(concat("PROU_",pp," is a member.")))
else (push([Dwf,w[i]],Dws),push([Dwf,w[i]],DP)),
/* non-zero aix and ai */
tellrat(DA),algebraic:on,
T:cs1[i+1],for g0 in cs1[i] while member(g:g0,T) do 0,
LRx:makelist(apply("*",makelist(x-assoc(t,GGRS),t,T:map(lambda([u],mul(g,u)),T))),j,1,pp),
printn(LRx),
algebraic:off,untellrat(A),
ai:0,for cnt:1 while ai=0 do
(CNT:cnt,aix:rrem(sum((w[i]^mod(cnt*j,pp)*LRx[j]),j,1,pp),append(Dws,DP)),
 ai:if aix=0 then 0 else coeff(aix,x,hipow(aix,x)),
 printn([[cnt],ai])),
/* w[i] with a[1],w[1],...,a[i-1],w[i-1] */
if degDwf=1 then w[i]:rrem(w[i],DP)
else (pop(DP),push([Dwf:poly_normalize(rrem(Dwf,DP),[w[i]]),w[i]],DP)),
/* Lagrange Resolvent */
LR[1]:aix,for cnt:2 thru pp do 
(LR[cnt]:rrem(sum(w[i]^mod(CNT*cnt*j,pp)*LRx[j],j,1,pp),DP),
 printn(concat("LR_",cnt)),printn(LR[cnt])),
F:rrem(sum(a[i]^(j-1)*ai^(pp-j)*LR[j],j,1,pp),DP),
F:subst(x=A,F),printn("F:"),printn(F),
/* def. of a[i] */
aipp:rrem(ai^pp,DP)-a[i]^pp,
/* reduction for F */
DP:delete([DA0,A],DP),fDP:map(first,DP),sDP:map(second,DP),
GR:poly_grobner(append([F,aipp],fDP),append([A,a[i]],sDP)),
printn("GR:"),printn(GR),
/* monic rather than primitive */
a1n:num(rat(aipp)),
 na1n:ifactors(abs(poly_content(subst(a[i]=0,a1n),sDP))),
 na1n:apply("*",map(lambda([s],s[1]^floor(s[2]/cs3[i])),na1n)),
 da1n:ifactors(abs(coeff(a1n,a[i],cs3[i]))),
 da1n:apply("*",map(lambda([s],s[1]^ceiling(s[2]/cs3[i])),da1n)),
printn(na1n/da1n),
a1n:poly_normalize(subst(a[i]=na1n/da1n*a[i],a1n),[a[i]]),
push([a1n,a[i]],DP),
push([a[i]-da1n/na1n*ai,a[i]],aiA),
DA0:sublist(GR,lambda([s],numberp(coeff(s,A,deg:hipow(s,A))) and deg<hipow(DA0,A)))[1],
DA0:poly_normalize(subst(a[i]=na1n/da1n*a[i],DA0),[A]),
printn([3,DA0,DP])
),
rr2:[solve(DA0,A)[1],DP],
append([[DA0,A]],DP))$

実行例.p の各根のリスト RA は長いのでプリントしていません.

(%i24) for i:1 thru 15 do print([p:PL[i],mp(p),RR6(cs(nGG(DA)))])$

[x^2-2,A^2-8,[[2*a1+A,A],[a1^2-2,a1]]] 
[x^3-3*x-1,A^3-3*A-1,[[a1^2*w1+a1+A,A],[w1+a1^3+1,a1],[w1^2+w1+1,w1]]] 
[x^4-2,A^8+28*A^4+2500,[[a3+2*a2+A,A],[a3^2-a1,a3],[a2^2+a1,a2],[a1^2-2,a1]]] 
[x^4+x^2-1,A^8+10*A^6+47*A^4+110*A^2+841,
 [[(a3+2*a2)/2+A,A],[a3^2-2*a1+2,a3],[a2^2+2*a1+2,a2],[a1^2-5,a1]]]
  
[x^4-2*x^3+2*x^2+2,A^12+4*A^10+24*A^8+48*A^6-560*A^4+3136,
 [[a3/21+A,A],[(-126*a1^2*w1)+a3^2+42*a2-84*a1^2+294*a1+294,a3],
  [(42*a1^2+588*a1)*w1+a2^2-168*a1^2+294*a1+1323,a2],[(-21*w1)+a1^3-7,a1],
  [w1^2+w1+1,w1]]]
  
[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,
 [[(a4+585)/585+A,A],
  [a1*(570*a2^2-2100*a2^2*w2)+1620*a2^2*w2+a4^2+2340*a3+2385*a2^2+114075*a2
                             +1711125,a4],
  [a1*((30*a2^2-29250*a2)*w2-660*a2^2-40950*a2)
    +2700*a2^2*w2+a3^2+1440*a2^2+1368900,a3],
  [4860*w2+a1*((-6300*w2)-8010)+a2^3-2295,a2],[w2^2+w2+1,w2],[a1^2-3,a1]]]
  
[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,
 [[a4/624+A,A],
  [a1*(38*a2^2-140*a2^2*w2)-648*a2^2*w2+a4^2+9984*a3-954*a2^2+64896*a2,a4],
  [a1*((-65*a2*w2)-91*a2)+(60*a2^2+351*a2)*w2+a3^2+32*a2^2-117*a2+32448,a3],
  [(-3888*w2)+a1*((-840*w2)-1068)+a2^3+1836,a2],[w2^2+w2+1,w2],[a1^2-229,a1]]]
  
[x^5-2,A^20+2500*A^10+50000,
 [[a3+A,A],[a3^5-5*a2,a3],[a2^2+22*a1+50,a2],[a1^2-5,a1]]]
  
[x^5-5*x+12,A^10-10*A^8-75*A^6+1500*A^4-5500*A^2+16000,
 [[A-(a1*((3*a2^4-5*a2^3+25*a2^2)*w2^3+(3*a2^4+15*a2^3+25*a2^2)*w2^2
                                      +10*a2^3*w2+9*a2^4+5*a2^3-50*a2^2)
     +((-20*a2^4)-25*a2^3)*w2^3+(10*a2^4-25*a2^3-250*a2^2)*w2^2
     +((-10*a2^4)-250*a2^2)*w2-5*a2^4-25*a2^3-125*a2^2-625*a2)
     /3125,A],
  [a1*(20625*w2^3+20625*w2^2+33750)+90625*w2^3-21875*w2^2+68750*w2+a2^5+34375,
   a2],[w2^4+w2^3+w2^2+w2+1,w2],[a1^2+10,a1]]]
  
[x^5+20*x+32,A^10-20*A^8+100*A^6+2000*A^4-32000*A^2+128000,
 [[(a1*((6*a2^4-10*a2^3+100*a2^2)*w2^3+(6*a2^4-70*a2^3+100*a2^2)*w2^2
                                      -80*a2^3*w2-7*a2^4-40*a2^3+50*a2^2)
    +(10*a2^4+100*a2^3)*w2^3+(20*a2^4+100*a2^3+500*a2^2)*w2^2
    +(30*a2^4+500*a2^2)*w2+15*a2^4-50*a2^3+250*a2^2+2500*a2)
    /12500
    +A,A],
  [a1*(22500*w2^3+22500*w2^2+42500)+87500*w2^3-12500*w2^2+75000*w2+a2^5+37500,
   a2],[w2^4+w2^3+w2^2+w2+1,w2],[a1^2+5,a1]]]
  
[x^5+11*x+44,A^10-22*A^8+77*A^6+4356*A^4-53724*A^2+189728,
 [[(a1*((23*a2^4-55*a2^3+605*a2^2)*w2^3+(23*a2^4-935*a2^3+605*a2^2)*w2^2
                                       -990*a2^3*w2-106*a2^4-495*a2^3
                                       +1815*a2^2)
    +(100*a2^4+715*a2^3)*w2^3+(50*a2^4+715*a2^3+6050*a2^2)*w2^2
    +(150*a2^4+6050*a2^2)*w2+75*a2^4-165*a2^3+3025*a2^2+33275*a2)
    /166375
    +A,A],
  [a1*(74800*w2^3+74800*w2^2+126775)+158125*w2^3-34375*w2^2+123750*w2+a2^5
                                    +61875,a2],[w2^4+w2^3+w2^2+w2+1,w2],
  [a1^2+2,a1]]]
  
[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,
 [[((10*a1^4+99*a1^3+121*a1^2)*w1^3+((-15*a1^4)+33*a1^3+242*a1^2)*w1^2
                                   +(20*a1^4+77*a1^3-242*a1^2)*w1+26*a1^4
                                   -33*a1^3+1331*a1+1331)
    /6655
    +A,A],[385*w1^3+110*w1^2+220*w1+a1^5-66,a1],[w1^4+w1^3+w1^2+w1+1,w1]]]
  
[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,
 [[(a1*(a2*(a3^4+2*a3^2)-6*a3^3)+a2*(a3^4+10*a3^2)+10*a3^3+40*a3)/40+A,A],
  [a3^5-50*a1*a2-110*a2,a3],[a2^2-22*a1+50,a2],[a1^2-5,a1]]]
  
[x^6+x^3+1,A^6+A^3+1,[[a2/2+A,A],[a2^3-4*a1-4,a2],[a1^2+3,a1]]] 
[x^6-2,A^12-1012*A^6+19307236,
 [[a3+A,A],[a3^3-18*a2-35*a1,a3],[a2^2+6,a2],[a1^2-2,a1]]]
  
Evaluation took 9.3890 seconds (9.4080 elapsed) using 4571.050 MB.

3.組成列

組成列の計算については,http://ehito.hatenablog.com/entry/2019/02/13/200937 で述べましたが,mnsg に積閉の必要条件などを加え,やや速くなっています.cs の引数は nGG の出力,出力は先述の通りです.

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])))$
isG(GG):=is(setify(apply('append,makelist(makelist(mul(i,j),i,GG),j,GG)))=setify(GG))$
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))$
mnsg(GG):=block([S:0,lenGG:length(GG)],
C0:sort(CC(GG)),C01:C0[1],
PS:map(lambda([s],apply('append,listify(s))),listify(powerset(setify(rest(C0))))),
PS:map(lambda([s],append(C0[1],s)),PS), printn(length(PS)),
for d in rest(reverse(listify(divisors(lenGG)))) while S=0 do (
PS:sublist(PS,lambda([s],length(s)<=d and member(lst2:mul(lst:last(s),lst),s) and member(mul(lst2,lst),s))), printn(length(PS)),
for g in PS unless length(S:g)=d and isG(g) do S:0),S)$
cs(GG):=block([S:[GG],m:GG],unless length(m)=1 do push(m:mnsg(m),S),
csGG:[S:reverse(S),S:map(length,S),makelist(S[i-1]/S[i],i,2,length(S))])$

実行例.

(%i7) cs(nGG(mp(x^12-3)));
Evaluation took 3.5550 seconds (3.5900 elapsed) using 679.381 MB.
(%o7) [[[[1,2,3,4,5,6,7,8,9,10,11,12],[12,9,8,7,5,6,4,3,2,11,10,1],
          [10,9,7,8,3,4,6,5,11,2,12,1],[11,2,4,3,8,7,6,5,10,9,1,12],
          [7,5,11,10,12,2,9,1,8,3,4,6],[4,5,10,11,1,9,2,12,3,8,7,6],
          [3,8,9,1,11,10,2,12,4,5,6,7],[8,3,2,12,10,11,9,1,7,5,6,4],
          [3,6,11,10,9,1,12,2,4,7,8,5],[8,6,10,11,2,12,1,9,7,4,3,5],
          [7,4,12,2,11,10,1,9,8,6,5,3],[4,7,1,9,10,11,12,2,3,6,5,8],
          [2,1,7,8,6,5,3,4,12,10,11,9],[9,12,4,3,6,5,8,7,1,11,10,2],
          [10,12,3,4,7,8,5,6,11,1,9,2],[11,1,8,7,4,3,5,6,10,12,2,9],
          [6,3,12,2,1,9,11,10,5,7,8,4],[6,8,1,9,12,2,10,11,5,4,3,7],
          [1,11,5,6,3,4,8,7,9,12,2,10],[12,10,5,6,8,7,3,4,2,1,9,11],
          [2,11,6,5,7,8,4,3,12,9,1,10],[9,10,6,5,4,3,7,8,1,2,12,11],
          [5,7,9,1,2,12,11,10,6,3,4,8],[5,4,2,12,9,1,10,11,6,8,7,3]],
         [[1,2,3,4,5,6,7,8,9,10,11,12],[1,11,5,6,3,4,8,7,9,12,2,10],
          [2,1,7,8,6,5,3,4,12,10,11,9],[9,10,6,5,4,3,7,8,1,2,12,11],
          [10,12,3,4,7,8,5,6,11,1,9,2],[11,2,4,3,8,7,6,5,10,9,1,12],
          [12,9,8,7,5,6,4,3,2,11,10,1],[2,11,6,5,7,8,4,3,12,9,1,10],
          [11,1,8,7,4,3,5,6,10,12,2,9],[9,12,4,3,6,5,8,7,1,11,10,2],
          [10,9,7,8,3,4,6,5,11,2,12,1],[12,10,5,6,8,7,3,4,2,1,9,11]],
         [[1,2,3,4,5,6,7,8,9,10,11,12],[1,11,5,6,3,4,8,7,9,12,2,10],
          [2,1,7,8,6,5,3,4,12,10,11,9],[11,2,4,3,8,7,6,5,10,9,1,12],
          [2,11,6,5,7,8,4,3,12,9,1,10],[11,1,8,7,4,3,5,6,10,12,2,9]],
         [[1,2,3,4,5,6,7,8,9,10,11,12],[2,11,6,5,7,8,4,3,12,9,1,10],
          [11,1,8,7,4,3,5,6,10,12,2,9]],[[1,2,3,4,5,6,7,8,9,10,11,12]]],
        [24,12,6,3,1],[2,2,2,3]]

2.Galois 群

現在のオリジナルプログラムでは,Galois 群の元と分解体の絶対定義式の根との対応の把握に若干の検索を行っています.

これに対して,下記の nGG では,1.で得た RA における A が絶対定義式 DA の任意の根であるという点に着目して,それらの根と「入力 p の根の置換」との対応を直接求めています.すなわち,1.で得た DA の数値根 B_k\ (k=1,\ldots,m) (絶対誤差は p の根差の半分より小とする)を RA=[s_1(A),\ldots,s_n(A)] の A に代入したものは [s_1(B_1),\ldots,s_n(B_1)] の順列  [s_{{\sigma_k}(1)}(B_1) ,\ldots, s_{{\sigma_k}(n)}(B_1)] であり,この \sigma_k が絶対定義式の根 A_k に対応する Galois 群の元,そして,A_k 自体も1.で得た KK により,A_k=\sum_{i=1}^{|KK|} KK[i] s_{{\sigma_k}(i)}(A) と表せます.

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]]

1.分解体の絶対定義多項式

今回のオリジナルプログラム https://github.com/YasuakiHonda/GaloisGroupSolver が行う主な処理は
1.分解体の絶対定義式
2.Galois 群
3.組成列
4.添加する元の定義式
5.分解体の相対定義式
の算出です.種々の方法の比較の意味も含め,私家版を作成しましたので,順次公開させて頂きます.

まず,オリジナルプログラムの1.では,入力の多項式 p の次数 n が低い場合を想定して私が気軽に提案した(数値根と n! 次式の因数分解を用いる)方法を採用して頂いたのですが,これは明らかにハイコストであり,例えば,x^5-2 の処理に数分を要します.他にも,primitive element と p の根との関係の調整はユーザーが行わねばなりません.

これに対し,下記の mp は有理数\mathbb{Q} に p の根を1つずつ添加し,その都度,拡大体の \mathbb{Q} 上の定義式(絶対定義式)の次数を上げてゆく(代数体上の 2 次以上の既約因子から primitive element を生成する),分解体の定義式の計算としては一般的なものです.これによれば,primitive element と p の根との関係(*)は自動的に得られ,また,n=5 の場合には,絶対定義式の次数が 20 を超えた段階で,非可解と判定できます.

mp の入力 p は,x を変数とする \mathbb{Z} 上の 2 次以上の monic な既約多項式であり,出力は p の分解体(単純拡大)の primitive element A の \mathbb{Z} 上の最小多項式の変数に 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]]

のように指摘されています.

一方,共役類に分割された群  G=\bigcup_{i\in A}C_i とその部分群 S について,次の2つは等価でした.
(1)SG正規部分群である.
(2) S=\bigcup_{i\in B}C_i となる空でない B \subseteq A がある.

従って,有限群 G を共役類に分割すれば,G の位数 d正規部分群の全体は,元の個数の合計が d,かつ,和が積閉である共役類の和の全体に一致します.

以下,この性質による組成列計算のコードを記します.なお,mnsg では powerset を用いましたが,入力の群のサイズが程々(笑)なら共役類も多くはないので,爆発しないと思われます(母関数を用いると,元の個数の合計が d になる組のみを抽出できますが,…).

/* 置換の積,逆元 */
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);