多項式の Galois 群計算 - ATPとCASのこと

本プログラムの Galois 群計算の仕組みは http://ehito.hatenablog.com/entry/2019/02/24/200919 で述べましたが，今回は具体例を用いて解説します．この方法は，単に入力多項式の Galois 群 $G$ を特定するのではなく，本プログラムの目標である入力多項式の根（実質的には分解体の primitive element ）の冪根表示に必要な情報を取得する中で，根の置換と同型写像の像を構成するものであり，他の目的での利用価値は不明ながら，同じ結果を得るにあたり代数体上の因数分解を繰り返す方法よりも高速な処理が望めます（以下は基礎体が有理数体 $\mathbb{Q}$ の場合ですが，円分体でも方針は同じです）．

計算は次の５つのステップからなります．
（ステップ１）入力多項式 $f$ の $\mathbb{Q}$ 上の分解体の primitive element $\alpha$ の $\mathbb{Q}$ 上の最小多項式 $g$ を取得．
（ステップ２） $f$ の根を $\alpha$ に対応する変数 $A$ の $\mathbb{Q}$ 上の多項式で表したリスト $L$ を取得．
（ステップ３） $g$ の数値根（多倍長浮動小数点数）を $L$ の $A$ に代入し，根の置換としての $G$ の各元のリスト $G_1$ を取得．
（ステップ４）内積 $L\cdot K$ が $A$ となる整数のリスト $K$ を取得．
（ステップ５） $L\cdot K=A$ の $L$ に（ステップ３）で得た各置換を施し， $G$ の各元による $A$ の像のリスト $G_2$ を取得．

それでは， $f=x^{32}-2$ を例に各ステップを実行してみましょう．

（ステップ１）PARI/GP において

? d=poldegree(ga=subst(g=nfsplitting(f=x^32-2),x,A));g

と入力すると，約 10 秒（Ubuntu 14.04, Core i5-3210M 2.50GHz）で次のように出力されます（環境によっては default コマンドで parisizemax を増やす必要があるかも分かりません）．

%1 = x^256+9617286208*x^224+16939743610198607616*x^192+2420550419437751206410678272*x^160+5415888922248827561943116410019840*x^128+21997636389176028314448861495081828352*x^96+6678882252223940045087005013533786112*x^64+370736334465544641104295493632*x^32+16777216

なお，Chebotarev's density theorem ( https://en.wikipedia.org/wiki/Chebotarev%27s_density_theorem ) を用いて

