有限群の不変式環の生成元を計算するアルゴリズム

Asirによる有限群の不変式環の生成元の計算
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42593
を、そのままぱくっただけ。実装が見当たらなかったので、自分で実装した。同じ悲しみに遭遇した人の役に立つこともあるかもしれない


そして、上の論文も、以下の本のほぼ"パクリ"であるらしい。この業界(計算不変式論?)の基本図書らしい。わたしは読んでないけど
Algorithms in Invariant Theory
http://www.amazon.co.jp/dp/3211774165


有限群の不変式の計算は、SingularやMagamaにも実装されてるらしいので、計算するだけなら、それ使うのが一番確実かもしれない
Online Manual - finvar_lib
http://www.singular.uni-kl.de/Manual/3-0-4/sing_1086.htm#SEC1145

Invariant Rings of Finite Groups
http://magma.maths.usyd.edu.au/magma/handbook/text/1219


アルゴリズムの流れは、大体以下の事実に依存している。どれも証明はそんなに難しくない(Hilbert級数の計算は、無駄な計算を減らすためのもので、必須ではない)
(1)有限群の不変式環は有限生成(by Hilbert?)
(2)Noetherの正規化定理
(3)有限群の不変式環はCohen-Macaulay環(Hochster-Eagon)


不変式環のNoether正規化を与える代数的に独立な不変式をprimary invariantsと呼んでいて(primary invariantsの生成する不変式環の部分環をプライマリ不変式環と呼ぶことにすると)、不変式環のCohen-Macaulay性から、不変式環は、プライマリ不変式環上の有限生成自由加群となる。その加群の基底をsecondary invariantsと呼んでいる。プライマリ不変式環は元の多項式環のNoether正規化でもあり、逆に不変式から生成される多項式環のNoether正規化は、プライマリ不変式環を与える。というわけで、primary invariantsの計算が、不変式から生成される多項式環のNoether正規化を探す問題に帰着する。primary invariantsが分かれば、secondary invariantsは、プライマリ不変式環の自由加群としての基底なので、容易に計算できる。あと、primary invariantsもsecondary invariantsも同次にとれるので、Reynolds作用素を単項式に作用させたものから探せばいい



例1:位数2の巡回群の二次元作用f(x,y)=f(-x,-y)による不変式の場合
primary invariantsとして、(例えば)x^2とy^2が取れる。そして、基礎体をKとして、不変式環は、K[x^2,y^2]上、1とxyで生成される自由加群となる。結局、X=x^2,Z=xy,Y=y^2を不変式の生成元として取って、関係式はXY-Z^2=0


例2:対称群の置換による不変式を考えると、primary invariantsとして、基本対称式が取れる。よく知られている通り、これで不変式環が生成されるので、secondary invariantsは1しかない。つまり、この場合は、不変式環=ring of primary invariantsとなる。Chevalley-Shephard-Toddの定理によって、任意の鏡映群に対して、適当なprimary invariantsを選べば、secondary invariantsは1のみになる


[future work]
(1)で、別に不変式が計算したかったわけではなく、本題は、自明表現以外の既約表現成分を計算したかった。一次元既約表現の場合は、相対不変式を考えるのに相当する。これらは不変式環上の加群で、商空間(一般に特異点を持ちうる)上の同変ベクトル束と見ることもできる。基本的には、各既約表現成分もプライマリ不変式環上の有限生成自由加群となり、Reynolds作用素を一般化した射影作用素も、既約指標を使って書くことができる(Molienの公式も任意の表現に容易に拡張できる)。なので、既約表現成分の計算は、secondary invariantsの計算とほぼ同じになる(はず)。


