===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v
retrieving revision 1.9
retrieving revision 1.15
diff -u -p -r1.9 -r1.15
--- OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2009/11/12 01:39:54	1.9
+++ OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2010/06/19 08:32:37	1.15
@@ -1,3 +1,4 @@
+/* $OpenXM$ */
 /* requires 'primdec' */
 
 #define TMP_H hhhhhhhh
@@ -12,6 +13,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 +23,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$
@@ -43,9 +45,13 @@ localf ideal_intersection$
 
 def bfunction(F)
 {
+	/* 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;
+	if ( type(Op=getopt(op)) == -1 ) Op = 0;
 	L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
 	Indata = L[0]; AllData = L[1]; VData = L[2];
 	GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4];
@@ -53,7 +59,34 @@ def bfunction(F)
 	dp_set_weight(W);
 	B = weyl_minipoly(GIN,VDV,0,WVDV);
 	dp_set_weight(0);	
-	return subst(B,s,-s-1);
+	if ( !Op ) return subst(B,s,-s-1);
+
+	V0 = VData[0]; DV0 = VData[1]; T = VData[2]; DT = VData[3];
+	BPT = weyl_subst(B,T*DT,VDV);
+
+	/* computation using G0,GIN0,VDV0 */
+	G0 = AllData[0]; GIN0 = AllData[1]; VDV0 = AllData[2]; WtV0 = AllData[5];
+	dp_set_weight(WtV0); dp_ord(0);
+	PS = map(dp_ptod,GIN0,VDV0); Len = length(PS);
+	for ( I = Len-1, Ind = []; I >= 0; I-- ) Ind = cons(I,Ind);
+	/* QR = [D,M,Coef] */
+	Ax = 1;
+	AxBPT = dp_ptod(Ax*BPT,VDV0);
+	QR = weyl_nf_quo(Ind,AxBPT,1,PS);
+	if ( !weyl_nf_quo_check(AxBPT,PS,QR) ) 
+		error("bfunction : invalid quotient");
+	if ( QR[0] ) error("bfunction : invalid quotient");
+	Den = QR[1]; Coef = QR[2];
+	for ( I = 0, R = Den*AxBPT; I < Len; I++ )
+		R -= dp_weyl_mul(Coef[I],dp_ptod(G0[I],VDV0));
+	R = dp_dtop(R,VDV0);
+	CR = conv_tdt(R,F,V0,DV0,T,DT);
+
+	dp_set_weight(0);
+	Cont = cont(CR); CR /= Cont; 
+	Cont *= dn(Fcont); Den *= nm(Fcont);
+	Gcd = igcd(Den,Cont);
+	return [subst(B,s,-s-1),(Cont*CR)/(Den*Ax)];
 }
 
 /*
@@ -469,6 +502,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]);
@@ -1182,6 +1261,9 @@ def replace_var(V,X,Y)
 
 def action_on_gfs(P,V,GFS)
 {
+	for ( T = V, DV = []; T != []; T = cdr(T) )
+		DV = cons(strtov("d"+rtostr(car(T))),DV);
+	V = append(append(V,[s]),reverse(cons(ds,DV)));
 	DP = dp_ptod(P,V);	
 	N = length(V)/2;
 	for ( I = N-1, V0 = []; I >= 0; I-- )
@@ -1323,7 +1405,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)