? N=0;P=10^4;for(n=1,P,if(#factor(Mod(f,prime(n)))[,2]==poldegree(f),N=N+1));
time = 1,929 ms.
? nfsplitting(f,round(P/N));
time = 1,049 ms.

のように高速化できる場合もあります．
（ステップ２）次に

? L=nfisincl(f,ga);

と入力すると，やはり約 10 秒で結果が得られ，次数のみを拾うと

? apply(poldegree,L)
%3 = [241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241,241]

となっています．もし， $L$ を因数分解で得ようと

? nffactor(ga,f)

とすれば．．．time = 11min, 6,323 ms. を要します．
（ステップ３）続いて適当な精度を設定し， $g$ の数値根を求めます（精度不足ならば後の Map 生成時に一対一でない旨が指摘されます）．

? default(realprecision,max(200,floor(1.5*d)));ng=polroots(g);

それらを $L$ の $A$ に代入すると， $G$ の元の準同型性から，置換後の $f$ の数値根が得られるので，誤差を含む部分をカットします（精度不足の場合には，やはり Map 生成時に指摘されます）．

? nL=round(10^5*[subst(L,A,n)|n<-ng]);

この先頭のリストを基準に番号を割り振れば，根の置換群としての Galois 群 $G_1$ が得られます．

? M=Map(Mat([n1=nL[1]~,[1..#n1]~]));G1=apply(L->Vecsmall(apply(s->mapget(M,s),L)),nL);
time = 11 ms.

結果の中身を覗くと

? G1[256]
%10 = Vecsmall([32,3,2,28,29,7,6,24,25,11,10,20,21,15,14,16,17,19,18,12,13,23,22,8,9,27,26,4,5,31,30,1])
? G1[1]
%11 = Vecsmall([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32])

となっています．
（ステップ４） $L\cdot K=A$ を満たす $K$ の例を求めると

? K=matinverseimage(matconcat(vector(d,i,subst(L,A,i))~),[1..d]~)
%13 = [-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]~

となり，確かに

? L*K
%14 = A

となっています．
（ステップ５）（この例では不具合となりませんが， $G$ の位数の大きな問題では $L$ のサイズも大きくなるので）係数を逆に置換し， $A$ の像のリスト $G_2$ を計算します．

? G2=apply(p->L*p~,apply(p->apply(s->K[s],p),apply(p->Vec(Vecsmall(p)^(-1)),G1)));
time = 245 ms.

結果の中身を覗くと

? G2[1]
%16 = A
? G2[256]
%17 = -140056475292934525788582404664588050206670069515460087095141106743/35076811263842343718440191339547566716049104354229383482035976026041335790325202944*A^241-220739581799915126746789240596502331038977678547961/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^225-45346189340689199256448759270211222876248385931559781648808637129058905/1180878375432343917265021254361283554943748463312327749866549152506104759975936*A^209-2122915735604011586670819590599825988863612107748621607679979/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^193-37070637225162918656856440873437887413938164543428519572466260441003191806210881907/548075175997536620600627989680430729938267255534834116906812125406895871723831296*A^177-233704495019564084166805443788958628996577558892455075732292184313125/34254698499846038787539249355026920621141703470927132306675757837930991982739456*A^161-662136250029854892894980319171183536258815554834684143009196788040044982406610465898717709/68509396999692077575078498710053841242283406941854264613351515675861983965478912*A^145-8348613864254667748980481589124398754157263605209916044263166480616892223077/8563674624961509696884812338756730155285425867731783076668939459482747995684864*A^129-185188064704155447844565893173075231389171389924625227439295053172443812433789217646889315355429/8563674624961509696884812338756730155285425867731783076668939459482747995684864*A^113-583740749888698264024996895890648074047860624270764984232490179816787253679769389/267614832030047178027650385586147817352669558366618221145904358108835874865152*A^97-94021954878039241476007349158712536023216697732706026258402809031488357341671233165947721189077867/1070459328120188712110601542344591269410678233466472884583617432435343499460608*A^81-1185485609681620689614824041651980277946739037704685809931592407246459157786629070661/133807416015023589013825192793073908676334779183309110572952179054417937432576*A^65-3568346813747491847830869204262761414148872467371423771007232604647490123832556602374866005560817/133807416015023589013825192793073908676334779183309110572952179054417937432576*A^49-22495943806238745534925077213822835558283166332763585419821490732033000276871382989/8362963500938974313364074549567119292270923698956819410809511190901121089536*A^33-24759304369138635226898064742818791415078768965335903993868383116707204581031081330150511/16725927001877948626728149099134238584541847397913638821619022381802242179072*A^17-4141308750802347261186343968762150510786111593016503005587014426837918040267/2090740875234743578341018637391779823067730924739204852702377797725280272384*A

となっています．勿論， $G_2$ を因数分解で得ようと

? nffactor(ga,g)

とすれば．．．time = 1h, 38min, 46,267 ms. を要し，galoisinit や nfgaloisconj でも同様です．

なお，この $f$ の場合，冪根拡大を含む全体の処理時間は約 30 秒，そのうち約 25 秒が上記５つのステップに費やされます．一方，online の Magma V2.24-5 ( http://magma.maths.usyd.edu.au/calc/ ) の SolveByRadicals ( https://magma.maths.usyd.edu.au/magma/handbook/text/401 https://carma.newcastle.edu.au/tdlc/sins/190503_nicole_sutherland.pdf http://www.fachgruppe-computeralgebra.de/data/CA-Rundbrief/car62.pdf ) で

P<x> := PolynomialRing(IntegerRing());
SolveByRadicals(x^32 - 2);

のように処理すると

Runtime error: too many transformation, cannot find a primitive element

となります．