load(sym)$ load("sym/compile")$ load(gcdex)$ /* 以下は逆元を求めるサブルーチンプログラム 使い方 A: 逆元を求めたい元 Const: 代数体を規定しているリスト 例えば　Const:[ [v,gv[0]],[a[1],B[1]]] など gv[0]: 最小多項式　変数v B[1]: 添加される数a[1]が満たす冪根方程式 具体的には　Const:[[v,v^8-v^6+v^4-v^2+1],[a[1],a[1]^2-5/4]]$ */ Inverse(A,Const):=block( [Nconst,vn,Tn,rwn,rpn], Nconst:length(Const), vn:Nconst, varL:[], maxpL:[], Tn:1, for i:1 thru Nconst do( con:Const[i], y[i]:con[1], varL:endcons(y[i],varL), D[i]:con[2], pw[i]:hipow(D[i],y[i]), maxpL:endcons(pw[i],maxpL),Tn:Tn*pw[i] ), /* 　基底の次数の組み合わせ Tn: pw[1]*pw[2]*pw[3] 2*3*2=12 pwL: 基底次数の組み合わせのリスト [[0,0,0],[0,0,1],[0,1,0],[0,1,1],[0,2,0],[0,2,1], 　[1,0,0],[1,0,1],[1,1,0],[1,1,1],[1,2,0],[1,2,1]] rwn: 繰り返し同じ数を書く数(re-wright number) rpn: (repetition number) 例えば上の例ではリスト[x1,x2,x3]のx3の例でいうとx2だけをみると 　　　　　 [0,0,1,1,2,2,0,0,1,1,2,2]となっているが、この [0,0,1,1,2,2]の繰り返しが2回繰り返されている。 この数を指定する数 Tn pw[i] rwn rpr 12　÷ 2 = 6 1=Tn/(pw[1]*rwn) 6　 ÷ 3 = 2 2=Tn/(pw[2]*rwn) 2　÷ 2 = 1 6=Tn/(pw[3]*rwn) という関係式があるので下記のプログラムとなる 　*/ pwL:[], for i:1 thru Tn do(pwL:endcons([],pwL)), rwn:Tn, for i:1 thru Nconst do( rwn:rwn/pw[i], rpn:Tn/(pw[i]*rwn), n:1, for j:1 thru rpn do( for k:1 thru pw[i] do( for m:1 thru rwn do ( pwL[n]:endcons(k-1,pwL[n]), n:n+1 ) ) ) ), baseL:[], for i:1 thru Tn do( b:1, for j:1 thru vn do( b:b*varL[j]^pwL[i][j] ), baseL:endcons(b,baseL) ), coefL:makelist(c[i],i,0,Tn-1), F:sum(coefL[i]*baseL[i],i,1,Tn), W:F*A, for i:1 thru Nconst do( W:remainder(W,Const[i][2],varL[i]) ), W:expand(W), CeL:[], for i:1 thru Tn do( ce:W, for j:1 thru vn do( ce:coeff(ce,varL[j],pwL[i][j]) ), CeL:endcons(ce,CeL) ), CeL[1]:CeL[1]-1, SLV:solve(CeL,coefL), IA:sum(coefL[i]*baseL[i],i,1,Tn), IA:ev(IA,SLV) )$ /*　 以下は、代数体の次数を低減させるサブルーチンプログラム Const: 代数体を規定しているリスト 例えば　Const:[ [v,gv],[a[1],B[1]]] など gv[i]: 現在の体で使用されている v　の最小多項式 B[1]: 添加される数a[1]が満たす冪根方程式 具体的には　Const:[[v,v^8-v^6+v^4-v^2+1],[a[1],a[1]^2-5/4]]$ 　*/ Reduce(A,Const):=block( crn:length(Const), for i:1 thru crn do( A:remainder(A,Const[i][2],Const[i][1]) ), A )$ /* ------------------------------------------------------------------------------ */ rtL:[α,β,γ,δ,ε,λ,ρ]; cL:[1,2,3,4,5,6,7]; /* 多項式の最高次数と対称式としての値を返す面白い命令　polynome2ele */ fx:x^3-3*x+1; N:hipow(fx,x); NN:N!; rtL:rest(rtL,N-7); cL:rest(cL,N-7); /* v[i]の計算の時に必要*/ /* v[i]の計算の時に必要*/ V(cL,rtL):=sum(cL[i]*rtL[i],i,1,N); V(cL,rtL); prtL:listify( permutations(rtL)); viL:makelist(v[i],i,1,NN); pVrtL:map(lambda([z],V(cL,z)),prtL); virtL:map( lambda([x,y],x=y),viL,pVrtL); /* 第1ステップ　 fxを(x-α), (x-β), (x-γ)で順番に割り算してゆく。 divideという命令を使う その余りが全てゼロと言う条件を課すと良い。 fx/(x-a)の余り　　　　　　　　　　　　　α^3-3*α-1 その商を(x-b)で割ったときの余り　　　　β^2+α*β+α^2-3 更にその商を(x-c)で割ったときの余り　　γ+β+α を使う */ fx1:divide(fx,x-α,x); fx2:divide(fx1[1],x-β,x); fx3:divide(fx2[1],x-γ,x); /*第2ステップ　其の二 終結式を使うのだが　eliminateでなくて resultantを使った方が使いやすい */ v0:V(cL,rtL); r1:resultant(v-v0,fx3[2],γ); r2:resultant(r1,fx2[2],β); Gv:resultant(r2,fx1[2],α); /* vの最小多項式の候補Gvは求まったが 既約かどうか判定する 既約でないので、因数分解されたGvの第一因子を vの最小多項式とする */ Gv:factor(Gv); Gpw:hipow(Gv,v)$ if Gpw=NN then gv[0]:Gv else gv[0]:part(Gv,1)$ gv[0]; gx[0]:subst(x,v,gv[0]); /*以上でvの最小多項式gx[0]が決定した*/ fxg0:factor(fx,gv[0]); SoL:solve(fxg0,x)$ SoL:expand(SoL); an:[]$ for i:1 thru 3 do( an:endcons(rhs(SoL[i]),an))$ an; kiL:makelist(k[i],i,1,N); sokiL:map(lambda([y,x],y=x),kiL,an); pkiL:listify(permutations(kiL)); VpkiL:map(lambda([z],V(cL,z)),pkiL); /* さて代数体上での因数分解で求めた[α,β,γ]の候補 SoLの中のリストをv=α+2β+3γを満たす様に 対応付けなければならない。以下はその計算 */ fix:0$ for i:1 thru NN do( vkev:expand(ev(VpkiL[i],sokiL)), if vkev=v then fix:i)$ fix; fixki:map( lambda([x,y],x=y),rtL,pkiL[fix]); SoL:ev(fixki,sokiL); VivL:expand(ev(virtL,SoL)); /* chekL:[]$ for i:1 thru 6 do( rot:[], gvx:subst(v[i],x,gx[0]), gvx1:expand(ev(gvx,VivL)), gvx2:remainder(gvx1,gv[0]), rot:endcons(gvx1,rot), rot:endcons(gvx2,rot), chekL:endcons(rot,chekL) )$ for i:1 thru 6 do(print(i,"=",chekL[i])); */ check:[]$ GgrL:[]$ for i:1 thru NN do( z:remainder(subst(rhs(VivL[i]),v,gv[0]),gv[0]), check:endcons(z,check), if z=0 then GgrL:endcons(i,GgrL) )$ check; GgrL; /* 第2ステップ　 3次多項式のGvが因数分解がされ3次の多項式がg[0]となった。 　　　　　従って ガロア群はGal0= [1,4,5]　のA3がガロア群G0と 　　　　　なることが分かった　　　　 　　　　　以下の計算は　Gal0=[1,4,5]から　Gal1=[1] への 　　　　　　計算である Gal0/Gal1=3 である 　*/ Ω:ω^2+ω+1$ Const:[]$ Const:cons([ω,Ω],Const); CurConst:cons([v,gv[0]],Const); h0:ev(x-v[1],VivL); h1:ev(x-v[4],VivL); h2:ev(x-v[5],VivL); t0:expand((h0+h1+h2)/3)$ t0:Reduce(t0,CurConst); t1:expand((h0+ω*h1+ω^2*h2)/3)$ t1:Reduce(t1,CurConst); t2:expand((h0+ω^2*h1+ω*h2)/3)$ t2:Reduce(t2,CurConst); T1:expand(t1^3)$ T1:Reduce(T1,CurConst); T2:expand(t2^3)$ T2:Reduce(T2,CurConst); t12:t1*t2$ t12:expand(t12); t12:Reduce(t12,CurConst); iT1:Inverse(T1,Const); B[1]:T1-a[1]^3; t1:a[1]$ a2:(a[1]^2)*t12*iT1$ a2:rat(a2); t2:a2$ h0:t0+t1+t2; gx[1]:h0$ gx[1]:expand(gx[1]); gv[1]:subst(v,x,gx[1]); gx[0]; gx[1]; B[1]; Const:cons([a[1],B[1]],Const); CurConst:cons([v,gv[1]],Const); W:solve(gx[1],x); /* 以下はg[1]が以下に分解されたかのチェック*/ h1:t0+t1*ω^2+t2*ω; h2:t0+t1*ω+t2*ω^2; h012:expand(h0*h1*h2); h012:Reduce(h012,CurConst); /* 以上で　 g[0]=h0*h1:h2　と因数分解できていることが分かった　*/ SoL; SoL:expand(SoL); Rt:subst(v,x,W)$ Rt:expand(Rt); AnsL:[]$ for i:1 thru 3 do( z:subst(rhs(Rt[1]),v,rhs(SoL[i])), z:Reduce(z,CurConst), AnsL:endcons(z,AnsL) )$ AnsL:expand(AnsL); AnsL:ratsimp(AnsL); /*　実際に数値的に正しいか計算してみる。　maximaで数値で解を求める命令はallroots()である*/ root:[]$ ωr:allroots(Ω)$ root:endcons(ω=rhs(ωr[1]),root); b1:ev(B[1],root)$ b1r:allroots(b1)$ root:endcons(a[1]=rhs(b1r[1]),root); FansL:ev(AnsL,root)$ FansL:expand(FansL); allroots(fx); /*FansLとallroota(fx)の実数値で一致している事が判ったので、AnsLの数式は正しい事が判った*/