=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.18 retrieving revision 1.24 diff -u -p -r1.18 -r1.24 --- OpenXM_contrib2/asir2000/lib/bfct 2002/01/28 02:42:27 1.18 +++ OpenXM_contrib2/asir2000/lib/bfct 2003/04/28 02:15:30 1.24 @@ -45,14 +45,47 @@ * 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.17 2002/01/28 01:02:03 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.23 2003/04/20 11:59:57 noro Exp $ */ /* requires 'primdec' */ +#define TMP_S ssssssss +#define TMP_DS dssssssss +#define TMP_T dtttttttt +#define TMP_DT tttttttt +#define TMP_Y1 yyyyyyyy1 +#define TMP_DY1 dyyyyyyyy1 +#define TMP_Y2 yyyyyyyy2 +#define TMP_DY2 dyyyyyyyy2 + +extern LIBRARY_GR_LOADED$ +extern LIBRARY_PRIMDEC_LOADED$ + +if(!LIBRARY_GR_LOADED) load("gr"); else ; LIBRARY_GR_LOADED = 1$ +if(!LIBRARY_PRIMDEC_LOADED) load("primdec"); else ; LIBRARY_PRIMDEC_LOADED = 1$ + +/* toplevel */ + +def bfunction(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); + return bfct_via_gbfct_weight(F,V); +} + /* annihilating ideal of F^s */ def ann(F) { + if ( member(s,vars(F)) ) + error("ann : the variable 's' is reserved."); V = vars(F); N = length(V); D = newvect(N); @@ -66,12 +99,12 @@ def ann(F) 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); + W = append([TMP_Y1,TMP_Y2,TMP_T],V); + DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV); - B = [1-y1*y2,t-y1*F]; + B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F]; for ( I = 0; I < N; I++ ) { - B = cons(DV[I]+y1*diff(F,V[I])*dt,B); + B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); } /* homogenized (heuristics) */ @@ -80,11 +113,11 @@ def ann(F) G1 = []; for ( T = G0; T != []; T = cdr(T) ) { E = car(T); VL = vars(E); - if ( !member(y1,VL) && !member(y2,VL) ) + if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) ) G1 = cons(E,G1); } - G2 = map(psi,G1,t,dt); - G3 = map(subst,G2,t,-1-s); + G2 = map(psi,G1,TMP_T,TMP_DT); + G3 = map(subst,G2,TMP_T,-1-s); return G3; } @@ -109,13 +142,13 @@ def ann0(F) DV = cons(strtov("d"+rtostr(V[I])),DV); /* XXX : heuristics */ - W = append([y1,y2,t],reverse(V)); - DW = append([dy1,dy2,dt],reverse(DV)); + W = append([TMP_Y1,TMP_Y2,TMP_T],reverse(V)); + DW = append([TMP_DY1,TMP_DY2,TMP_DT],reverse(DV)); WDW = append(W,DW); - B = [1-y1*y2,t-y1*F]; + B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F]; for ( I = 0; I < N; I++ ) { - B = cons(DV[I]+y1*diff(F,V[I])*dt,B); + B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); } /* homogenized (heuristics) */ @@ -124,17 +157,17 @@ def ann0(F) G1 = []; for ( T = G0; T != []; T = cdr(T) ) { E = car(T); VL = vars(E); - if ( !member(y1,VL) && !member(y2,VL) ) + if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) ) G1 = cons(E,G1); } - G2 = map(psi,G1,t,dt); - G3 = map(subst,G2,t,-1-s); + G2 = map(psi,G1,TMP_T,TMP_DT); + G3 = map(subst,G2,TMP_T,-1-TMP_S); /* G3 = J_f(s) */ - V1 = cons(s,V); DV1 = cons(ds,DV); V1DV1 = append(V1,DV1); + V1 = cons(TMP_S,V); DV1 = cons(TMP_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); + Bf = weyl_minipoly(G4,V1DV1,0,TMP_S); FList = cdr(fctr(Bf)); for ( T = FList, Min = 0; T != []; T = cdr(T) ) { @@ -143,30 +176,9 @@ def ann0(F) if ( dn(Root) == 1 && Root < Min ) Min = Root; } - return [Min,map(subst,G3,s,Min)]; + return [Min,map(subst,G3,TMP_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]); @@ -367,6 +379,9 @@ def initial_part(F,V,MW,W) def bfct(F) { + /* XXX */ + F = replace_vars_f(F); + V = vars(F); N = length(V); D = newvect(N); @@ -386,6 +401,29 @@ def bfct(F) return Minipoly; } +/* called from bfct() only */ + +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; +} + /* b-function computation via generic_bfct() (experimental) */ def bfct_via_gbfct(F) @@ -403,14 +441,14 @@ def bfct_via_gbfct(F) for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); - B = [t-F]; + B = [TMP_T-F]; for ( I = 0; I < N; I++ ) { - B = cons(DV[I]+diff(F,V[I])*dt,B); + B = cons(DV[I]+diff(F,V[I])*TMP_DT,B); } - V1 = cons(t,V); DV1 = cons(dt,DV); + V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV); W = newvect(N+1); W[0] = 1; - R = generic_bfct_1(B,V1,DV1,W); + R = generic_bfct(B,V1,DV1,W); return subst(R,s,-s-1); } @@ -440,11 +478,11 @@ def bfct_via_gbfct_weight(F,V) for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); - B = [t-F]; + B = [TMP_T-F]; for ( I = 0; I < N; I++ ) { - B = cons(DV[I]+diff(F,V[I])*dt,B); + B = cons(DV[I]+diff(F,V[I])*TMP_DT,B); } - V1 = cons(t,V); DV1 = cons(dt,DV); + V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV); W = newvect(N+1); W[0] = 1; R = generic_bfct_1(B,V1,DV1,W); @@ -477,18 +515,114 @@ def bfct_via_gbfct_weight_1(F,V) for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); - B = [t-F]; + B = [TMP_T-F]; for ( I = 0; I < N; I++ ) { - B = cons(DV[I]+diff(F,V[I])*dt,B); + B = cons(DV[I]+diff(F,V[I])*TMP_DT,B); } - V1 = append(V,[t]); DV1 = append(DV,[dt]); + V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]); W = newvect(N+1); W[N] = 1; - R = generic_bfct(B,V1,DV1,W); + R = generic_bfct_1(B,V1,DV1,W); dp_set_weight(0); return subst(R,s,-s-1); } +def bfct_via_gbfct_weight_2(F,V) +{ + N = length(V); + D = newvect(N); + Wt = getopt(weight); + if ( type(Wt) != 4 ) { + for ( I = 0, Wt = []; I < N; I++ ) + Wt = cons(1,Wt); + } + Tdeg = w_tdeg(F,V,Wt); + + /* a weight for the first GB computation */ + /* [t,x1,...,xn,dt,dx1,...,dxn,h] */ + WtV = newvect(2*(N+1)+1); + WtV[0] = Tdeg; + WtV[N+1] = 1; + WtV[2*(N+1)] = 1; + /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */ + for ( I = 1; I <= N; I++ ) { + WtV[I] = Wt[I-1]; + WtV[N+1+I] = Tdeg-Wt[I-1]+1; + } + dp_set_weight(WtV); + + /* a weight for the second GB computation */ + /* [x1,...,xn,t,dx1,...,dxn,dt,h] */ + WtV2 = newvect(2*(N+1)+1); + WtV2[N] = Tdeg; + WtV2[2*N+1] = 1; + WtV2[2*(N+1)] = 1; + for ( I = 0; I < N; I++ ) { + WtV2[I] = Wt[I]; + WtV2[N+1+I] = Tdeg-Wt[I]+1; + } + + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + B = [TMP_T-F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+diff(F,V[I])*TMP_DT,B); + } + V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV); + V2 = append(V,[TMP_T]); DV2 = append(DV,[TMP_DT]); + W = newvect(N+1,[1]); + dp_weyl_set_weight(W); + + VDV = append(V1,DV1); + N1 = length(V1); + N2 = N1*2; + + /* create a non-term order MW in D */ + MW = newmat(N2+1,N2); + for ( J = 0; J < N1; J++ ) { + MW[0][J] = -W[J]; MW[0][N1+J] = W[J]; + } + for ( J = 0; J < N2; J++ ) MW[1][J] = 1; + for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1; + + /* homogenize F */ + VDVH = append(VDV,[h]); + FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH); + + /* compute a groebner basis of FH w.r.t. MWH */ + GH = dp_weyl_gr_main(FH,VDVH,0,1,11); + + /* 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 < N1; I++ ) + T += W[I]*V1[I]*DV1[I]; + + /* change of ordering from VDV to VDV2 */ + VDV2 = append(V2,DV2); + dp_set_weight(WtV2); + for ( Pind = 0; ; Pind++ ) { + Prime = lprime(Pind); + GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-Prime,0); + if ( GIN2 ) break; + } + + R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */ + dp_set_weight(0); + return subst(R,s,-s-1); +} + +/* minimal polynomial of s; modular computation */ + def weyl_minipolym(G,V,O,M,V0) { N = length(V); @@ -530,6 +664,8 @@ def weyl_minipolym(G,V,O,M,V0) } } +/* minimal polynomial of s over Q */ + def weyl_minipoly(G0,V0,O0,P) { HM = hmlist(G0,V0,O0); @@ -542,20 +678,28 @@ def weyl_minipoly(G0,V0,O0,P) PS[I] = dp_ptod(car(T),V0); for ( I = Len - 1, GI = []; I >= 0; I-- ) GI = cons(I,GI); + PSM = newvect(Len); DP = dp_ptod(P,V0); - for ( I = 0; ; I++ ) { - Prime = lprime(I); + for ( Pind = 0; ; Pind++ ) { + Prime = lprime(Pind); if ( !valid_modulus(HM,Prime) ) continue; - MP = weyl_minipolym(G0,V0,O0,Prime,P); - D = deg(MP,var(MP)); + setmod(Prime); + for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ ) + PSM[I] = dp_mod(dp_ptod(car(T),V0),Prime,[]); NFP = weyl_nf(GI,DP,1,PS); + NFPM = dp_mod(NFP[0],Prime,[])/ptomp(NFP[1],Prime); + NF = [[dp_ptod(1,V0),1]]; LCM = 1; - for ( J = 1; J <= D; J++ ) { + TM = dp_mod(<<0>>,Prime,[]); + TTM = dp_mod(dp_ptod(1,V0),Prime,[]); + GM = NFM = [[TTM,TM]]; + + for ( D = 1; ; D++ ) { if ( dp_gr_print() ) print(".",2); NFPrev = car(NF); @@ -564,7 +708,18 @@ def weyl_minipoly(G0,V0,O0,P) NFJ = remove_cont(NFJ); NF = cons(NFJ,NF); LCM = ilcm(LCM,NFJ[1]); + + /* modular computation */ + TM = dp_mod(<>,Prime,[]); + TTM = dp_mod(NFJ[0],Prime,[])/ptomp(NFJ[1],Prime); + NFM = cons([TTM,TM],NFM); + LM = dp_lnf_mod([TTM,TM],GM,Prime); + if ( !LM[0] ) + break; + else + GM = insert(GM,LM); } + if ( dp_gr_print() ) print(""); U = NF[0][0]*idiv(LCM,NF[0][1]); @@ -656,13 +811,6 @@ def flatmf(L) { return S; } -def member(A,L) { - for ( ; L != []; L = cdr(L) ) - if ( A == car(L) ) - return 1; - return 0; -} - def intersection(A,B) { for ( L = []; A != []; A = cdr(A) ) @@ -697,6 +845,33 @@ def w_tdeg(F,V,W) for ( R = 0; T; T = cdr(T) ) { D = dp_td(T); if ( D > R ) R = D; + } + return R; +} + +def replace_vars_f(F) +{ + return subst(F,s,TMP_S,t,TMP_T,y1,TMP_Y1,y2,TMP_Y2); +} + +def replace_vars_v(V) +{ + V = replace_var(V,s,TMP_S); + V = replace_var(V,t,TMP_T); + V = replace_var(V,y1,TMP_Y1); + V = replace_var(V,y2,TMP_Y2); + return V; +} + +def replace_var(V,X,Y) +{ + V = reverse(V); + for ( R = []; V != []; V = cdr(V) ) { + Z = car(V); + if ( Z == X ) + R = cons(Y,R); + else + R = cons(Z,R); } return R; }