===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v
retrieving revision 1.18
retrieving revision 1.21
diff -u -p -r1.18 -r1.21
--- OpenXM_contrib2/asir2000/lib/bfct	2002/01/28 02:42:27	1.18
+++ OpenXM_contrib2/asir2000/lib/bfct	2002/01/30 02:12:58	1.21
@@ -45,7 +45,7 @@
  * 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.20 2002/01/29 05:37:12 noro Exp $
  */
 /* requires 'primdec' */
 
@@ -410,7 +410,7 @@ def bfct_via_gbfct(F)
 	V1 = cons(t,V); DV1 = cons(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);
 }
@@ -484,11 +484,105 @@ def bfct_via_gbfct_weight_1(F,V)
 	V1 = append(V,[t]); DV1 = append(DV,[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 = [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<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);
+}
+
 def weyl_minipolym(G,V,O,M,V0)
 {
 	N = length(V);
@@ -542,20 +636,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 +666,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]);