load(sym)$
	load("sym/compile")$
	load(gcdex)$
	rtL:[α,β,γ,δ,ε,λ,ρ];
	cL:[1,2,3,4,5,6,7];
	eiL:[e1,e2,e3,e4,e5,e6,e7];
	/* 多項式の最高次数と対称式としての値を返す面白い命令　polynome2ele */
	fx:x^3+3*x+1;
	fxinfoL:polynome2ele(fx,x);
	N:fxinfoL[1]; NN:N!;
	feL:rest(fxinfoL,1); 
	eiL:rest(eiL,N-7);
	fxeL:map( lambda([x,y],x=y),eiL,feL);
	/*本当は eiL:makelist(e[i],i,1,Nfpwr);   としたい
	
	　のだが e[i]のようなサフィックスだと根[α,β,γ]の対称式から方程式の
	係数に変換すると他の関数　" tcontract(...),  elem(...)  "
	の出力がe1,e2,e3を使っておりと一致しないので、やめにした　　*/
	/*　
	3次方程式 fx=0の根を ( α, β, γ　)  とするとき
	(α,β,γ)の全ての置換の組をリスト化する命令
	
	permutations([α,β,γ])
	
	*/
	rtL:rest(rtL,N-7);  cL:rest(cL,N-7);  zL:[u,w,y];
	/*  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:makelist(V(cL,prtL[i]),i,1,NN);
	virtL:map( lambda([x,y],x=y),viL,pVrtL);
	awL:[1,0,0]; bwL:[0,1,0]; cwL:[0,0,1];
	pawL:makelist(V(awL,prtL[i]),i,1,NN); 
	pbwL:makelist(V(bwL,prtL[i]),i,1,NN); 
	pcwL:makelist(V(cwL,prtL[i]),i,1,NN);
	vxL:makelist(x-v[i],i,1,NN);
	awxL:makelist(pawL[i]/(x-v[i]),i,1,NN);
	bwxL:makelist(pbwL[i]/(x-v[i]),i,1,NN);
	cwxL:makelist(pcwL[i]/(x-v[i]),i,1,NN);
	Vx:apply("*",vxL);
	/*　α, β, γ  の計算　の準備　*/
	Pax:0$
	for i:1 thru NN do( Pax:Pax+Vx*awxL[i])$
	Pax;
	Pbx:0$
	for i:1 thru NN do( Pbx:Pbx+Vx*bwxL[i])$
	Pbx$
	Pcx:0$
	for i:1 thru NN do( Pcx:Pcx+Vx*cwxL[i])$
	Pcx$
	Vrt:ev(Vx,virtL);  Vrt:expand(Vrt)$
	Part:ev(Pax,virtL);  Pbrt:ev(Pbx,virtL)$  Pcrt:ev(Pcx,virtL)$
	for i:0 thru NN do(c[i]:coeff(Vrt,x,i), 
	                            ct[i]:tcontract(c[i],[α,β,γ]), 
	                            ce[i]:elem([3],ct[i],[α,β,γ]),
	                            ce[i]:expand(ce[i])
	                            )$
	for i:0 thru NN do (print(NN-i," ",c[NN-i]));
	for i:0 thru NN do (print(NN-i," ",ce[NN-i]));
	Vei:0$
	for i:0 thru NN do (Vei:Vei+ce[i]*x^i)$  Vei:expand(Vei);
	Vx:ev(Vei,fxeL);
	/*  Vxは最小多項式gx[0]の候補となる。 Vx が因数分解されたかどうかチェックする。
	因数分解されないとき　gx[0]=Vx　 とする
	Vxが因数分解されたとき　　gx[0]=part(Vx,1)　として因数分解の第１因子をgx[0]とする
	以上の判断をして以下の議論をする。
	　　　参考：式" f "の1番目(n=1)の因数を取り出す命令　part(f,n)
	
	hipow(Vx,x)
	と言うコマンドも面白い
	
	*/
	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]);
	/*  P(x)の計算　*/
	/*　α の計算　　最終的にαをvの多項式で表現する計算*/
	Pw:expand(Part)$
	for i:0 thru (NN-1) do(c[i]:coeff(Pw,x,i), 
	                            ct[i]:tcontract(c[i],[α,β,γ]), 
	                            ce[i]:elem([3],ct[i],[α,β,γ]),
	                            ce[i]:expand(ce[i])
	                            )$
	for i:0 thru (NN-1) do (print(NN-1-i," ",c[NN-1-i]));
	for i:0 thru (NN-1) do (print(NN-1-i," ",ce[NN-1-i]));
	Pei:0$
	for i:0 thru (NN-1) do (Pei:Pei+ce[i]*x^i)$  Pei:expand(Pei);
	Pax:ev(Pei,fxeL);
	/*  βの計算　*/
	Pw:expand(Pbrt)$
	for i:0 thru (NN-1) do(c[i]:coeff(Pw,x,i), 
	                            ct[i]:tcontract(c[i],[α,β,γ]), 
	                            ce[i]:elem([3],ct[i],[α,β,γ]),
	                            ce[i]:expand(ce[i])
	                            )$
	for i:0 thru (NN-1) do (print(NN-1-i," ",ce[NN-1-i]))$
	Pei:0$
	for i:0 thru (NN-1) do (Pei:Pei+ce[i]*x^i)$  Pei:expand(Pei)$
	Pbx:ev(Pei,fxeL);
	/*  γの計算　*/
	Pw:expand(Pcrt)$
	for i:0 thru (NN-1) do(c[i]:coeff(Pw,x,i), 
	                            ct[i]:tcontract(c[i],[α,β,γ]), 
	                            ce[i]:elem([3],ct[i],[α,β,γ]),
	                            ce[i]:expand(ce[i])
	                            )$
	for i:0 thru (NN-1) do (print(NN-1-i," ",ce[NN-1-i]))$
	Pei:0$
	for i:0 thru (NN-1) do (Pei:Pei+ce[i]*x^i)$  Pei:expand(Pei)$
	Pcx:ev(Pei,fxeL);
	/* 
	
	Vxの微分した dV の逆元　IdV を計算する
	これらを使いα,β,γを計算する
	
	逆元を求める際には、ここではmaximaが持っているgcdexという命令を
	使ってみる
	
	　*/
	dV:diff(Vx,x);
	gc:gcdex(Vx,dV,x);
	IdV:gc[2]$  check:remainder(dV*IdV,Vx);
	/*
	上記check=1となったのでIdVがdVの逆元である事が
	確認された
	*/
	/*
	逆元がvの多項式としても止まったので、α,β,γ の計算が出来るようになった
	*/
	SoL:[]$
	αx:remainder(Pax*IdV,Vx)$   αv:subst(v,x,αx)$ SoL:endcons(α=αv,SoL)$
	βx:remainder(Pbx*IdV,Vx)$   βv:subst(v,x,βx)$ SoL:endcons(β=βv,SoL)$
	γx:remainder(Pcx*IdV,Vx)$   γv:subst(v,x,γx)$ SoL:endcons(γ=γv,SoL)$
	SoL;
	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;
	/* 第１ステップ　
	
	        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の数式は正しい事が判った*/