# 代数体上の線型方程式系の求解

c5^4+c5^3+c5^2+c5+1,
e1^2+10,
e2^5+((-3750*c5^3)-3750*c5^2-1875)*e1-3125*c5^3+9375*c5^2+6250*c5+3125

である体 Q(c5,e1,e2) 上で

(((10*a4-2*a3+30*a2-70*a1+360)*e1-10*a4-20*a3-150*a2-100*a1-600)*c5^3
+((10*a4-14*a3+30*a2-90*a1+360)*e1-30*a4-20*a3-50*a2-100*a1-1000)
*c5^2+(((-16*a3)-160*a1)*e1-40*a4-200*a2-1600)*c5
+(5*a4-8*a3+15*a2-80*a1+180)*e1-240*e2-20*a4-10*a3-100*a2+250*a1-800)
/240,
(((2*a4-20*a3-10*a2-140*a1+40)*e1-20*a4-60*a3-100*a2+100*a1-800)*c5^3
+((4*a4-20*a3-20*a2-140*a1+80)*e1-20*a4-100*a2+200*a1-800)*c5^2
+((6*a4-30*a2+120)*e1-60*a3+300*a1)*c5
+(3*a4-40*a3-15*a2-120*a1+60)*e1-48*e2^2-60*a4-30*a3-100*a2+150*a1
-2000)
/48,
-(((125*a4+10*a3+575*a2+750*a1+4900)*e1
+150*a4-450*a3+250*a2-2750*a1+5000)
*c5^3
+((125*a4-30*a3+575*a2-250*a1+4900)*e1
-250*a4-450*a3-1750*a2-2750*a1-11000)
*c5^2+(((-20*a3)+500*a1)*e1-100*a4-1500*a2-6000)*c5
+(100*a4-10*a3+100*a2+250*a1+3200)*e1+48*e2^3-50*a4-600*a3-750*a2
-4000*a1-3000)
/48,
-(((450*a4+1000*a3+3750*a2+3000*a1+21000)*e1
-1500*a3-5000*a2-7500*a1-10000)
*c5^3
+((400*a4+1000*a3+3000*a1+12000)*e1-5000*a3-5000*a2-10000*a1-10000)
*c5^2+((850*a4+3750*a2+33000)*e1-6500*a3-17500*a1)*c5
+(425*a4-1000*a3+1875*a2-1000*a1+16500)*e1+48*e2^4-2500*a4-3250*a3
-12500*a2-8750*a1-100000)
/48

