M1 セミナー 2000 年 1 月 21 日配布 M1171023 森島聖

1. saturation を計算するプログラム

Saturation[ideal_,f_]:=Saturation[ideal,f]=Module[{a,b,aa},
              a=Groebner[IdealQuotient[ideal,{f}],grlex,
              Variables[IdealQuotient[ideal,{f}]]];
              b=Groebner[IdealQuotient[ideal,{f^2}],grlex,
              Variables[IdealQuotient[ideal,{f^2}]]];n=1;aa=a;
              While[Length[a]=!=Length[b]||
              Do[aa=DeleteCases[aa,b[[i]]],{i,1,Length[b]}];aa=!={},a=b;
              b=Groebner[IdealQuotient[ideal,{f^(n+2)}],grlex,
              Variables[IdealQuotient[ideal,{f^(n+2)}]]];n=n+1;aa=a];{n,a}]

Saturation2[ideal_,f_]:=Saturation2[ideal,f]=
              GroebnerBasis[Append[ideal,1-f*yyy],{},{yyy}]

但し,Saturation のプログラム中の Groebner, IdealQuotient のプログラムは, 注を参照の

こと。

2. [プログラム説明] I = ⟨f_{1}, ⋯, f_{s}⟩ ⊂ k [x_{1}, ⋯, x_{n}] , f ∈ k [x_{1}, ⋯, x_{n}]とするとき, Saturation[{f_{1}, ⋯, f_{s}},f], Saturation2[{f_{1}, ⋯, f_{s}},f] と入力すると,I の f に関しての saturation, I : f^{∞} の生成元 g_{1}, ⋯, g_{t} を出力する。但し,Saturation の方は,I : f^{∞} = I : f^{m} を満たす最小の m もあわせて, {m,{g_{1}, ⋯, g_{t}}} と出力する。

ちなみに,Saturation は,EXERCISE.8-e, Saturation2 は,EXERCISE.9-b の algorithm を用いてプログラムを作成した。

3. [Expample.1] I = ⟨x^{10}⟩ ⊂ k [x], f = x ∈ k [x] とする。

   In[1]:= Saturation[{x^10},x]//Timing

Out[1]= {1.27995 Second, {10, {1}}}

   In[2]:= Saturation2[{x^10},x]//Timing

Out[2]= {0.079997 Second, {1}}

よって, I : f^{∞} = k [x], I : f^{∞} = I : f^{m} を満たす最小の m は 10。

4. [Expample.2] I = ⟨(x + y)^{2}(x − y)(x + z^{2})⟩ ⊂ k [x, y, z], f = (x + z^{2})^{3}(x − y)(z + y) ∈ k [x, y, z] とする。

   In[3]:= Saturation[{(x+y)^2*(x-y)*(x+z^2)},(x+z^2)^3*(x-y)*(z+y)]//Timing

Out[3]= {1.66993 Second, {1, {x^{2} + 2xy + y^{2}}}}

   In[4]:= Saturation2[{(x+y)^2*(x-y)*(x+z^2)},(x+z^2)^3*(x-y)*(z+y)]//Timing

Out[4]= {0.18999 Second, {x^{2} + 2xy + y^{2}}}

よって, I : f^{∞} = ⟨x^{2} + 2xy + y^{2}⟩, I : f^{∞} = I : f^{m} を満たす最小の m は 1。

5. [Expample.3] I = ⟨xz − y^{2}, x^{3} − yz⟩ ⊂ k [x, y, z], f = x + y ∈ k [x, y, z] とする。

   In[5]:= a=Saturation[{x*z-y^2,x^3-y*z},x+y]//Timing

Out[5]= {2.3799 Second, {1, {x^{3} − yz, x^{2}y − z^{2}, −y^{2} + xz, y^{5} − z^{4}, xy^{3} − z^{3}}}}

   In[6]:= b=Saturation2[{x*z-y^2,x^3-y*z},x+y]//Timing

Out[6]= {0.069996 Second, {x^{5} − z^{3}, −x^{3} + yz, x^{2}y − z^{2}, y^{2} − xz}}

   In[7]:= GroebnerBasis[a[[2]][[2]],{x,y,z}]===GroebnerBasis[b[[2]],{x,y,z}]
   Out[7]= True

よって, I : f^{∞} = ⟨x^{5} − z^{3}, −x^{3} + yz, x^{2}y − z^{2}, y^{2} − xz⟩, I : f^{∞} = I : f^{m} を満たす最小の m は 1。

6. [考察] Mathematica で 2 つの algorithm を用いて suturation のプログラムを作成した。 2 つの algorithm の計算スピードを比較したかったが,Saturation の方のプログラムの中の与えられた 2 つのイデアルが等しいかどうかを判定するところで自分の作成した被約 Gröbner 基底を計算するプログラム Groebner を用いたので, そこで極端に計算スピードを費やしてしまい, 比較の対象にならなかった。また, 単項生成でないイデアルに対し,m が 2 以上でかつ saturation が多項式環全体にならないような例が見つけられず Example としてつまらないものばかりになってしまった。

7. [注]

