File: [local] / OpenXM / src / k097 / debug / ahg2.k (download)
Revision 1.3, Mon Jan 8 05:26:51 2001 UTC (23 years, 8 months ago) by takayama
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_3_1_13b, RELEASE_1_2_3_12, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, KNOPPIX_2006, HEAD, DEB_REL_1_2_3-9 Changes since 1.2: +1 -1
lines
* Cleaned unnecessary files and functions.
* Completed the new help system. Type in make in the directory
OpenXM/src/k097/Doc
|
/* SSWork/yacc/debug/ahg.k
/* cf. debug/toric0.k */
/* toric の generator を求める関数.
A-hypergeometric の indicial ideal を求める関数.
*/
load("indexed.k");
ShimomuraSpecial = true ; OnePath= true;
Vvv = false;
SetRingVariables_Verbose = false;
def void QuietKan() {
sm1(" [(KanGBmessage) 0] system_variable ");
}
def testhg1() {
a = [[1,1,1,1,1,1],
[0,0,0,1,1,1],
[0,1,0,0,1,0],
[0,0,1,0,0,1]];
return(idhg(a));
}
def testhg2() {
a = [[1,1,1,1,1],
[0,2,3,4,3],
[0,1,1,0,2]];
return(idhg(a));
}
def idhg(a) {
local a,ans,rd,i,ans2,ans3,n,ff,d,zlist;
ans = toric(a);
if (ShimomuraSpecial) {
/* 先に, toric の initial part をとってしまう. */
/* 本当は, FW の initial part をとるべきなのかも? */
if (Vvv) {Println("-------- S-special ---------");}
ans = Map(ans,"Init");
}
ans = Map(ans,"ToString");
if (Vvv) {Println(ans);}
rd = RingDonIndexedVariables("z",Length(a[0])+1+Length(a));
/* 4 秒程度かかる. */
ans = Map(ans,"Poly");
n = Length(a[0]); d = Length(a);
ans2 = NewArray(Length(ans)); /* ans2 には, toric */
PSfor (i=0; i< Length(ans); i++) {
ans2[i] = ztoDz(ans[i],n);
}
if (Vvv) {Println(ans2);}
ans3 = atolin(a); /* ans3 には, 一次式 */
if (Vvv) {Println(ans3);}
ff = Map(Join(ans2,ans3),"ToString");
ans = zindicial(ff,n,d);
zlist = [ ];
PSfor(i= n; i<n+d+1; i++) {
zlist = Append(zlist,Indexed("z",i));
}
return([ans,zlist]);
}
def toric0_toMonom(aa,i,offset, ring)
{
local j,ans,m;
m = Length(aa);
ans = PolyR("1",ring);
for (j=0; j<m; j++) {
ans = ans*(z[offset+j]^(aa[j,i]));
}
return(ans);
}
def toric(aa) {
local i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
d = Length(aa); n = Length(aa[0]);
if (Vvv) {Println(aa);}
weight = [ ]; elim = [ ];
PSfor (i= n; i< n+d; i++) {
weight = Join(weight,[Indexed("z",i), 1]);
elim = Append(elim, Indexed("z",i));
}
weight = Append([weight],[Indexed("z",n-1),1]);
if (Vvv) {Println(weight); Println(elim);}
rz = RingPonIndexedVariables("z",n+d, weight); /* 4 秒程度かかる. */
/* z[0], ..., z[n-1], ... , z[n+d]*/
ideal = [ ];
PSfor (i=0; i< n; i++) {
ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,n,rz));
}
if (Vvv) {
Println(" --------- input ideal -------------");
Print(" z["); Print( n ); Print( "] --- z["); Print( n+d-1);
Println("] should be eliminated.");
Println(ideal);
}
ans = Groebner(ideal);
if (Vvv) {Println(" -------------- gb is ----------------- "); Println(ans);}
ans = Eliminatev(ans,elim);
if (Vvv) {Println(" ------------ eliminated -------------- "); Println(ans);}
rule = [[h, PolyR("1",rz) ] ];
nn = Length(ans); univ = [ ];
PSfor (i=0; i<nn; i++) {
univ = Append(univ,Replace(ans[i],rule));
}
ans = ReducedBase(univ);
if (Vvv) {
Println(" ----------- removed redundant elements ----------- ");
Println(" ---------- generators of the toric ideal are ----- ");
Println(ans);
Println(" ");
}
return(ans);
}
def zindicial0(input,n,m) {
local rz,weight, ww,i,rule,zinverse,m,d,ans,elim,tmp;
if (!OnePath) {
ww = [ ]; elim = [ ];
weight = [[Indexed("z",n),1]];
if (Vvv) {Print("-------- weight ---------: "); Println(weight);}
rz = RingDonIndexedVariables("z",n+1+m, weight);
z = NewArray(n+1+m); Dz = NewArray(n+1+m);
PSfor(i=0; i<n+1+m; i++) {
z[i] = PolyR(Indexed("z",i),rz);
Dz[i] = PolyR(Indexed("Dz",i),rz);
}
input = Mapto(input,rz);
if (Vvv) {Println("------------ input ------------"); Println(input);}
/* F-homogenization. z[0], ..., z[n-1],
z[n] is the homogenization variable*/
/* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
input = Replace(input,rule);
m = Length(input);
PSfor (i=0; i<m; i++) {
d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
if (d < 0) {
input[i] = z[n]^(-d)*input[i];
}
}
if (Vvv) {Print("------ input : "); Println(input);}
ans = GroebnerTime(input);
m = Length(ans);
PSfor (i=0; i<m; i++) {
/* FW principal part をとる. */
/* ans[i] = Init_w(ans[i],[z[n]],[1]); この関数は遅い. */
tmp = Coefficients(ans[i],z[n]);
ans[i] = tmp[1,0];
}
if (Vvv) {Print("--------FW principal parts : ");Println(ans);}
/* もう一回, GB の計算. */
input = Map(ans,"ToString");
}
/* もう一回, GB の計算. */
ww = [ ]; elim = [ ];
/* z[n+1], ..., z[n+m] がパラメータ変数 */
/* Dz[0] --- Dz[n-2], z[0] --- z[n-2] を消去する. */
PSfor (i=0; i<n-1; i++) {
ww = Join(ww,[Indexed("Dz",i), 1]);
/* ww = Join(ww,[Indexed("z",i), 1]); */
if (i != n-1) {
elim = Append(elim,Indexed("Dz",i));
/* elim = Append(elim,Indexed("z",i)); */
}
}
weight = [[Indexed("z",n),1] , ww];
if (Vvv) {Print("-------- weight ---------: "); Println(weight);}
rz = RingDonIndexedVariables("z",n+1+m, weight);
z = NewArray(n+1+m); Dz = NewArray(n+1+m);
PSfor(i=0; i<n+1+m; i++) { /* これをもう一度やらないと, mklib で sm1 に
したときのみ z が undefined
になる. どうしてか? */
z[i] = PolyR(Indexed("z",i),rz);
Dz[i] = PolyR(Indexed("Dz",i),rz);
}
input = Mapto(input,rz);
if (OnePath) {
/* F-homogenization. z[0], ..., z[n-1],
z[n] is the homogenization variable*/
/* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
input = Replace(input,rule);
m = Length(input);
PSfor (i=0; i<m; i++) {
d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
if (d < 0) {
input[i] = z[n]^(-d)*input[i];
}
}
}
if (Vvv) {Print("------ input : "); Println(input);}
ans = GroebnerTime(input);
m = Length(ans);
PSfor (i=0; i<m; i++) {
/* FW principal part をとる. */
/* ans[i] = Init_w(ans[i],[PolyR(Indexed("z",n),rz)],[1]); */
tmp = Coefficients(ans[i],z[n]);
ans[i] = tmp[1,0];
}
if (Vvv) {Print("--------FW principal parts : ");Println(ans);}
ans = Eliminatev(ans,elim);
m = Length(ans);
/* h, z[n] を 1 にする. */
for (i=0; i<m; i++) {
ans[i] = Replace(ans[i],[[h,PolyR("1",rz)],[PolyR(Indexed("z",n),rz),PolyR("1",rz)]]);
}
if (Vvv) {Println(" "); Println(" ");}
return(ans);
}
def zrho(f,n) {
local ans,i,top,w,rz;
ans = 0;
rz = GetRing(f);
while(true) {
if ( f == Poly("0")) sm1(" exit ");
top = Init(f);
f = f-top;
w = Exponent(top,[PolyR(Indexed("Dz",n-1),rz)]);
top = Replace(top,[[PolyR(Indexed("Dz",n-1),rz),PolyR("1",rz)]])*zipoch(z[n],w[0]);
ans = ans + top;
}
return(ans);
}
def zipoch(f,w) {
local ans,i;
ans = 1;
PSfor (i=0; i<w; i++) {
ans = ans*(f-i);
}
return(ans);
}
def zindicial(fff,n,mm) {
local ans,n,i,m,r,tmp;
ans = zindicial0(fff,n,mm);
if (Vvv) {Println(ans);}
m = Length(ans);
r = [ ];
if (Vvv) {
Println(AddString(["------ The generic indicial polynomial along z[",
ToString(n-1),
"] = 0 is the minimal degree polynomial of the following",
"polynomials."]));
Println(AddString(["z[",ToString(n),"] is equal to s."]));
}
PSfor (i=0; i<m; i++) {
tmp = ans[i];
tmp = Replace(tmp,[[Poly(Indexed("z",n-1)),Poly("1")]]);
tmp = zrho(tmp,n);
if (Vvv) {Print(i); Print(" : ");Println(tmp);}
r = Append(r,tmp);
}
if (Vvv) {Println(" ");}
return(r);
}
def ztoDz(f,n) {
local rule,i;
rule = NewArray(n);
PSfor(i=0; i<n; i++) {
rule[i] = [z[i],Dz[i]];
}
return(Replace(f,rule));
}
/* 行列より A-HG の線形方程式をだす. */
def atolin(a) {
local d,n,eqs,ans,i,j;
d = Length(a);
n = Length(a[0]);
eqs = NewArray(d);
PSfor (i=0; i<d; i++) {
ans = 0;
PSfor (j=0; j<n; j++) {
ans = ans + a[i,j]*z[j]*Dz[j];
}
ans = ans - z[n+1+i];
eqs[i] = ans;
}
return(eqs);
}
sm1(" ; ");