の零点 [a4,a3,a2,a1] を（少なくとも一つ）求めるといった処理であり，ツールとしては，Nemo/Julia シリーズ AbstractAlgebra.jl の solve_left ( https://nemocas.github.io/AbstractAlgebra.jl/latest/matrix/#AbstractAlgebra.Generic.solve_left-Union{Tuple{T},%20Tuple{MatElem{T},MatElem{T}}}%20where%20T%3C:RingElem )

Welcome to Nemo version 0.14.1

Nemo comes with absolutely no warranty whatsoever

Welcome to

_    _           _
| |  | |         | |
| |__| | ___  ___| | _____
|  __  |/ _ \/ __| |/ / _ \
| |  | |  __/ (__|   <  __/
|_|  |_|\___|\___|_|\_\___|

Version 0.6.2 ...
... which comes with absolutely no warranty whatsoever
(c) 2015-2019 by Claus Fieker, Tommy Hofmann and Carlo Sircana

julia> MXl2=MXl;

julia> MX=matrix(MXl,MXl2,MXl2,MXl);

julia> cx=matrix(MXl,1,MXl2,MXl);

julia> @time xx=solve_left(MX,cx);
1.407416 seconds (611.78 k allocations: 29.782 MiB, 1.75% gc time)

および，maxima の algebraic モードでの linsolve と fast_linsolve

(%i11) (map(lambda([s],tellrat(s)),DP),algebraic:on,linsolve(sAL1,AL0))$Evaluation took 0.0040 seconds (0.0040 elapsed) using 1.843 MB. (%i12) (map(lambda([s],tellrat(s)),DP),algebraic:on,fast_linsolve(sAL1,AL0))$
Assuming entries of type ANY-MACSYMA
Starting to solve.  There are 4 equations with 4 unknowns occurring.
The value of (SP-TYPE-OF-ENTRIES SP-MAT) is ANY-MACSYMA
The dimension of the solution space is 0
Evaluation took 0.0040 seconds (0.0040 elapsed) using 1.186 MB.

があります．

これらのうち solve_left は最近公開されたもので，高速 arithmetic で知られる Nemo シリーズゆえ相応のパフォーマンスが期待され，従来の hnf_with_transform と can_solve_left_reduced_triu の組み合わせよりは高速化していますが，サイズの大きな問題に対する処理速度は maxima の fast_linsolve が優位であり，例えば，47 分多項式を入力とする処理に現れる，文字列としてのサイズが 405801 バイト（上記の例は 1087 バイト）の方程式系の場合，

julia> @time xx=solve_left(MX,cx);
80.273864 seconds (65.62 M allocations: 31.598 GiB, 1.47% gc time)

に対し，無印の linsolve では

(%i14) (map(lambda([s],tellrat(s)),DP),algebraic:on,linsolve(sAL1,AL0))$Evaluation took 458.1150 seconds (472.1600 elapsed) using 92401.623 MB. そして， (%i15) (map(lambda([s],tellrat(s)),DP),algebraic:on,fast_linsolve(sAL1,AL0))$
Assuming entries of type ANY-MACSYMA
Starting to solve.  There are 22 equations with 22 unknowns occurring.
The value of (SP-TYPE-OF-ENTRIES SP-MAT) is ANY-MACSYMA
The dimension of the solution space is 0
Evaluation took 54.8670 seconds (56.4100 elapsed) using 20716.948 MB.

となります．

# 位数の取得

? for(i=1,30,nfsplitting2(x^i-2));
[0,1,1]
[0.03125000000,2,1]
[0.004801097394,6,1]
[0.001226218787,8,1]
[0.0008579881657,20,1]
[0.002494223903,12,1]
[0.0007927650227,42,1]
[0.002837500000,16,1]
[0.0004409526648,54,1]
[0.003434499314,40,1]
[7.133182736E-5,110,1]
[0.001070169231,48,1]
[0.0004472513280,156,1]
[0.002071109694,84,1]
[0.0006910009183,120,1]
[0.001790521626,64,1]
[0.001964397853,272,1]
[0.001106438745,108,1]
[0.0002551165675,342,1]
[0.001759224398,160,1]
[0.001657411051,252,1]
[0.0006244174148,220,1]
[0.001948888552,506,1]
[0.0008731821745,500,1]
[0.0008779652962,312,1]
[0.0003889834061,486,1]
[0.002439757082,336,1]
[0.0008391018743,812,1]
*** nfsplitting: Warning: increasing stack size to 16000000.
time = 10,490 ms.

? for(i=1,30,nfsplitting(x^i-2));
*** nfsplitting: Warning: increasing stack size to 16000000.
*** nfsplitting: Warning: increasing stack size to 32000000.
*** nfsplitting: Warning: increasing stack size to 64000000.
time = 3min, 24,135 ms.

GAP の組み込み関数 ProbabilityShapes は同じ原理で作動し，例えば

gap> ProbabilityShapes(x^16-x^15-x^14+4*x^13-3*x^12-x^11+5*x^10-7*x^9+2*x^8+4*x^7-5*x^6+5*x^5-x^4-2*x^3+2*x^2-2*x+1);time;
[ 1947 ]
111832
gap> TransitiveGroup(16,1947);
t16n1947
gap> Size(last);
7962624

のように位数を得ることも可能ですが，一般に

? galoisord(x^16-x^15-x^14+4*x^13-3*x^12-x^11+5*x^10-7*x^9+2*x^8+4*x^7-5*x^6+5*x^5-x^4-2*x^3+2*x^2-2*x+1)
time = 90 ms.
%2 = [0.006932446267,7962624,1]

のような高速処理は望めません．

# 多項式の Galois 群計算

（ステップ１）入力多項式  上の分解体の primitive element  上の最小多項式 を取得．
（ステップ２） の根を に対応する変数  上の多項式で表したリスト を取得．
（ステップ３） の数値根（多倍長浮動小数点数）を  に代入し，根の置換としての の各元のリスト を取得．
（ステップ４）内積  となる整数のリスト を取得．
（ステップ５）  に（ステップ３）で得た各置換を施し， の各元による の像のリスト を取得．

それでは， を例に各ステップを実行してみましょう．

（ステップ１）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]

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

? nffactor(ga,f)

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

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

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

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

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

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

? G1
%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
%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])