(2)それで、既約成分のプライマリ不変式環上の加群としての基底が分かったら、不変式環上の加群としての生成系と構造を調べないといけない。例えば、aを1の原始3乗根として、f(x,y)=f(ax , y/a)という不変式を考えると、x^3,xy,y^3が生成元となり、x^3とy^3をprimary invariantsとできる。それで、{x,y^2,x^2*y}を基底とするプライマリ不変式環上の加群は、既約成分の一つとなる。これはプライマリ不変式環上自由であるけども、不変式環上では
x^2*y=(不変式)*x + (不変式)*y^2
となる不変式が存在して、不変式環上の加群としての独立な生成系は、{x,y^2}のみ。この計算を実装したい(どうやって?)

load("gr");
load("sp"); /* 代数的数の操作 */

/* 補助関数 */
def choose1(CUR , LS){
   RET = [];
   if(length(CUR)==0){
      for(I=0;I<length(LS);I++){
         RET = cons( [LS[I]] , RET);
      }
      return RET;
   }
   for(I = 0 ; I<length(CUR) ; I++){
       TMP = CUR[I];
       NewElems = [];
       for(K=0;K<length(LS);K++){
          OK = 1;
          for(J=0;J<length(TMP);J++){
              if(TMP[J]==LS[K])OK = 0;
          }
          if(OK)NewElems = cons(LS[K] , NewElems);
       }
       for(L=0;L<length(NewElems);L++){
          RET = cons(cons(NewElems[L],TMP) , RET);
       }
   }
   return RET;
}

/* LSの中からN個の要素を抽出したリストのリスト */
def choose(LS , N){
   RET = [];
   for(I=0;I<N;I++){
      RET = choose1(RET , LS);
   }
   return RET;
}

/* ObjがLSの要素かどうか */
def memq(Obj , Ls){
  for(K = 0; K < length(Ls) ; K++){
     if(Obj==Ls[K])return 1;
  }
  return 0;
}

/* 行列群の生成元から全ての要素の集合を返す */
def grp(Gens){
  Ret = [];
  if(length(Gens)==0)return Ret;
  Dim = size( Gens[0] )[0];
  Id = newmat(Dim , Dim);
  for(I = 0 ; I < Dim ; I++)Id[I][I] = 1;
  Zero = newmat(Dim , Dim);
  Ret = Gens;
  for(;;){
     Cnt = 0;
     for(I = 0 ; I <length(Gens) ; I++){
        for(J = 0 ; J < length(Ret) ; J++){
            E = map(simpalg,Gens[I] * Ret[J]);
            if(E==Id)continue;
            /* 新しい要素かどうか */
            Check = 1;
            for(K=0 ; K<length(Ret) ; K++){
               /* 代数的数同士の比較では、simpalg(a)==simpalg(b)では不十分 */
               if(map(simpalg,E-Ret[K])==Zero){
                   Check = 0;
                   break;
               }
            }
            if(Check==1){
                 Ret = cons(E , Ret);
                 Cnt += 1;
            }
        }
     }
     if(Cnt == 0)return cons(Id,Ret);
  }
}


/* 
0次元イデアルかどうか

注意:grライブラリに同名のzero_dim関数があるけど、バグっている気がするので再定義している

例1:grライブラリのzero_dimが正しくない結果を返すケース
zero_dim([a*b*c*d-1,a*b*c+b*c*d+c*d*a+d*a*b,a*b+b*c+c*d+d*a,a+b+c+d],[a,b,c,d],0);

例2:以下は、primedec(diff(F,x),diff(F,y),diff(F,z),F] ,[x,y,z])によって、x=y=z=0のみが解と分かる
F=x^2 + y^2 - z^5;
zero_dim([diff(F,x),diff(F,y),diff(F,z),F] ,[x,y,z] , 0);

*/
def zero_dim(I , Vars , Order){
      G = gr(I , Vars , Order);
      if(length(G) < length(Vars))return 0;
      IT = [];
      for(J = 0 ; J < length(G) ; J++){
          IT = cons(p_terms(G[J] , Vars , Order) ,IT);
      }
      for(J = 0 ; J < length(Vars) ; J++){
         Check = 0;
         for(K = 0 ; K < length(IT) ; K++){
            Deg = deg(IT[K],Vars[J]);
            if(Deg!=0 && IT[K] != Vars[J]^Deg){
                break;
            }else if(IT[K]==Vars[J]^Deg){
                Check = 1;
                break;
            }
         }
         if(Check!=0){return 0;}
      }
      return 1;
}

