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