円分体の利用

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