/*
根基イデアルの所属判定
例:
radmember(x^2+y^2+z^2 , [x+y+z,x*y+y*z+z*x] , [x,y,z]);
*/

def radmember(Expr , Ideal , Vars){
  T = uc();
  K = cons(Expr*T -1 , Ideal);
  return ( gr(K , cons(T,Vars) , 0)==[1] );
}

/*
イデアルIが、零次元である時、K[X]/Iの単項式基底を計算する

例:
p_mbase([x^2,x*y,y^2],[x,y]);
p_mbase([x^3+y^3+z^3 , x^2+y^2+z^2 , x+y+z] , [x,y,z]);
p_mbase([x^8+14*x^4*y^4+y^8 , 
         1025*(x^24+y^24)+10626*(x^4*y^20+x^20*y^4)+735471*(x^8*y^16+x^16*y^8)+2704156*x^12*y^12],
        [x,y]);
*/
def p_mbase(I , Vars){
  GBasis = gr(I , Vars , 0);
  H = map(dp_ptod , GBasis , Vars);
  return map(dp_dtop , dp_mbase(H) , Vars);
}


/* 
Molienの公式によるHilbert級数の計算(有理式を返すけど、redしてないので注意)
Elementsには、群の要素が行列で入っているとする
(つまり、有限群の線形表現の像の集合)

例:C_4の3次元表現
S = hilbertSeries([newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]),
               newmat(3,3,[[0,1,0],[-1,0,0],[0,0,-1]]),
               newmat(3,3,[[-1,0,0],[0,-1,0],[0,0,1]]),
               newmat(3,3,[[0,-1,0],[1,0,0],[0,0,-1]])]);
*/
def hilbertSeries( Elements ){
  T = uc();  /* Hilbert級数の変数 */
  G = length(Elements); /* 群の位数 */
  RET = 0;
  Dim = size(Elements[0])[0];
  Id = newmat(Dim , Dim); /* 恒等行列 */
  for (I = 0; I < Dim ; I++)Id[I][I] = 1;
  for (I = 0 ; I < G; I++){
     RET = simpalg(RET + 1/det(Id - T*Elements[I]));
  }
  return RET/G;
}


/* 次数Degの単項式の生成 */
def gen_monomial(Vars , Deg){
  RET = [];
  if( length(Vars)<=1 ){
     return [Vars[0]^Deg];
  }else{
     for(I= 0 ; I <= Deg ; I++){
        Cur = Vars[0]^I;
        Tmp = gen_monomial(cdr(Vars) , Deg-I);
        for(J = 0 ; J < length(Tmp) ; J++){
          RET = cons(Cur*Tmp[J] , RET);
        }
     }
  }
  return RET;
}


/* 同次多項式の次数 */
def hdeg(Expr , Vars){
  if(type(Expr)==1)return 0;
  RET = 0;
  T = p_terms(Expr , Vars, 0)[0];
  for (J = 0 ; J < length(Vars) ; J++){
      RET += deg(T , Vars[J]);
  }
  return RET;
}


