===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v
retrieving revision 1.3
retrieving revision 1.11
diff -u -p -r1.3 -r1.11
--- OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2009/10/13 01:28:00	1.3
+++ OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2010/04/28 15:27:40	1.11
@@ -12,6 +12,7 @@
 
 if (!module_definedp("gr")) load("gr")$ else{ }$
 if (!module_definedp("primdec")) load("primdec")$ else{ }$
+if (!module_definedp("newsyz")) load("noro_module_syz.rr")$ else{ }$
   /* Empty for now. It will be used in a future. */
 
 /* toplevel */
@@ -21,7 +22,7 @@ module ndbf$
 /* bfunction */
 
 localf bfunction, in_ww, in_ww_main, ann, ann_n$
-localf ann0, psi, ww_weight, compare_first, generic_bfct$
+localf ann0, ann_fa, psi, ww_weight, compare_first, generic_bfct$
 localf generic_bfct_1, initial_part, bfct, indicial1, bfct_via_gbfct$
 localf bfct_via_gbfct_weight, bfct_via_gbfct_weight_1, bfct_via_gbfct_weight_2$
 localf weyl_minipolym, weyl_minipoly, weyl_nf, weyl_nf_quo_check$
@@ -29,6 +30,7 @@ localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial,
 localf replace_vars_f, replace_vars_v, replace_var$
 localf action_on_gfs, action_on_gfs_1$
 localf nd_gb_candidate$
+localf in_gb_oaku$
 
 /* stratification */
 
@@ -263,15 +265,30 @@ def ann(F)
 {
 	if ( member(s,vars(F)) )
 		error("ann : the variable 's' is reserved.");
+	F = ptozp(F);
 	V = vars(F);
 	N = length(V);
 	D = newvect(N);
+	if ( type(Wt=getopt(weight)) == -1 )
+		for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt);
 
-	for ( I = 0; I < N; I++ )
-		D[I] = [deg(F,V[I]),V[I]];
+	Wt1 = vector(N);
+	for ( I = 0, F1 =F; I < N; I++ ) {
+		VI = Wt[2*I]; WI = Wt[2*I+1];
+		for ( J = 0; J < N; J++ )
+			if ( VI == V[J] ) break;
+		F1 = subst(F1,VI,VI^WI);
+	}
+	for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]];
 	qsort(D,compare_first);
-	for ( V = [], I = N-1; I >= 0; I-- )
-		V = cons(D[I][1],V);
+	for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V);
+	V = reverse(V); 
+	for ( I = 0; I < N; I++ ) {
+		VI = Wt[2*I]; WI = Wt[2*I+1];
+		for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break;
+		Wt1[J] = WI;
+	}
+	Wt = vtol(Wt1);
 
 	for ( I = N-1, DV = []; I >= 0; I-- )
 		DV = cons(strtov("d"+rtostr(V[I])),DV);
@@ -284,9 +301,24 @@ def ann(F)
 		B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
 	}
 
-	/* homogenized (heuristics) */
-	dp_nelim(2);
-	G0 = nd_weyl_gr(B,append(W,DW),0,[[0,2],[0,length(W)*2-2]]);
+	Tdeg = w_tdeg(F,V,Wt);
+	/* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */
+	/* weight for [y1,y2,t,   x1,...,xn,  dy1,dy2, dt,dx1,...,dxn,   h]   */
+	/*              0  1 2    3      N3-1 N3  N3+1 N3+2              2*N3 */
+	/*              1  1 D+1  w1     wn    1   1    1  D       D      1    */
+	N3 = N+3;
+	WtV = newvect(2*N3+1);
+	WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1; 
+	for ( I = 3; I < N3; I++ ) WtV[I] = Wt[I-3];
+	for ( ; I <= N3+2; I++ ) WtV[I] = 1;
+	for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg;
+	WtV[2*N3] = 1;
+
+	/* B is already a GB => modular change of ordering can be applied */
+	/* any prime is available => HC=[1] */
+	dp_set_weight(WtV);
+	G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1);
+	dp_set_weight(0);
 	G1 = [];
 	for ( T = G0; T != []; T = cdr(T) ) {
 		E = car(T); VL = vars(E);
@@ -298,6 +330,67 @@ def ann(F)
 	return G3;
 }
 
+def in_gb_oaku(F)
+{
+	if ( member(s,vars(F)) )
+		error("ann : the variable 's' is reserved.");
+	F = ptozp(F);
+	V = vars(F);
+	N = length(V);
+	D = newvect(N);
+	if ( type(Wt=getopt(weight)) == -1 )
+		for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt);
+
+	Wt1 = vector(N);
+	for ( I = 0, F1 =F; I < N; I++ ) {
+		VI = Wt[2*I]; WI = Wt[2*I+1];
+		for ( J = 0; J < N; J++ )
+			if ( VI == V[J] ) break;
+		F1 = subst(F1,VI,VI^WI);
+	}
+	for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]];
+	qsort(D,compare_first);
+	for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V);
+	V = reverse(V); 
+	for ( I = 0; I < N; I++ ) {
+		VI = Wt[2*I]; WI = Wt[2*I+1];
+		for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break;
+		Wt1[J] = WI;
+	}
+	Wt = vtol(Wt1);
+
+	for ( I = N-1, DV = []; I >= 0; I-- )
+		DV = cons(strtov("d"+rtostr(V[I])),DV);
+
+	W = append([TMP_Y1,TMP_Y2,TMP_T],V);
+	DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV);
+
+	B = [TMP_T-TMP_Y1*F];
+	for ( I = 0; I < N; I++ ) {
+		B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
+	}
+
+	Tdeg = w_tdeg(F,V,Wt);
+	/* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */
+	/* weight for [y1,y2,t,   x1,...,xn,  dy1,dy2, dt,dx1,...,dxn,   h]   */
+	/*              0  1 2    3      N3-1 N3  N3+1 N3+2              2*N3 */
+	/*              1  1 D+1  1      1    1   1    1  D       D      1    */
+	N3 = N+3;
+	WtV = newvect(2*N3+1);
+	WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1; 
+	for ( I = 3; I <= N3+2; I++ ) WtV[I] = 1;
+	for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg;
+	WtV[2*N3] = 1;
+
+	/* B is already a GB => modular change of ordering can be applied */
+	/* any prime is available => HC=[1] */
+	dp_set_weight(WtV);
+	G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1);
+	dp_set_weight(0);
+	G1 = map(subst,G0,TMP_Y1,1);
+	return [G1,append(V,DV)];
+}
+
 /* F = [F0,F1,...] */
 
 def ann_n(F)
