=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.9 retrieving revision 1.16 diff -u -p -r1.9 -r1.16 --- OpenXM_contrib2/asir2000/lib/bfct 2000/12/14 09:36:17 1.9 +++ OpenXM_contrib2/asir2000/lib/bfct 2001/01/18 00:52:32 1.16 @@ -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.15 2001/01/11 08:43:23 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,100 @@ 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; + + /* If W is a list, convert it to a vector */ + if ( type(W) == 4 ) + W = newvect(length(W),W); + dp_weyl_set_weight(W); + + /* create a term order M in 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 */ + 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 */ + 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 */ + dp_gr_flags(["Top",1,"NoRA",1]); + GH = dp_weyl_gr_main(FH,VDVH,0,1,11); + dp_gr_flags(["Top",0,"NoRA",0]); + + /* 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) = 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 +322,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 +375,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(<>,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 +493,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)