/* 
primary invariantsの計算 
Elements:群の線形表現の像の集合
Vars:不変式の変数の集合

例1:
priminv([newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]),
         newmat(3,3,[[0,1,0],[-1,0,0],[0,0,-1]]),
         newmat(3,3,[[-1,0,0],[0,-1,0],[0,0,1]]),
         newmat(3,3,[[0,-1,0],[1,0,0],[0,0,-1]])] , [x,y,z]);

例2:3次対称群の3次元表現
priminv([newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]),
         newmat(3,3,[[0,1,0],[1,0,0],[0,0,1]]),
         newmat(3,3,[[1,0,0],[0,0,1],[0,1,0]]),
         newmat(3,3,[[0,1,0],[0,0,1],[1,0,0]]),
         newmat(3,3,[[0,0,1],[1,0,0],[0,1,0]]),
         newmat(3,3,[[0,0,1],[0,1,0],[1,0,0]])] , [x,y,z]);

例3:
Omega = newalg(x^2+x+1);
priminv( grp([newmat(2,2,[[Omega,0],[0,1/Omega]])]) , [x,y]);

例4:Shephard-Toddの複素鏡映群No.4
Omega = (-1-Root3I)/2;
Root3I = newalg(x^2+3);
G = grp( [newmat(2,2,[[1,2],[1,-1]])/(Root3I) , newmat(2,2,[[1,0],[0,Omega]])]);
priminv(G , [x,y]);

例5:Shephard-Toddの複素鏡映群No.8
Root2 = newalg(x^2-2);
Imaginary = newalg(x^2+1);
G = grp( [(1-Imaginary)*newmat(2,2,[[1,-1],[1,1]])/2 , newmat(2,2,[[-Imaginary,0],[0,1]])]);
priminv(G , [x,y]);

*/
def priminv( Elements ,Vars ){
  RET = [];
  TmpVars = newvect(length(Vars));
  for(I = 0 ; I < length(Vars) ; I++)TmpVars[I] = uc();
  /* assert(length(Elements) > 0); */
  /* length(Vars)==Sz[0]==Sz[1]==表現空間の次元 */
  Dim = size(Elements[0])[0];
  /* 群作用の行先 */
  GMap = newvect(length(Elements));
  for (K = 0 ; K < length(Elements) ; K++){
      Pt = newvect(Dim);
      for (I=0 ; I < length(Vars) ; I++){
         for (J = 0 ; J < Dim ; J++){
             Pt[I] += Elements[K][I][J] * Vars[J];
         }
      }
      GMap[K] = Pt;
   }
  /* main loop */
  H = hilbertSeries(Elements);
  GRet = [];
  for(D=1 ; ; D++){
      /* Hilbert級数を、次数Dまで形式的べき級数に展開 */
      S = uinv_as_power_series(dn(H) , D) * nm(H);
      C = coef(S , D);
      if(C==0)continue;
      Tmp = [];  /* 次数Dの不変式 */
      Monomials = gen_monomial(Vars , D);
      for(J = 0 ; J < length(Monomials) ; J++){
         /* Monomials[J]にReynolds作用素をapply */
         RM = 0;
         for(K = 0 ; K < length(Elements) ; K++){
            M = Monomials[J];
            for(I=0 ; I<length(Vars) ; I++)M=subst(M , Vars[I] , TmpVars[I]);
            for(I=0 ; I<length(Vars) ; I++)M=subst(M , TmpVars[I] , GMap[K][I]);
            RM += M;
         }
         RM = simpalg(RM);
         if(RM==0)continue;
         RM = ptozp( RM/length(Elements) );
         /* 既に同じ要素がTmpに含まれていないかチェック */
         if( memq(RM , Tmp)== 1)continue;
         /* RMがRETから独立であれば追加 */
         if( radmember(RM , GRet , Vars)==0 ){
            RET = cons(RM , RET);
            GRet = gr(RET , Vars , 0);
         }
         Tmp = cons(RM , Tmp);
         if( length(Tmp)==C )break;
      }
      /* RETの生成するイデアルの次元が0なら終了 */
      if( length(RET) >= Dim && zero_dim(GRet , Vars , 0)==1 ){
         break;
      }
  }
  Cand = choose(RET , Dim);
  for(I = 0 ; I < length(Cand) ; I++){
     if( zero_dim(Cand[I] , Vars , 0)==1 )return Cand[I];
  }
}


