===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/dp-supp.c,v
retrieving revision 1.44
retrieving revision 1.49
diff -u -p -r1.44 -r1.49
--- OpenXM_contrib2/asir2000/builtin/dp-supp.c	2007/09/15 10:17:08	1.44
+++ OpenXM_contrib2/asir2000/builtin/dp-supp.c	2009/01/04 05:44:51	1.49
@@ -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/builtin/dp-supp.c,v 1.43 2007/09/07 00:45:50 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.48 2007/10/21 07:47:59 noro Exp $ 
 */
 #include "ca.h"
 #include "base.h"
@@ -71,6 +71,30 @@ void print_composite_order_spec(struct order_spec *spe
  *
  */
 
+static NODE RatDenomList;
+
+void init_denomlist()
+{
+	RatDenomList = 0;
+}
+
+void add_denomlist(P f)
+{
+	NODE n;
+
+	if ( OID(f)==O_P ) {
+		MKNODE(n,f,RatDenomList); RatDenomList = n;
+	}
+}
+
+LIST get_denomlist()
+{
+	LIST l;
+
+	MKLIST(l,RatDenomList); RatDenomList = 0;
+	return l;
+}
+
 void dp_ptozp(DP p,DP *rp)
 {
 	MP m,mr,mr0;
@@ -448,6 +472,7 @@ void dp_prim(DP p,DP *rp)
 	P *w;
 	Q *c;
 	Q dvr;
+	NODE tn;
 
 	if ( !p )
 		*rp = 0;
@@ -502,6 +527,7 @@ void dp_prim(DP p,DP *rp)
 				NEXTMP(mr0,mr); divsp(CO,m->c,g,&mr->c); mr->dl = m->dl;
 			}
 			NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar;
+			add_denomlist(g);
 		}
 	}
 }
@@ -805,6 +831,7 @@ void dp_red(DP p0,DP p1,DP p2,DP *head,DP *rest,P *dnp
 	} else {
 		ezgcdpz(CO,(P)c1,(P)c2,&g);
 		divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
+		add_denomlist(g);
 	}
 	NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
 	*multp = s;
@@ -853,9 +880,9 @@ void dp_red_marked(DP p0,DP p1,DP p2,DP hp2,DP *head,D
 		ezgcdpz(CO,(P)c1,(P)c2,&g);
 		divsp(CO,(P)c1,g,&a); c1 = (Q)a; divsp(CO,(P)c2,g,&a); c2 = (Q)a;
 	}
-	NEWMP(m); m->dl = d; chsgnp((P)c1,&m->c); NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
+	NEWMP(m); m->dl = d; m->c = (P)c1; NEXT(m) = 0; MKDP(n,m,s); s->sugar = d->td;
 	*multp = s;
-	muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); addd(CO,s,t,&r);
+	muld(CO,s,p2,&t); muldc(CO,p1,(P)c2,&s); subd(CO,s,t,&r);
 	muldc(CO,p0,(P)c2,&h);
 	*head = h; *rest = r; *dnp = (P)c2;
 }
@@ -1204,6 +1231,63 @@ void dp_true_nf_marked_mod(NODE b,DP g,DP *ps,DP *hps,
 	*rp = d; *dnp = dn;
 }
 
