load(gcdex)$  load(grobner)$ load(sym)$ load("sym/compile")$
	rtL:[α,β,γ,δ,ε,λ,ρ];
	cL:[1,2,3,4,5,6,7];
	/* 多項式の最高次数と対称式としての値を返す面白い命令　polynome2ele */
	fx:x^3+3*x+1;
	N:hipow(fx,x); NN:N!;
	/*　
	3次方程式 fx=0の根を ( α, β, γ　)  とするとき
	(α,β,γ)の全ての置換の組をリスト化する命令
	
	permutations([α,β,γ])
	
	*/
	rtL:rest(rtL,N-7);  cL:rest(cL,N-7); 
	/*  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);
	vxL:makelist(x-v[i],i,1,NN);
	Vx:apply("*",vxL);
	Vrt:ev(Vx,virtL);  Vrt:expand(Vrt)$
	Vrt;
	fx1:divide(fx,x-α,x);
	fx2:divide(fx1[1],x-β,x);
	fx3:divide(fx2[1],x-γ,x);
	/*
	ここからの計算が秀逸
	Vrtを fx3[2]->fx2[2]->fx1[2]と順番に
	剰余を取ると自然に[α,β,γ]が消去されて
	基礎体での多項式となる
	*/
	Vrt:remainder(Vrt,fx3[2],γ);
	Vrt:remainder(Vrt,fx2[2],β);
	Vrt:remainder(Vrt,fx1[2],α);
	Vx:Vrt;
	Vx:factor(Vx);
	Vpw:hipow(Vx,x)$
	if Vpw=NN then gx[0]:Vx else gx[0]:part(Vx,1)$
	gx[0]; 
	gv[0]:subst(v,x,gx[0]);
	/*　
	
	以上でvの最小多項式　gx[0]が求まった　
	次に[α,β,γ] と　v[1]...v[6]をvの
	多項式表現として求める計算
	
	*/
	pev:[v=α+2*β+3*γ]$
	pvL:[]$
	for i:0 thru (NN-1) do(
	    z:v^i, z:ev(z,pev),  z:remainder(z,fx3[2],γ), z:remainder(z,fx2[2],β), z:remainder(z,fx1[2],α),
	    z:expand(z), pvL:endcons(z,pvL)
	)$ 
	pvL;
	BV:[1,β,α,α*β,α^2,(α^2)*β]$
	Mab:[]$
	for i:1 thru NN do(
	    z:pvL[i], ab:[], 
	for j:0 thru 2 do(     
	           zj:coeff(z,α,j),
	for k:0 thru 1 do(   
	            zjk:coeff(zj,β,k),       
	            ab:endcons(zjk,ab)
	         )
	    ),
	     Mab:endcons(ab,Mab)
	)$
	Mab;
	Mav:apply('matrix, Mab);
	Nav:invert(Mav);
	Nav18:18*Nav;
	VL:[1,v,v^2,v^3,v^4,v^5]$
	Ansabc:Nav.VL;
	SoL:[α=Ansabc[3][1],β=Ansabc[2][1],γ=-Ansabc[3][1]-Ansabc[2][1]]; 
	VivL:expand(ev(virtL,SoL));
	check:[]$ GgrL:[]$
	for i:1 thru NN do(
	    z:remainder(subst(rhs(VivL[i]),x,Vx),gv[0]),
	    check:endcons(z,check), if z=0 then GgrL:endcons(i,GgrL)
	                          )$
	check; GgrL;
	/* 
	これで準備は終わった　RT2はここまで
	
	以下はRT4の計算　　
	*/
	/*  RT4の第１ステップ　
	
	        6次式のg[0]を分解して3次のg[1]を計算する
	        ガロア群はGal0= [1,2,3,4,5,6]からGal1=[1,4,5]に縮小
	　　　　　Gal0/Gail1=2　の計算である　　　　　　　　　　　　　　　　
	*/
	h0:(x-v[1])*(x-v[4])*(x-v[5]);
	h0:ev(h0,VivL);
	h0:expand(h0); 
	h0:remainder(h0,gv[0]);
	h1:(x-v[2])*(x-v[3])*(x-v[6]);
	h1:ev(h1,VivL);
	h1:expand(h1);
	h1:remainder(h1,gv[0]);
	t0:(h0+h1)/2;
	t1:(h0-h1)/2;
	T1:expand(t1^2);
	T1:remainder(T1,gv[0]);
	A[1]:T1;  B[1]:a[1]^2-A[1];   t1:a[1]$
	h0:expand(t0+t1);
	gx[1]:h0; gv[1]:subst(v,x,h0);
	/* 以下は念のためのチェック　*/
	h1:t0-t1;
	h01:expand(h0*h1);  h01:remainder(h01,B[1],a[1]);
	/*   第2ステップ　
	        3次多項式のg[1]をさらに分解してゆき1次多項式のg[2]を求める
	        ガロア群はGal0= [1,2,3,4,5,6]からGal1=[1,4,5]になり
	　　　　以下の計算は　Gal1=[1,4,5]から　Gal2=[1] への最後の
	　　　　計算である
	        Gal1/Gal2=3  である
	　*/
	Ω:ω^2+ω+1$
	h0:ev(x-v[1],VivL);  
	h1:ev(x-v[4],VivL)$ h1:remainder(h1,gv[1])$ h1:expand(h1);
	h2:ev(x-v[5],VivL)$ h2:remainder(h2,gv[1])$  h2:expand(h2);
	
	t0:expand((h0+h1+h2)/3)$ 
	      t0:remainder(t0,g[1]); t0:remainder(t0,B[1]); t0:remainder(t0,Ω);
	
	t1:(h0+ω*h1+ω^2*h2)/3$  
	        t1:remainder(t1,g[1]); t1:remainder(t1,B[1]); 
	        t1:remainder(t1,Ω);  t1:expand(t1);
	
	t2:(h0+ω^2*h1+ω*h2)/3$  
	        t2:remainder(t2,g[1]); t2:remainder(t2,B[1]);
	        t2:remainder(t2,Ω);   t2:expand(t2);
	
	/* T1=t1^3 を求める計算*/
	T1:expand(t1^3);  T1:remainder(T1,gv[1]); T1:remainder(T1,Ω);   
	  T1:remainder(T1,B[1]); T1:expand(T1); 
	t1; t2; T12:t1*t2; T12:expand(T12); T12:remainder(T12,gv[1]); T12:remainder(T12,Ω);   
	  T12:remainder(T12,B[1]); T12:expand(T12);
	t1; t2; T12:t1*t2; T12:expand(T12); T12:remainder(T12,gv[1]);  
	  T12:remainder(T12,B[1]); 
	T12:remainder(T12,Ω);  
	T12:expand(T12);
	/* ここまでのまとめが以下の式*/
	A[2]:T1;  B[2]:a[2]^3-A[2];   t1:a[2]$
	/*　
	T1の逆数を求める計算 
	IAと言う逆元多項式を仮定して係数を決定する方式　個の逆数を求めるやり方は重要
	*/
	IA:d[0]+d[1]*ω+d[2]*a[1]+d[3]*ω*a[1];
	W:IA*A[2];   W:expand(W);   W:remainder(W,Ω);   
	 W:remainder(W,B[1]);   W:expand(W);
	D[0]:coeff(W,ω,0)$ D[0]:coeff(D[0],a[1],0);
	D[1]:coeff(W,ω,1)$ D[1]:coeff(D[1],a[1],0);
	D[2]:coeff(W,ω,0)$ D[2]:coeff(D[2],a[1],1);
	D[3]:coeff(W,ω,1)$ D[3]:coeff(D[3],a[1],1);
	Dsol:solve([D[0]=1,D[1]=0,D[2]=0,D[3]=0],[d[0],d[1],d[2],d[3]]);
	IA:ev(IA,Dsol[1]);
	t2:-3*a[2]^2*IA; t2:expand(t2); t2:factor(t2);
	t0; t1; t2; h0:t0+t1+t2;
	A[1]; B[1];  A[2]; B[2];
	h0:t0+t1+t2$
	gx[2]:h0; gv[2]:subst(v,x,h0);
	gv[0];  gv[1]; gv[2];  B[1]; B[2]; Ω; 
	W:solve(gx[2],x);  W:expand(W);
	/* 以下はg[1]が以下に分解されたかのチェック*/
	h1:t0+t1*ω^2+t2*ω;  h2:t0+t1*ω+t2*ω^2;
	h012:expand(h0*h1*h2)$
	h012:remainder(h012,B[2],a[2])$  h012:remainder(h012,Ω,ω)$
	h012:remainder(h012,B[1],a[1]);
	
	
	SoL; Rt:subst(v,x,W);
	/*
	以下はgx[2]=0で求まるvの値を、vの多項式で表現された[α,β,γ]に
	代入して、添加数a[1],a[2]で3根を表す計算
	*/
	SoL:ev(SoL,Rt);
	AnsL:[]$
	for i:1 thru 3 do(
	        z:rhs(SoL[i]),
	        z:remainder(z,B[2],a[2]), 
	        z:remainder(z,Ω,ω),
	        z:remainder(z,B[1],a[1]), 
	        AnsL:endcons(z,AnsL)
	)$
	AnsL$
	/*　　以下の式が添加数a[1],a[2],ωで記述されたf(x)の3根の数式である　　*/
	AnsL:map( lambda([x,y],x=y),rtL,AnsL);
	/*　実際に数値的に正しいか計算してみる。　maximaで数値で解を求める命令はallroots()である*/
	root:[]$
	ωr:allroots(Ω)$  root:endcons(ω=rhs(ωr[1]),root);
	a1r:allroots(B[1])$  root:endcons(a[1]=rhs(a1r[1]),root);
	b2:ev(B[2],root)$  b2r:allroots(b2)$  root:endcons(a[2]=rhs(b2r[1]),root);
	
	FansL:ev(AnsL,root)$  FansL:expand(FansL);
	allroots(fx);
	/*FansLとallroota(fx)の実数値で一致している事が判ったので、AnsLの数式は正しい事が判った*/