となっています．
（ステップ４） を満たす の例を求めると

? 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

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

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

? G2
%16 = A
? G2
%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

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

? nffactor(ga,g)

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

なお，この の場合，冪根拡大を含む全体の処理時間は約 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

となります．

# 円分体上の Galois 群

・サンプル番号
・入力多項式 上の分解体の primitive element の 上の最小多項式の次数
・同じく円分体上の最小多項式の（主係数についての）次数
・各中間体の拡大次数のリスト
・組成列までの処理時間（秒）
・冪根拡大の処理時間（秒）
のリストです．

[1,x^2-2,2,2,,0.152,0.002]
[2,x^3-3*x-1,3,3,,0.137,0.003]
[3,x^4-2,8,8,[2,2,2],0.184,0.006]
[4,x^4+x^2-1,8,8,[2,2,2],0.163,0.007]
[5,x^4-2*x^3+2*x^2+2,12,12,[3,2,2],0.181,0.03]
[6,x^4+2*x^3+3*x^2+4*x+5,24,24,[2,3,2,2],0.275,0.125]
[7,x^4+x+1,24,24,[2,3,2,2],0.234,0.089]
[8,x^5-2,20,5,,0.209,0.002]
[9,x^5-5*x+12,10,10,[2,5],0.205,0.022]
[10,x^5+20*x+32,10,10,[2,5],0.192,0.027]
[11,x^5+11*x+44,10,10,[2,5],0.166,0.026]
[12,x^5+x^4-4*x^3-3*x^2+3*x+1,5,5,,0.172,0.009]
[13,x^5+100*x^2+1000,20,5,,0.21,0.012]
[14,x^6+x^3+1,6,3,,0.144,0.003]
[15,x^6-2,12,6,[2,3],0.184,0.004]
[16,x^5-5*x^3+5*x-5,20,10,[2,5],0.19,0.022]
[17,x^6+x^3+7,18,9,[3,3],0.166,0.011]
[18,x^6-3*x^4+1,24,24,[2,3,2,2],0.266,0.061]
[19,x^6+x^4-9,24,24,[2,3,2,2],0.219,0.083]
[20,x^6+x^4-8,48,48,[2,2,3,2,2],0.611,0.703]
[21,x^7+7*x^3+7*x^2+7*x-1,14,7,,0.184,0.039]
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,,0.19,0.033]
[23,x^7-2,42,7,,0.21,0.004]
[24,x^8-2*x^6-x^4+7*x^2-5*x+1,16,16,[2,2,2,2],0.28,0.034]
[25,x^9+x^8+3*x^6+3*x^3-3*x^2+5*x-1,18,18,[2,3,3],0.253,0.049]
[26,x^12-3,24,12,[2,2,3],0.516,0.009]
[27,x^7-14*x^5+56*x^3-56*x+22,21,7,,0.192,0.029]
[28,
x^15-x^14-14*x^13+13*x^12+78*x^11-66*x^10-220*x^9+165*x^8+330*x^7-210*x^6
-252*x^5+126*x^4+84*x^3-28*x^2-8*x+1,15,15,[3,5],0.606,0.182]

