===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v
retrieving revision 1.9
retrieving revision 1.14
diff -u -p -r1.9 -r1.14
--- OpenXM_contrib2/asir2000/lib/bfct 2000/12/14 09:36:17 1.9
+++ OpenXM_contrib2/asir2000/lib/bfct 2001/01/10 04:30:35 1.14
@@ -45,8 +45,8 @@
* DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
* PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
*
- * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.8 2000/12/14 09:13:37 noro Exp $
-*/
+ * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.13 2000/12/27 07:17:39 noro Exp $
+ */
/* requires 'primdec' */
/* annihilating ideal of F^s */
@@ -73,20 +73,79 @@ def ann(F)
for ( I = 0; I < N; I++ ) {
B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
}
+
+ /* homogenized (heuristics) */
dp_nelim(2);
- G0 = dp_weyl_gr_main(B,append(W,DW),0,0,6);
+ G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
G1 = [];
for ( T = G0; T != []; T = cdr(T) ) {
E = car(T); VL = vars(E);
if ( !member(y1,VL) && !member(y2,VL) )
G1 = cons(E,G1);
}
- G2 = map(subst,G1,dt,1);
- G3 = map(b_subst,G2,t);
- G4 = map(subst,G3,t,-1-s);
- return G4;
+ G2 = map(psi,G1,t,dt);
+ G3 = map(subst,G2,t,-1-s);
+ return G3;
}
+/*
+ * compute J_f|s=r, where r = the minimal integral root of global b_f(s)
+ * ann0(F) returns [MinRoot,Ideal]
+ */
+
+def ann0(F)
+{
+ V = vars(F);
+ N = length(V);
+ D = newvect(N);
+
+ for ( I = 0; I < N; I++ )
+ D[I] = [deg(F,V[I]),V[I]];
+ qsort(D,compare_first);
+ for ( V = [], I = 0; I < N; I++ )
+ V = cons(D[I][1],V);
+
+ for ( I = N-1, DV = []; I >= 0; I-- )
+ DV = cons(strtov("d"+rtostr(V[I])),DV);
+
+ /* XXX : heuristics */
+ W = append([y1,y2,t],reverse(V));
+ DW = append([dy1,dy2,dt],reverse(DV));
+ WDW = append(W,DW);
+
+ B = [1-y1*y2,t-y1*F];
+ for ( I = 0; I < N; I++ ) {
+ B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
+ }
+
+ /* homogenized (heuristics) */
+ dp_nelim(2);
+ G0 = dp_weyl_gr_main(B,WDW,1,0,6);
+ G1 = [];
+ for ( T = G0; T != []; T = cdr(T) ) {
+ E = car(T); VL = vars(E);
+ if ( !member(y1,VL) && !member(y2,VL) )
+ G1 = cons(E,G1);
+ }
+ G2 = map(psi,G1,t,dt);
+ G3 = map(subst,G2,t,-1-s);
+
+ /* G3 = J_f(s) */
+
+ V1 = cons(s,V); DV1 = cons(ds,DV); V1DV1 = append(V1,DV1);
+ G4 = dp_weyl_gr_main(cons(F,G3),V1DV1,0,1,0);
+ Bf = weyl_minipoly(G4,V1DV1,0,s);
+
+ FList = cdr(fctr(Bf));
+ for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
+ LF = car(car(T));
+ Root = -coef(LF,0)/coef(LF,1);
+ if ( dn(Root) == 1 && Root < Min )
+ Min = Root;
+ }
+ return [Min,map(subst,G3,s,Min)];
+}
+
def indicial1(F,V)
{
W = append([y1,t],V);
@@ -99,10 +158,10 @@ def indicial1(F,V)
B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
}
dp_nelim(1);
- /* we use homogenization (heuristically determined) */
+
+ /* homogenized (heuristics) */
G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
G1 = map(subst,G0,y1,1);
- Mat = newmat(2,2,[[-1,1],[0,1]]);
G2 = map(psi,G1,t,dt);
G3 = map(subst,G2,t,-s-1);
return G3;
@@ -146,6 +205,93 @@ def compare_first(A,B)
return 0;
}
+/* generic b-function w.r.t. weight vector W */
+
+def generic_bfct(F,V,DV,W)
+{
+ N = length(V);
+ N2 = N*2;
+
+ /* create a term order M in D<x,d> (DRL) */
+ M = newmat(N2,N2);
+ for ( J = 0; J < N2; J++ )
+ M[0][J] = 1;
+ for ( I = 1; I < N2; I++ )
+ M[I][N2-I] = -1;
+
+ VDV = append(V,DV);
+
+ /* create a non-term order MW in D<x,d> */
+ MW = newmat(N2+1,N2);
+ for ( J = 0; J < N; J++ )
+ MW[0][J] = -W[J];
+ for ( ; J < N2; J++ )
+ MW[0][J] = W[J-N];
+ for ( I = 1; I <= N2; I++ )
+ for ( J = 0; J < N2; J++ )
+ MW[I][J] = M[I-1][J];
+
+ /* create a homogenized term order MWH in D<x,d,h> */
+ MWH = newmat(N2+2,N2+1);
+ for ( J = 0; J <= N2; J++ )
+ MWH[0][J] = 1;
+ for ( I = 1; I <= N2+1; I++ )
+ for ( J = 0; J < N2; J++ )
+ MWH[I][J] = MW[I-1][J];
+
+ /* homogenize F */
+ VDVH = append(VDV,[h]);
+ FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
+
+ /* compute a groebner basis of FH w.r.t. MWH */
+ GH = dp_weyl_gr_main(FH,VDVH,0,1,MWH);
+
+ /* dehomigenize GH */
+ G = map(subst,GH,h,1);
+
+ /* G is a groebner basis w.r.t. a non term order MW */
+ /* take the initial part w.r.t. (-W,W) */
+ GIN = map(initial_part,G,VDV,MW,W);
+
+ /* GIN is a groebner basis w.r.t. a term order M */
+ /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
+
+ /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
+ for ( I = 0, T = 0; I < N; I++ )
+ T += W[I]*V[I]*DV[I];
+ B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
+ return B;
+}
+
+def initial_part(F,V,MW,W)
+{
+ N2 = length(V);
+ N = N2/2;
+ dp_ord(MW);
+ DF = dp_ptod(F,V);
+ R = dp_hm(DF);
+ DF = dp_rest(DF);
+
+ E = dp_etov(R);
+ for ( I = 0, TW = 0; I < N; I++ )
+ TW += W[I]*(-E[I]+E[N+I]);
+ RW = TW;
+
+ for ( ; DF; DF = dp_rest(DF) ) {
+ E = dp_etov(DF);
+ for ( I = 0, TW = 0; I < N; I++ )
+ TW += W[I]*(-E[I]+E[N+I]);
+ if ( TW == RW )
+ R += dp_hm(DF);
+ else if ( TW < RW )
+ break;
+ else
+ error("initial_part : cannot happen");
+ }
+ return dp_dtop(R,V);
+
+}
+
/* b-function of F ? */
def bfct(F)
@@ -169,6 +315,35 @@ def bfct(F)
return Minipoly;
}
+/* b-function computation via generic_bfct() (experimental) */
+
+def bfct_via_gbfct(F)
+{
+ V = vars(F);
+ N = length(V);
+ D = newvect(N);
+
+ for ( I = 0; I < N; I++ )
+ D[I] = [deg(F,V[I]),V[I]];
+ qsort(D,compare_first);
+ for ( V = [], I = 0; I < N; I++ )
+ V = cons(D[I][1],V);
+ V = reverse(V);
+ for ( I = N-1, DV = []; I >= 0; I-- )
+ DV = cons(strtov("d"+rtostr(V[I])),DV);
+
+ B = [t-F];
+ for ( I = 0; I < N; I++ ) {
+ B = cons(DV[I]+diff(F,V[I])*dt,B);
+ }
+ V1 = cons(t,V); DV1 = cons(dt,DV);
+ W = newvect(N+1);
+ W[0] = 1;
+ R = generic_bfct(B,V1,DV1,W);
+
+ return subst(R,s,-s-1);
+}
+
def weyl_minipolym(G,V,O,M,V0)
{
N = length(V);
@@ -193,70 +368,82 @@ def weyl_minipolym(G,V,O,M,V0)
G = H = [[TT,T]];
for ( I = 1; ; I++ ) {
+ if ( dp_gr_print() )
+ print(".",2);
T = dp_mod(<<I>>,M,[]);
TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M);
H = cons([TT,T],H);
L = dp_lnf_mod([TT,T],G,M);
- if ( !L[0] )
- return dp_dtop(L[1],[V0]);
- else
+ if ( !L[0] ) {
+ if ( dp_gr_print() )
+ print("");
+ return dp_dtop(L[1],[t]); /* XXX */
+ } else
G = insert(G,L);
}
}
-def weyl_minipoly(G0,V0,O0,V)
+def weyl_minipoly(G0,V0,O0,P)
{
+ HM = hmlist(G0,V0,O0);
+
+ N = length(V0);
+ Len = length(G0);
+ dp_ord(O0);
+ PS = newvect(Len);
+ for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
+ PS[I] = dp_ptod(car(T),V0);
+ for ( I = Len - 1, GI = []; I >= 0; I-- )
+ GI = cons(I,GI);
+ DP = dp_ptod(P,V0);
+
for ( I = 0; ; I++ ) {
Prime = lprime(I);
- MP = weyl_minipolym(G0,V0,O0,Prime,V);
- for ( D = deg(MP,V), TL = [], J = 0; J <= D; J++ )
- TL = cons(V^J,TL);
- dp_ord(O0);
- NF = weyl_gennf(G0,TL,V0,O0)[0];
+ if ( !valid_modulus(HM,Prime) )
+ continue;
+ MP = weyl_minipolym(G0,V0,O0,Prime,P);
+ D = deg(MP,var(MP));
- LHS = weyl_nf_tab(-car(TL),NF,V0);
- B = weyl_hen_ttob(cdr(TL),NF,LHS,V0,Prime);
+ NFP = weyl_nf(GI,DP,1,PS);
+ NF = [[dp_ptod(1,V0),1]];
+ LCM = 1;
+
+ for ( J = 1; J <= D; J++ ) {
+ if ( dp_gr_print() )
+ print(".",2);
+ NFPrev = car(NF);
+ NFJ = weyl_nf(GI,
+ dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS);
+ NFJ = remove_cont(NFJ);
+ NF = cons(NFJ,NF);
+ LCM = ilcm(LCM,NFJ[1]);
+ }
+ if ( dp_gr_print() )
+ print("");
+ U = NF[0][0]*idiv(LCM,NF[0][1]);
+ Coef = [];
+ for ( J = D-1; J >= 0; J-- ) {
+ Coef = cons(strtov("u"+rtostr(J)),Coef);
+ U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]);
+ }
+
+ for ( UU = U, Eq = []; UU; UU = dp_rest(UU) )
+ Eq = cons(dp_hc(UU),Eq);
+ M = etom([Eq,Coef]);
+ B = henleq(M,Prime);
+ if ( dp_gr_print() )
+ print("");
if ( B ) {
- R = ptozp(B[1]*car(TL)+B[0]);
+ R = 0;
+ for ( I = 0; I < D; I++ )
+ R += B[0][I]*s^I;
+ R += B[1]*s^D;
return R;
}
}
}
-def weyl_gennf(G,TL,V,O)
-{
- N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
- for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
- PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
- }
- for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
- DTL = cons(dp_ptod(car(TL),V),DTL);
- for ( I = Len - 1, GI = []; I >= 0; I-- )
- GI = cons(I,GI);
- T = car(DTL); DTL = cdr(DTL);
- H = [weyl_nf(GI,T,T,PS)];
-
- T0 = time()[0];
- while ( DTL != [] ) {
- T = car(DTL); DTL = cdr(DTL);
- if ( dp_gr_print() )
- print(".",2);
- if ( L = search_redble(T,H) ) {
- DD = dp_subd(T,L[1]);
- NF = weyl_nf(GI,dp_weyl_mul(L[0],dp_subd(T,L[1])),dp_hc(L[1])*T,PS);
- } else
- NF = weyl_nf(GI,T,T,PS);
- NF = remove_cont(NF);
- H = cons(NF,H);
- }
- if ( dp_gr_print() ) print("");
- TNF = time()[0]-T0;
- if ( dp_gr_print() )
- print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
- return [H,PS,GI];
-}
-
def weyl_nf(B,G,M,PS)
{
for ( D = 0; G; ) {
@@ -299,78 +486,6 @@ def weyl_nf_mod(B,G,PS,Mod)
}
}
return D;
-}
-
-def weyl_hen_ttob(T,NF,LHS,V,MOD)
-{
- T0 = time()[0]; M = etom(weyl_leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
- T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
- if ( dp_gr_print() ) {
- print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
- }
- return U ? vtop(T,U,LHS) : 0;
-}
-
-def weyl_leq_nf(TL,NF,LHS,V)
-{
- TLen = length(NF);
- T = newvect(TLen); M = newvect(TLen);
- for ( I = 0; I < TLen; I++ ) {
- T[I] = dp_ht(NF[I][1]);
- M[I] = dp_hc(NF[I][1]);
- }
- Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
- for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
- D = dp_ptod(car(L),V);
- for ( I = 0; I < TLen; I++ )
- if ( D == T[I] )
- break;
- INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
- }
- if ( !LHS ) {
- COEF[0] = 1; NM = 0; DN = 1;
- } else {
- NM = LHS[0]; DN = LHS[1];
- }
- for ( J = 0, S = -NM; J < Len; J++ ) {
- DNJ = M[INDEX[J]];
- GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
- S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
- DN *= CS;
- }
- for ( D = S, E = []; D; D = dp_rest(D) )
- E = cons(dp_hc(D),E);
- BOUND = LHS ? 0 : 1;
- for ( I = Len - 1, W = []; I >= BOUND; I-- )
- W = cons(COEF[I],W);
- return [E,W];
-}
-
-def weyl_nf_tab(A,NF,V)
-{
- TLen = length(NF);
- T = newvect(TLen); M = newvect(TLen);
- for ( I = 0; I < TLen; I++ ) {
- T[I] = dp_ht(NF[I][1]);
- M[I] = dp_hc(NF[I][1]);
- }
- A = dp_ptod(A,V);
- for ( Z = A, Len = 0; Z; Z = dp_rest(Z), Len++ );
- INDEX = newvect(Len); COEF = newvect(Len);
- for ( Z = A, J = 0; Z; Z = dp_rest(Z), J++ ) {
- D = dp_ht(Z);
- for ( I = 0; I < TLen; I++ )
- if ( D == T[I] )
- break;
- INDEX[J] = I; COEF[J] = dp_hc(Z);
- }
- for ( J = 0, S = 0, DN = 1; J < Len; J++ ) {
- DNJ = M[INDEX[J]];
- GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
- S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
- DN *= CS;
- }
- return [S,DN];
}
def remove_zero(L)