GGRR on SageMath - ATPとCASのこと

今回の「可解な方程式の全ての根を 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()