<permutation group with 72 generators>
Welcome to Nemo version 0.13.5
0.634780 seconds (798.89 k allocations: 42.184 MiB, 1.88% gc time)
3.561398 seconds (3.36 M allocations: 1014.657 MiB, 2.55% gc time)
[29,x^6+x^4-x^2+5*x-5,72,72,[2,2,2,3,3],4.496,15.317]
*** nfisincl: Warning: increasing stack size to 16000000.
*** write: Warning: increasing stack size to 16000000.
<permutation group with 13 generators>
Welcome to Nemo version 0.13.5
0.729594 seconds (834.13 k allocations: 43.711 MiB, 8.42% gc time)
0.795107 seconds (1.01 M allocations: 50.953 MiB, 1.18% gc time)
[30,2*(x+1)^13-x^13,156,13,,12.036,12.616]
*** nfisincl: Warning: increasing stack size to 16000000.
<permutation group with 49 generators>
Welcome to Nemo version 0.13.5
0.632685 seconds (811.09 k allocations: 42.324 MiB, 9.44% gc time)
0.994494 seconds (1.16 M allocations: 75.809 MiB, 1.63% gc time)
Welcome to Nemo version 0.13.5
0.639144 seconds (811.09 k allocations: 42.326 MiB, 9.76% gc time)
0.808816 seconds (1.01 M allocations: 51.011 MiB, 1.18% gc time)
[31,
x^14+28*x^11+28*x^10-28*x^9+140*x^8+360*x^7+147*x^6+196*x^5+336*x^4-546*x^3
-532*x^2+896*x+823,98,49,[7,7],5.772,33.549]

*** nfisincl: Warning: increasing stack size to 16000000.
*** write: Warning: increasing stack size to 16000000.
<permutation group with 75 generators>
Welcome to Nemo version 0.13.5
0.635679 seconds (802.47 k allocations: 41.997 MiB, 1.84% gc time)
3.301811 seconds (2.28 M allocations: 812.627 MiB, 2.59% gc time)
Welcome to Nemo version 0.13.5
0.688764 seconds (802.47 k allocations: 41.997 MiB, 8.95% gc time)
1.009562 seconds (1.17 M allocations: 82.609 MiB, 1.53% gc time)
Welcome to Nemo version 0.13.5
0.690487 seconds (802.47 k allocations: 41.997 MiB, 8.90% gc time)
0.809244 seconds (1.01 M allocations: 50.949 MiB, 1.17% gc time)
[32,
x^15-470*x^13-305*x^12+71840*x^11+85357*x^10-4292700*x^9-3714805*x^8
+119761820*x^7+25284495*x^6-1542190154*x^5+717324725*x^4+7178878600*x^3
-5452953875*x^2-7998223215*x+4461221029,75,75,[3,5,5],15.793,132.7]

Evaluation took 180.3810 seconds (240.6020 elapsed) using 74460.010 MB.