+/* true nf by a marked GB and collect quotients */
+
+DP *dp_true_nf_and_quotient_marked (NODE b,DP g,DP *ps,DP *hps,DP *rp,P *dnp)
+{
+	DP u,p,d,s,t,dmy,hp,mult;
+	DP *q;
+	NODE l;
+	MP m,mr;
+	int i,n,j;
+	int *wb;
+	int sugar,psugar,multiple;
+	P nm,tnm1,dn,tdn,tdn1;
+	Q cont;
+
+	dn = (P)ONE;
+	if ( !g ) {
+		*rp = 0; *dnp = dn; return;
+	}
+	for ( n = 0, l = b; l; l = NEXT(l), n++ );
+	wb = (int *)ALLOCA(n*sizeof(int));
+	for ( i = 0, l = b; i < n; l = NEXT(l), i++ )
+		wb[i] = QTOS((Q)BDY(l));
+	q = (DP *)MALLOC(n*sizeof(DP));
+	for ( i = 0; i < n; i++ ) q[i] = 0;
+	sugar = g->sugar;
+	for ( d = 0; g; ) {
+		for ( u = 0, i = 0; i < n; i++ ) {
+			if ( dp_redble(g,hp = hps[wb[i]]) ) {
+				p = ps[wb[i]];
+				dp_red_marked(d,g,p,hp,&t,&u,&tdn,&mult);
+				psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
+				sugar = MAX(sugar,psugar);
+				for ( j = 0; j < n; j++ ) {
+					muldc(CO,q[j],(P)tdn,&dmy); q[j] = dmy;
+				}
+				addd(CO,q[wb[i]],mult,&dmy); q[wb[i]] = dmy;
+				mulp(CO,dn,tdn,&tdn1); dn = tdn1;
+				d = t;
+				if ( !u ) goto last;
+				break;
+			}
+		}
+		if ( u ) {
+			g = u;
+		} else {
+			m = BDY(g); NEWMP(mr); mr->dl = m->dl; mr->c = m->c;
+			NEXT(mr) = 0; MKDP(g->nv,mr,t); t->sugar = mr->dl->td;
+			addd(CO,d,t,&s); d = s;
+			dp_rest(g,&t); g = t;
+		}
+	}
+last:
+	if ( d ) d->sugar = sugar;
+	*rp = d; *dnp = dn;
+	return q;
+}
+
 /* nf computation over Z */
 
 void dp_nf_z(NODE b,DP g,DP *ps,int full,int multiple,DP *rp)
