===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v
retrieving revision 1.14
retrieving revision 1.25
diff -u -p -r1.14 -r1.25
--- OpenXM_contrib2/asir2000/lib/bfct 2001/01/10 04:30:35 1.14
+++ OpenXM_contrib2/asir2000/lib/bfct 2003/04/28 03:02:52 1.25
@@ -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.13 2000/12/27 07:17:39 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.24 2003/04/28 02:15:30 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;
}
@@ -95,47 +128,10 @@ def ann(F)
def ann0(F)
{
- V = vars(F);
- N = length(V);
- D = newvect(N);
+ F = subst(F,s,TMP_S);
+ Ann = ann(F);
+ Bf = bfunction(F);
- 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));
@@ -143,30 +139,9 @@ def ann0(F)
if ( dn(Root) == 1 && Root < Min )
Min = Root;
}
- return [Min,map(subst,G3,s,Min)];
+ return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)];
}
-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]);
@@ -212,6 +187,11 @@ 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<x,d> (DRL) */
M = newmat(N2,N2);
for ( J = 0; J < N2; J++ )
@@ -244,7 +224,9 @@ def generic_bfct(F,V,DV,W)
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);
+ 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);
@@ -263,6 +245,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<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 */
+/* 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<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);
@@ -296,6 +342,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);
@@ -315,6 +364,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)
@@ -332,11 +404,11 @@ 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(B,V1,DV1,W);
@@ -344,6 +416,176 @@ 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 = [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);
+ 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 = [TMP_T-F];
+ for ( I = 0; I < N; I++ ) {
+ B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
+ }
+ V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_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 = [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<x,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<x,d>) = D<x,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);
@@ -362,6 +604,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,[]);
@@ -384,6 +627,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);
@@ -396,20 +641,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);
@@ -418,7 +671,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(<<D>>,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]);
@@ -510,13 +774,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) )
@@ -541,6 +798,44 @@ 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;
+}
+
+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;
}
end$