===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/dp-supp.c,v
retrieving revision 1.46
retrieving revision 1.50
diff -u -p -r1.46 -r1.50
--- OpenXM_contrib2/asir2000/builtin/dp-supp.c	2007/09/19 05:42:59	1.46
+++ OpenXM_contrib2/asir2000/builtin/dp-supp.c	2009/01/04 10:02:00	1.50
@@ -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.45 2007/09/16 09:08:25 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/dp-supp.c,v 1.49 2009/01/04 05:44:51 noro Exp $ 
 */
 #include "ca.h"
 #include "base.h"
@@ -880,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;
 }
@@ -1231,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)
@@ -1626,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 ) {
@@ -1664,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;
 }
@@ -2127,17 +2193,18 @@ void homogenize_order(struct order_spec *old,int n,str
 					error("homogenize_order : invalid input");
 			}
 			break;
-		case 1:
+		case 1: case 257:
 			length = old->ord.block.length;
 			l = (struct order_pair *)
 				MALLOC_ATOMIC((length+1)*sizeof(struct order_pair));
 			bcopy((char *)old->ord.block.order_pair,(char *)l,length*sizeof(struct order_pair));
 			l[length].order = 2; l[length].length = 1;
-			new->id = 1; new->nv = n+1;
+			new->id = old->id; new->nv = n+1;
 			new->ord.block.order_pair = l;
 			new->ord.block.length = length+1;
+			new->istop = old->istop;
 			break;
-		case 2:
+		case 2: case 258:
 			nv = old->nv; row = old->ord.matrix.row;
 			oldm = old->ord.matrix.matrix; newm = almat(row+1,nv+1);
 			for ( i = 0; i <= nv; i++ )
@@ -2147,10 +2214,11 @@ void homogenize_order(struct order_spec *old,int n,str
 					newm[i+1][j] = oldm[i][j];
 				newm[i+1][j] = 0;
 			}
-			new->id = 2; new->nv = nv+1;
+			new->id = old->id; new->nv = nv+1;
 			new->ord.matrix.row = row+1; new->ord.matrix.matrix = newm;
+			new->istop = old->istop;
 			break;
-		case 3:
+		case 3: case 259:
 			onv = old->nv;
 			nnv = onv+1;
 			olen = old->ord.composite.length;
@@ -2184,12 +2252,33 @@ void homogenize_order(struct order_spec *old,int n,str
 				(struct sparse_weight *)MALLOC(sizeof(struct sparse_weight));
 			nwb[i].body.sparse_weight[0].pos = onv;
 			nwb[i].body.sparse_weight[0].value = 1;
-			new->id = 3;
+			new->id = old->id;
 			new->nv = nnv;
 			new->ord.composite.length = nlen;
 			new->ord.composite.w_or_b = nwb;
+			new->istop = old->istop;
 			print_composite_order_spec(new);
 			break;
+		case 256: /* simple module order */
+			switch ( old->ord.simple ) {
+				case 0:
+					new->id = 256; new->ord.simple = 0; break;
+				case 1:
+					l = (struct order_pair *)
+						MALLOC_ATOMIC(2*sizeof(struct order_pair));
+					l[0].length = n; l[0].order = old->ord.simple;
+					l[1].length = 1; l[1].order = 2;
+					new->id = 257;
+					new->ord.block.order_pair = l;
+					new->ord.block.length = 2; new->nv = n+1;
+					break;
+				case 2:
+					new->id = 256; new->ord.simple = 1; break;
+				default:
+					error("homogenize_order : invalid input");
+			}
+			new->istop = old->istop;
+			break;
 		default:
 			error("homogenize_order : invalid input");
 	}
@@ -2669,6 +2758,95 @@ int compare_facet_preorder(int n,int *u,int *v,
 		else if ( t < 0 ) return 0;
 	}
 	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 */