[1,x^2-2,2,2,,0.133,0.001]
[2,x^3-3*x-1,3,3,,0.169,0.003]
[3,x^4-2,8,8,[2,2,2],0.168,0.005]
[4,x^4+x^2-1,8,8,[2,2,2],0.16,0.007]
[5,x^4-2*x^3+2*x^2+2,12,12,[3,2,2],0.174,0.029]
[6,x^4+2*x^3+3*x^2+4*x+5,24,24,[2,3,2,2],0.273,0.126]
[7,x^4+x+1,24,24,[2,3,2,2],0.224,0.087]
[8,x^5-2,20,5,,0.182,0.002]
[9,x^5-5*x+12,10,10,[2,5],0.18,0.023]
[10,x^5+20*x+32,10,10,[2,5],0.163,0.032]
[11,x^5+11*x+44,10,10,[2,5],0.146,0.025]
[12,x^5+x^4-4*x^3-3*x^2+3*x+1,5,5,,0.141,0.008]
[13,x^5+100*x^2+1000,20,5,,0.161,0.012]
[14,x^6+x^3+1,6,3,,0.139,0.002]
[15,x^6-2,12,6,[2,3],0.167,0.004]
[16,x^5-5*x^3+5*x-5,20,10,[2,5],0.171,0.023]
[17,x^6+x^3+7,18,9,[3,3],0.162,0.009]
[18,x^6-3*x^4+1,24,24,[2,3,2,2],0.243,0.061]
[19,x^6+x^4-9,24,24,[2,3,2,2],0.218,0.083]
[20,x^6+x^4-8,48,48,[2,2,3,2,2],0.612,0.702]
[21,x^7+7*x^3+7*x^2+7*x-1,14,7,,0.17,0.04]
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,,0.175,0.028]
[23,x^7-2,42,7,,0.185,0.004]
[24,x^8-2*x^6-x^4+7*x^2-5*x+1,16,16,[2,2,2,2],0.274,0.033]
[25,x^9+x^8+3*x^6+3*x^3-3*x^2+5*x-1,18,18,[2,3,3],0.241,0.048]
[26,x^12-3,24,12,[2,2,3],0.52,0.009]
[27,x^7-14*x^5+56*x^3-56*x+22,21,7,,0.177,0.029]
[28,
x^15-x^14-14*x^13+13*x^12+78*x^11-66*x^10-220*x^9+165*x^8+330*x^7-210*x^6
-252*x^5+126*x^4+84*x^3-28*x^2-8*x+1,15,15,[3,5],0.605,0.189]

<permutation group with 72 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2

Evaluation time: 8.8
"Done","Done","Done","Done"
// Time 8.8
// Total time 8.8
[29,x^6+x^4-x^2+5*x-5,72,72,[2,2,2,3,3],4.448,15.771]
*** nfisincl: Warning: increasing stack size to 16000000.
*** write: Warning: increasing stack size to 16000000.
<permutation group with 13 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
"Done","Done","Done","Done"
// Time 0
// Total time 0
[30,2*(x+1)^13-x^13,156,13,,12.067,8.419]
*** nfisincl: Warning: increasing stack size to 16000000.
<permutation group with 49 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2

Evaluation time: 0.43
"Done","Done","Done","Done"
// Time 0.43
// Total time 0.43
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
"Done","Done","Done","Done"
// Time 0
// Total time 0
[31,
x^14+28*x^11+28*x^10-28*x^9+140*x^8+360*x^7+147*x^6+196*x^5+336*x^4-546*x^3
-532*x^2+896*x+823,98,49,[7,7],5.794,25.028]

*** nfisincl: Warning: increasing stack size to 16000000.
*** write: Warning: increasing stack size to 16000000.
<permutation group with 75 generators>
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2

Evaluation time: 8.13
"Done","Done","Done","Done"
// Time 8.13
// Total time 8.13
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
"Done","Done","Done","Done"
// Time 0.35
// Total time 0.35
// Using locale /usr/share/locale/
// ja_JP.UTF-8
// /usr/share/locale/
// giac
// UTF-8
// Maximum number of parallel threads 2
"Done","Done","Done","Done"
// Time 0.01
// Total time 0.01
[32,
x^15-470*x^13-305*x^12+71840*x^11+85357*x^10-4292700*x^9-3714805*x^8
+119761820*x^7+25284495*x^6-1542190154*x^5+717324725*x^4+7178878600*x^3
-5452953875*x^2-7998223215*x+4461221029,75,75,[3,5,5],15.858,123.12]

Evaluation took 179.0620 seconds (218.4700 elapsed) using 74669.995 MB.

