=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.4 retrieving revision 1.16 diff -u -p -r1.4 -r1.16 --- OpenXM_contrib2/asir2000/lib/bfct 2000/12/08 08:26:09 1.4 +++ OpenXM_contrib2/asir2000/lib/bfct 2001/01/18 00:52:32 1.16 @@ -45,88 +45,454 @@ * 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.3 2000/08/22 05:04:20 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 ? */ +/* annihilating ideal of F^s */ def ann(F) { V = vars(F); - W = append([y1,y2,t],V); N = length(V); - B = [1-y1*y2,t-y1*F]; + D = newvect(N); + + for ( I = 0; I < N; I++ ) + D[I] = [deg(F,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = N-1; I >= 0; I-- ) + V = cons(D[I][1],V); + for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); + + W = append([y1,y2,t],V); DW = append([dy1,dy2,dt],DV); + + 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,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; } -/* b-function of F ? */ +/* + * compute J_f|s=r, where r = the minimal integral root of global b_f(s) + * ann0(F) returns [MinRoot,Ideal] + */ -def bfct(F) +def ann0(F) { - G4 = ann(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); - N1 = 2*(N+1); + /* XXX : heuristics */ + W = append([y1,y2,t],reverse(V)); + DW = append([dy1,dy2,dt],reverse(DV)); + WDW = append(W,DW); - M = newmat(N1+1,N1); - for ( J = N+1; J < N1; J++ ) + 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); + N = length(V); + B = [t-y1*F]; + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + DW = append([dy1,dt],DV); + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+y1*diff(F,V[I])*dt,B); + } + dp_nelim(1); + + /* homogenized (heuristics) */ + G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6); + G1 = map(subst,G0,y1,1); + G2 = map(psi,G1,t,dt); + G3 = map(subst,G2,t,-s-1); + return G3; +} + +def psi(F,T,DT) +{ + D = dp_ptod(F,[T,DT]); + Wmax = weight(D); + D1 = dp_rest(D); + for ( ; D1; D1 = dp_rest(D1) ) + if ( weight(D1) > Wmax ) + Wmax = weight(D1); + for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) ) + if ( weight(D1) == Wmax ) + Dmax += dp_hm(D1); + if ( Wmax >= 0 ) + Dmax = dp_weyl_mul(<>,Dmax); + else + Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax); + Rmax = dp_dtop(Dmax,[T,DT]); + R = b_subst(subst(Rmax,DT,1),T); + return R; +} + +def weight(D) +{ + V = dp_etov(D); + return V[1]-V[0]; +} + +def compare_first(A,B) +{ + A0 = car(A); + B0 = car(B); + if ( A0 > B0 ) + return 1; + else if ( A0 < B0 ) + return -1; + else + 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 ( J = 0; J < N+1; J++ ) - M[1][J] = 1; -#if 0 - for ( I = 0; I < N+1; I++ ) - M[I+2][N-I] = -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) +{ + V = vars(F); + N = length(V); + D = newvect(N); + for ( I = 0; I < N; I++ ) - M[I+2+N+1][N1-1-I] = -1; -#elif 1 - for ( I = 0; I < N1-1; I++ ) - M[I+2][N1-I-1] = 1; -#else - for ( I = 0; I < N1-1; I++ ) - M[I+2][I] = 1; -#endif + 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); V1 = cons(s,V); DV1 = cons(ds,DV); - G5 = dp_weyl_gr_main(cons(F,G4),append(V1,DV1),0,0,M); - for ( T = G5, G6 = []; T != []; T = cdr(T) ) { - E = car(T); - if ( intersection(vars(E),DV1) == [] ) - G6 = cons(E,G6); + + G0 = indicial1(F,reverse(V)); + G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0); + Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s); + 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); } - G6_0 = remove_zero(map(z_subst,G6,V)); - F0 = flatmf(cdr(fctr(dp_gr_main(G6_0,[s],0,0,0)[0]))); - for ( T = F0, BF = []; T != []; T = cdr(T) ) { - FI = car(T); - for ( J = 1; ; J++ ) { - S = map(srem,map(z_subst,idealquo(G6,[FI^J],V1,0),V),FI); - for ( ; S != [] && !car(S); S = cdr(S) ); - if ( S != [] ) + 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); + Len = length(G); + dp_ord(O); + setmod(M); + PS = newvect(Len); + PS0 = newvect(Len); + + for ( I = 0, T = G; T != []; T = cdr(T), I++ ) + PS0[I] = dp_ptod(car(T),V); + for ( I = 0, T = G; T != []; T = cdr(T), I++ ) + PS[I] = dp_mod(dp_ptod(car(T),V),M,[]); + + for ( I = Len - 1, GI = []; I >= 0; I-- ) + GI = cons(I,GI); + + U = dp_mod(dp_ptod(V0,V),M,[]); + + T = dp_mod(<<0>>,M,[]); + TT = dp_mod(dp_ptod(1,V),M,[]); + 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] ) { + if ( dp_gr_print() ) + print(""); + return dp_dtop(L[1],[t]); /* XXX */ + } else + G = insert(G,L); + } +} + +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); + if ( !valid_modulus(HM,Prime) ) + continue; + MP = weyl_minipolym(G0,V0,O0,Prime,P); + D = deg(MP,var(MP)); + + 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 = 0; + for ( I = 0; I < D; I++ ) + R += B[0][I]*s^I; + R += B[1]*s^D; + return R; + } + } +} + +def weyl_nf(B,G,M,PS) +{ + for ( D = 0; G; ) { + for ( U = 0, L = B; L != []; L = cdr(L) ) { + if ( dp_redble(G,R=PS[car(L)]) > 0 ) { + GCD = igcd(dp_hc(G),dp_hc(R)); + CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD); + U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R); + if ( !U ) + return [D,M]; + D *= CG; M *= CG; break; + } } - BF = cons([FI,J],BF); + if ( U ) + G = U; + else { + D += dp_hm(G); G = dp_rest(G); + } } - return BF; + return [D,M]; +} + +def weyl_nf_mod(B,G,PS,Mod) +{ + for ( D = 0; G; ) { + for ( U = 0, L = B; L != []; L = cdr(L) ) { + if ( dp_redble(G,R=PS[car(L)]) > 0 ) { + CR = dp_hc(G)/dp_hc(R); + U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod); + if ( !U ) + return D; + break; + } + } + if ( U ) + G = U; + else { + D += dp_hm(G); G = dp_rest(G); + } + } + return D; } def remove_zero(L)