===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v
retrieving revision 1.19
retrieving revision 1.22
diff -u -p -r1.19 -r1.22
--- OpenXM_contrib2/asir2000/lib/bfct	2002/01/29 02:03:41	1.19
+++ 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.18 2002/01/28 02:42:27 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)
@@ -410,7 +432,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,7 +506,7 @@ 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);
 }
@@ -572,7 +594,11 @@ def bfct_via_gbfct_weight_2(F,V)
 	/* change of ordering from VDV to VDV2 */
 	VDV2 = append(V2,DV2);
 	dp_set_weight(WtV2);
-	GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-1,0);
+	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);
@@ -632,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);
@@ -654,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(<<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]);