# GGRR on SageMath def GGRR():
L1=gp('px='+Px+';DA=subst(nfsplitting(px),x,A);ls3=select(s->s>2,cs3=factor(poldegree(DA,A))[,1]~);PAgs=[[sb=nfsubfields(PA=polcyclo(j,ci=eval(concat(\"c\",j)))),PA,lift(nffactor(nfinit(PA),subst(lift(nffactor(nfinit((sf=sb[#sb-1])),DA)[,1]~),ci,sf))[,1]~)][2..3]|j<-ls3];');
L2=maxima('(kill(all),linsolvewarn:off,gcd:subres,[DA,PAgs,ls3]:'+gp.eval('[DA,PAgs,ls3]')+',map(lambda([s],tellrat(s)),PAs:map(first,PAgs)),algebraic:on,DAC:DA,map(lambda([s],DAC:gcd(DAC,s,A)),map(last,PAgs)),algebraic:off);');
L3=gp('[PAs,DAC]='+maxima.eval('[PAs,DAC]')+';default(realprecision,max(200,floor(1.5*poldegree(DA,A))));nRS=polroots(substvec(DAC,apply(s->variable(s),PAs),[polroots(j)|j<-PAs]));lRA=apply(s->[s],RA=nfisincl(px,DA));nRAS=round(10^10*[subst(lRA,A,nRSi)|nRSi<-nRS]);M=Map(Mat([n1=nRAS~,vector(#n1,i,i)~]));GG2=apply(L->(apply(s->mapget(M,s),L)),nRAS);');
L4=maxima('(load(\"ratpow\"),[RA,GG2]:'+gp.eval('[RA,GG2]')+',if not(DA=DAC) then print(concat(\"using a factor \",string(DAC),\" over the cyclotomic_\",apply(\"*\",ls3))),lRA:length(RA),solc:subst(linsolve(ratp_dense_coeffs((cis:makelist(concat(c,i),i,1,lRA)).RA-A,A),cis),cis),KK:subst(map(lambda([s],s=0),listofvars(solc)),solc),RS:rat(subst(map(\"=\",cis,RA),map(lambda([s],KK.s),subst(map(\"=\",GG2,cis),GG2)))),GGRS:map(\"[\",GG2,RS),GG2);');
L5=gap.eval('GG2:='+gp.eval('GG2')+';;GG:=Group(List(GG2,s->PermList(s)));LL:=[GG];;while Size(GG)>1 do GG:=MaximalNormalSubgroups(GG);Append(LL,[GG]);od;GGG:=List(List(LL,Elements),L->List(L{[2..Length(L)]},s->ListPerm(s,Length(GG2))));;');
L6=maxima('(S:map(lambda([s],append(['+gp.eval('GG2')+'],s)),'+gap.eval('GGG').replace('\n','')+'),[cs1,cs2,cs3]:csGG:[S,S:map(length,S),makelist(S[i-1]/S[i],i,2,length(S))]);');
L7=maxima('redLR():=block([],linsolve_params:off,AL1:reverse(makelist(A^j=concat(A,j),j,1,hipow(DA0,A)-1)),AL0:map(rhs,AL1),TIM0:elapsed_real_time(),prd:1,sAL1:(makelist(remainder(prd:(prd*da1n*ai),DA0,A)-(na1n*a[i])^j,j,1,pp-1)),sAL1:rat(subst(AL1,sAL1)),ALsol:linsolve(sAL1,AL0),print([i],\"~\",elapsed_real_time()-TIM0),ans:subst(x=A,LRx1A:subst(ALsol,subst(AL1,LRx1))),ans);mul(A,B):=map(lambda([s],A[s]),B);');
L8=maxima('load(\"grobner\");(aiA:DAs:DA0s:[],GGRS0:GGRS,kill(cpd),map(lambda([s],cpd[s]:cpd0[s]:s),GGRS),w:if member(2,ls3:listify(setify(cs3))) then -1 else 1,DP:map(lambda([s],Ri:rat(((ci:concat(c,s))^s-1)/(ci-1)),w:w*ci, [Ri,ci]),delete(2,ls3)),ppp:apply(\"*\",ls3),map(lambda([s],tellrat(s)),DP),algebraic:on,DA0:DAC,for i:1 thru length(cs3) do (push(DA0,DA0s),pp:cs3[i],map(lambda([s],tellrat(s)),DP),T0:T:cs1[i+1],for g in (Ti:cs1[i]) while member(h:g,T) do 0,glim:1400,if not(DA=DA0) and (slength(string(DA0))<=glim or i>1) then (map(lambda([s],cpd[s]:remainder(cpd[s],DA0,A)),Ti)),DA00:DA0, LRx1:LRy1:1,if slength(string(DA0))>glim*100 then (  LRy1:(remainder(expand_giac(),DA0,A)) )else for t in T do LRy1:(remainder(LRy1*(1+y*rat(cpd[t])),DA0,A)), LRx1:(ratnum((subst(y=-1/x,LRy1)))),if member(A,listofvars(LRx1)) then( LRy1:rat(LRy1), degy:0,unless hipow(coeff(LRy1,y,degy),A)>0 do degy:degy+1, tellrat(y^(degy+1)), LRy:makelist((prd:1,for t in (T:map(lambda([u],mul(h,u)),T)) do prd:(remainder(prd*(1+y*cpd[t]),DA0,A)),prd),j,2,pp), untellrat(y), push(LRy1,LRy), LRA:coeff(rat(LRy),y,degy),ai:0,for cnt:1 while ai=0 do(CNT:cnt,aix:rat(sum(w^mod(ppp/pp*cnt*j,ppp)*LRA[1+j],j,0,pp-1)), ai:if aix=0 then 0 else aix),if member(A,listofvars(ai)) then(a[i]:concat(e,i),prd:1,for i:1 thru pp do prd:remainder(prd*ai,DA0,A),aipp:rat(prd-a[i]^pp),a1n:num(aipp), na1n:ifactors(abs(poly_content(coeff(a1n,a[i],0),map(second,DP)))), na1n:apply(\"*\",map(lambda([s],s^floor(s/cs3[i])),na1n)), da1n:ifactors(abs(coeff(a1n,a[i],pp))), da1n:apply(\"*\",map(lambda([s],s^ceiling(s/cs3[i])),da1n)),a1n:poly_normalize(subst(a[i]=na1n/da1n*a[i],a1n),[a[i]]),tellrat(a1n),push([a1n,a[i]],DP),push([aiai:da1n/na1n*ai-a[i],a[i]],aiA),LRx:[LRx1],if hipow(rla:remainder(LRx1,aiai,A),A)=0 then (DA0:gcd(r1:subst(x=A,rla),DA0,A)) else (DA0:r1:redLR()),DA0:poly_normalize(DA0,[A])))),algebraic:off,map(lambda([s],untellrat(s)),DP),rr:append([[DA0,A]],DP));');
return maxima('rr')