/*
有限群の不変式を計算する

例1:
finvar([newmat(2,2,[[1,0],[0,1]]),newmat(2,2,[[-1,0],[0,-1]])] , [x,y]);

例2:
finvar([newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]),
        newmat(3,3,[[0,1,0],[-1,0,0],[0,0,-1]]),
        newmat(3,3,[[-1,0,0],[0,-1,0],[0,0,1]]),
        newmat(3,3,[[0,-1,0],[1,0,0],[0,0,-1]])] , [x,y,z]);

例3:三次対称群
finvar([newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]),
        newmat(3,3,[[0,1,0],[1,0,0],[0,0,1]]),
        newmat(3,3,[[1,0,0],[0,0,1],[0,1,0]]),
        newmat(3,3,[[0,1,0],[0,0,1],[1,0,0]]),
        newmat(3,3,[[0,0,1],[1,0,0],[0,1,0]]),
        newmat(3,3,[[0,0,1],[0,1,0],[1,0,0]])] , [x,y,z]);

例4:
finvar([newmat(2,2,[[1,0],[0,1]]) , newmat(2,2,[[-1,0],[0,1]]) , 
        newmat(2,2,[[1,0],[0,-1]]) , newmat(2,2,[[-1,0],[0,-1]])] , [x,y]);

例5:4次巡回群 generated by f(x,y)=f(-y,x)
finvar( grp([newmat(2,2,[[0,-1],[1,0]])]) , [x,y]);

例6:Shephard?Toddの複素鏡映群No.9(少し時間がかかる)
Root2 = newalg(x^2-2);
Z2 = newalg(x^2+1);
G = grp([ newmat(2,2,[[1,1],[1,-1]])/Root2 , newmat(2,2,[[1,0],[0,Z2]])]);
finvar( G , [x,y]);

例7:binary tetrahedral group
Z2 = newalg(x^2+1);
Id = newmat(2,2,[[1,0],[0,1]]);
I = newmat(2,2,[[Z2 , 0], [0,-Z2]]);
J = newmat(2,2,[[0,1],[-1,0]]);
K = newmat(2,2,[[0 , Z2],[Z2,0]]);
G= [Id,-Id ,I,-I,J ,-J,K,-K , (Id+I+J+K)/2 , (-Id+I+J+K)/2 ,
    (Id-I+J+K)/2 , (Id+I-J+K)/2 ,(Id+I+J-K)/2,(-Id-I+J+K)/2 , (-Id+I-J+K)/2 ,
    (-Id+I+J-K)/2 , (Id-I-J+K)/2 , (Id-I+J-K)/2 , (Id+I-J-K)/2 ,
    (-Id-I-J-K)/2 , (Id-I-J-K)/2 , (-Id+I-J-K)/2 , (-Id-I+J-K)/2 , (-Id-I-J+K)/2];
finvar( G , [x,y]);

例8:yet another binary tetrahedral group(例7と同じ結果)
Root2 = newalg(x^2-2);
Z2 = newalg(x^2+1);
Z8 = (1+Z2)/Root2;
G = grp([newmat(2,2,[[Z8,Z8^3],[Z8,Z8^7]])/Root2 , newmat(2,2,[[Z8^2,0],[0,1/Z8^2]])]);
finvar( G , [x,y]);

例9:binary octahedral group
Root2 = newalg(x^2-2);
Z2 = newalg(x^2+1);
Z8 = (1+Z2)/Root2;
G = grp([newmat(2,2,[[Z8^3,0],[0,Z8^5]]) , newmat(2,2,[[Z8^7,Z8^7],[Z8^5,Z8]])/Root2]);
finvar( G , [x,y]);

例10:"二項二十面体群"(Gの位数が240になってるので、1の5乗根をきちんと書くこと)<-遅い
Root5 = newalg(x^2-5);
Z2 = newalg(x^2+1);
Z5 = newalg(1+x+x^2+x^3+x^4);
G = grp([newmat(2,2,[[-Z5+Z5^4,Z5^2-Z5^3],[Z5^2-Z5^3,Z5-Z5^4]])/Root5,
        -newmat(2,2,[[Z5^3,0],[0,Z5^2]])]);
finvar( G , [x,y]);

*/
def finvar(Elements , Vars){
   Ps = priminv(Elements , Vars);
   GPs = gr(Ps , Vars , 0);
   Ss = []; /*secondary invariants */
   TmpVars = newvect(length(Vars));
   for(I = 0 ; I < length(Vars) ; I++)TmpVars[I] = uc();
   /* 群作用の行先 */
   Dim = length(Vars);
   GMap = [];
   for (K = 0 ; K < length(Elements) ; K++){
      Pt = newvect(Dim);
      for (I=0 ; I < length(Vars) ; I++){
         for (J = 0 ; J < Dim ; J++){
             Pt[I] += Elements[K][I][J] * Vars[J];
         }
      }
      GMap = cons(Pt, GMap);
   }
   /* ================= */
   H = hilbertSeries(Elements);
   for(I = 0 ; I < length(Ps) ; I++){
       H = H*(1-var(H)^hdeg(Ps[I],Vars));
   }
   H = red(H);
   if(H==1)return [Ps,[1]];
   /*
   Monomials = [];
   for(D=0;D<=deg(H,var(H));D++){
         Monomials = append(gen_monomial(Vars,D) , Monomials);
   }
   */
   Monomials = p_mbase(Ps , Vars);
   for (L = 0 ; L < length(Monomials) ; L++){
      Deg = hdeg(Monomials[L] , Vars);
      if( coef(H,Deg)==0 )continue;
      /* Monomials[L]にReynolds作用素をapply */
      RM = 0;
      for(K = 0 ; K < length(Elements) ; K++){
         M = Monomials[L];
         for(I=0 ; I<length(Vars) ; I++)M=subst(M , Vars[I] , TmpVars[I]);
         for(I=0 ; I<length(Vars) ; I++)M=subst(M , TmpVars[I] , GMap[K][I]);
         RM += M;
      }
      RM = simpalg(RM);
      if(RM==0)continue;
      RM = ptozp(RM/length(Elements));
      /* 既出かどうか */
      if( memq(RM , Ss)==1 )continue;
      /* --------------------------- */
      TMP = [];
      for(K = 0 ; K < length(Ss) ; K++){
          if (hdeg(Ss[K] , Vars)==Deg)TMP=cons(Ss[K] , TMP);
      }
      G = gr(append(TMP,GPs) , Vars , 0);
      if( p_nf(RM , G , Vars , 0)!=0 ){
         Ss = cons(RM , Ss);
      }
   }
   return [Ps,Ss];
}


