2019-07-18

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

今回の冪根拡大の構成では，中間体上の線型方程式系（連立一次方程式）を解くことが一つの柱となっています（ http://ehito.hatenablog.com/entry/2019/04/11/193241 ）．それは例えば，定義多項式が

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[2];

julia> MX=matrix(MXl[4],MXl2,MXl2,MXl[1]);

julia> cx=matrix(MXl[4],1,MXl2,MXl[3]);

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[1])),DP),algebraic:on,linsolve(sAL1,AL0))$
Evaluation took 0.0040 seconds (0.0040 elapsed) using 1.843 MB.

(%i12) (map(lambda([s],tellrat(s[1])),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[1])),DP),algebraic:on,linsolve(sAL1,AL0))$
Evaluation took 458.1150 seconds (472.1600 elapsed) using 92401.623 MB.

そして，

(%i15) (map(lambda([s],tellrat(s[1])),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.

となります．

2019-07-08

位数の取得

前回も触れたように，nfsplitting の処理時間は結果の多項式の次数を第二引数として与えることで一般に短縮されます．この性質を利用するため本プログラムでは，その次数，つまり，入力多項式 $f$ の $\mathbb{Q}$ 上の Galois 群 $G$ の位数を，予め用意したある範囲内の全ての cycle index polynomial (https://en.wikipedia.org/wiki/Cycle_index) の中から $G$ の cycle index polynomial $Z(G)$ を特定することで取得しています．

置換群の cycle index polynomial とは，その群に属する元の cycle type ( https://groupprops.subwiki.org/wiki/Cycle_type_of_a_permutation ) の頻度を変数名と次数，そして係数により，一つの多変数多項式で表したもので，群の位数の逆数は（1個しかない）単位元の cycle type の頻度（対応する単項式の係数）として現れます．参照する cycle index polynomial のデータは，次数が $37$ 以下の $f$ の Galois 群を全て含む，GAP の Transitive Groups Library transgrp ( https://www.gap-system.org/Packages/transgrp.html ) から抽出したもので，その制約から次数が $38$ 以上の $f$ についてはノーヒントで nfsplitting を実行しています．なお，同型でない群でも同じ cycle index polynomial を有するものが存在するため，群自体の特定には至らないことはよく目にする注意です．

実際の $f$ から $Z(G)$ を特定するには，前回も引用した「正の数 $M$ 以下で $f$ の判別式の因数ではない素数 $p$ を法とする $f$ の既約因子の次数型（どのような次数の因子がいくつずつあるか）の密度が， $M\to+\infty$ のとき， $G$ の同じ型の cycle type の頻度に収束する」という Chebotarev's density theorem を踏まえ， $p$ が値の小さいほうから 10^3 個目の素数となる，または，単位元の cycle type が 3 回現れるまで，既約分解を繰り返し，各次数型と密度から cycle index polynomial と同じ書式の多項式を作り，用意したデータのうち，この多項式との差の係数ベクトルの $l_2$ ノルムが最小となるものを $Z(G)$ としています．

以下，この処理を組み込んだ関数 nfsplitting2 と標準の nfsplitting の実行例です．

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

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

2019-06-23

多項式の Galois 群計算

本プログラムの 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

となります．

2019-05-22

円分体上の Galois 群

前回までの冪根拡大では，円分体を基礎体とする一方，中間体の生成は $\mathbb{Q}$ 上の Galois 群の組成列に従っており， $M_i(X,0)\in F_{i-1}[X]$ となる $i$ の発生原因となっていました．今回はこれを「まず円分体とその上の分解体の Galois 群を構成し，その組成列を用いて中間体を生成する手順」に，また，「最小多項式の文字列としてのサイズに応じて， $M_i(X,0)$ の展開方法を切り替える方式」（下の表示に出力があるものは各外部ツール，その他は組み込み関数を使用）に改めましたので，サンプルの実行結果を公開します．

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

展開の外部ツールとして Nemo を用いた場合．

[1,x^2-2,2,2,[2],0.152,0.002] 
[2,x^3-3*x-1,3,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,[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,[5],0.172,0.009] 
[13,x^5+100*x^2+1000,20,5,[5],0.21,0.012] 
[14,x^6+x^3+1,6,3,[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,[7],0.184,0.039] 
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.19,0.033] 
[23,x^7-2,42,7,[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,[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,[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.

展開の外部ツールとして Giac を用いた場合．

[1,x^2-2,2,2,[2],0.133,0.001] 
[2,x^3-3*x-1,3,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,[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,[5],0.141,0.008] 
[13,x^5+100*x^2+1000,20,5,[5],0.161,0.012] 
[14,x^6+x^3+1,6,3,[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,[7],0.17,0.04] 
[22,x^7-14*x^5+56*x^3-56*x+22,21,7,[7],0.175,0.028] 
[23,x^7-2,42,7,[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,[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
Added 0 synonyms

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
Added 0 synonyms
"Done","Done","Done","Done"
// Time 0
// Total time 0
[30,2*(x+1)^13-x^13,156,13,[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
Added 0 synonyms

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
Added 0 synonyms
"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
Added 0 synonyms

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
Added 0 synonyms
"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
Added 0 synonyms
"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.

2019-05-15

GGRR on SageMath

今回の「可解な方程式の全ての根を Galois 理論に基づく方法により冪根で表すプログラム」を SageMath ( http://www.sagemath.org/ ) に移植しました．

f:id:ehito:20190515191939j:plain

SageMath は https://cocalc.com/?utm_source=sagemath.org&utm_medium=icon から，Facebook，Github，Google，Twitter のアカウントでも利用できます．

適当な名前の project と，その中に Sage work sheet タイプのファイルを作成して，次の定義をコピペし，Run ボタンを押せば準備完了です．

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])[1]),DA)[,1]~[1]),ci,sf[2]))[,1]~[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)[1]|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[1]~,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[1],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)[1];Append(LL,[GG]);od;GGG:=List(List(LL,Elements),L->List(L{[2..Length(L)]},s->ListPerm(s,Length(GG2[1]))));;');
    L6=maxima('(S:map(lambda([s],append(['+gp.eval('GG2[1]')+'],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[1]]:cpd0[s[1]]:s[2]),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[1])),DP),algebraic:on,DA0:DAC,for i:1 thru length(cs3) do (push(DA0,DA0s),pp:cs3[i],map(lambda([s],tellrat(s[1])),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[1]^floor(s[2]/cs3[i])),na1n)), da1n:ifactors(abs(coeff(a1n,a[i],pp))), da1n:apply(\"*\",map(lambda([s],s[1]^ceiling(s[2]/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[2])),DP),rr:append([[DA0,A]],DP));');
    return maxima('rr')

以下のサンプルを一行ずつコピペ，Run ボタンを押せば，それぞれが実行されます（最初だけ外部プログラムの起動に時間が掛かります）．

%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()

2019-05-05

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

一連の冪根拡大の構成では， $\{\sigma(A)\,|\,\sigma\in G_i\}$ を根全体とする多項式 $M_i(X,0)$ の係数 $C$ から冪根を生成し，その過程で得られる情報を $M_i(X,0)$ に還元して分解体の primitive element の最小多項式 $g_i$ を得ていますが， $C$ については $\sigma(A)$ 達の基本対称式を漸化式から得る方法，最小多項式についても幾度か述べたように例えば $g_i=\text{gcd}(g_{i-1},e_i-r_i)$ といった位置付けがあり，規模の小さな問題に対しては実用的でもある一方，こうした工夫はサイズの大きな問題に対しては実装上の問題により逆効果となります．つまり， $M_i(X,0)$ をそのまま展開することとなり，そのうちとくに負担となりえる $M_1(X,0)$ を如何に整理，展開するべきかという視点から今回の表題に至ります．

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

例１　 $x^6+x^4-x^2+5 x-5$
Hirokazu Anai,Kazuhiro Yokoyama, Radical Representation of Polynomial Roots
http://www.jssac.org/Editor/Suushiki/V09/No1/V9N1_108.pdf より．
例２　 $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$
Andreas Distler, Ein Algorithmus zum Lösen einer Polynomgleichung durch Radikale
http://www.icm.tu-bs.de/ag_algebra/software/distler/Diplom.pdf より．
例３　 $2(1+x)^{13}-x^{13}$
例４　 $x^{13}-2$

例１，２は各先行研究において処理時間が最大，或いは結果が得られなかったもの（なお，これらは冪根拡大に限らない中間体の列を生成した後，それらの primitive element を冪根に置き換える方法によるもので，例１，２はそれぞれ逐次拡大を推奨，単純拡大を採用しています）．例４は標準的（？）なサンプルで，分解体は例３と同型ですが，生成されるデータサイズは全く異なります．

Galois 群の位数は，順に 72，98，156，156．
$\mathbb{Q}$ 上の最小分解体 $\operatorname{sp}(f)$ の $\mathbb{Q}$ 上の定義多項式の文字列としてのサイズは，順に 2208，1367，3638，159 バイト．
$\operatorname{sp}(f)$ の $\mathbb{Q}$ 上の primitive element $A$ の円分体上の最小多項式 $g_0$ の文字列としてのサイズは，順に 2208，2284，1413，132 バイト．
$\sigma(A)$ 全体の文字列としてのサイズは，順に 3276035，4452826，20461351，88263 バイト．

実際の $\sigma(A)$ の利用には
（Ａ） $\sigma(A)$ のまま用いる方法
（Ｂ） $\sigma(A)$ の $g_0$ による剰余を用いる方法
（ $f_{\mathbb{Q}}=g_0$ の場合，（Ａ），（Ｂ）は同じ）があり，展開には
（Ｃ）maxima プロパーの（乗算毎に $g_0$ による剰余をとる）方法
（Ｄ）外部ツール（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 秒

となります．

2019-04-19

改善策と問題点

前々回の記事（ http://ehito.hatenablog.com/entry/2019/04/11/193241 ）で述べた方法では， $M_i(X,1)$ 以降をただ $1$ つの $r_i$ を見付けるために生成していますが，実は $M_i(X,0)$ を得た段階で，それが $F_{i-1}[X]$ に属する場合には，先述の通り拡大を行わず，属さない場合，つまり， $M_i(X,0)$ が $C\in F_{i}\setminus F_{i-1}$ となる係数 $C$ をもつ場合には，その Lagrange Resolvent

$\begin{align} R_i(k)&=\sum_{j=0}^{p_i-1} \zeta^{j k P/p_i} \tau^j(C)\ \ (k=1,\ldots,p_i-1) \end{align}$

の中に $0$ でないものがあるので，それを $r_i$ とすることができます．

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

ところが，この $\tau^j(C)\ (j=1,\ldots,p_i-1)$ を定義通りに計算すると， $C$ は $\sigma(A)\, (\sigma\in G_i)$ の基本対称式の $1$ つゆえ， $C=S(A)$ となる $F_{i-1}$ 上の多項式 $S$ があり， $\tau$ の準同型性から， $\tau^j(S(A))=S(\tau^j(A))$ なので， $\tau^j$ に対応する $f_{\mathbb{Q}}$ の根を $S(A)$ の $A$ に代入したもの，つまり，最も悪い場合には， $f_{\mathbb{Q}}$ の次数のおよそ $2$ 乗の次数の多項式（しかも係数は分子分母の桁数が大変なことになっている既約分数）の簡約整理が必要となります．

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