%time Px='1+x+x^2+x^3+x^4+x^5+x^6+x^7+x^8+x^9+x^10';GGRR()
maxima('DA');maxima('RA');maxima('GGRS0');maxima('csGG')
%time Px='x^4+x+1';GGRR()
%time Px='x^6+x^4-x^2+5*x-5';GGRR() #CPU time: 0.04 s, Wall time: 45.74 s
%time x=var('x');for i in range(2,31):Px=str(x^i-2);print(Px);GGRR() #CPU time: 1.00 s, Wall time: 1033.85 s
maxima('DA');maxima('RA');maxima('GGRS0');maxima('csGG')
PLS=["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-5*x^3+5*x-5","x^6+x^3+7","x^6-3*x^4+1","x^6+x^4-9","x^6+x^4-8","x^7+7*x^3+7*x^2+7*x-1","x^7-14*x^5+56*x^3-56*x+22","x^7-2","x^8-2*x^6-x^4+7*x^2-5*x+1","x^9+x^8+3*x^6+3*x^3-3*x^2+5*x-1","x^12-3","x^7-14*x^5+56*x^3-56*x+22","x^15-x^14-14*x^13+13*x^12+78*x^11-66*x^10-220*x^9+165*x^8+330*x^7-210*x^6-252*x^5+126*x^4+84*x^3-28*x^2-8*x+1","x^14+28*x^11+28*x^10-28*x^9+140*x^8+360*x^7+147*x^6+196*x^5+336*x^4-546*x^3-532*x^2+896*x+823","x^15-470*x^13-305*x^12+71840*x^11+85357*x^10-4292700*x^9-3714805*x^8+119761820*x^7+25284495*x^6-1542190154*x^5+717324725*x^4+7178878600*x^3-5452953875*x^2-7998223215*x+4461221029"]
%time for Px in PLS:print(Px);GGRR() 

