留数を計算するアルゴリズム
Higher residueについて眺めていて、多変数留数を計算するアルゴリズムを調べようと思った
Higher Residuesとその応用 (微分方程式の超局所解析)
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/102684
#Higher residueを理解したい動機は、以下のあたり
Descendants constructed from matter fields in topological Landau-Ginzburg theories coupled to topological gravity
http://arxiv.org/abs/hep-th/9211090
Topological Strings, Flat Coordinates and Gravitational Descendants
http://arxiv.org/abs/hep-th/9302048
[一変数留数の場合]
以下、基本的には、有理数関数の場合を想定する。一変数の留数の場合は、部分分数分解してやれば、原理的には気合いで何とかなるはずではある。そうはいっても、コンピュータ上では、普通、有理数係数でやるから例えば、h(z)/(z^3+z+1)とかなってると、因数分解がめんどいとか、細かい問題はある
一応、一変数の場合のアルゴリズムの解説は、以下にあって、
有理関数のローラン展開アルゴリズムと代数的局所コホモロジ.
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/25926
高速留数計算アルゴリズム
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/47838
一変数留数計算アルゴリズムについて
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/58572
D-加群を用いた留数値計算アルゴリズムの局所化
http://www.jssac.org/Editor/Suushiki/V07/No3/V7N3_102.pdf
現在のRisa/Asirには、これを実装したのだと思われるtaji_alc.rrがデフォルトで配布されている
1変数代数的局所コホモロジー類用パッケージ taji_alc
http://www.math.kobe-u.ac.jp/OpenXM/Current/doc/asir-contrib/ja/taji_alc-html/taji_alc-ja.html
[多変数留数のアウトライン]
多変数の場合は、$n$変数多項式の組$f_1, \cdots , f_n$に対して
という形の積分で、多変数留数が定義される。但し、f_1= ... = f_n=0の解集合は0次元で、hは、その解の適当な近傍で正則とする(積分は、あるn次元サイクル上で行う)。留数(local residue)は、f_1, ... , f_nの各共通零点に対して定義される。代数的に定義されたGrothendieck留数というのもあって、上の解析的な留数と一致することが証明されているっぽい(けど、わたしは、証明を知らない)。留数値は、f_1= ... = f_n=0の解集合の各点で定義されることから、f_1, ... , f_nの積が等しくても、留数は異なる値を取りうる
多変数の場合でも、基本になるのは、f_i = (z_i - a_i)の形で、一般に、$f_{i} = (z_{i} - a_{i})^{n_i}$の時、Cauchyの留数公式の多変数版が成立する。f_1= ... = f_n=0の解集合が0次元であるかどうかは、グレブナー基底を使って容易に判定できるし、解がどこにあるかも、イデアルの準素分解によって探すことができる(Risa/Asirでは、primdec.rrをロードして、primadecとかを使えばいい)。けど、多変数の場合は、一変数の場合より厄介なことが起こる。例えば、2次元で、f_1(x,y)=xy , f_2(x,y)=x+yという形の場合、f_1=f_2=0の解は、(x,y)=(0,0)のみであるけど、留数値の計算は、Cauchyの留数公式に頼ることはできない
これを補うために、transformation lawと呼ばれる便利な定理が知られていて、
という形の関係式($A_{ij}$は多項式)で結ばれている多項式の組$f_1 , \cdots , f_n$と$g_1 , \cdots , g_n$が共に、共通零点が0次元になってるとき、
が成り立つ。但し、Resは、global residueで、全てのlocal residueの和。local residueが分からなくても、例えば、完全交叉Artin代数に、Gorenstein環/Frobenius代数の構造を入れる時なんかは、global residueを使うので、local residueは分からなくてもいいという状況はある。上のhigher residueも、その一般化なので、local residueは要らんのだと思う。transformation lawの証明は、やっぱり、わたしは知らない
transformation lawの例として、f_1(x,y) = xy , f_2(x,y)=x+yの場合を考えると
x^2 = -1 * f_1 + x * f_2
y^2 = -1 * f_2 + y * f_2
とかやると、Cauchyの留数公式が使える形に持ちこめる。とはいえ、こういう変換を直接見つけるのは、通常結構しんどい
[多変数のglobal residueを計算するアルゴリズム]
多変数のglobal residueを計算するアルゴリズムは、以下の論文がある
Computing Multidimensional Residues
http://arxiv.org/abs/alg-geom/9404011
アルゴリズムそのものは、17ページに書いてあって、Algorithm(3.1)というのは、$f_1 , \cdots f_n$がある重み付き順序に関するグレブナー基底になってる時のみ使えるアルゴリズム。例えば、f_1=xy , f_2=x+yは、どんな項順序を持ってきても、グレブナー基底にはなってない。グレブナー基底になってない(一般の場合)を扱うのが、Algorithm(3.2)で、これの言ってることは、Buchbergerアルゴリズムの一回のループで、イデアルの生成元$p_1 \cdots , p_n$から、新しい生成元$q_1 , \cdots , q_m$を得た時、注意深く見ると$q_j = \sum A_{ij} p_i$という関係を満たす多項式$A_{ij}$が存在している。これを繰り返してグレブナー基底に至るわけだから、transformation lawを適用して留数の計算ができる(グレブナー基底の個数が、丁度n個になる理由は、イデアルの0次元性判定とか思い出せば分かる)
Algorithm(3.2)は、Buchbergerのアルゴリズムを実装する必要があるのだけど、Risa/AsirでBuchbergerって車輪の再発明ってレベルじゃない感じなので、とりあえず放置して、Algorithm(3.1)だけ実装した。(3.2)は、気が向いたら、そのうちやろうと思う。まぁ、論文にも書いてあるけど、このアルゴリズムは場合によっては、次数が大きいとこの計算が急速にしんどくなるので、実装した甲斐がないなぁという感。アルゴリズム(3.2)を実装してないから、使いにくくてしょうがないし
/*
例1:
residue(y , [x+y , y^2] , [x,y], [2,1]);
例2:
Computing Multidimensional Residues
http://arxiv.org/abs/alg-geom/9404011
の末尾にある例
residue(x^2*y^3*z^2 , [x^5+y^3+z^2-1 , x^2+y^2+z-1 , x^6+y^5+z^3-1] , [x,y,z] , [3,4,7]);
*/
def residue(H , Polys , Vars , Weights){
if(length(Polys)!=length(Vars)){
error("mismatch number of polynomials and variables@residue");
}
if(length(Vars)!=length(Weights)){
error("mismatch number of variables and weights@residue");
}
/* step 1 */
T = uc();
HPolys = newvect(length(Polys));
Degs = newvect(length(Weights));
for(I = 0 ; I < length(Weights) ; I++){
Degs[I] = deg(Polys[I] , Vars[I]) - 1;
HPolys[I] = Polys[I];
for(J = 0 ; J < length(Vars) ; J++){
HPolys[I] = subst(HPolys[I] , Vars[J] , Vars[J]*T^(-Weights[J]));
}
HPolys[I] *= T^(Weights[I]*(Degs[I]+1));
}
/* step 2 */
D_min = 0;
D_max = 0;
for(I = 0 ; I < length(Weights) ; I++){
D_max += Weights[I] * (Degs[I]+1);
D_min += Weights[I] * Degs[I];
}
DS = [];
for(Terms = p_terms(H,Vars,0) ; Terms!=[] ; Terms=cdr(Terms)){
Term = car(Terms);
D = 0;
for(I = 0; I < length(Vars) ; I++){
D += deg(Term , Vars[I]) * Weights[I];
}
DS = cons(D,DS);
}
D = D_min;
for(I = 0 ; I < length(DS) ; I++){
if(DS[I] > D_max)continue;
if(DS[I] > D)D=DS[I];
}
/* step 3 */
B = 1;
for(I = 0 ; I < length(HPolys) ; I++)B *= 1/HPolys[I];
Bt = newvect(D-D_min+1);
for(C=red(B) , N=1 , I=0 ; I<=D - D_min ; I++){
Bt[I] = subst(C , T, 0)/N;
if(I == D - D_min)break;
N *= (I+1);
C = red(diff(C , T));
}
/* step 4 */
Ret = 0;
for(Terms = p_terms(H,Vars,0) ; Terms!=[] ; Terms=cdr(Terms)){
Term = car(Terms);
Indices = newvect(length(Vars));
D = 0;
Coef = H;
for(I = 0; I < length(Vars) ; I++){
Indices[I] = deg(Term , Vars[I]);
D += Indices[I] * Weights[I];
Coef = coef(Coef , Indices[I] , Vars[I]);
}
if (D > D_max || D < D_min)continue;
LP = Bt[D-D_min];
DN = dn(LP);
V = nm(LP);
for(I = 0 ; I < length(Vars) ; I++){
V = coef(V , deg(DN , Vars[I])-Indices[I]-1 , Vars[I]);
DN = coef(DN , deg(DN , Vars[I]) , Vars[I]);
if(V==0)break;
}
Ret += Coef*V/DN;
}
return Ret;
}
[多変数のlocal residueを計算するアルゴリズム]
local residueを計算するアルゴリズムもあるっぽい。
多変数有理関数の留数計算について
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/62806
零次元代数的局所コホモロジー類に付随するホロノミック系の構成アルゴリズム (超局所解析の展望)
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/24911
多変数留数の計算代数解析とホロノミーD加群
http://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/58958
完全なアルゴリズムの説明は、2つ目にある模様。けれど、どう実装すればいいのか、あんまり検討つかないのと、差し当たって必要ないので、とりあえず保留
