===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v
retrieving revision 1.11
retrieving revision 1.19
diff -u -p -r1.11 -r1.19
--- OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2010/04/28 15:27:40	1.11
+++ OpenXM/src/asir-contrib/testing/noro/ndbf.rr	2011/03/30 05:07:01	1.19
@@ -1,10 +1,11 @@
+/* $OpenXM: OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v 1.18 2011/01/16 08:46:10 noro Exp $ */
 /* requires 'primdec' */
 
 #define TMP_H hhhhhhhh
 #define TMP_S ssssssss
 #define TMP_DS dssssssss
-#define TMP_T t
-#define TMP_DT dt
+#define TMP_T tttttttt
+#define TMP_DT dtttttttt
 #define TMP_Y1 yyyyyyyy1
 #define TMP_DY1 dyyyyyyyy1
 #define TMP_Y2 yyyyyyyy2
@@ -30,7 +31,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$
+localf in_gb_oaku, homogenize_oaku$
 
 /* stratification */
 
@@ -44,9 +45,15 @@ localf ideal_intersection$
 
 def bfunction(F)
 {
+	if ( member(s,vars(F)) )
+		error("ann : the variable 's' is reserved.");
+	/* 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];
@@ -54,7 +61,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)];
 }
 
 /*
@@ -265,8 +299,13 @@ def ann(F)
 {
 	if ( member(s,vars(F)) )
 		error("ann : the variable 's' is reserved.");
+	if ( type(Vord=getopt(vord)) == -1 ) Vord = 0;
 	F = ptozp(F);
 	V = vars(F);
+	if ( Vord ) {
+		Param = setminus(V,Vord);
+		V = Vord;
+	}
 	N = length(V);
 	D = newvect(N);
 	if ( type(Wt=getopt(weight)) == -1 )
@@ -391,6 +430,31 @@ def in_gb_oaku(F)
 	return [G1,append(V,DV)];
 }
 
+/* homogenization w.r.t. (-W,W)-weight */
+/* VDV = [x1,...,xn,dx1,...,dxn] */
+/* homogenize F w.r.t. (W,-W,1) for (x,dx,y) */
+
+def homogenize_oaku(F,VDV,W,Y)
+{
+	N = length(VDV);
+	if ( N%2 ) error("invalid variable list");
+	N2 = N/2;
+	if ( length(W) != N2 ) error("inconsistent weight vector");
+	W0 = dp_set_weight();
+	Wt = append(W,append(vtol(-ltov(W)),[1]));
+	dp_set_weight(Wt);
+	H = homogenize(F,VDV,Y);
+	dp_set_weight(W0);
+	if ( type(Vars=getopt(vars)) != -1 && Vars ) {
+		DY = strtov("d"+rtostr(Y));
+		for ( I = 0, T = VDV, V = []; I < N2; I++, T = cdr(T) )
+			V = cons(car(T),V);
+		T = cons(Y,append(T,[DY]));
+		for ( ; V != []; V = cdr(V) ) T = cons(car(V),T);
+		return [H,T];
+	} else return H;
+}
+
 /* F = [F0,F1,...] */
 
 def ann_n(F)
@@ -1229,6 +1293,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-- )
@@ -1375,11 +1442,14 @@ def weyl_ideal_quotient(B,F,VDV)
 
 def bf_strat(F)
 {
+	if ( member(s,vars(F)) )
+		error("ann : the variable 's' is reserved.");
 	dp_ord(0);
 	T0 = time();
 	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(ElimIdeal=getopt(elimideal)) == -1 ) ElimIdeal = 0;
 	L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
 	T1 = time();
 	print(["in_ww",(T1[0]+T1[1])-(T0[0]+T0[1])]);
@@ -1418,6 +1488,7 @@ def bf_strat(F)
 	}
 
 	L2 = bf_strat_stage2(L);
+	if ( ElimIdeal ) return L2;
 	S = bf_strat_stage3(L2);
 	R = [];
 	for ( T = S; T != []; T = cdr(T) ) {
@@ -1528,6 +1599,8 @@ def bf_strat_stage3(L)
 
 def bf_local(F,P)
 {
+	if ( member(s,vars(F)) )
+		error("ann : the variable 's' is reserved.");
 	/* F -> F/Fcont */
 	F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1;
 	if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;