# 代数体上の多項式の高速乗算

まず， の計算に現れるデータのサイズがどの程度となりえるのかを例示しましょう．入力のサンプルは次の４つです．

Hirokazu Anai,Kazuhiro Yokoyama, Radical Representation of Polynomial Roots
http://www.jssac.org/Editor/Suushiki/V09/No1/V9N1_108.pdf より．

Andreas Distler, Ein Algorithmus zum Lösen einer Polynomgleichung durch Radikale
http://www.icm.tu-bs.de/ag_algebra/software/distler/Diplom.pdf より．

Galois 群の位数は，順に 72，98，156，156． 上の最小分解体  上の定義多項式の文字列としてのサイズは，順に 2208，1367，3638，159 バイト．  上の primitive element の円分体上の最小多項式 の文字列としてのサイズは，順に 2208，2284，1413，132 バイト． 全体の文字列としてのサイズは，順に 3276035，4452826，20461351，88263 バイト．

（Ａ） のまま用いる方法
（Ｂ）  による剰余を用いる方法 の場合，（Ａ），（Ｂ）は同じ）があり，展開には
（Ｃ）maxima プロパーの（乗算毎に による剰余をとる）方法
（Ｄ）外部ツール（Julia の Nemo パッケージ http://nemocas.org/ の ResidueRing）を呼び出す方法
があります．各場合の入力から出力までの時間は，次の通りです．

（Ａ），（Ｃ）の場合，順に 23.89，127.47，643.56，7.15 秒．
（Ａ），（Ｄ）の場合，順に 19.88，50.35，175.85，7.43 秒．
（Ｂ），（Ｃ）の場合，順に 24.05，86.62，128.49，3.82 秒．
（Ｂ），（Ｄ）の場合，順に 19.80，50.75，229.23，3.86 秒．

なお，上記の Nemo の ResidueRing は同パッケージの NumberField や Hecke パッケージ http://thofma.github.io/Hecke.jl/latest/ の代数体関連の関数よりも高速です．また，giac https://www-fourier.ujf-grenoble.fr/~parisse/giac.html で乗算毎に rem を用いる方法では

（Ｂ），（Ｄ）の場合，順に 21.83，51.03，397.46，3.94 秒

となります．

# 改善策と問題点 の中に でないものがあるので，それを とすることができます．

この方法では，全ての係数を生成するのは に直結した のみで，後は つの係数 に対する処理となり，相応の高速化が期待できるように思えます．

ところが，この を定義通りに計算すると，  の基本対称式の つゆえ， となる 上の多項式 があり， の準同型性から， なので， に対応する の根を  に代入したもの，つまり，最も悪い場合には， の次数のおよそ 乗の次数の多項式（しかも係数は分子分母の桁数が大変なことになっている既約分数）の簡約整理が必要となります．

これまでのところ
【１】前々回の記事のままの実装
【２】上記の改善策による実装
そして
【３】上記の に対する次数よりも大きな次数の項を消去しながら を簡約整理する実装
の処理時間は，サイズの小さな問題では【２】が短く，大きな問題では【１】が短く，【３】はその中間といった実験結果を得ています．