QE による検証
Mathematica の QE のデフォルトの議論領域は
Reduce[Exists[x, a*x^2 == 1]]
に対する出力
からも判るように複素数体であり，前回公開した Mathematica での実装の目的の 1 つは，この QE による結果の検証でした．
Mathematica で前回の各関数の定義コードを実行した後，例えば，
(p1=PL[[1]])@x
を
mp[p1];nGG[];cs[];RR[];InputForm@DP
のように処理すると
といった出力が得られます．
このとき，A についての 2 つの条件
(1) A が分解体の primitive element の定義多項式 p2 の零となる
(2) A の多項式で表された根全体 r1 /. pe -> A が入力の多項式 p1 の根全体となる
が等価であるという式
ForAll[A, Equivalent[p2@A == 0, ForAll[x, Times @@ (x - r1 /. pe -> A) == p1@x]]]
すなわち
を複素数体上で
Reduce[%,Complexes]
のように QE すると
と答えてくれます．
また，A についての 2 つの条件
(1) A が分解体の primitive element の定義多項式 p2 の零となる
(3) 出力 DP に属する全ての多項式の共通の零点 (w,e[1]) が存在する
が等価であるという式
ForAll[A, Equivalent[p2@A == 0, Exists[Evaluate@Rest@vs, And @@ (# == 0 & /@ DP)]]]
すなわち
を複素数体上で
Reduce[%,Complexes]
のように QE すると
と答えてくれます．
Mathematica 推参
オリジナルプログラムの元になったコードは Mathematica で書かれていますが，汎用性を重視した筆致なので，今回は
https://reference.wolfram.com/language/tutorial/AlgebraicNumberFields.html.ja
にある代数体関連の組み込み関数などを用いて書いてみました（取り敢えず動く程度）．
なお，https://develop.open.wolframcloud.com/app/ の Create a New Notebook ボタンを押すと Mathematica がオンラインで利用できます．
(* mp *)
mp[p1_]:=Block[{c,x},
n1=Range@Exponent[p1@x,x];
r1=ToNumberField[AlgebraicNumber[Root[p1,#],{0,1}]&/@n1,All];
p2=(pe=r1[[1,1]])[[1]];(* pe is not a "AlgebraicNumber" *)
c1=Map[c,n1];
c1=c1/.SolveAlways[x==r1.c1/.pe->x,x][[1]]/.c[_]->0];
(* nGG *)
nGG[]:=Block[{x},
r12=Map[ToString,Chop[N[r1/.pe->x/.NSolve[p2@x,WorkingPrecision->50],10]],{2}];
gg=PermutationCycles/@(r12/.Thread[r12[[pe[[2]]]]->n1]);
r2=r1.Permute[c1,#]&/@gg];
(* csc and cs *)
csc[gg0_]:=Block[{},
cc=Rest@GroupOrbits[PermutationGroup@gg0,gg0];
Last@SortBy[Select[Append[#,Cycles[{}]]&/@(Union@@#&/@Subsets@cc),2*Length@#<=Length@gg0&&
GroupOrder@PermutationGroup@#==Length@#&],Length]];
cs[]:=Block[{},
cs1=NestWhileList[csc,gg,Length@#>1&];
cs2=Length/@cs1;cs3=Table[cs2[[i]]/cs2[[i+1]],{i,Length@cs2-1}]];
(* RR *)
Off[ToNumberField::nalg];
RR[]:=Block[{},
w1=ToNumberField@Exp[2*Pi*I/(ppp=LCM@@cs3)];
Print["using Cyclotomic_",ppp];
DP={Cyclotomic[ppp,w],
If[ppp==2,p2@A,First@Last@FactorList[p2@A,Extension->w1]/.p_AlgebraicNumber->ToNumberField[p,w1]/.w1[[1]]->w]};
For[i=1,i<Length@cs2,i++,LRx=Sum[w^Mod[ppp/pp*j,ppp]*Product[x-r2[[Position[gg,s][[1,1]]]],{s,PermutationProduct[PermutationPower[Complement[cs1[[i]],cs1[[i+1]]][[1]],j],#]&/@cs1[[i+1]]}],{j,pp=cs3[[i]]}];
ai=0;j=1;While[ai==0,ai=Simplify[LRx/.{w->w^j,pe->A},And@#==0&/@DP];j++];
DP=GroebnerBasis[Append[DP,Last@CoefficientList[ai,x]-e[i]],vs=Join[{A},Map[e,Reverse@Range@i],{w}]];(*Print@nputForm@DP*)]];
実行例．例によって，入力多項式の根の表示は長いので略します．
(* expamles *)
PL={#^3-3*#-1&,#^4-2&,#^4+#^2-1&,#^4-2*#^3+2*#^2+2&,#^4+2*#^3+3*#^2+4*#+5&,#^4+#+1&,#^5-2&,#^5-5*#+12&,#^5+20*#+32&,#^5+11*#+44&,#^5+#^4-4*#^3-3*#^2+3*#+1&,#^5+100*#^2+1000&,#^6+#^3+1&,#^6-2&};
(mp[#];nGG[];cs[];RR[];Print@InputForm@DP)&/@PL;//Timing
円分体の利用
RR シリーズでは冪根が代数的整数となるよう変数の導入時に係数を調整しています．これは algebraic モードへの適用を念頭に置いたものですが，1 の原始冪乗根の中間体上の定義式は必ずしも monic でないため，RR6m では，tellrat への入力の前に線形変換を施し，拡大体上の gcd 計算（それは係数に含まれる変数の線形変換と可換）の後，逆変換を施していました．
こうした局所的な往復を中間体の生成以前に 1 の原始冪乗根を導入すること，つまり，分解体の基礎体を然るべき円分体とすることで解消したのが RR7 です．
利用する円分体 Q(c_n) の n は中間体の拡大次数の最小公倍数でよいのですが，RR7 では，因数分解のコストを鑑み，各奇素因数についての円分体上での定義式の gcd を分解体の定義式としています．
今回は必要なコードを一括掲載しましたので，以下を例えば RR7.mac としてファイルに保存，load("RR7.mac")$ と読み込めば，下記の例のように実行できます．
/* examples */
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]$
PL3:[x^5-5*x^3+5*x+5,x^7+7*x^3+7*x^2+7*x-1,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]$
/* verbose mode on, off */
pon():=(printn(s):=print(s),printx(s):=print(s))$
poff():=(printn(s):=0,printx(s):=0)$
poff()$
/* mp */
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"),
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), printn(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)
)$
/* nGG */
nGG(DA):=block([i],
er:sort(map(lambda([s],abs(s[1]-s[2])),listify(map(listify,powerset(setify(map('rhs,allroots(p))),2)))))[1]/10,
nRS:allroots(%i*DA),
GG2:nRAS:expand(map(lambda([s],subst(s,RA)),nRS)),
i:0,for s in nRAS[1] do (
i:i+1,map(lambda([t],j:1,while(j>0) do if not(integerp(t[j])) and abs(t[j]-s)<er then (t[j]:i,j:0) else j:j+1,t),GG2)),
RS:rat(map(lambda([s],KK.s),fullmapl(lambda([s],RA[s]),map(lambda([s],firstn(s,length(KK))),GG2)))),
GGRS:map("[",GG2,RS),
GG2)$
/* cs */
mul(A,B):=map(lambda([s],A[s]),B)$
inv(A):=block([B,i],B:A*(i:0),for s in A do B[s]:i:i+1,B)$
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))])$
/* tools */
factor2list(f0):=block([f:num(f0)],if op(f)="-" then args(-f) elseif op(f)="*" then args(f) else [f])$
load("grobner")$
/* Root Of Unity */
ROU(x,N):=rat(apply("*",makelist((x^(N/s)-1)^moebius(s),s,listify(divisors(N)))))$
/* RR7 */
RR7(csGG):=block([],[cs1,cs2,cs3]:csGG,aiA:[],DAs:[],
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,
push(factor2list(factor(DA,Ri))[1],DAs),
[Ri,ci]),delete(2,ls3)),
printn(DP),ppp:apply("*",ls3),
map(lambda([s],tellrat(s[1])),DP),algebraic:on,
DA0:DA,map(lambda([s],DA0:gcd(DA0,s,A)),DAs),
if hipow(DA0,A)<hipow(DA,A) then
print(concat("using a factor ",string(DA0)," over the cyclotomic_",ppp)),
for i:1 thru length(cs3) do (
pp:cs3[i],map(lambda([s],tellrat(s[1])),DP),
printn("phase:"),printn([i,pp]),
/* non-zero aix and ai */
T:cs1[i+1],for g in cs1[i] while member(h:g,T) do 0,
LRx:makelist((prd:1,map(lambda([t],prd:remainder(prd*(x-assoc(t,GGRS)),DA0,A)),T:map(lambda([u],mul(h,u)),T)),prd),j,1,pp),
printn("LRx:"),printn(LRx),
ai:0,for cnt:1 while ai=0 do
(aix:sum(w^mod(ppp/pp*cnt*j,ppp)*LRx[j],j,1,pp),
ai:if aix=0 then 0 else coeff(rat(aix),x,hipow(aix,x)),
printn(concat("ai_",cnt)),printn(ai)),
if member(A,listofvars(ai)) then
(
/* def. of a[i] */
a[i]:concat(e,i),
aipp:rat(remainder(ai^pp,DA0,A)-a[i]^pp),
printn("aipp:"),printn(aipp),
/* monic rather than primitive */
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)),
printn(na1n/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),
/* rel-def. is gcd */
printn(LD:[last(LRx),pp,ai,a1n,aiai,DA0,DP,tellrat()]),
DA0:if hipow(aiai,A)=1 then aiai else gcd(aiai,DA0,A),
DA0:poly_normalize(DA0,[A])
),
printn([[i],DA0,DP])
),algebraic:off,map(lambda([s],untellrat(s[2])),DP),
rr5:rr3:append([[DA0,A]],DP),
for i in rest(rr3) do
if not(member(i[2],listofvars(rr4:delete(i,rr3)))) then rr3:rr4,
rr3)$
実行例．
/* オリジナルプログラム付属の例 */ for i:1 thru 15 do (print([i]),print([p:PL[i],mp(p),RR7(cs(nGG(DA)))]))$ /* 15 次，21 次拡大などの例 */ for i:1 thru 7 do (print([i]),print([p:PL3[i],mp(p),RR7(cs(nGG(DA)))]))$
それぞれの timer_info()．
[[total,0.001112600037223153*sec,5373,5.978*sec,0],
[mp,0.2166666666666667*sec,15,3.25*sec,0],
[RR7,0.1295333333333333*sec,15,1.943*sec,0],
[nGG,0.03886666666666666*sec,15,0.583*sec,0],
[mul,2.706185567010309e-5*sec,3880,0.105*sec,0],
[inv,2.832415420928404e-5*sec,1271,0.036*sec,0],
[mnsg,6.216216216216216e-4*sec,37,0.023*sec,0],
[CC,5.135135135135136e-4*sec,37,0.019*sec,0],
[cs,5.333333333333334e-4*sec,15,0.008*sec,0],
[isG,1.891891891891892e-4*sec,37,0.007*sec,0],
[ROU,2.5e-4*sec,12,0.003*sec,0],
[factor2list,8.333333333333333e-5*sec,12,0.001*sec,0]]
[[total,0.003578925872983459*sec,9794,35.052*sec,0],
[mp,4.366*sec,7,30.562*sec,0],
[RR7,0.4417142857142857*sec,7,3.092*sec,0],
[nGG,0.1021428571428571*sec,7,0.715*sec,0],
[mnsg,0.0173*sec,20,0.346*sec,0],
[mul,2.711284807034685e-5*sec,8188,0.222*sec,0],
[inv,4.290657439446366e-5*sec,1445,0.062*sec,0],
[isG,4.736842105263157e-4*sec,57,0.027*sec,0],
[CC,0.00115*sec,20,0.023*sec,0],[ROU,2.5e-4*sec,8,0.002*sec,0],
[factor2list,1.25e-4*sec,8,0.001*sec,0]]
inv 再考
置換の表現には，リストタイプと巡回タイプとがあり，それぞれ積，逆元の実装に適していますが，今回のコードは根の置換という出自からリストタイプを利用しており，逆元の実装に工夫の余地がありました．現在公開中の inv は
inv(A):=map(second,sort(map("[",A,sort(A)),lambda([s,t],s[1]<t[1])))$
という sort ベースのもので
(%i3) L:random_permutation(makelist(i,i,1,10)); (%o3) [4,9,8,3,5,10,7,1,2,6] (%i5) for j:1 thru 10^5 do inv(L)$ Evaluation took 19.6490 seconds (19.6760 elapsed) using 2073.668 MB.
といった具合でした．これを積（左作用）
mul(A,B):=map(lambda([s],A[s]),B)$
と同じく成分をインデックスと見做し，第 i 成分が s である入力のリストに対して，第 s 成分が i であるリストを出力とするもの
inv(A):=block([B,i],B:A*(i:0),for s in A do B[s]:i:i+1,B)$
に変更すると
(%i7) for i:1 thru 10^5 do inv(L)$ Evaluation took 4.1150 seconds (4.1320 elapsed) using 482.180 MB.
という結果となりました．
GCD の利用
現在のオリジナルプログラムや RR6 では，分解体の primitive element A の中間体 K(w,a) 上の定義式が Lagrange resolvents の和となることを利用していますが，http://ehito.hatenablog.com/entry/2019/03/02/163023 で述べたように，A の K(w,a) 上の定義式と a の K(w) 上の定義式は，それぞれ A の K 上の定義式と冪根 a の一次の定義式との A に関する GCD と終結式（を因数分解したもの）なので，代数体での計算が高速なら，これに基づく実装の価値があります．
残念ながら maxima の因数分解は逐次代数拡大に対応していませんが，GCD は対応（algebraic モード）しており，例えば
(%i4) tellrat(a2^3+21*a1+14,a1^2+a1+1)$algebraic:on$gcd:subres$ (%i7) gcd(A^12+4*A^10+24*A^8+48*A^6-560*A^4+3136,485968*a2+A^4*(18236*a1-9524)+A^10*(83*a1-563)+655424*a1+A^8*((-739*a1)-2672)+A^6*((-2012*a1)-12700)+A^2*(474768-28560*a1)+390432,A); Evaluation took 0.0020 seconds (0.0020 elapsed) using 255.969 KB. (%o7) A^2*((12*a1+4)*a2^2+28*a2+28)-28*a2^2+((-84*a1)-28)*a2+21*A^4+392 (%i8) algebraic:off$untellrat(a1,a2)$
といった具合に，tellrat の定義式が monic（つまり，扱えるのは代数的整数のみ）という制約さえクリアすれば，．．．という訳で，次の RR6m を作りました（m は monic の m）．
(%i15) fundef(RR6m);
(%o15) RR6m(csGG):=block([],[cs1,cs2,cs3]:csGG,kill(w),DA0:DA,Dws:[[DA,A]],DP:aiA:[],
for i thru length(cs3) do
(push([DA0,A],DP),w[i]:concat(a,2*i-1),a[i]:concat(a,2*i),
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)),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 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])),
printn(Dwf),
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)),
printn(Dwf),printn([ai,pp,DP]),aipp:rrem(ai^pp,DP)-a[i]^pp,a1n:num(rat(aipp)),
na1n:ifactors(abs(poly_content(subst(a[i] = 0,a1n),map(second,DP)))),
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:"),printn(na1n/da1n),
a1n:poly_normalize(subst(a[i] = (na1n*a[i])/da1n,a1n),[a[i]]),
DP:delete([DA0,A],DP),push([a1n,a[i]],DP),push([a[i]-(da1n*ai)/na1n,a[i]],aiA),
DA0:if hipow(ai,A) = 1 then a[i]-(da1n*ai)/na1n
else (DPr:reverse(DP),sDPr:map(second,DPr),
for j thru (lDPr:length(DPr)) do
DPr:subst((vj2:sDPr[j])
= vj2
/mc[j]
:coeff(num(rat(fj1:DPr[j][1])),vj2,hipow(fj1,vj2)),
DPr),DPm:map("[",fDPr:map(first,DPr),sDPr),
[DA0m,aiaim]:subst(makelist((vj2:sDPr[j]) = vj2/mc[j],j,1,lDPr),
[DA0,a[i]-(da1n*ai)/na1n]),
print([num(rat(DA0m)),num(rat(aiaim)),A,reverse(DPm)]),
map(tellrat,fDPr),algebraic:on,
gcdA:subst(makelist((vj2:sDPr[j]) = vj2*mc[j],j,1,lDPr),
gcd(aiaim,DA0m,A)),algebraic:off,map(untellrat,sDPr),
gcdA),printn([3,DA0,DP])),rr2:[solve(DA0,A)[1],DP],
append([[DA0,A]],DP))
実行例 for i:1 thru 15 do print([p:PL[i],mp(p),RR6m(cs(nGG(DA)))]) の timer_info() では，[RR6,0.1656*sec,15,2.484*sec,0] の約 3/4 の処理時間になり，反復除算の簡約 rrem の回数も 210 から 112 になっています．
[[total,0.001120116300199891*sec,5503,6.164*sec,0],
[mp,0.2176*sec,15,3.264*sec,0],[RR6m,0.1252*sec,15,1.878*sec,0],
[nGG,0.0378*sec,15,0.5670000000000001*sec,0],
[rrem,0.002366071428571429*sec,112,0.265*sec,0],
[inv,6.372934697088908e-5*sec,1271,0.081*sec,0],
[mul,1.597938144329897e-5*sec,3880,0.062*sec,0],
[CC,6.486486486486487e-4*sec,37,0.024*sec,0],
[mnsg,4.594594594594595e-4*sec,37,0.017*sec,0],
[isG,1.351351351351351e-4*sec,37,0.005*sec,0],
[factor2list,2.702702702702703e-5*sec,37,0.001*sec,0]]
微調整
http://ehito.hatenablog.com/entry/2019/02/24/234939 の実行例の timer_info() を見ると
[[total,0.001663981382026495*sec,5586,9.295*sec,0],
[mp,0.214*sec,15,3.21*sec,0],
[nGG,0.1895333333333333*sec,15,2.843*sec,0],
[RR6,0.167*sec,15,2.505*sec,0],
[rrem,0.002457142857142857*sec,210,0.516*sec,0],
[mul,2.216494845360825e-5*sec,3880,0.08600000000000001*sec,0],
[inv,6.136900078678206e-5*sec,1271,0.078*sec,0],
[CC,8.378378378378379e-4*sec,37,0.031*sec,0],
[mnsg,4.054054054054054e-4*sec,37,0.015*sec,0],
[isG,2.162162162162163e-4*sec,37,0.008*sec,0],
[factor2list,5.405405405405407e-5*sec,37,0.002*sec,0],
[cs,6.666666666666667e-5*sec,15,0.001*sec,0]]
のように nGG が重く，その原因は例によって bfloat であろう，という訳で，書き換えてみました．
nGG(DA):=block([i],
er:sort(map(lambda([s],abs(s[1]-s[2])),listify(map(listify,powerset(setify(map('rhs,allroots(p))),2)))))[1]/10,
nRS:allroots(%i*DA),
GG2:nRAS:expand(map(lambda([s],subst(s,RA)),nRS)),
i:0,for s in nRAS[1] do (
i:i+1,map(lambda([t],j:1,while(j>0) do if not(integerp(t[j])) and abs(t[j]-s)<er then (t[j]:i,j:0) else j:j+1,t),GG2)),
RS:rat(map(lambda([s],KK.s),fullmapl(lambda([s],RA[s]),map(lambda([s],firstn(s,length(KK))),GG2)))),
GGRS:map("[",GG2,RS),
GG2)$
（なお，こうした root separation の評価は計算機代数の古典的なテーマのひとつで，例えば https://core.ac.uk/download/pdf/82808521.pdf ）
結果は，．．．
[[total,0.001227491531467285*sec,5609,6.885*sec,0],
[mp,0.2095333333333333*sec,15,3.143*sec,0],
[RR6,0.1656*sec,15,2.484*sec,0],
[nGG,0.03653333333333333*sec,15,0.548*sec,0],
[rrem,0.002447619047619048*sec,210,0.514*sec,0],
[inv,7.002360346184106e-5*sec,1271,0.089*sec,0],
[mul,1.314432989690722e-5*sec,3880,0.051*sec,0],
[CC,8.378378378378379e-4*sec,37,0.031*sec,0],
[mnsg,4.054054054054054e-4*sec,37,0.015*sec,0],
[isG,2.432432432432433e-4*sec,37,0.009000000000000001*sec,0],
[factor2list,2.702702702702703e-5*sec,37,0.001*sec,0]]
実行例
円分多項式の例です．
(%i9) ROU(x,N):=rat(apply("*",makelist((x^(N/s)-1)^moebius(s),s,listify(divisors(N)))))$
Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes.
(%i10) for i:3 thru 22 do print([p:ROU(x,i),mp(p),RR6(cs(nGG(DA)))])$
[x^2+x+1,A^2+3,[[a1+A,A],[a1^2+3,a1]]]
[x^2+1,A^2+4,[[2*a1+A,A],[a1^2+1,a1]]]
[x^4+x^3+x^2+x+1,A^4+A^3+A^2+A+1,
[[(a2+a1+1)/4+A,A],[a2^2-2*a1+10,a2],[a1^2-5,a1]]]
[x^2-x+1,A^2+3,[[a1+A,A],[a1^2+3,a1]]]
[x^6+x^5+x^4+x^3+x^2+x+1,A^6+A^5+A^4+A^3+A^2+A+1,
[[A-(a1*(6*a2^2*w2+2*a2^2-7)-14*a2^2-14*a2-7)/42,A],[a1*(3*w2+2)+a2^3+7,a2],
[w2^2+w2+1,w2],[a1^2+7,a1]]]
[x^4+1,A^4+1,[[(a2+a1)/2+A,A],[a2^2+2,a2],[a1^2-2,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^4-x^3+x^2-x+1,A^4-A^3+A^2-A+1,
[[(a2+a1-1)/4+A,A],[a2^2+2*a1+10,a2],[a1^2-5,a1]]]
[x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1,
A^10+A^9+A^8+A^7+A^6+A^5+A^4+A^3+A^2+A+1,
[[(a1*((58*a2^4+198*a2^3-11*a2^2)*w2^3+((-322*a2^4)-33*a2^3-121*a2^2)*w2^2
+(240*a2^4+132*a2^3+55*a2^2)*w2
-290*a2^4+99*a2^3-110*a2^2+121)
+((-1892*a2^4)-374*a2^3-605*a2^2)*w2^3
+((-462*a2^4)-473*a2^3-121*a2^2)*w2^2+((-880*a2^4)+88*a2^3-363*a2^2)*w2
-1650*a2^4-671*a2^3-484*a2^2+242*a2+121)
/1210
+A,A],[a1*(40*w2^3+20*w2^2+30*w2+47)-55*w2^3+55*w2^2-55*w2+a2^5+77,a2],
[w2^4+w2^3+w2^2+w2+1,w2],[a1^2+11,a1]]]
[x^4-x^2+1,A^4-A^2+1,[[(a2+a1)/2+A,A],[a2^2+1,a2],[a1^2-3,a1]]]
[x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1,
A^12+A^11+A^10+A^9+A^8+A^7+A^6+A^5+A^4+A^3+A^2+A+1,
[[(a1*(a2*(4*a3^2*w3-25*a3^2)-26*a3^2*w3-52*a3^2+312)
+a2*((-26*a3^2*w3)+65*a3^2+312)+286*a3^2*w3+104*a3^2+624*a3+312)
/3744
+A,A],[a2*(18*w3+7)+a1*(a2*((-6*w3)-3)-8)+a3^3+52,a3],[w3^2+w3+1,w3],
[a2^2+6*a1+26,a2],[a1^2-13,a1]]]
[x^6-x^5+x^4-x^3+x^2-x+1,A^6-A^5+A^4-A^3+A^2-A+1,
[[A-(a1*(6*a2^2*w2+2*a2^2-7)+14*a2^2-14*a2+7)/42,A],[a1*(3*w2+2)+a2^3-7,a2],
[w2^2+w2+1,w2],[a1^2+7,a1]]]
[x^8-x^7+x^5-x^4+x^3-x+1,A^8-A^7+A^5-A^4+A^3-A+1,
[[(a3+a2+a1-1)/8+A,A],[a3^2+2*a2+a1*((-2*a2)-4)+28,a3],[a2^2-6*a1-30,a2],
[a1^2-5,a1]]]
[x^8+1,A^8+1,[[(a3+a2)/2+A,A],[a3^2+a1+2,a3],[a2^2+a1-2,a2],[a1^2-2,a1]]]
[x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1,
A^16+A^15+A^14+A^13+A^12+A^11+A^10+A^9+A^8+A^7+A^6+A^5+A^4+A^3+A^2+A+1,
[[(a4+a3+a2+a1+1)/16+A,A],[a4^2-2*a3+a1*(8-2*a3)+a2*((-2*a3)-8)+136,a4],
[a3^2+a1*(2*a2+12)-6*a2-68,a3],[a2^2-2*a1-34,a2],[a1^2-17,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^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2
+x+1,
A^18+A^17+A^16+A^15+A^14+A^13+A^12+A^11+A^10+A^9+A^8+A^7+A^6+A^5+A^4+A^3+A^2
+A+1,
[[(a1*(w2*(a2*(3672*a3^2*w3+2112*a3^2)+a2^2*(231*a3^2*w3+448*a3^2+22924944))
+a2*(2340*a3^2*w3+1784*a3^2)+a2^2*((-324*a3^2*w3)-293*a3^2+33113808)
+15504*a3^2*w3-10944*a3^2+774353664)
+w2*(a2^2*(399*a3^2*w3+266*a3^2-145191312)
+a2*((-7296*a3^2*w3)-17632*a3^2))
+a2*(4332*a3^2*w3-3496*a3^2+774353664)
+a2^2*((-342*a3^2*w3)+1501*a3^2-48397104)-55632*a3^2*w3-41344*a3^2
+20377728*a3+774353664)
/13938365952
+A,A],
[a1*(w2*(a2*(1481544*w3-493848)+a2^2*((-136458*w3)-227430))
+11852352*w3+a2^2*((-155952*w3)-12996)-493848*a2+7901568)
+w2*(a2*(7407720*w3+6420024)+a2^2*(1111158*w3+370386))
+a2*(2963088*w3-493848)+a2^2*((-740772*w3)-493848)+a3^3+75064896,a3],
[w3^2+w3+1,w3],[228*w2+a1*(16-36*w2)+a2^3+152,a2],[w2^2+w2+1,w2],
[a1^2+19,a1]]]
[x^8-x^6+x^4-x^2+1,A^8-A^6+A^4-A^2+1,
[[(a3+a2)/4+A,A],[a3^2+2*a1+6,a3],[a2^2+2*a1-10,a2],[a1^2-5,a1]]]
[x^12-x^11+x^9-x^8+x^6-x^4+x^3-x+1,A^12-A^11+A^9-A^8+A^6-A^4+A^3-A+1,
[[A-(a1*(3*a2*a3^2+8*a3^2-56)+a2*(7*a3^2-56)+28*a3^2-112*a3+56)/672,A],
[a3^3-7*a2+a1*(a2-12)+56,a3],[a2^2+2*a1+10,a2],[a1^2-21,a1]]]
[x^10-x^9+x^8-x^7+x^6-x^5+x^4-x^3+x^2-x+1,
A^10-A^9+A^8-A^7+A^6-A^5+A^4-A^3+A^2-A+1,
[[(a1*((58*a2^4-198*a2^3-11*a2^2)*w2^3+((-322*a2^4)+33*a2^3-121*a2^2)*w2^2
+(240*a2^4-132*a2^3+55*a2^2)*w2
-290*a2^4-99*a2^3-110*a2^2+121)
+(1892*a2^4-374*a2^3+605*a2^2)*w2^3+(462*a2^4-473*a2^3+121*a2^2)*w2^2
+(880*a2^4+88*a2^3+363*a2^2)*w2+1650*a2^4-671*a2^3+484*a2^2+242*a2-121)
/1210
+A,A],[a1*(40*w2^3+20*w2^2+30*w2+47)+55*w2^3-55*w2^2+55*w2+a2^5-77,a2],
[w2^4+w2^3+w2^2+w2+1,w2],[a1^2+11,a1]]]
Evaluation took 41.3960 seconds (41.4410 elapsed) using 10160.376 MB.
次は 11 次拡大を含むので時間が掛かります．
(%i12) print([p:ROU(x,23),mp(p),RR6(cs(nGG(DA)))])$
[x^22+x^21+x^20+x^19+x^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7
+x^6+x^5+x^4+x^3+x^2+x+1,
A^22+A^21+A^20+A^19+A^18+A^17+A^16+A^15+A^14+A^13+A^12+A^11+A^10+A^9+A^8+A^7
+A^6+A^5+A^4+A^3+A^2+A+1,
[[(a1*((22711467475534198*a2^10+847656939034980*a2^9+328864770731640*a2^8
+19604741278064*a2^7-10713350489280*a2^6
+899469127648*a2^5+655331361792*a2^4
+12250319616*a2^3-6160979456*a2^2)
*w2^9
+(39210376235165047*a2^10+2041710859169370*a2^9+733881801397984*a2^8
+14975939529168*a2^7-26094666163296*a2^6
+3085175093696*a2^5+1220699341568*a2^4
+77155521792*a2^3-22064903168*a2^2)
*w2^8
+(44258439323364695*a2^10+3203058077583282*a2^9+1086461046496968*a2^8
-12416794104704*a2^7-41260490775696*a2^6
+5863155824672*a2^5+1516563404544*a2^4
+174298407168*a2^3-42697020416*a2^2)
*w2^7
+(36252932377559731*a2^10+3962979058122594*a2^9+1274661086788456*a2^8
-53876459982320*a2^7-51395783099968*a2^6
+8351433082784*a2^5+1449040253312*a2^4
+272730799872*a2^3-61179958784*a2^2)
*w2^6
+(17735547288835748*a2^10+4080204258656724*a2^9+1238729739093416*a2^8
-96239906925832*a2^7-53282659141264*a2^6
+9759997523200*a2^5+1039568725888*a2^4
+341074688256*a2^3-72069131776*a2^2)
*w2^5
+((-5414577014723690*a2^10)+3517515506329206*a2^9+990074952435016*a2^8
-126057041160920*a2^7-46322048846112*a2^6
+9641627602944*a2^5+418118078976*a2^4
+357838283520*a2^3-71925853184*a2^2)
*w2^4
+((-25847439637650301*a2^10)+2453562504314970*a2^9+607642864832560*a2^8
-133861132247864*a2^7-32723893820416*a2^6
+8033919936032*a2^5-217984361344*a2^4
+317541179520*a2^3-60177008640*a2^2)
*w2^3
+((-37075751068922942*a2^10)+1226142813215268*a2^9+212852962075632*a2^8
-117174436025296*a2^7-16805512936224*a2^6
+5447303727776*a2^5-666801436032*a2^4
+233078449536*a2^3-41264234496*a2^2)
*w2^2
+((-35534601764595846*a2^10)+224953512486384*a2^9-68951751229648*a2^8
-81294861335224*a2^7-3620882635952*a2^6
+2703003978208*a2^5-785831683712*a2^4
+131314829568*a2^3-20775395840*a2^2)
*w2-21713295739522928*a2^10-232134869059746*a2^9-148300269904712*a2^8
-37613918803656*a2^7+2643970591392*a2^6+672332846240*a2^5
-537243326592*a2^4+44380543872*a2^3-5301307904*a2^2+6590815232)
+(23127481912300174*a2^10-4266420106368780*a2^9-1138621329910408*a2^8
+187360039390496*a2^7+56494249686592*a2^6
-12676273242976*a2^5-124724014336*a2^4
-485284591104*a2^3+95566820864*a2^2)
*w2^9
+(101470356548088363*a2^10-5657737299806442*a2^9-1243802641976576*a2^8
+395808948442240*a2^7+76242423677632*a2^6
-21008077937344*a2^5+1469393600256*a2^4
-861390895104*a2^3+158179565568*a2^2)
*w2^8
+(210155314754852863*a2^10-3732218201923602*a2^9-282149613567896*a2^8
+559165668119584*a2^7+52974608932720*a2^6
-22350123734048*a2^5+4276364496128*a2^4
-1009039484160*a2^3+171361196032*a2^2)
*w2^7
+(314675650233351123*a2^10+898798478208330*a2^9+1441019713076728*a2^8
+625565597409600*a2^7-5921818285152*a2^6
-16276263112800*a2^5+7405103289984*a2^4
-882237930240*a2^3+128520897024*a2^2)
*w2^6
+(381846894974485910*a2^10+6764997665935512*a2^9+3378611252436072*a2^8
+573927223073160*a2^7-81747656536464*a2^6
-4714980758016*a2^5+9862008765952*a2^4
-520316206848*a2^3+42840299008*a2^2)
*w2^5
+(390342653382438482*a2^10+12003902568405402*a2^9+4915453381952560*a2^8
+420645366298040*a2^7-150428730646176*a2^6
+8663142278528*a2^5+10867108088832*a2^4
-38327023360*a2^3-56021929472*a2^2)
*w2^4
+(337465582190696883*a2^10+14952197895233958*a2^9+5563609584583344*a2^8
+214385938621864*a2^7-190159286564736*a2^6
+19610665087776*a2^5+10101319834240*a2^4
+410385707136*a2^3-138407119872*a2^2)
*w2^3
+(240003777890885208*a2^10+14673820711808060*a2^9+5117294845027984*a2^8
+20634845593008*a2^7-188325170689696*a2^6
+24651756486752*a2^5+7807716133504*a2^4
+683976178560*a2^3-177952011264*a2^2)
*w2^2
+(128900674716232452*a2^10+11257153806707688*a2^9+3718210938875040*a2^8
-99093307914712*a2^7-145508697393104*a2^6
+22185977861024*a2^5+4714689528704*a2^4
+694829531904*a2^3-161474973184*a2^2)
*w2+39430722907678478*a2^10+5786964780828642*a2^9+1810557120845464*a2^8
-106785679762504*a2^7-75303787322272*a2^6+12996173068448*a2^5
+1804128210816*a2^4+439614539904*a2^3-95566820864*a2^2+6590815232*a2
+6590815232)
/144997935104
+A,A],
[a1*(981871616*w2^9+573213696*w2^8-345849856*w2^7+1610481664*w2^6
+2609913856*w2^5+2784956416*w2^4+3007206400*w2^3
+863858688*w2^2+549716992*w2+958533632)
+3252648960*w2^9+13344021504*w2^8+9279959040*w2^7+3458611200*w2^6
+6119798784*w2^5+2020761600*w2^4+3886339072*w2^3+1821276160*w2^2
-6719032320*w2+a2^11-2574940160,a2],
[w2^10+w2^9+w2^8+w2^7+w2^6+w2^5+w2^4+w2^3+w2^2+w2+1,w2],[a1^2+23,a1]]]
Evaluation took 691.5360 seconds (692.1830 elapsed) using 92077.766 MB.
その後は，．．．
(%i13) for i:24 thru 28 do print([p:ROU(x,i),mp(p),RR6(cs(nGG(DA)))])$
[x^8-x^4+1,A^8-A^4+1,
[[(a3+a2+a1)/4+A,A],[a3^2-2*a1*a2+8,a3],[a2^2-2,a2],[a1^2-6,a1]]]
[x^20+x^15+x^10+x^5+1,A^20+A^15+A^10+A^5+1,
[[a3/2+A,A],[a3^5-8*a2-8*a1-8,a3],[a2^2-2*a1+10,a2],[a1^2-5,a1]]]
[x^12-x^11+x^10-x^9+x^8-x^7+x^6-x^5+x^4-x^3+x^2-x+1,
A^12-A^11+A^10-A^9+A^8-A^7+A^6-A^5+A^4-A^3+A^2-A+1,
[[A-(a1*(a2*(4*a3^2*w3-25*a3^2)+26*a3^2*w3+52*a3^2-312)
+a2*(26*a3^2*w3-65*a3^2-312)+286*a3^2*w3+104*a3^2-624*a3+312)
/3744,A],[a1*(a2*(6*w3+3)-8)+a2*(18*w3+7)+a3^3-52,a3],[w3^2+w3+1,w3],
[a2^2-6*a1+26,a2],[a1^2-13,a1]]]
[x^18+x^9+1,A^18+A^9+1,
[[a3/2+A,A],[a3^3-4*a2,a3],[a2^3-4*a1-4,a2],[a1^2+3,a1]]]
[x^12-x^10+x^8-x^6+x^4-x^2+1,A^12-A^10+A^8-A^6+A^4-A^2+1,
[[(a1*(6*a3^2*w3+2*a3^2+7)+a2*(7-14*a3^2)+14*a3)/42+A,A],
[a1*((-3*w3)-2)+a3^3-7*a2,a3],[w3^2+w3+1,w3],[a2^2+1,a2],[a1^2-7,a1]]]
Evaluation took 77.5920 seconds (77.6710 elapsed) using 9047.339 MB.
