計算不変式論の使い方:合同変換群の多項式不変量の計算事例 -

簡約群の(線形な作用に関する)不変式環(の生成元)を計算するアルゴリズムは、既に知られている(標数0の場合も正標数の場合も)(といっても、発見されたのは、ここ10年ほどのこと)。一方、簡約とは限らない一般の(アフィン)代数群については、不変式環が有限生成となることが保証されてないので、難しい。簡約でない最も簡単な群は、加法群で、並進移動に関する不変性を考えたい場合は多いので、実用上は面倒なことになる。最近(2013年10月)、Kemperによって、簡約とは限らない一般の代数群の不変式環を計算するアルゴリズムが与えられた(有限生成でない場合は停止しないので、正確には半アルゴリズムであるけども)。
Using Extended Derksen Ideals in Computational Invariant Theory
http://arxiv.org/abs/1310.6851

#どーでもいいけど、単に簡約群というと、普通、簡約代数群を指していると思うけど、簡約Lie群というものもあって、これらはオーバーラップするにも関わらず、違う概念であるという、紛らわしい用語の使い方になっている。簡約代数群は代数的閉体上でのみ定義される概念で、特に、複素数体上の簡約代数群は、コンパクトLie群の複素化として得られる代数群という特徴づけができて(コンパクトLie群には実代数群の構造が一意に入るというのも有名な事実)、全ての有理表現は完全可約である(この性質を線形簡約という。標数0の体で考える時は、簡約と線形簡約は同じことであるけど、正標数の体では簡約でも線形簡約とは限らない)。一方、簡約Lie群は、随伴表現が完全可約であるという特徴づけを持つものの、他の表現が完全可約であることは全然保証しない。またLie群のLie代数が簡約Lie代数であれば、それは簡約Lie群である(これは簡約Lie群の標準的な定義)けど、(アフィン)代数群の場合はLie代数が簡約Lie代数でも簡約代数群とは限らない(e.g.加法群のLie代数は可換なので簡約Lie代数だけど、簡約代数群ではない)

アルゴリズム自体を理解するには、多少の代数幾何と可換環の知識が必要だけど、アルゴリズムが正しいと信用してしまえば、あれこれ面倒なことを考えず(簡約群とは何かとか、自分が扱っているのは線形簡約群なのかとか)、一般の代数群の不変式をコンピュータに計算させられる時代がやってきた(多くのケースで、線形な作用を扱えれば事足りるので、ユーザとしては行列と多項式の扱いが理解できてれば十分)。適用事例として
(1)三角形の合同条件の導出
(2)平面二次曲線の合同変換による分類
(3)Huモーメント不変量の計算
について書いてみる(Kemperのアルゴリズムを、Risa/Asirに実装したものを、最後に書いた)

例0:

合同変換群の前に最も簡単な例として、次のような回転群の作用による不変式(不変式と言えば、有理関数を考える場合もあるけど、以下、多項式のみを想定する)を考える
$\left( \begin{matrix} x' \\ y' \end{matrix} \right) = \left( \begin{matrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{matrix} \right) \left( \begin{matrix} x \\ y \end{matrix} \right)$
という作用の下で、
$f(x',y')=f(x,y)$
を満たす多項式fを探す。任意の一変数多項式Fに対して
$f(x,y)=F(x^2+y^2)$
は条件を満たし、条件を満たすfはこれで尽きる。対象となる群は
$G= \{ \left( \begin{matrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{matrix} \right) \mid \theta \in \mathbf{R} \}$
であるけども、これは以下のように代数群として定義できる
$G = \{ \left( \begin{matrix} c & s \\ -s & c \end{matrix} \right) \mid c^2+s^2=1 , c \in \mathbf{R} , s \in \mathbf{R} \}$

Risa/Asirで書いたKemperのアルゴリズムを使って、不変式を計算するには

linvar([c*c+s*s-1] , newmat(2,2,[[c,s],[-s,c]]) , [x,y]);

とかやればいい。最初の2つの引数は、代数群の定義を与えていて、[x,y]は座標変数の名前なので任意(Risa/Asirでは小文字から開始していると不定元と見なされる)。返り値は

[x^2+y^2]

となる

例1:三角形の合同条件

三角形の合同条件と不変式論
http://d.hatena.ne.jp/m-a-o/20130104#p1
に書いたのと同じことだけど、当時は、まだアルゴリズムが存在してなかったのだった。

三角形の三頂点の座標を(x1,y1),(x2,y2),(x3,y3)とする。これだけでは、平行移動が表現できないので、ダミーの変数zを用意して
$\left(\begin{matrix} x_1' \\ y_1' \\ x_2' \\ y_2' \\ x_3' \\ y_3' \\ z' \end{matrix}\right) = \left( \begin{matrix} c & s & 0 & 0 & 0 & 0 & a \\ -s & c & 0 & 0 & 0 & 0& b \\ 0 & 0 & c & s& 0& 0 & a \\ 0 & 0 & -s & c & 0 & 0 & b \\ 0 & 0 & 0 & 0 & c & s & a \\ 0 & 0 & 0 & 0 & -s & c & b \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} \right) \left( \begin{matrix} x_1 \\ y_1 \\ x_2 \\ y_2 \\ x_3 \\ y_3 \\ z \end{matrix} \right)$
というような変換を考える。但し
$c^2+s^2=1$
とする。こうすると、例えば
$x_1 = cx_1 + sy_1 + az$
なので、z=1と置けば、平行移動を表現できる。なので、不変式を計算したあと、z=1とおく。不変式の計算は

linvar([c*c+s*s-1] ,
        newmat(7,7,[[c,s,0,0,0,0,a],
                    [-s,c,0,0,0,0,b],
                    [0,0,c,s,0,0,a],
                    [0,0,-s,c,0,0,b],
                    [0,0,0,0,c,s,a],
                    [0,0,0,0,-s,c,b],
                    [0,0,0,0,0,0,1]]) ,
        [x1,y1,x2,y2,x3,y3,z]);

とかやればいい。結果は

[x1^2-2*x3*x1+x3^2+y1^2-2*y3*y1+y3^2,
x2^2-2*x3*x2+x3^2+y2^2-2*y3*y2+y3^2,
(-y2+y3)*x1+(y1-y3)*x2+(-y1+y2)*x3,
(x2-x3)*x1-x3*x2+x3^2+(y2-y3)*y1-y3*y2+y3^2,
z]

みたいな感じで、ぱっと見では、ちょっと分かりにくいけど、整理すると
$P=(x_1-x_3)^2 + (y_1-y_3)^2$
$Q=(x_2-x_3)^2 + (y_2-y_3)^2$
$R=(y_1-y_2)(x_2-x_3) - (y_2-y_3)(x_1-x_2)$
$S=(x_1-x_3)(x_2-x_3) + (y_1-y_3)(y_2-y_3)$
$T=z$
という5つの不変式が出る。z=1とすると、Tは定数になるので、自明となる。残るP,Q,R,Sは高校生なら、よく知っている式で、以前計算したのと同じような結果を再現する

実は、本来のKemperのアルゴリズムは、ここで扱ってるように、行列で書けないような作用でも計算することができる。これには、論文のalgorithm6.1とalgorithm6.2を使えばいい。algorithm6.1は全ての計算で必要なので実装してあるけども、algorithm6.2は、面倒なので実装してない。これを、ちゃんと実装すれば、ダミーの変数zを導入するというトリックは要らなくなるはず。一応、algorithm6.1の適用結果だけ示しておくと

localizedInvariants([] , [x1,y1,x2,y2,x3,y3] , [c^2+s^2-1] , [c,s,a,b] , 
                    [c*x1+s*y1+a , -s*x1+c*y1+b,
                     c*x2+s*y2+a , -s*x2+c*y2+b,
                     c*x3+s*y3+a , -s*x3+c*y3+b]);

から、

[-x1^2+2*x2*x1-y1^2+2*y2*y1-x2^2-y2^2,
  [-x1^2+(x2+x3)*x1-y1^2+(y2+y3)*y1-x3*x2-y3*y2,
   (-y2+y3)*x1+(x2-x3)*y1-y3*x2+x3*y2,
   (x2-x3)*x1+(y2-y3)*y1-x2^2+x3*x2-y2^2+y3*y2,
   -x1^2+2*x2*x1-y1^2+2*y2*y1-x2^2-y2^2]]

みたいな結果を返す。ここに出てくる5つの式は、整理すると
$-(x_1-x_2)^2-(y_1-y_2)^2$
$-(x_1-x_3)(x_1-x_2) - (y_1-y_3)(y_1-y_2)$
$(-y_2+y_3)(x_1-x_2) - (y_1-y_2)(-x_2+x_3)$
$(x_2-x_3)(x_1-x_2)+(y_2-y_3)(y_1-y_2)$
$-(x_1-x_2)^2-(y_1-y_2)^2$
で、最初のと最後のは、同じ式であるけども、全部不変式となっている。

例2:平面二次曲線の分類

これはよくある話題で、
幾何学と不変量
http://www.amazon.co.jp/dp/4535784639
とかにも書いてあるっぽい(目次だけからの推測で未読なので、書いてないかもしれない。多分4章あたり)。岩波から出ていた『モジュライ理論I』という本にも書いてあった気がする(もう本は手元にないので、記憶が確かならだが)

平面二次曲線は、一般に
$\{(x,y) \in \mathbf{R}^2 | px^2+qxy+ry^2+tx+uy+v=0 \}$
という形で定義される。一般に、6つのパラメータ(p,q,r,t,u,v)を指定すれば、一つ曲線が定まる。p=q=r=0であれば、直線であるけども、そういう場合も含めて考える。

それで、例えば、円周
$x^2+y^2=1$
のx座標をaだけ平行移動することを考えると、その曲線の定義式は
$(x-a)^2+y^2=x^2+2ax+a^2+y^2=1$
という形となる。一般に、x座標をaだけ平行移動すると(p,q,r,t,u,v)で定義されていた曲線は
$(p,q,r,-2ap+t,-aq+u,a^2p-at+v)$
で定義される曲線に変換される。

というように、(x,y)平面上の合同変換で移りあう曲線を区別することを考える。このような目的には、不変式を計算するといい。一応、(p,q,r,t,u,v)を一斉に(0でない)定数倍しても、定義される曲線は同じだけど、このような変換の下での不変式は、有理関数しか存在しない。なので、合同変換による不変式を計算する。それらの不変式が同次式であれば、スケール変換も含めた不変有理式は、合同変換による不変式から作ることができる

計算用のRisa/Asirコードは(monocoefやmsubstは末尾のコード参照)

P = p*x^2+q*x*y+r*y^2+t*x+u*y+v;
P = msubst(P , [x,y] , [c*x+s*y+a , -s*x+c*y+b]);  /* 代入処理 */
Monomials = [x^2,x*y,y^2,x,y,1];
GVars = [p,q,r,t,u,v];
M=newmat(6,6);
for(I = 0 ; I < length(Monomials) ; I++){
   for(J = 0 ; J < length(GVars) ; J++){
      M[I][J] = monocoef(monocoef(P,Monomials[I],[x,y]),GVars[J] , GVars);
   }
}

linvar([c*c+s*s-1],M,GVars);

計算結果は

[p*u^2-t*q*u+(-4*r*p+q^2)*v+t^2*r,-4*r*p+q^2,p+r]

となる。2次の不変式
$q^2-4pr$
は、二次曲線の判別式というやつで、(曲線が直線や一点でないなら)これが正の時、双曲線、0の時、放物線、負の時、楕円となる。3次の不変式は、退化条件に対応するもので、これが0の時、その時に限り、曲線は、直線や空集合などになる。(p,q,r,t,u,v)の定数倍についても、不変な有理式は
$S=(pu^2-qtu-4prv+q^2v+t^2r)/(p+r)^3$
$T=(-4pr+q^2)/(p+r)^2$
の2つを取ることができる。つまり、平面二次曲線の分類空間は2次元であるということで、高校で習う楕円や双曲線の標準形が2つの独立なパラメータを含むという事実と合致している

例3:Huモーメント不変量の計算

Huモーメント不変量は、画像処理などで使われるものであるけど、ここでは数学的な話だけ扱う。Huモーメント不変量で、どういうことができるかは、
Huモーメント不変量による類似画像検索
http://d.hatena.ne.jp/m-a-o/20131104#p2

平面上の関数fに対して
$M_{pq}(f) = \int x^p y^q f(x,y) dxdy$
で定義される量がモーメントである。積分領域は全平面で、積分が収束することは一般に保証されないが、細かいことは考えず、全ての自然数p,qに対して、モーメントが有限であるような関数を対象にする。以下、p+qをモーメントの次数と呼ぶことにする

関数fを合同変換や拡大・縮小によって変換すると、別の関数に移る。この時、新しい関数のモーメントは、元の関数fのモーメントを使って書ける。例えば、x座標をaだけ平行移動すると
$\int x^p y^q f(x+a,y) dx dy = \int (x-a)^p y^q f(x,y) dxdy$
となる。このモーメントの変換のされ方は、関数fには依存していない。というわけで、今度は、合同変換群は、モーメント全体の空間に作用している。今までと違って、(p,q)が異なうモーメント間には何の関係もないので、独立なモーメントは無限個ある。けれども、あるモーメントの変換のされ方は、そのモーメントの次数以下のモーメントのみによって記述できる。従って、適当な自然数nに対して、n次以下のモーメントの空間への合同変換群の作用が定義でき、このようにして、高々有限個の独立なモーメントを考えればよくなる。

Huモーメント不変量は、拡大・縮小変換によっても不変でないといけないが、スケール変換な不変量は、例2の時と同様、多項式の範疇で発見できなくなる。モーメントの拡大・縮小による変換のされ方は
$\int x^p y^q f(cx,cy) dxdy = c^{-2-p-q} \int x^p y^q f(x,y) dxdy$
で、0次のモーメント(p=q=0)は合同変換で不変なので、合同変換による不変式を、0次のモーメントのべきで割ればいい。Huモーメント不変量は実際に、そのようにして定義されている

Huモーメント不変量というのは、一般に3次までのモーメントの組み合わせによって得られる不変式から得られるけど、例として、2次以下のモーメントで作られる合同変換の不変量を計算する(まぁ、3次までやろうとすると、計算が全然終わらないからであるけど)。2次以下のモーメントは、(p,q)=(0,0),(1,0),(0,1),(2,0),(1,1),(0,2)の6つがある(0次のモーメントは自明なので、実質的には5次元)。計算用のコードは

Monomials = [1,x,y,x^2,x*y,y^2];
M=newmat(6,6);
for(I = 0 ; I < length(Monomials) ; I++){
   P = msubst(Monomials[I] , [x,y] , [c*x+s*y+a , -s*x+c*y+b]);
   for(J = 0 ; J < length(Monomials) ; J++){
      M[I][J] = monocoef(P , Monomials[J] , [x,y]);
   }
}

linvar([c*c+s*s-1], M , [m00,m10,m01,m20,m11,m02]);

で、計算結果は

[(-m02*m20+m11^2)*m00+m02*m10^2-2*m11*m01*m10+m20*m01^2,
 (-m20-m02)*m00+m10^2+m01^2,
 m00]

の3つが返ってくる。このうち、m00は0次のモーメントなので、残りの2つが本質的にHuモーメント不変量

WikipediaのImage moment
http://en.wikipedia.org/wiki/Image_moment
を参照すると、scale invariant momentというので、式が書かれているけど、例えば
$\mu_{00}^2 I_1 = \mu_{20}+\mu_{02} = ( (M_{20} M_{00} - M_{10}^2)+(M_{02}M_{00}-M_{01}^2) )/M_{00}$
なので、M_{00}の冪を除いて、2番目の不変式と一致している。

Wikipediaによる、もう一つの2次のHuモーメント不変量
$\mu_{00}^4 M_{00}^2 I_2 = (M_{00}M_{20}-M_{10}^2-M_{00}M_{02}+M_{01}^2)^2+4(M_{11}M_{00}-M_{10}M_{01})^2$
については
$(M_{00}^5 I_2 - M_{00}^3 I_1^2)/4$
が1番目の不変式と等しい。2次以下の不変式については、正しいものが返ってきている

3次の場合もやりたいのだけど、うちの環境では、24時間たっても計算が終わらなかった。localizedInvariantsの最初のグレブナー基底の計算から進めない

一応実装

load("noro_pd.rr");
load("noro_matrix.rr");

def memq(Obj , Ls){
  for(C=Ls;C!=[];C=cdr(C)){
     if(Obj==car(C))return 1;
  }
  return 0;
}

def union(Ls , Rs){
   Ret = [];
   for(C=append(Ls,Rs) ; C!=[] ; C=cdr(C)){
      It = car(C);
      if( !memq(It,Ret) )Ret=cons(It , Ret);
   }
   return Ret;
}

def intersect(Ls , Rs){
    Ret = [];
    for(C=Ls;C!=[];C=cdr(C)){
       if(memq(car(C) , Rs))Ret=cons(car(C) , Ret);
    }
    return Ret;
}

def diffset(Ls , Rs){
   Ret = [];
   for(C=Ls;C!=[];C=cdr(C)){
      if(memq(car(C) , Rs))continue;
      Ret = cons(car(C) , Ret);
   }
   return Ret;
}

/*
多項式P中の単項式Mの係数を計算する
monocoef(y^3+x^3*y^2+2*x*y^2+2*x*y , x*y^2 , [x,y]);
なら、2を返す
monocoef(y^3+x^3*y^2+2*y^2+2*x*y , y^2 , [y]);
なら、x^3+2
monocoef(y^3+x^3*y^2+2*y^2+2*x*y , x*y^2 , [y]);
でも、x^3+2
*/
def monocoef(P,M,V){
   if( nmono(M)!=1 )return 0;
   HT = dp_ht(dp_ptod(M , V));
   for(T = dp_ptod(P,V) ; T ; T = dp_rest(T)){
      if(dp_ht(T)==HT)return dp_hc(T);
   }
}


def gen_monomials(Deg , V){
   if(Deg==0)return [1];
   Ret = [];
   for(C=gen_monomials(Deg-1 , V) ; C!=[] ; C=cdr(C)){
      for(I=0; I<length(V) ; I++){
         M = V[I]*car(C);
         if(!memq(M,Ret))Ret = cons(M , Ret);
      }
   }
   return Ret;
}

/*
msubst(t1 , [t1,t2,t3,t4],[t1*x^5,t2*x^3,t3*y^2,t4*y^4]);

*/
def msubst(P , Vars , Exprs){
   N = length(Vars);
   TV = [];
   TP = P;
   for(I = 0 ; I < N ; I++){
      NV = uc();
      TV = cons(NV , TV);
      TP = subst(TP , Vars[I] , NV);
   }
   TV=reverse(TV);
   for(I = 0 ; I < N ; I++){
      TP = subst(TP , TV[I] , Exprs[I]);
   }
   return TP;
}

/*
多項式係数の行列の定数係数のKernelを計算する
*/
def p_compute_kernel(PMat){
   Sz = size(PMat);
   Vars = vars(PMat);
   N = Sz[0];
   M = Sz[1];
   RPol = 0;
   for(I = 0 ; I < N ; I++){
      for(J = 0 ; J < M ; J++){
          for(V = dp_ptod(PMat[I][J] , Vars) ; V ; V=dp_rest(V)){
             RPol = RPol + dp_ht(V);
          }
      }
   }
   Terms = [];
   for(V=RPol;V;V=dp_rest(V)){Terms = cons(dp_ht(V) , Terms);}
   CMat = newmat(length(Terms)*N , M);
   for(Idx = 0 , Terms = Terms ; Terms!=[] ; Terms=cdr(Terms)){
      Term = dp_dtop(car(Terms) , Vars);
      for(W = 0 ; W < N ; W++ , Idx++){
           for(H = 0 ; H < M ; H++){
               CMat[Idx][H] = monocoef(PMat[W][H] , Term , Vars);
           }
      }
   }
   return linalg.compute_kernel(CMat);
}


/*
Algorithm 2.7 in
"Computing invariants of reductive groups in positive characteristic" by G.Kemper

IG:代数群の座標環の定義イデアル
RM:代数群の表現行列
Xs:代数群が作用する空間の座標(RMの行数/列数と等しくないといけない)
D:degree of homogeneous invariants(non-negative integer)


例:
homogeneousInvariants([a*d-b*c-1] , newmat(3,3,[[a^2 , 2*a*b , b^2],[a*c , a*d+b*c , b*d],[c^2 ,2*c*d ,d^2]]) , [x,y,z] , 2);

return [-z*x+y^2];

*/
def homogeneousInvariants(IG , RM , Xs , D){
    GVars = append(vars(IG) , vars(RM));
    AllVars = append(Xs , GVars);
    G = nd_gr(IG , GVars , 0 , 0);
    Monomials = gen_monomials(D , Xs);
    HPols = [];
    Ys = newvect(length(Xs));
    for(I = 0 ; I < length(Xs) ; I++){
       for(J = 0 ; J < length(Xs) ; J++){
          Ys[I] += RM[I][J]*Xs[J];
       }
    }
    for(C = Monomials ; C!=[] ; C=cdr(C)){
        P = car(C);
        Q = msubst(P , Xs , Ys) - P;
        HPols = cons(p_nf(Q , G , GVars , 0) , HPols);
    }
    Kernel = p_compute_kernel(newmat(1 , length(HPols) , [reverse(HPols)]));
    Ret = [];
    for(I = 0 ; I < length(Kernel) ; I++){
       Evec = Kernel[I][0];
       P = 0;
       for(J = 0 ; J < length(Evec) ; J++){
          P = P + Evec[J]*Monomials[J];
       }
       if(P!=0)Ret = cons(P , Ret);
    }
    return Ret;
}



/*
step(5)(6)(7) of Algorithm 6.1 (Computation of a localization of an invariant ring)
*/
def rationalInvariants(IY , Ys , IX , Xs , Ord){
    Ret = [];
    CurOrd = dp_ord();
    dp_ord(Ord);
    for(RG = nd_gr_trace(IY , Ys , 1 , 1 , Ord) ; RG!=[] ; RG=cdr(RG)){
       F = p_nf(car(RG) , IX , Xs , 0);   /* IX should be Groebner basis */
       if(F==0)continue;
       V = dp_ptod(F , Ys);
       LC_F = dp_hc(V);
       for(V=dp_rest(V) ; V ; V=dp_rest(V)){
          C_F = dp_hc(V);
          Q = red(C_F/LC_F);
          if(length(vars(Q))==0)continue;
          Check = 0;
          for(K=0;K<length(Ret);K++){
             R=red(Ret[K]/Q);
             if(length(vars(R))==0){
               Check=1;
               break;
             }
          }
          if(Check==0){
            Ret = cons(Q , Ret);
          }
       }
    }
    dp_ord(CurOrd);
    return Ret;
}


/*
algorithm6.1

IX:代数群が作用する空間の定義イデアル
Xs:代数群が作用する空間の座標変数
IG:代数群の座標環の定義イデアル
Zs:代数群の座標環の変数
GReps:代数群の作用を表す多項式
*/
def localizedInvariants(IX , Xs , IG , Zs , GReps){
   N = length(Xs);
   /* step(1) is omitted */
   /* step(2) */
   Ys = [];
   GG = nd_gr(IG , Zs , 0 , 0);
   Ehat = append(IX , GG);
   for(I = 0 ; I < N ; I++){
     YI = uc();
     Ys = cons(YI , Ys);
     Ehat = cons(YI - GReps[I] , Ehat);
   }
   /* step(3) */
   AllVars = append(Zs , append(Ys, Xs));
   Ghat = nd_gr_trace(Ehat , AllVars , 1 , 1 , [[0,length(Zs)],[0,length(Ys)],[0,length(Xs)]]);
   G = noro_pd.elimination(Ghat , append(Xs , Ys));
   Gx = noro_pd.elimination(Ghat , Xs);
   Gy = diffset(G , Gx);
   /* step(4) is skipped and proceed step(5)(6)(7) */
   QI = rationalInvariants(Gy , Ys , Gx , Xs , 0);
   QRats = [];
   QPols = [];
   for(I=0;I<length(QI);I++){
      if(length(vars(dn(QI[I])))!=0){
         QRats = cons(QI[I] , QRats);
      }else{
         QPols = cons(QI[I] , QPols);
      }
   }
   if(QRats==[]){
      return [1,QPols];
   }
   Hs = map(dn , QRats);
   Fs = map(nm , QRats);
   /* step(8)(9)(10) */
   Monomials = [];
   K = length(QRats);
   for(R = 0 ;; R++){
      for(C=gen_monomials(R , Xs) ; C!=[] ; C=cdr(C)){
         Check = 0;
         for(Ms = Monomials ; Ms!=[] ; Ms=cdr(Ms)){
            if(p_nf(car(Ms)-car(C) , Gx , Xs , 0)==0){
               Check=1;
               break;
            }
         }
         if(Check==0)Monomials = cons(car(C) , Monomials);
      }
      L = length(Monomials);
      NF_fm = newmat(L,K);
      NF_hm = newmat(L,K);
      NF_G = [];
      for(I = 0 ; I < L ; I++){
         P = msubst(Monomials[I] , Xs , GReps) - Monomials[I];
         P = p_nf(P , Gx , Xs , 0);
         P = p_nf(P , GG , Zs , 0);
         NF_G = cons(P , NF_G);
         for(J = 0 ; J < K ; J++){
            NF_fm[I][J] = p_nf(Fs[J]*Monomials[I] , Gx , Xs , 0);
            NF_hm[I][J] = p_nf(Hs[J]*Monomials[I] , Gx , Xs , 0);
         }
      }
      NF_G = reverse(NF_G);
      /* solve equations in step(9) */
      PMat = newmat(K+1 , L*(K+1));
      for(J = 0 ; J < K ; J++){
          for(I = 0 ; I < L ; I++){
              PMat[J][I] = NF_fm[I][J];
          }
          for(JI = 0 ; JI < L ; JI++){
             for(JJ = 0 ; JJ < K ; JJ++){
                if(J==JJ)PMat[J][L+JI*K+JJ] = NF_hm[JI][JJ];
             }
          }
          for(I = 0 ; I < L ; I++){
             PMat[K][I] = NF_G[I];
          }
      }
      Kernel = p_compute_kernel(PMat);
      if(length(Kernel)!=0){
         Vs = Kernel[0][0];
         AP = 0;
         for(I = 0 ; I < L ; I++)AP = AP+Vs[I]*Monomials[I];
         if(AP!=0)break;
      }
   }
   Bs = [];
   for(J = 0 ; J < K ; J++){
      BP = 0;
      for(I = 0 ; I < L ; I++){
          BP = BP+Vs[L+I*K+J]*Monomials[I];
      }
      if(BP!=0)Bs = cons(BP , Bs);
   }
   Bs = union(QPols , Bs);
   Check = 0;
   for(I = 0 ; I < length(Bs) ; I++){
      if( length(vars(red(Bs[I]/AP)))==0 ){Check=1;break;}
   }
   if(Check==0){Bs=cons(AP,Bs);}
   return [AP,Bs];
}


def linvar(IG , RM , Xs){
   GVars = append(vars(IG) , vars(RM));
   NM = size(RM);
   Ts = [];
   for(I = 0 ; I < length(Xs) ; I++){
      NP = 0;
      for(J = 0 ; J < length(Xs) ; J++){
          NP = NP + RM[I][J]*Xs[J];
      }
      Ts = cons(NP , Ts);
   }
   AB = localizedInvariants([] , Xs , IG , GVars , reverse(Ts));
   AP = AB[0];
   BPols = AB[1];
   /* start algorithm 7.1 */
   Ret = [];
   for(Deg = 1 ; ; Deg++){
      G = [];
      Ys = newvect(length(Ret));
      for(I = 0 ; I < length(Ret) ; I++){
         Ys[I] = uc();
         G = cons(Ys[I] - Ret[I] , G);
      }
      /* algorithm7.1 step(2) */
      Jhat = nd_gr(G , append(Xs , vtol(Ys)) , 0 , [[0,length(Xs)],[0,length(Ys)]]);
      J = noro_pd.elimination(Jhat , vtol(Ys));
      Check1 = 1;  /* $B \subset A$? */
      Check2 = 0;  /* $K[X] \dot a \cap A= A \dot a$? */
      /* Check B \subset A */
      /* 
         algorithm for subring contaiment:
         * subroutine 2.5.4 in "Algorithms in Invariant Theory" by Sturmfels
         * subroutine 2.5 in 
         "Efficient computation offundamental invariants - An approach using Buchberger's Grobner bases method" 
         (Sturmfels and White)
      */
      for(I = 0 ; I < length(BPols) ; I++){
         QPol = p_nf(BPols[I] , Jhat , append(Xs , vtol(Ys)) , 0);
         if( intersect(vars(QPol) , Xs)!=[] ){
             Check1 = 0;
             break;
         }
      }
      if(AP==1){
         Check2 = Check1;
      }else if(Check1==1){
         Check2 = 1;
         for(Hs = J ; Hs!=[] ; Hs=cdr(Hs)){
             HPol = red(msubst(car(Hs) , vtol(Ys) , Ret)/AP);
             if(dn(HPol)!=1){Check2=0;break;}
             QPol = p_nf(HPol , Jhat , append(Xs , vtol(Ys)) , 0);
             if( intersect(vars(QPol) , Xs)!=[] ){
                Check2 = 0;
                break;
             }
         }
      }
      if(Check2==1)return Ret;
      /* algorithm 7.2 step(3)(4) */
      for(S=homogeneousInvariants(IG,RM,Xs,Deg) ; S!=[] ; S=cdr(S)){
         SPol = p_nf(car(S) , Jhat , append(Xs , vtol(Ys)) , 0);
         if( intersect(vars(SPol) , Xs)!=[] ){
            Ret = cons(car(S), Ret);
         }
      }
   }
}