@@ -335,7 +428,11 @@ def ann_n(F)
 	VA = append(U,append(W,V));
 	DVA = append(DU,append(DW,DV));
 	VDV = append(VA,DVA);
+#if 0
 	G0 = nd_weyl_gr(B,VDV,0,[[0,2*L],[0,length(VDV)-2*L]]);
+#else
+	G0 = nd_gb_candidate(B,VDV,[[0,2*L],[0,length(VDV)-2*L]],0,[1],1);
+#endif
 	G1 = [];
 	for ( T = G0; T != []; T = cdr(T) ) {
 		E = car(T); VL = vars(E);
@@ -373,6 +470,52 @@ def ann0(F)
 	return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)];
 }
 
+/*
+ * For a polynomial F and a scalar A,
+ * compute generators of Ann(F^A).
+ */
+
+def ann_fa(F,A)
+{
+	if ( type(Syz=getopt(syz)) == -1 ) Syz = 0;
+	
+	F = subst(F,s,TMP_S);
+	Ann = ann(F);
+	Bf = bfunction(F);
+
+	FList = cdr(fctr(Bf));
+	for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
+		LF = car(car(T));
+		Root = -coef(LF,0)/coef(LF,1);
+		if ( dn(Root) == 1 && Root < Min )
+			Min = Root;
+	}
+	D = A-Min;
+	if ( dn(D) != 1 || D <= 0 )
+		return map(ptozp,map(subst,Ann,s,A,TMP_S,s,TMP_DS,ds));
+	
+	V = vars(F);
+	for ( I = length(V)-1, DV = []; I >= 0; I-- )
+		DV = cons(strtov("d"+rtostr(V[I])),DV);
+	VDV = append(V,DV);
+	R = map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds);
+	F = ptozp(F);
+	
+	if ( Syz ) {
+		/* syzygy method */
+		S = newsyz.module_syz(cons(F^D,R),VDV,1,0|weyl=1);
+		B = car(S);
+		for ( R = []; B != []; B = cdr(B) )
+			if ( H = car(car(B)) )
+				R = cons(ptozp(H),R);
+	} else {
+		/* colon method */
+		for ( I = 0; I < D; I++ )
+			R = weyl_ideal_quotient(R,F,VDV);
+	}
+	return R;
+}
+
 def psi0(F,T,DT)
 {
 	D = dp_ptod(F,[T,DT]);
@@ -1227,7 +1370,7 @@ def weyl_ideal_quotient(B,F,VDV)
 	O1 = [[0,1],[0,N+1]];
 	GJ = nd_weyl_gr(J,VDV1,0,O1);
 	R = elimination(GJ,VDV);
-	return R;
+	return map(weyl_divide_by_right,R,F,VDV,0);
 }
 
 def bf_strat(F)
@@ -1354,9 +1497,15 @@ def bf_strat_stage3(L)
 		}
 		dp_set_weight(0);
 		if ( K < J ) {
-			for ( L = 0, T = []; L < NDI; L++ )
+			for ( L = 0, T = []; L < NDI; L++ ) {
+#if 0
+				NewId = DK[L][1];
+#else
+				NewId = ideal_intersection(DK[L][1],DI[L][1],V0,0);
+#endif
 				T = cons([[DK[L][0][0]*DI[L][0][0],DK[L][0][1]],
-					DK[L][1],DK[L][2]],T);
+					NewId,DK[L][2]],T);
+			}
 			Data[K] = reverse(T);
 		} else
 			Data[J++] = DI;
@@ -1379,6 +1528,8 @@ def bf_strat_stage3(L)
 
 def bf_local(F,P)
 {
+	/* F -> F/Fcont */
+	F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1;
 	if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
 	if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
 	if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
@@ -1443,7 +1594,10 @@ def bf_local(F,P)
 	CR = conv_tdt(R,F,V0,DV0,T,DT);
 
 	dp_set_weight(0);
-	return [LB,Ax,CR];
+	Cont = cont(CR); CR /= Cont; 
+	Cont *= dn(Fcont); Den *= nm(Fcont);
+	Gcd = igcd(Den,Cont);
+	return [LB,(Den/Gcd)*Ax,(Cont/Gcd)*CR];
 }
 
 /* t^(l+k)*dt^l (k>l) -> (s-k)(s-k-1)...(s-(k+l-1))t^k */