===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v
retrieving revision 1.20
retrieving revision 1.26
diff -u -p -r1.20 -r1.26
--- OpenXM_contrib2/asir2000/lib/bfct	2002/01/29 05:37:12	1.20
+++ OpenXM_contrib2/asir2000/lib/bfct	2003/10/20 00:58:47	1.26
@@ -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.19 2002/01/29 02:03:41 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.25 2003/04/28 03:02:52 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
+
+if (!module_definedp("gr")) load("gr") $$
+if (!module_definedp("primdec")) load("primdec") $$
+module bfct $
+  /* Empty for now. It will be used in a future. */
+endmodule $
+
+/* 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]);
@@ -249,9 +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 */
-	dp_gr_flags(["Top",1]);
+	dp_gr_flags(["Top",1,"NoRA",1]);
 	GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
-	dp_gr_flags(["Top",0]);
+	dp_gr_flags(["Top",0,"NoRA",0]);
 
 	/* dehomigenize GH */
 	G = map(subst,GH,h,1);
@@ -367,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);
@@ -386,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)
@@ -403,14 +404,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 +441,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,14 +478,14 @@ 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);
 }
@@ -527,12 +528,12 @@ def bfct_via_gbfct_weight_2(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);
-	V2 = append(V,[t]); DV2 = append(DV,[dt]);
+	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);
 
@@ -583,6 +584,8 @@ def bfct_via_gbfct_weight_2(F,V)
 	return subst(R,s,-s-1);
 }
 
+/* minimal polynomial of s; modular computation */
+
 def weyl_minipolym(G,V,O,M,V0)
 {
 	N = length(V);
@@ -624,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);
@@ -769,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) )
@@ -810,6 +808,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;
 }