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];
	zz:[z1,z2,z3,z4,z5,z6,z7];
	/* 多項式の最高次数と対称式としての値を返す面白い命令　polynome2ele */
	fx:x^4+4*x+2;
	N:hipow(fx,x); NN:N!;
	rtL:rest(rtL,N-7);  cL:rest(cL,N-7);  zz:rest(zz,N-7);
	/*  v[i]の計算の時に必要*/
	V(x,y):=sum(x[i]*y[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-γ),(x-δ)で順番に割り算してゆく。
	divideという命令を使う
	
	その余りが全てゼロと言う条件を課すと良い。
	
	 */
	fx1:divide(fx,x-α,x);
	fx2:divide(fx1[1],x-β,x);
	fx3:divide(fx2[1],x-γ,x);
	fx4:divide(fx3[1],x-δ,x);
	/*第2ステップ　其の二
	終結式を使うのだが　eliminateでなくて
	resultantを使った方が使いやすい
	*/
	v0:V(cL,rtL);
	r1:resultant(v-v0,fx4[2],δ);
	r2:resultant(r1,fx3[2],γ);
	r3:resultant(r2,fx2[2],β);
	Gv:resultant(r3,fx1[2],α);
	/* 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]);
	/*　　第3ステップ　*/
	/*ここ重要　拘束条件のリスト作成*/
	
	Const:[]$
	CurConst:[]$
	CurConst:cons([v,gv[0]],Const);
	
	/*現時点での単拡大体F(v)上での拘束条件作成は以下の通り*/
	/*
	gv{0}が生成する代数体上でfxを因数分解して
	fxの根[α,β,γ,δ]の候補を計算する
	*/
	fxg0:factor(fx,gv[0])$
	SoL:solve(fxg0,x);
	/*
	以上でfxの根[α,β,γ,δ]の候補が計算された。少し煩雑な式ですが
	
	以下は、v=α+2β+3γ+4δ を満足する順列を探す 
	
	*/
	an:[]$
	for i:1 thru N 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)$
	/*HTML出力用*/
	VpkiL[1]; VpkiL[10]; VpkiL[24];
	/*HTML出力完了*/
	vkev:expand(ev(VpkiL[1],sokiL));
	vkev:expand(ev(VpkiL[10],sokiL));
	vkev:expand(ev(VpkiL[24],sokiL));
	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)$ SoL:expand(SoL);
	VivL:expand(ev(virtL,SoL))$
	/*HTML　用出力*/
	expand(VivL[1]);  expand(VivL[2]); expand(VivL[23]);   expand(VivL[24]);
	/*HTML　終了*/
	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;
	/*
	これで準備が整った。後は組成列に従って巡回拡大に分解して
	gx[0] - > gx{1} -. gx[2] -> gx[3] -> gx[4]　と
	最小多項式の次数を提言してゆく計算
	*/
	/*   第４ステップ　 P=2　*/
	Gal11:[1,4,5,8,9,12,13,16,17,20,21,24]$
	Gal12:[2,3,6,7,10,11,14,15,18,19,22,23]$
	N1:length(Gal11);
	vx1L:makelist(x-v[Gal11[i]],i,1,N1);
	vx2L:makelist(x-v[Gal12[i]],i,1,N1);
	h0:apply("*",vx1L);
	h1:apply("*",vx2L);
	h0:ev(h0,VivL)$   h0:Reduce(h0,CurConst);
	h1:ev(h1,VivL)$   h1:Reduce(h1,CurConst);
	t0:ratexpand((h0+h1)/2)$  t0:Reduce(t0,CurConst)$
	t1:ratexpand((h0-h1)/2)$  t1:Reduce(t1,CurConst)$
	ratvars(x)$ t0; t1:ratsimp(t1);
	/* t1:expand(t1)$ t1:ratsimp(t1); */
	t1:expand(t1)$
	cv:coeff(t1,x,hipow(t1,x));
	cv2:cv^2$ cv2:Reduce(cv2,CurConst);
	A[1]:cv2$  B[1]:a[1]^2-A[1];
	icv:Inverse(cv,CurConst);
	q1:Reduce(icv*t1,CurConst)$ q1:expand(q1);
	t1:a[1]*q1;
	h0:t0+t1$  gx[1]:h0;  gv[1]:subst(v,x,h0);
	/*
	新たな制約条件B[1]と新たな最小多項式gv[1]が生成されたので、
	Const  CurConst　を後進する必要があるので以下の命令をする
	*/
	
	cn:[a[1],B[1]]$  Const:cons(cn,Const);
	Ω:ω^2+ω+1$
	cn:[ω,Ω]$  Const:cons(cn,Const)$
	CurConst:cons([v,gv[1]],Const);
	
	/*拘束条件作成は終*/
	
	
	/*   第５ステップ　　P=3  */
	/*   Gal11:[1,4,5,8,9,12,13,16,17,20,21,24]$  */
	Gal21:[1,8,17,24]$
	Gal22:[4,12,13,21]$
	Gal23:[5,9,16,20]$
	N2:length(Gal21);
	vx1L:makelist(x-v[Gal21[i]],i,1,N2);
	vx2L:makelist(x-v[Gal22[i]],i,1,N2);
	vx3L:makelist(x-v[Gal23[i]],i,1,N2);
	h0:apply("*",vx1L);
	h1:apply("*",vx2L);
	h2:apply("*",vx3L);
	h0:ev(h0,VivL)$  h0:Reduce(h0,CurConst); 
	h1:ev(h1,VivL)$  h1:Reduce(h1,CurConst);
	h2:ev(h2,VivL)$  h2:Reduce(h2,CurConst)$ 
	
	t0:(h0+h1+h2)/3$              t0:Reduce(t0,CurConst);
	t1:(h0+ω*h1+ω^2*h2)/3$  t1:Reduce(t1,CurConst)$
	t2:(h0+ω^2*h1+ω*h2)/3$   t2:Reduce(t2,CurConst)$
	/*HTML用出力*/ 
	t1:expand(t1)$ t2:expand(t2)$
	mp:hipow(t1,x);
	a2:coeff(t1,x,mp);
	b2:coeff(t2,x,mp);
	/*出力終了*/
	t1:expand(t1)$  a2:coeff(t1,x,hipow(t1,x))$
	t2:expand(t2)$  b2:coeff(t2,x,hipow(t2,x))$
	
	ia2:Inverse(a2,CurConst)$
	ib2:Inverse(b2,CurConst)$
	q1:Reduce(ia2*t1,CurConst)$ q1:expand(q1);
	q2:Reduce(ib2*t2,CurConst)$ q2:expand(q2);
	/*HTML 出力用*/
	ia2:Inverse(a2,CurConst);
	ib2:Inverse(b2,CurConst);
	d2:a2^3$  d2:Reduce(d2,CurConst)$ d2:expand(d2);
	A[2]:d2$  B[2]:a[2]^3-A[2];
	
	id2:Inverse(d2,CurConst);
	ev2:a2*b2$  ev2:Reduce(ev2,CurConst);
	ev2:Reduce(ev2*id2,CurConst)$ ev2:expand(ev2)*a[2]^2;
	
	t1:a[2]*q1;  t2:ev2*q2;
	t2:Reduce(t2,CurConst)$ t2:expand(t2);
	ratvars(a[2])$  ratsimp(t1); ratsimp(t2);
	gx[2]:t0+t1+t2; gv[2]:subst(v,x,gx[2]);
	/*　　以上で3乗根多項式が求められた　*/
	hipow(gx[2],x);
	c4:coeff(gx[2],x,4);c3:coeff(gx[2],x,3);  c2:coeff(gx[2],x,2);  c1:coeff(gx[2],x,1);  c0:coeff(gx[2],x,0);
	/*　拘束条件の更新　*/
	
	cn:[a[2],B[2]]$  Const:cons(cn,Const);
	
	CurConst:cons([v,gv[2]],Const);
	
	
	/*拘束条件の更新完了*/
	/*   第6ステップ　　P=2 */
	ratvars(x)$
	Gal31:[1,8]$
	Gal32:[17,24]$
	N3:length(Gal31);
	vx1L:makelist(x-v[Gal31[i]],i,1,N3);
	vx2L:makelist(x-v[Gal32[i]],i,1,N3);
	h0:apply("*",vx1L);
	h1:apply("*",vx2L);
	h0:ev(h0,VivL)$  h0:Reduce(h0,CurConst);
	h1:ev(h1,VivL)$  h1:Reduce(h1,CurConst)$
	t0:expand((h0+h1)/2);  t1:expand((h0-h1)/2);
	pw:hipow(t1,x);
	t1:expand(t1)$  cv1:coeff(t1,x,hipow(t1,x)); cv0:coeff(t1,x,0);
	t1:expand(t1)$  cv1:coeff(t1,x,hipow(t1,x))$
	
	dv1:cv1^2$ dv1:Reduce(dv1,CurConst)$ dv1:expand(dv1);
	A[3]:dv1$  B[3]:a[3]^2-A[3];
	
	icv1:Inverse(cv1,CurConst);
	q1:Reduce(icv1*t1,CurConst)$ q1:expand(q1);
	/*HTML出力用*/
	
	dv1:expand(dv1); ratvars(a[2])$ dv1:ratsimp(dv1);
	
	/*HTML　終了*/
	t1:a[3]*q1;
	gx[3]:t0+t1;   gv[3]:subst(v,x,gx[3]);
	/*　拘束条件の更新　*/
	
	cn:[a[3],B[3]]$  Const:cons(cn,Const);
	
	CurConst:cons([v,gv[3]],Const);
	
	
	/*拘束条件の更新完了*/
	/*  第7ステップ　*/
	Gal41:[1]$
	Gal42:[8]$
	N4:length(Gal41);
	vx1L:makelist(x-v[Gal41[i]],i,1,N4);
	vx2L:makelist(x-v[Gal42[i]],i,1,N4);
	h0:apply("*",vx1L);
	h1:apply("*",vx2L);
	h0:ev(h0,VivL)$ h0:Reduce(h0,CurConst);
	h1:ev(h1,VivL)$ h1:Reduce(h1,CurConst);
	
	t0:expand((h0+h1)/2);  t1:expand((h0-h1)/2);
	t1:expand(t1)$  cv1:coeff(t1,x,hipow(t1,x))$
	
	dv1:cv1^2$ dv1:Reduce(dv1,CurConst)$ dv1:expand(dv1);
	A[4]:dv1$  B[4]:a[4]^2-A[4];
	
	icv1:Inverse(cv1,CurConst)$
	q1:Reduce(icv1*t1,CurConst)$ q1:expand(q1);
	t1:a[4]*q1$
	gx[4]:expand(t0+t1); gv[4]:subst(v,x,gx[4]);
	/*　拘束条件の更新　*/
	
	cn:[a[4],B[4]]$  Const:cons(cn,Const);
	
	CurConst:cons([v,gv[4]],Const);
	
	
	/*拘束条件の更新完了*/
	W:solve(gx[4],x); 
	B[1]; B[2]; B[3]; B[4]; Ω;
	s:rhs(SoL[1]);
	s:remainder(s,gv[4],v)$
	s:remainder(s,B[4],a[4])$
	s:remainder(s,B[3],a[3]);
	 s:remainder(s,B[2],a[2]);
	s:remainder(s,B[1],a[1]);
	s:remainder(s,Ω,ω);
	rsoL:[]$
	for i:1 thru N do(
	
	    s:rhs(SoL[i]),  
	    s:remainder(s,gv[1],v),  s:remainder(s,gv[2],v), s:remainder(s,gv[3],v),
	    s:remainder(s,gv[4],v),   s:remainder(s,Ω,ω),
	    s:remainder(s,B[4],a[4]),  s:remainder(s,B[3],a[3]),  s:remainder(s,B[2],a[2]), 
	    s:remainder(s,B[1],a[1]),  s:remainder(s,Ω,ω), 
	    rsoL:endcons(s,rsoL)
	)$
	rsoL;
	rω:allroots(Ω)$  fω:rhs(rω[1]);
	
	for i:1 thru 4  do(  B[i]:subst(fω,ω,B[i]) )$ 
	
	ra1:allroots(B[1])$  fa1:rhs(ra1[1]);
	
	for i:2 thru 4 do( B[i]:subst(fa1,a[1],B[i]))$
	
	ra2:allroots(B[2])$  fa2:rhs(ra2[1]);
	
	for i:3 thru 4 do( B[i]:subst(fa2,a[2],B[i]))$
	
	ra3:allroots(B[3])$  fa3:rhs(ra3[1]);
	
	for i:4 thru 4 do( B[i]:subst(fa3,a[3],B[i]))$
	
	ra4:allroots(B[4])$  fa4:rhs(ra4[1]);
	AnsLf:[]$
	for i:1 thru 4 do(
	     z:rsoL[i], z:subst(fω,ω,z), z:subst(fa1,a[1],z), z:subst(fa2,a[2],z), z:subst(fa3,a[3],z), z:subst(fa4,a[4],z),
	    z:expand(z),
	    AnsLf:endcons(z,AnsLf)
	)$
	AnsLf;
	allroots(fx);