Groebner[generator_,monoorder_,valueorder_]:=
  Groebner[generator,monoorder,valueorder]=
    Module[{LT2,ListPoly,PolyList,Polymono2,PolyDiv2,S2,Groeb3,MinGroeb,
        RedGroeb},
      If[monoorder==lex,
        LT2[f_]:=
          LT2[f]=Module[{s,lt,a},s=Length[f];lt=f[[1]];
              Do[If[DeleteCases[lt[[2]]-f[[a]][[2]],0][[1]]<0,lt=f[[a]]],{a,
                  2,
                  s}];lt]];
      If[monoorder==grlex,
        LT2[f_]:=
          LT2[f]=Module[{s,lt,t,sum1,sum2,a},s=Length[f];lt=f[[1]];
              t=Length[lt[[2]]];
              Do[sum1=Sum[lt[[2]][[i]],{i,t}];

sum2=Sum[f[[a]][[2]][[i]],{i,t}];
If[sum1<sum2,lt=f[[a]],
  If[sum1==sum2,
    If[DeleteCases[lt[[2]]-f[[a]][[2]],0][[1]]<0,
      lt=f[[a]]]]],{a,2,s}];lt]];
If[monoorder==grevlex,
LT2[f_]:=
LT2[f]=Module[{s,lt,t,sum1,sum2,a},s=Length[f];lt=f[[1]];
t=Length[lt[[2]]];
Do[sum1=Sum[lt[[2]][[i]],{i,t}];
sum2=Sum[f[[a]][[2]][[i]],{i,t}];
If[sum1<sum2,lt=f[[a]],
  If[sum1==sum2,
    If[DeleteCases[lt[[2]]-f[[a]][[2]],0][[-1]]>0,
      lt=f[[a]]]]],{a,2,s}];lt]];
If[monoorder=!=lex&&monoorder=!=grlex&&monoorder=!=grevlex,
LT2[f_]:=
LT2[f]=Module[{s,lt,t,c,j},s=Length[f];lt=f[[1]];
t=Length[monoorder];
Do[c=False;j=1;
While[c==False,
  If[monoorder[[j]].lt[[2]]=!=monoorder[[j]].f[[i]][[2]],
    c=True;If[
      monoorder[[j]].lt[[2]]<monoorder[[j]].f[[i]][[2]],
      lt=f[[i]]]];j=j+1],{i,2,s}];lt]];
PolyList[F_,val_]:=
PolyList[F,val]=
Module[{s,PolyListone},
PolyListone[f_,v_]:=
PolyListone[f,v]=
Module[{getcoeff,aa,bb},
  getcoeff[m_,b_]:=
    getcoeff[m,b]=
      Module[{a,c},
        If[Variables[m]=!={},
          a=Exponent[m,b];{Coefficient[m,Times@@(b^a)],a},{m,
            Table[0,{Length[b]}]}]];aa=MonomialList[f];
  bb=Length[aa];Table[getcoeff[aa[[i]],v],{i,1,bb}]];
s=Length[F];Table[PolyListone[F[[i]],val],{i,s}]];
ListPoly[F_,val_]:=
ListPoly[F,val]=
Module[{ListPolyone,s,FF},
ListPolyone[f_,v_]:=
ListPolyone[f,v]=
Module[{s,ss},s=Length[f];ss=Length[v];
  Sum[f[[i]][[1]]*
      Product[v[[j]]^(f[[i]][[2]][[j]]),{j,1,ss}],{i,1,s}]];
s=Length[F];FF={};

Do[FF=Append[FF,ListPolyone[F[[i]],val]],{i,s}];
FF];Polymono2[f_,m_]:=
Polymono2[f,m]=
Module[{ff,s,c,a},ff=f;s=Length[ff];
If[s==0,ff={m},c=True;a=1;
  While[c,If[DeleteCases[m[[2]]-ff[[a]][[2]],0]=={},
      ff=ReplacePart[ff,{ff[[a]][[1]]+m[[1]],m[[2]]},a];
      If[ff[[a]][[1]]==0,ff=Delete[ff,a]];c=False,a=a+1];
    If[a==s+1,ff=Append[ff,m];c=False]]];ff];
PolyDiv2[f_,g_]:=
PolyDiv2[f,g]=
Module[{s,a,r,p,i,c,h1,h2},s=Length[g];a={};Do[a=Append[a,{}],{s}];
r={};p=f;
While[p=!={},i=1;c=True;
  While[i<=s&&c,h1=LT2[p][[2]]-LT2[g[[i]]][[2]];
    If[Select[h1,#1<0&]=={},h2=LT2[p][[1]]/(LT2[g[[i]]][[1]]);
      a=ReplacePart[a,Polymono2[a[[i]],{h2,h1}],i];gg=g[[i]];
      While[gg=!={},
        p=Polymono2[p,{-h2*LT2[gg][[1]],h1+LT2[gg][[2]]}];
        gg=Polymono2[gg,{-LT2[gg][[1]],LT2[gg][[2]]}]];c=False,
      i=i+1]];
  If[c,r=Polymono2[r,LT2[p]];
    p=Polymono2[p,{-LT2[p][[1]],LT2[p][[2]]}]]];{a,r}];
S2[f_,g_]:=
S2[f,g]=Module[{mf,mg,s,c,a,b,z,aa,ff,gg},mf=LT2[f][[2]];
mg=LT2[g][[2]];s=Length[mf];c={};
Do[c=Append[c,Max[{mf[[i]],mg[[i]]}]],{i,1,s}];
a={1/(LT2[f][[1]]),c-mf};b={1/(LT2[g][[1]]),c-mg};z={};aa={};
ff=f;
While[ff=!=z,
  aa=Polymono2[aa,{a[[1]]*LT2[ff][[1]],a[[2]]+LT2[ff][[2]]}];
  ff=Polymono2[ff,{-LT2[ff][[1]],LT2[ff][[2]]}]];gg=g;
While[gg=!=z,
  aa=Polymono2[aa,{-b[[1]]*LT2[gg][[1]],b[[2]]+LT2[gg][[2]]}];
  gg=Polymono2[gg,{-LT2[gg][[1]],LT2[gg][[2]]}]];aa];
Groeb3[F_]:=
Groeb3[F]=
Module[{Criterion,s,B,G,t,ss,k,S},
Criterion[i_,j_,SS_]:=
  Criterion[i,j,SS]=
    Module[{sss,c,cc,ii},sss=Length[G];c=False;cc=True;ii=1;
      While[cc,
        If[(ii=!=i)&&(ii=!=j)&&(
              MemberQ[SS,If[ii>i,{i,ii},{ii,i}]]==False)&&(
              MemberQ[SS,If[ii>j,{j,ii},{ii,j}]]==False)&&(
              Select[
                  Table[Max[{LT2[G[[i]]][[2]][[kk]],
                          LT2[G[[j]]][[2]][[kk]]}],{kk,ss}]-

                    LT2[G[[ii]]][[2]],#1<0&]=={}),c=True;
          cc=False];If[ii==sss,cc=False];ii=ii+1];c];s=Length[F];
B={};
Do[Do[B=Append[B,{i,j}],{j,i+1,s}],{i,1,s-1}];G=F;t=s;
ss=Length[F[[1]][[1]][[2]]];
While[B=!={},k=B[[1]];
  If[(Table[Max[{LT2[G[[k[[1]]]]][[2]][[i]],
                LT2[G[[k[[2]]]]][[2]][[i]]}],{i,ss}]=!=
          LT2[G[[k[[1]]]]][[2]]+LT2[G[[k[[2]]]]][[2]])&&(
        Criterion[k[[1]],k[[2]],B]==False),
    S=PolyDiv2[S2[G[[k[[1]]]],G[[k[[2]]]]],G][[2]];
    If[S=!={},t=t+1;G=Append[G,S];
      Do[B=Append[B,{tt,t}],{tt,1,t-1}]]];B=Delete[B,1]];G];
MinGroeb[F_]:=
MinGroeb[F]=
Module[{G,s,a,h,ss,t,sss,aa},G=Groeb3[F];s=Length[G];
Do[Do[If[i=!=j&&G[[j]]=!=a&&G[[i]]=!=a,
      h=LT2[G[[i]]][[2]]-LT2[G[[j]]][[2]];
      If[Select[h,#1<0&]=={},G=ReplacePart[G,a,i]]];j=j+1,{j,s}];
  i=i+1,{i,s}];G=DeleteCases[G,a];ss=Length[G];
Do[t=LT2[G[[k]]][[1]];
  If[t=!=1,sss=Length[G[[k]]];aa=G[[k]];
    Do[aa=ReplacePart[aa,{aa[[m]][[1]]/t,aa[[m]][[2]]},m],{m,
        sss}];G=ReplacePart[G,aa,k]],{k,ss}];G];
RedGroeb[F_]:=
RedGroeb[F]=
Module[{G,s,k},G=MinGroeb[F];s=Length[G];
Do[k=PolyDiv2[G[[j]],Delete[G,j]][[2]];
  G=ReplacePart[G,k,j],{j,s}];G];
ListPoly[RedGroeb[PolyList[generator,valueorder]],valueorder]]

IdealQuotient[f_,g_]:=
IdealQuotient[f,g]=
Module[{IdealIntersection,s,a},
IdealIntersection[ff_,gg_]:=
IdealIntersection[ff,gg]=
Module[{kk},kk=GroebnerBasis[Join[ttt*ff,(1-ttt)*gg],{},{ttt}];
If[kk=={},kk={0}];kk];s=Length[g];
a=If[g[[1]]===0,{1},
Expand[Simplify[IdealIntersection[f,{g[[1]]}]/g[[1]]]]];
Do[a=IdealIntersection[a,
If[g[[i]]===0,{1},
  Expand[Simplify[IdealIntersection[f,{g[[i]]}]/g[[i]]]]]],{i,2,
s}];Expand[Simplify[a]]]

Groebner に関しては,

http://www.ipc.hiroshima-u.ac.jp/ morishi/groeb3.txt

参照のこと。 IdealQuotient に関しては, 前回の配布プリント参照のこと。