=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.15 retrieving revision 1.22 diff -u -p -r1.15 -r1.22 --- OpenXM_contrib2/asir2000/lib/bfct 2001/01/11 08:43:23 1.15 +++ OpenXM_contrib2/asir2000/lib/bfct 2003/04/20 08:54:28 1.22 @@ -45,10 +45,32 @@ * 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.14 2001/01/10 04:30:35 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.21 2002/01/30 02:12:58 noro Exp $ */ /* requires 'primdec' */ +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) @@ -212,6 +234,9 @@ 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) */ @@ -267,6 +292,70 @@ def generic_bfct(F,V,DV,W) return B; } +/* all term reduction + interreduce */ +def generic_bfct_1(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); @@ -348,6 +437,174 @@ def bfct_via_gbfct(F) return subst(R,s,-s-1); } +/* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */ + +def bfct_via_gbfct_weight(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); + WtV = newvect(2*(N+1)+1); + WtV[0] = Tdeg; + WtV[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; + } + WtV[2*(N+1)] = 1; + dp_set_weight(WtV); + 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_1(B,V1,DV1,W); + dp_set_weight(0); + return subst(R,s,-s-1); +} + +/* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */ + +def bfct_via_gbfct_weight_1(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); + WtV = newvect(2*(N+1)+1); + /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */ + for ( I = 0; I < N; I++ ) { + WtV[I] = Wt[I]; + WtV[N+1+I] = Tdeg-Wt[I]+1; + } + WtV[N] = Tdeg; + WtV[2*N+1] = 1; + WtV[2*(N+1)] = 1; + dp_set_weight(WtV); + 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 = append(V,[t]); DV1 = append(DV,[dt]); + W = newvect(N+1); + W[N] = 1; + 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 = [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); + V2 = append(V,[t]); DV2 = append(DV,[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); +} + def weyl_minipolym(G,V,O,M,V0) { N = length(V); @@ -366,6 +623,7 @@ def weyl_minipolym(G,V,O,M,V0) GI = cons(I,GI); U = dp_mod(dp_ptod(V0,V),M,[]); + U = dp_weyl_nf_mod(GI,U,PS,1,M); T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]); @@ -400,20 +658,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); @@ -422,7 +688,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]); @@ -545,6 +822,17 @@ def v_factorial(V,N) { for ( J = N-1, R = 1; J >= 0; J-- ) R *= V-J; + return R; +} + +def w_tdeg(F,V,W) +{ + dp_set_weight(newvect(length(W),W)); + T = dp_ptod(F,V); + for ( R = 0; T; T = cdr(T) ) { + D = dp_td(T); + if ( D > R ) R = D; + } return R; } end$