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