/*
[おまけ]多項式の間に成立する代数的な関係式を決定する

例1:
relation([t^3,t^4,t^5],[t],[x,y,z]);
x=t^3,y=t^4,z=t^5に対して
-z*x+y^2=0
-y*x^2+z^2=0
x^3-z*y=0

例2:
relation([x^3,x*y,y^3] , [x,y] , [s,t,u]);
X=x^3 , Y= x*y, Z=y^3に対して、
X*Y-Z^3=0

例3:
relation([a+b+c , a*b+b*c+c*a , a*b*c] , [a,b,c] , [x,y,z]);
対称多項式a+b+c , a*b+b*c+c*a , a*b*cは代数的に独立

*/
def relation(Exprs , Vars , NewVars){
  /* assert(length(Exprs)==length(NewVars)); */
  S = [];
  for(I = 0 ; I < length(Exprs) ; I++)S = cons(Exprs[I] - NewVars[I] , S);
  Basis = nd_gr(S , append(Vars , NewVars) , 0 , [[0, length(Vars)] , [1,length(NewVars)]]);
  Com = [];
  for(I = 0 ; I < length(Basis) ; I++){
    Check = 1;
    for(J= 0 ; J < length(Vars) ; J++){
       if( deg(Basis[I] , Vars[J]) > 0 ){
          Check = 0;
          break;
       }
    }
    if(Check == 1)Com = cons(Basis[I] , Com);
  }
  return gr(Com , NewVars , 0);
}