@@ -1599,7 +1683,8 @@ int create_order_spec(VL vl,Obj obj,struct order_spec 
 	struct order_pair *l;
 	NODE node,t,tn;
 	MAT m;
-	pointer **b;
+	VECT v;
+	pointer **b,*bv;
 	int **w;
 
 	if ( vl && obj && OID(obj) == O_LIST ) {
@@ -1637,6 +1722,14 @@ int create_order_spec(VL vl,Obj obj,struct order_spec 
 		spec->nv = col; spec->ord.matrix.row = row;
 		spec->ord.matrix.matrix = w;
 		return 1;
+	} else if ( OID(obj) == O_VECT ) {
+		v = (VECT)obj; bv = BDY(v);
+		if ( v->len < 2 ) error("create_order_spec : invalid argument");
+		create_order_spec(0,(Obj)bv[1],&spec);
+		spec->id += 256; spec->obj = obj;
+		spec->istop = bv[0]==0;
+		*specp = spec;
+		return 1;
 	} else
 		return 0;
 }
@@ -2644,6 +2737,95 @@ int compare_facet_preorder(int n,int *u,int *v,
 	return 1;
 }
 
+Q inner_product_with_small_vector(VECT w,int *v)
+{
+	int n,i;
+	Q q,s,t,u;
+
+	n = w->len;
+	s = 0;
+	for ( i = 0; i < n; i++ ) {
+		STOQ(v[i],q); mulq((Q)w->body[i],q,&t); addq(t,s,&u); s = u;
+	}
+	return s;
+}
+
+Q compute_last_t(NODE g,NODE gh,Q t,VECT w1,VECT w2,NODE *homo,VECT *wp)
+{
+	int n,i;
+	int *wt;
+	Q last,d1,d2,dn,nm,s,t1;
+	VECT wd,wt1,wt2,w;
+	NODE tg,tgh;
+	MP f;
+	int *h;
+	NODE r0,r;
+	MP m0,m;
+	DP d;
+
+	n = w1->len;
+	wt = W_ALLOC(n);
+	last = ONE;
+	/* t1 = 1-t */
+	for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
+		f = BDY((DP)BDY(tg));
+		h = BDY((DP)BDY(tgh))->dl->d;
+		for ( ; f; f = NEXT(f) ) {
+			for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
+			for ( i = 0; i < n && !wt[i]; i++ );
+			if ( i == n ) continue;
+			d1 = inner_product_with_small_vector(w1,wt);
+			d2 = inner_product_with_small_vector(w2,wt);
+			nm = d1; subq(d1,d2,&dn);
+			/* if d1=d2 then nothing happens */
+			if ( !dn ) continue;
+			/* s satisfies ds = 0*/
+			divq(nm,dn,&s);
+
+			if ( cmpq(s,t) > 0 && cmpq(s,last) < 0 ) 
+				last = s;
+			else if ( !cmpq(s,t) ) {
+				if ( cmpq(d2,0) < 0 ) {
+					last = t;
+					break;
+				}
+			}
+		}
+	}
+	if ( !last ) {
+		dn = ONE; nm = 0;
+	} else {
+		NTOQ(NM(last),1,nm);
+		if ( INT(last) ) dn = ONE;
+		else {
+			NTOQ(DN(last),1,dn);
+		}
+	}
+	/* (1-n/d)*w1+n/d*w2 -> w=(d-n)*w1+n*w2 */
+	subq(dn,nm,&t1); mulvect(CO,(Obj)w1,(Obj)t1,(Obj *)&wt1); 
+	mulvect(CO,(Obj)w2,(Obj)nm,(Obj *)&wt2); addvect(CO,wt1,wt2,&w);
+
+	r0 = 0;
+	for ( tg = g, tgh = gh; tg; tg = NEXT(tg), tgh = NEXT(tgh ) ) {
+		f = BDY((DP)BDY(tg));
+		h = BDY((DP)BDY(tgh))->dl->d;
+		for ( m0 = 0; f; f = NEXT(f) ) {
+			for ( i = 0; i < n; i++ ) wt[i] = h[i]-f->dl->d[i];
+			for ( i = 0; i < n && !wt[i]; i++ );
+			if ( !inner_product_with_small_vector(w,wt) ) {
+				NEXTMP(m0,m); m->c = f->c; m->dl = f->dl;
+			}
+		}
+		NEXT(m) = 0;
+		MKDP(((DP)BDY(tg))->nv,m0,d);  d->sugar = ((DP)BDY(tg))->sugar;
+		NEXTNODE(r0,r); BDY(r) = (pointer)d;
+	}
+	NEXT(r) = 0;
+	*homo = r0;
+	*wp = w;
+	return last;
+}
+
 /* return 0 if last_w = infty */
 
 NODE compute_last_w(NODE g,NODE gh,int n,int **w,
@@ -2708,86 +2890,67 @@ NODE compute_last_w(NODE g,NODE gh,int n,int **w,
 }
 
 /* compute a sufficient set of d(f)=u-v */
-static int comp_vector_lex_nv;
 
-int comp_vector_lex(int **a,int **b)
-{
-	int i;
-	int *pa,*pb;
-
-	pa = *a; pb = *b;
-	for ( i = 0; i < comp_vector_lex_nv; i++ )
-		if ( pa[i] < pb[i] ) return -1;
-		else if ( pa[i] > pb[i] ) return 1;
-	return 0;
-}
-
 NODE compute_essential_df(DP *g,DP *gh,int ng)
 {
-	VECT v;
-	Q q;
+	int nv,i,j,k,t,lj;
+	NODE r,r1,ri,rt,r0;
 	MP m;
-	NODE r,r1;
-	int nv,len,i,j,k;
-	int *p,*dm,*mi,*mj,*h;
-	int **mat;
+	MP *mj;
+	DL di,hj,dl,dlt;
+	int *d,*dt;
+	LIST l;
+	Q q;
 
-	nv = comp_vector_lex_nv = g[0]->nv;
-	for ( len = 0, j = 0; j < ng; j++ ) {
-		for ( m = BDY(g[j]); m; m = NEXT(m), len++ );
-	}
-	mat = almat(len,nv);
-	for ( i = 0, j = 0; j < ng; j++ ) {
-		h = BDY(gh[j])->dl->d;
-		for ( m = BDY(g[j]); m; m = NEXT(m) ) {
-			dm = m->dl->d;
-			for ( k = 0; k < nv; k++ )
-				if ( dm[k] ) break;
-			if ( k == nv ) continue;
-			else {
-				p = mat[i];
-				for ( k = 0; k < nv; k++ )
-					p[k] = h[k]-dm[k];
-				i++;
+	nv = g[0]->nv;
+	r = 0;
+	for ( j = 0; j < ng; j++ ) {
+		for ( m = BDY(g[j]), lj = 0; m; m = NEXT(m), lj++ );
+		mj = (MP *)ALLOCA(lj*sizeof(MP));
+		for ( m = BDY(g[j]), k = 0; m; m = NEXT(m), k++ )
+			mj[k] = m;
+		for ( i = 0; i < lj; i++ ) {
+			for ( di = mj[i]->dl, k = i+1; k < lj; k++ )
+				if ( _dl_redble(di,mj[k]->dl,nv) ) break;
+			if ( k < lj ) mj[i] = 0;
+		}
+		hj = BDY(gh[j])->dl;
+		_NEWDL(dl,nv); d = dl->d;
+		r0 = r;
+		for ( i = 0; i < lj; i++ ) {
+			if ( mj[i] && !dl_equal(nv,di=mj[i]->dl,hj) ) {
+				for ( k = 0, t = 0; k < nv; k++ ) {
+					d[k] = hj->d[k]-di->d[k];
+					t += d[k];
+				}
+				dl->td = t;	
+#if 1
+				for ( rt = r0; rt; rt = NEXT(rt) ) {
+					dlt = (DL)BDY(rt);
+					if ( dlt->td != dl->td ) continue;
+					for ( dt = dlt->d, k = 0; k < nv; k++ )
+						if ( d[k] != dt[k] ) break;
+					if ( k == nv ) break;
+				}
+#else
+				rt = 0;
+#endif
+				if ( !rt ) {
+					MKNODE(r1,dl,r); r = r1;
+					_NEWDL(dl,nv); d = dl->d;
+				}
 			}
 		}
 	}
-	len = i;
-	qsort(mat,len,sizeof(int *),
-		(int (*)(const void *,const void *))comp_vector_lex);
-	for ( i = 0; i < len; i++ ) {
-		for ( j = 0; j < nv; j++ )
-			printf("%d ",mat[i][j]);
-		printf("\n");
-	}
-	for ( i = 0; i < len; i++ ) {
-		mi = mat[i];
-		if ( !mi ) continue;
-		for ( j = i+1; j < len; j++ ) {
-			mj = mat[j];
-			if ( !mj ) continue;
-			for ( k = 0; k < nv; k++ )
-				if ( mi[k] > mj[k] )  break;
-			if ( k == nv ) mat[j] = 0;
+	for ( rt = r; rt; rt = NEXT(rt) ) {
+		dl = (DL)BDY(rt); d = dl->d;
+		ri = 0;
+		for ( k = nv-1; k >= 0; k-- ) {
+			STOQ(d[k],q);
+			MKNODE(r1,q,ri); ri = r1;
 		}
-	}
-	for ( i = 0; i < len; i++ ) {
-		if ( mat[i] ) {
-			for ( j = 0; j < nv; j++ )
-				printf("%d ",mat[i][j]);
-			printf("\n");
-		}
-	}
-	r = 0;
-	for ( i = 0; i < len; i++ ) {
-		if ( mi = mat[i] ) {
-			MKVECT(v,nv);
-			for ( k = 0; k < nv; k++ ) {
-				STOQ(mi[k],q);
-				v->body[k] = (pointer)q;
-			}
-			MKNODE(r1,v,r); r = r1;
-		}
+		MKNODE(r1,0,ri); MKLIST(l,r1);
+		BDY(rt) = (pointer)l;
 	}
 	return r;
 }